The Tabula Muris Senis compendium is a single-cell survey across 23 C57BL/6JN mice tissues. From its website (https://tabula-muris-senis.ds.czbiohub.org/all/droplet/), one of the two previously accepted pericyte markers, Cspg4 or Pdgfrb, was positively expressed mostly in lung, heart, kidney, and bladder (Supplementary Figure 1). We therefore moved forward to analyze scRNA-seq data from lung (n = 16 mice; 24,540 cells), heart (n = 11; 8,613), kidney (n = 16; 21,647), and bladder (n = 8; 8,945) obtained from the Senis (Figure 1A) (16). Unsupervised clustering revealed substantial cellular heterogeneity within each tissue, with the number of distinct cell clusters ranging from 15 in the bladder to 27 in the kidney. To identify pericyte-enriched clusters, we first considered the individual and joint expression of two previously accepted pericyte markers, Cspg4 and Pdgfrb. Consistent with prior observations (10, 11, 14), the preponderance of cells expressing either Cspg4 and Pdgfrb localized to a limited number of clusters (one each in lung and kidney, two in heart and bladder). In contrast to all other clusters, which either consisted of few Cspg4 or Pdgfrb-expressing cells or had low/undetectable expression of both markers, these clusters were tentatively deemed pericyte-enriched because of their large portions of cells expressing high levels of both Cspg4 and Pdgfrb (Figure 1A). However, within these pericyte-enriched clusters, substantial heterogeneity in the expression of Cspg4 and Pdgfrb was observed. Therefore, within these pericyte-enriched clusters, we defined a set of “stringent pericytes” by considering only those cells that expressed both Cspg4 and Pdgfrb (Figure 1B). A total of 393 cells met these criteria (28.4% of all cells), including 28 in the lung, 219 in the heart, 71 in the kidney, and 75 in the bladder, with the proportion of stringent pericytes within each cluster varying widely by tissue (highest in heart, 55.8%; lowest in the lung, 9.8%). These refined subsets were considered tissue-specific stringent pericytes in all downstream analyses.

We conducted a differential expression (DE) analysis between the stringent pericytes and all other cells based on the Wilcoxon rank-sum test to identify potential pericyte-specific markers. We identified transcripts that were enriched in the stringent pericytes in each tissue, prioritizing transcripts with high pericyte DE (adjusted p-value < 0.05) that were expressed (count > 0) in at least 80% of the stringent pericytes but <5% of all other cells. Based on these criteria, we identified 18, 4, 22, and 4 potential pericyte markers for the lung, heart, kidney, and bladder, respectively (Figure 2, left). No marker fulfilled these criteria in all four cell types. However, several markers were common to two or three tissues (Figure 2, right), including Cox4i2 (in lung, kidney, and bladder); Gucy1a3, Ndufa4i2, and Myl9 (in lung and kidney); Gm13889 and Notch3 (in heart and kidney); Pcp4l1 (in kidney and bladder); and Higd1b (in lung and heart). Also, for each cluster, we determined the percentage of cells expressing the identified markers (Figure 2, left and Supplementary Figure 6). The pct3 represents the maximum percentage found across the clusters, excluding the pericyte enriched cluster(s). When considering markers with pct3 below 10%, we identified Kcnk3 (in the lung); Rgs4 (in the heart); Myh11 and Kcna5 (in the kidney); Pcp4l1 (in the bladder); and Higd1b (in lung and heart), which showed that they were very specific to pericyte-enriched clusters.

To determine pathways associated with stringent pericytes across the lung, heart, kidney, and bladder, we conducted a pre-ranked Gene Set Enrichment Analysis (GSEA) (17) (Figure 3). First, genes were sorted based on their average log2(Fold Change) when comparing their expression in tissue-specific stringent pericytes to all cells in the corresponding tissue. GSEA was then performed on the resulting pre-ranked list using the Kyoto Encyclopedia of Genes and Genomes pathway database (18). We identified 7, 20, 36, and 8 pathways enriched (adjusted q-value < 0.05) in the stringent pericytes of the lung, heart, kidney, and bladder, respectively.

Vascular smooth muscle contraction dilated cardiomyopathy, and focal adhesion were top enriched pathways in the stringent pericytes of all four tissues, suggesting pathways involved in pericyte physiological and metabolic activity (Figure 3). This may be because pericytes share a similar lineage with smooth muscle cells (SMCs) and differentiate into SMCs under stress. Additionally, as broadly defined mural cells, their biological function can be identical to vital regulators of vessel structure support, permeability, and homeostasis (19, 20).

The other identified pathways in heart tissue were involved with neurodegenerative disorders and cardiac muscle contraction. Pericytes are one of the key regulators of the blood-brain barrier. Their involvement in the pathogenesis of neurological disorders may be associated with cardiovascular manifestations under the heart-brain axis (7, 21, 22). Pathways related to cancers are also enriched in the kidney-specific stringent pericytes. Pericytes contribute to the invasion and metastasis of cancer or tissue fibrosis through pericyte-fibroblast transition (PFT) mechanisms, which may be highlighted by the cell cycle, TGFb/Wnt/MAPK/Notch/Hedgehog pathways (19, 23). The urinary bladder is highly enriched with blood vessels, though the function of pericytes in this organ has been sorely under-investigated.

To determine whether the potential pericyte markers identified from the mouse tissue were conserved in humans, we analyzed single-cell data from the Human Lung Cell Atlas and adult human heart single-cell data collected by Litvinukova et al. (24, 25) (Figure 4). The Human Lung Cell Atlas data consisted of scRNA-seq on 75,000 human lung cells and circulating blood collected from three individuals aged 75 (male), 46 (male), and 51 years (female) (24). The human heart data consisted of scRNA-seq on 486,134 human cells from all heart compartments and collected 14 individuals of both sexes with ages ranging from 40 to 75 years (25). To unify our analysis of these data with the Tabula Muris Senis data analysis, we selected only droplet-based scRNA-seq data processed using the 10 × Genomics platform (10 × 3′v2), which left us with three individuals-worth of lung tissue, and seven individuals worth of heart tissue. The data across individuals were then integrated and processed using the Seurat pipeline, resulting in a final dataset consisting of 65,662 cells from the human lung tissue and 238,154 cells from the human heart tissue (26). Pericytes were found annotated based on CSPG4, PDGFRB, and TRPC6 within the Human Lung Cell Atlas, while within the human heart data, pericytes were found annotated based on RGS5, ABCC9, and KCNJ8. We mapped the annotated pericytes onto the UMAP plot. We compared it with the expression density of the human orthologs of the potential pericyte markers identified in mouse lung and heart.

CSPG4 or PDGFRB was expressed in the human lung and heart (Figures 4A,B; Supplementary Figures 7, 8). We identified high expression of GUCY1B1, HIGD1B, NDUFAL2, COX4I2, and KCNK3 in the annotated pericyte cluster from the Human Lung Cell Atlas. Except for GUCY1A1, PDE5A, GUCY1B1, and DES, all other markers were non-specific and expressed in other regions. The expression of HIGD1B, NDUFAL2, COX4I2, and KCNK3 was more specific and almost fully overlapped with CSPG4 or PDGFRB and mostly did not overlap with the expression of vascular smooth muscle cell (VSMC) markers ACTA2, MYH11, and TAGLN (Figure 4A and Supplementary Figure 7). For the human orthologs of the potential pericyte markers identified from the mouse heart tissue, we identified the expression of NOTCH3 and HIGD1B to be relatively high and specific to the annotated pericytes compared with the other markers and relatively did not overlap with the VSMC markers.

Most scRNA-seq approaches rely on dissociating cells from tissues, thereby losing the crucial spatial context of the cells' locations within the tissue. Techniques such as the image-based single-cell transcriptomics method, multiplexed RNA FISH, hybridized tissue sections to spatially barcoded microarrays, sequencing on intact tissues would simultaneously profile the expression of hundreds or thousands of genes within single cells whose spatial location is preserved (27). We, therefore, investigated the pericyte spatial localization within the human heart tissue. We compared it with regions of capillaries and arteries using the spatial transcriptomic data of human heart tissue from the Spatial Gene Expression Dataset by 10 × Genomics Space Ranger (10 × Genomics Space Ranger 1.1.0) (Figure 5).

Each spot in the dataset represented tissue regions where the transcriptomic reads were obtained. The dataset was projected onto a 2-dimensional map using UMAP, and unsupervised clustering of the spots was performed in which we identified 14 clusters (28, 29). In a sequential tissue section stained by hematoxylin and eosin (H&E), blood cells in dark pink color were visible within the hollow cavities where the same areas as vessels were mainly estimated to reside in. Integrating the transcriptomic and spatial data, the tissue section can be annotated as the same 14 clusters applied to the UMAP plot using the same color-coding (Figure 5A).

Co-expression of CSPG4 and PDGFRB was mainly found in Cluster 1/6, from which 1,167 spots were identified and designated as pericyte spots (Figures 5B,C, red dots). A total of 289 spots that expressed MYH11 were identified as smooth muscle coverage of arteries (orange dots), while 1,018 spots that expressed PECAM1 were endothelial cells and capillary lumens (dark pink dots). While capillary spots were randomly distributed, we observed the pattern and shape of arteries as continuous and vessel-like.

We further identified 1,349 spots expressing HIGD1B (dark blue dots) and 1,222 spots expressing NOTCH3 (cyan dots), in which 449 spots co-expressed both HIGD1B and NOTCH3 (Venn diagram, 176+273) (Figure 5C). After colocalization with Cluster 1/6, 37.87% overlapped spots were co-expressing HIGD1B (yellow dots on the upper-side flow), while 37.36% overlapped spots were co-expressing NOTCH3 (yellow dots on the lower-sided flow). A total of 176 dots (6.8%) co-expressed HIGD1B and NOTCH3 among the pericyte spots (sky blue dots on the middle right-sided flow), in concordance with their expected location independent from MYH11 and PECAM-1. Although these spots were not cells (but rather aggregates of transcripts) and would have lower resolution than immunofluorescence staining, our spatial data revealed that pericytes were exclusively located on distal vessels, not on large vessels. This method also opened a new avenue to evaluate quantitative and spatially resolved maps of gene expression at a transcriptomic level no longer dependent on the accessibility of their corresponding antibodies.

In this manuscript, we analyzed microfluidic droplet based scRNA-seq data processed using the 10 × Genomics platform (10 × 3′v2) of the (1) Tabula Muris Senis, (2) Human Lung Cell Atlas, and the (3) Cells of the human heart that were all obtained from the Cellxgene collections (https://cellxgene.cziscience.com) (16, 24, 25). The specific pre-processed Seurat objects for these datasets we obtained are as follows:

In addition, spatial transcriptomics for the human heart was obtained from the Spatial Gene Expression Dataset by 10x Genomics Space Ranger (10x Genomics Space Ranger 1.1.0) at: (https://www.10xgenomics.com/resources/datasets/human-heart-1-standard-1-1-0).

The obtained scRNA-seq data were analyzed based on the pipeline of the Seurat R package (26). The dataset of each subject and tissue was identified and processed through the Seurat pipeline in which the feature counts were log-normalized and subsequently scaled. The top 2,000 variable features were then selected based on the Variance Stabilizing Transformation method for integration (50). Within each tissue, the datasets of the according subjects were then integrated using the IntegrateData function of Seurat. The annotated cell types within the Seurat objects of the human lung dataset and the human heart datasets were used to identify pericytes.

For each integrated scRNA-seq dataset, principal component (PC) analysis was conducted on the highly variable genes to reduce the dimensionality. Based on the percentage of variance explained by the top PCs, we selected the top 45 PCs for the mouse kidney, the top 30 PCs for the mouse bladder, and the top 35 PCs for the remaining (lung and heart of both the mouse and human) datasets to conduct the downstream analysis. Based on the selected PCs, for each dataset, cells were projected into a 2-dimensional space using the Uniform Manifold Approximation and Projection (UMAP) algorithm (28). We conducted an unsupervised clustering of single cells using the FindCluster function in Seurat, setting the resolution measure to 0.5, which constructs a shared nearest neighbor graph by identifying the k-nearest neighbors for each cell in the PC space and weights the edges between the cells based on the number of the k-nearest neighbors or their Jaccard similarity (29). The densely connected cells are then considered as clusters in which the modularity is optimized using the Louvain community detection method (26, 29).

Differential expression analysis was conducted by comparing the expression of stringent pericytes (cells that co-expressed Cspg4 and Pdgfrb within the identified pericyte-enriched cluster) with all other cells according to the tissue for the Tabula Muris Senis dataset using a Wilcoxon rank sum test (26). Genes with an adjusted p-value below 0.05 that are expressed in more than 80% of the stringent pericytes and in <5% of all other cells were selected as potential pericyte markers according to the tissue.

The Gene Set Enrichment Analysis (GSEA) (17) was used to interpret the gene expression data using the software available at https://www.gsea-msigdb.org/gsea/index.jsp. A pre-ranked GSEA determines whether genes within an a priori defined specific gene tend to occur toward the top or the bottom of a pre-ranked gene list, which then identifies whether the corresponding gene set shows statistically significant differences between two biological states. We generated pre-ranked gene lists by sorting the DE genes based on their log2 fold change from the DE analysis in a descending order. GSEA was conducted using the resulting pre-ranked gene lists and compared with the gene sets defined by the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets. A permutation test estimates the statistical significance of the enrichment score. GSEA generates a version of the dataset with phenotype labels randomly scrambled, produces the corresponding ranked list, and then recomputes the enrichment score of the gene set for this permuted dataset. Pathways with an adjusted q-value below 0.05 in this analysis were considered statistically significant.