Single-cell ‘omic profiles of human aortic endothelial cells in vitro and human atherosclerotic lesions ex vivo reveal heterogeneity of endothelial subtype and response to activating perturbations
- The Department of Cellular and Molecular Medicine, The University of Arizona, United States
- The Clinical Translational Sciences Graduate Program, The University of Arizona, United States
- Asthma and Airway Disease Research Center, The University of Arizona, United States
Figures
Figure 1 with 2 supplements
Human aortic endothelial cell (HAEC) transcriptomic profiling discover a heterogeneous cell population.
(A) Schematic diagram of the experimental design. Endothelial cells (ECs) were isolated from six human heart transplant donor’s ascending aortic trimmings and treated with IL1B, TGFB2, or siERG (ERG siRNA) for 7 d. (B) Weighted nearest-neighbor UMAP (WNNUMAP) of aggregate cells from all perturbations and donors is shown. Each dot represents a cell, and the proximity between each cell approximates their similarity of both transcriptional and epigenetic profiles. Colors denote cluster membership. (C) The proportion of cells from each donor for each EC subtype. (D) Gene expression across top markers for each cluster including pan EC (ERG), EC1 (KDR), EC2 (TOP2A), and EC4 (COL1A1). (E) Top markers for pan EC (PECAM1, CDH5, ERG), EC1 (KDR, PGF), EC2 (CENPE, TOP2A), EC3 (SEMA3C, ACKR3), EC4 (COL1A1, COL6A1), and EC5 (LRRC17, LAMA2). The size of the dot represents the percentage of cells within each EC subtype that express the given gene, while the shade of the dot represents the level of average expression (‘Avg. Expn.’ in the legend). (F) Heatmap of pathway enrichment analysis (PEA) results from submitting top 200 differentially expressed genes (DEGs; by ascending p-value) between EC subtypes. Rows (pathways) and columns (EC subtypes) are clustered based on -Log10(P). (G) Violin plots of top Metascape pathway module scores across EC subtypes. Module scores are generated for each cell barcode with the Seurat function AddModuleScore().
Figure 1—figure supplement 1
Marker genes and pathways for endothelial clusters.
(A) Violin plots of representative marker genes for angiogenic EC1 (KDR, GAB1, PGF, NRP1), proliferative EC2 (CENPE, CENPF, KIF11, KIF4A), and mesenchymal EC4 (COL1A1, COL1A2, COL3A1, COL5A1) sub-phenotypes. (B) Top 20 pathway enrichment analysis (PEA) results from submitting top 200 differentially expressed genes (DEGs; by ascending p-value) regulated in EC3 versus EC1-2 and EC4-5. (C) Top 20 Gene Ontology (GO) PEA results from submitting top 200 DEGs (by ascending p-value) regulated in EC3 versus EC1-2 and EC4-5.
Figure 1—figure supplement 2
Violin plot of XIST showing expected expression in female in vitro donor cells (1 and 3) and lack of expression in male in vitro donor cells (2, 4, 5, and 6).
Figure 2
Endothelial cells (ECs) have epigenetically distinct cell states.
(A) Upset plot of differential peaks across EC subtypes. Intersection size represents the number of genes at each intersection, while set size represents the number of genes for each EC subtype. (B) Genomic annotation for the complete peak set. (C) Heatmap of top transcription factors (TFs) from motif enrichment analysis for marker peaks in each EC subtype. Top TFs for each EC subtype are selected based on ascending p-value. Rows (TFs) and columns (EC subtype) are clustered based on enrichment score (ES). (D) Feature plots and position weight matrices (PWMs) for top TF binding motifs for EC1 (TCF12), EC2 (ETV1), EC3 (GATA5), and EC4 (TEAD3). Per-cell motif activity scores are computed with chromVAR, and motif activities per cell are visualized using the Signac function FeaturePlot. (E), PWMs comparing Jaspar 2020 ETV1 motif to ERG motif reported in Hogan et al.
Figure 3 with 3 supplements
Endothelial cell (EC)-activating perturbations modestly shift cells into the EC3 subtype.
(A) The proportion of cells in 7-day control and 7-day IL1B treatment are shown per human aortic endothelial cell (HAEC) donor and cluster on the top and for 7-day control and 7-day TGFB2 on the bottom. (B) The proportion of cells in 7-day siSCR control and 7-day siERG knockdown are shown per HAEC donor and cluster. EC1 was omitted in (A) due to lack of cells in both conditions.
-
Figure 3—source data 1
Annotated western blots for Figure 3—figure supplement 3B where the leftmost six wells are shown from left to right as (1) the protein ladder (labeled in kD), (2) the lipofectamine transfected control, (3) the scrambled siRNA control, (4, 5) two lanes using different siRNAs against the TCF4 gene (not relevant to these studies), and lastly, (6) siRNA against ERG.
The ERG blot (middle) and the H3 blots (bottom) originated from the same gel and membrane, which was cut at about 30 kD, and each piece was blotted for either anti-ERG antibody or anti-H3 antibody.
- https://cdn.elifesciences.org/articles/91729/elife-91729-fig3-data1-v1.pdf
-
Figure 3—source data 2
Original western blots for ERG, annotated and uncropped.
Lanes are the same as described in Figure 3—source data 1.
- https://cdn.elifesciences.org/articles/91729/elife-91729-fig3-data2-v1.pdf
-
Figure 3—source data 3
Original western blots for H3, annotated and uncropped.
Lanes are the same as described in Figure 3—source data 1.
- https://cdn.elifesciences.org/articles/91729/elife-91729-fig3-data3-v1.pdf
-
Figure 3—source data 4
Original western blots for ERG, unannotated and uncropped.
Lanes are the same as described in Figure 3—source data 1.
- https://cdn.elifesciences.org/articles/91729/elife-91729-fig3-data4-v1.pdf
-
Figure 3—source data 5
Original western blots for H3, unannotated and uncropped.
Lanes are the same as described in Figure 3—source data 1.
- https://cdn.elifesciences.org/articles/91729/elife-91729-fig3-data5-v1.pdf
Figure 3—figure supplement 1
Principal component analysis using RNA-seq data.
(A) Principal component analysis (PCA) of EC1-4 snRNA-seq samples +/-siERG or control across donor replicates. (B) PCA of EC1-4 snRNA-seq samples +/-TGFB2 or control across donor replicates. (C) PCA of EC1-4 snRNA-seq samples +/-IL1B or control across donor replicates.
Figure 3—figure supplement 2
Receptor expression profiles across endothelial clusters.
(A) Feature plots of expression of TGFB pathway receptors: TGFBR1, TGFBR2, and ACVRL1. (B) Feature plots of IL1B pathway receptors: IL1R1, IL1R2, IL1RL1, and IL1RL2.
Figure 3—figure supplement 3
Validation of ERG knockdown.
(A) qPCR results for ERG knockdown across donors. (B) Western blot representing a typical knockdown of siERG with the siRNA pools used in this study. TCF4 samples are irrelevant for the purpose of this study.
Figure 4 with 1 supplement
Endothelial cell (EC)-activating perturbations in vitro elicit EC subtype-specific transcriptional responses.
(A) Upset plots of up- and downregulated differentially expressed genes (DEGs) across EC subtypes with siERG (gray), IL1B (pink), and TGFB2 (blue). Upset plots visualize intersections between sets in a matrix, where the columns of the matrix correspond to the sets, and the rows correspond to the intersections. Intersection size represents the number of genes at each intersection. (B) Pathway enrichment analysis (PEA) for EC3-4 up- and downregulated DEGs with TGFB2 compared to control media. (C) PEA for EC2-4 up- and downregulated DEGs with IL1B compared to control media. (D) PEA for EC1-4 up- and downregulated DEGs with siERG compared to siSCR. (E) PEA comparing up- and downregulated DEGs that are mutually exclusive and shared between IL1B and siERG in EC3.
Figure 4—figure supplement 1
Profiles of shared peaks and motif enrichments across endothelial subtypes.
(A) Upset plot of induced peaks for siERG (gray), IL1B (pink), and TGFB2 (blue) across EC1, EC2, EC3, and EC4. Upset plots visualize intersections between sets in a matrix, where the columns of the matrix correspond to the sets, and the rows correspond to the intersections. Intersection size represents the number of genes at each intersection. (B) Heatmap of top motifs enriched in IL1B-induced peaks. Top transcription factors (TFs) for each EC subtype are selected based on ascending p-value. Rows (TFs) and columns (EC subtype) are clustered based on enrichment score (ES). (C) Heatmap of top motifs enriched in siERG-induced peaks. (D), Heatmap of top motifs enriched in TGFB2-induced peaks.
Figure 5 with 5 supplements
Endothelial cells (ECs) from ex vivo human atherosclerotic plaques show two major populations.
(A) scRNA-seq UMAP of different cell subtypes across 17 samples of ex vivo human atherosclerotic plaques. (B) Dot plot of top markers for each cell type. (C) Heatmap of pathway enrichment analysis (PEA) results generated from submitting 200 differentially expressed genes (DEGs) between endothelial cells 1 (Endo1) and endothelial cells 2 (Endo2). Rows (pathways) and columns (cell subtypes) are clustered based on -Log10(P). (E) Heatmap displaying expression of genes belonging to ribosome cytoplasmic pathway for Endo1 and Endo2.
Figure 5—figure supplement 1
Characterization of RNA-seq profiles from ex vivo arterial samples.
(A) UMAP displaying original clusters formed from scRNA-seq data taken from 17 samples across four studies of human ex vivo atherosclerotic plaques. Colors denote different clusters. (B) UMAP from (A). Colors denote anatomical location from which cells derived. (C) Stacked bar graph showing the distribution of anatomic location (red denoting carotid, blue denoting coronary arteries) from which cells derived. (D) Heatmap of pathway enrichment analysis (PEA) results from submitting top 100 differentially expressed genes (DEGs) (by ascending p-value) between ex vivo cell types. Rows (pathways) and columns (cell subtypes) are clustered based on -Log10(P). (E) PEA of the top 100 DEGs (by ascending p-value) for plasmacytoid dendritic cells (PDCs). (F) PEA of the top 100 DEGs (by ascending p-value) for neurons. (G) PEA of the top 100 DEGs (by ascending p-value) for basophils. Adjusted p-value<0.05 for DEGs submitted in (D–G).
Figure 5—figure supplement 2
Violin plots displaying upregulation of several EndMT markers in Endo2, compared to Endo1, including FN1, BGN, COL8A1, ELN, CCN1, FBLN5.
Figure 5—figure supplement 3
Cross comparisons between in vitro and ex vivo RNA-based module scores.
(A) Feature plots for each ex vivo module score across in vitro cells. Briefly, ex vivo module scores are generated using top marker genes for each ex vivo cell subtype. The Seurat function AddModuleScore is used to score each cell for visualization. (B) Violin plots displaying endothelial cells 1 (Endo1), endothelial cells 2 (Endo2), and vascular smooth muscle cell 5 (VSMC5) module scores for each perturbation across in vitro EC1–4. (C) Feature plots displaying distribution of Endo1 (red) versus Endo2 (green) module scores across in vitro cells. (D) Feature plots displaying distribution of Endo1 (red) versus VSMC5 (green) module scores across in vitro cells. (E) Feature plots displaying distribution of Endo2 (red) versus VSMC5 (green) modules scores across in vitro cells.
Figure 5—figure supplement 4
Heatmap and pathway analysis for marker genes of the VSMC5 ex vivo cluster.
(A) Heatmap displaying average expression of VSMC5 maker genes (black arrow) across in vitro and ex vivo datasets. Rows (genes) and columns (cell subtypes) are clustered based on average expression for each given gene. (B) Pathway enrichment analysis (PEA) of the top 200 genes for VSMC5 (adjusted p-value<0.05).
Figure 5—figure supplement 5
Breakdown of ex vivo module scores across in vitro clusters and sample identity.
(A) Violin plots of ex vivo Endo1 module scores across EC1–4. (B), Violin plots of ex vivo Endo2 module scores across EC1–4. (C) Violin plots of ex vivo VSMC5 module scores across EC1–4. Adjusted p-value for (A–C) generated using Wilcoxon rank sum test with continuity correction by setting the alternative hypothesis to ‘two.sided’.
Figure 6 with 1 supplement
Endothelial cell (EC) subtype is a major determinant in the ability to recapitulate ‘omic profiles seen in atherosclerosis.
(A) Heatmap displaying average expression between in vitro perturbation-subtype combinations and ex vivo cell subtypes using all up- and downregulated genes between IL1B, TGFB2, or siERG versus respective controls. Spearman correlation was used as the distance metric. Rows (in vitro EC subtypes) and columns (ex vivo cell subtypes) are clustered using all significant genes (adjusted p-value<0.05) induced and attenuated across all in vitro EC subtypes for each perturbation versus its respective control. (B) Heatmap of coronary artery disease (CAD)-associated single-nucleotide polymorphism (SNP) enrichments across in vitro EC subtypes and perturbation–subtype combinations. Rows (EC subtypes and perturbation-subtype combinations) are clustered using -Log10(P) for enrichment in significant CAD-associated SNPs (p-value<5 × 10–8). Note that ‘diff’ represents peaks common to more than one EC subtype; it is found by subtracting EC1–5 subtype-specific peaks from the entire in vitro peak set (termed ‘panEC’). (C) Coverage plots displaying links for COL4A1/COL4A2 genes to EC4-specific peaks, including one overlapping with CAD-associated SNP rs9515203. (D) Coverage plot showing links for PECAM1 gene to EC4-specific peaks, including one overlapping with CAD-associated SNP rs1108591. (E) Coverage plot showing links for BMP6 gene to EC4-specific peaks, including one overlapping with CAD-associated SNP rs6597292.
Figure 6—figure supplement 1
Pathway enrichment analysis (PEA) of significant (p-value<0.05) EC4-linked genes that overlap with significant (p-value<10–8) coronary artery disease (CAD)-associated single-nucleotide polymorphism (SNP).
Additional files
-
Supplementary file 1
(a) The numbers of total reads and reads per nucleus in snRNA-seq and snATAC-seq data: total reads and reads per nucleus in snRNA-seq and snATAC-seq data are shown. (b) Cells and transcripts detected per nucleus for each snRNA-seq library after filtering during quality control. (c) Proportion of cells belonging to each donor for each EC cluster. (d) Markers for each cluster: discovered using Seurat’s FindAllMarkers function with latent.vars set to ‘Treatment’. Each gene is tested for differential expression between cells of each cluster and the cells outside of that cluster. p_val corresponds to the p-value of the Wilcoxon rank sum test, avg_log2FC is the log2-fold change in expression of the marker for the average cell in the cluster. Subtype is the designated subtype chosen for each cluster, cluster is the cluster of cells for which the marker was discovered, gene is the markers, pct.1 is the percentage of cells expressing the marker within the tested cluster, pct.2 is the percentage of cells expressing the marker outside the cluster, and p_val_adj is the Bonferroni adjusted p-value (based on the number of genes tested). (e) Top biological pathways from Metascape pathway enrichment analysis (and corresponding p-values) taken from top 100 significant genes for each cluster from FindAllMarkers. (f) Cells and peaks detected per nucleus for each snATAC-seq library after filtering during quality control. (g) List of chemical and genetic perturbations that each primary cell donor underwent: TGFB2=7-day TGFB2 (10 ng/ml); IL1B=7-day IL1B (10 ng/ml); IL1B6h=6 hr IL1B (1 ng/ml); siERG = siRNA-mediated knock-down of ERG (siERG), achieved with serial transfections on days 0 and 4. Ancestries were determined according to 1000 Genomes Project superpopulations. (h) Differential accessible regions (p<0.005) across EC subtypes. (i) Differentially expressed genes (DEGs) with TGFB2, IL1B, or siERG induction across EC subtypes: briefly, pseudobulk DEG analysis is performed to find genes that are significantly differentially expressed between perturbations and respective controls for pan ECs (EC1-4) and each EC subtype. DEGs with adjusted p-value<0.05 are considered significant and listed below. (j) Differentially accessible regions affected by TGFB2, IL1B, or siERG across EC subtypes: briefly, significant DARs were found using FindMarkers using logistic regression with ‘nCount_ATAC’ set as a latent variable. DARs with p<0.005 were considered significant and listed below. (k) Description of public data used in this study. (l) Markers for each cluster: discovered using Seurat’s FindAllMarkers function; briefly, each gene is tested for differential expression between cells of each cluster and the cells outside of that cluster. p_val corresponds to the p-value of the Wilcoxon rank sum test, avg_log2FC is the log2-fold change in expression of the marker for the average cell in the cell type, cell type is the collapsed cluster of cells for which the marker was discovered, gene is the markers, pct.1 is the percentage of cells expressing the marker within the tested cluster, pct.2 is the percentage of cells expressing the marker outside the cluster, and p_val_adj is the Bonferroni adjusted p-value (based on the number of genes tested). (m) Effect of EndMT perturbation on in vitro EC subtypes according to ex vivo cell type module scores: Adjusted p-values<0.05 are generated using Wilcoxon rank sum test with continuity correction setting alternative hypothesis to ‘two.sided’. Colored arrows represent significantly (adjusted p-value<0.05) upregulated (green) and downregulated (red) module scores for each EC sub-phenotype and perturbation combination. (n) Results table for effect of EndMT perturbation on in vitro EC subtypes according to ex vivo cell type module scores: Adjusted p-values<0.05 are generated using Wilcoxon rank sum test with continuity correction setting alternative hypothesis to ‘two.sided’. (o) SNP enrichment analysis results across clusters and perturbation-cluster combinations: briefly, SNP enrichment analysis was used with traseR with 'test.method' set to ‘fisher’ and 'alternative' set to ‘greater’. (p) Significantly enriched SNPs (adjusted p-value<0.05) across clusters and perturbation-cluster combinations: briefly, SNP enrichment analysis was used with traseR with 'test.method' set to ‘fisher’ and 'alternative' set to ‘greater’. (q) Significant links between genes and EC4-specific peaks which overlap with CAD-associated SNPs (p-value<0.05): Briefly, peak-to-gene linkage was performed using Signac ‘LinkPeaks’ function on a dataset consisting of EC4-specific peaks and genes. Results were filtered based on peak regions which overlap with significant (p<5e-8) CAD associated SNPs. (r) siRNAs and qPCR primers used in this study. siERG #1, #2, #4, and #5 were pooled together. Non-targeting siRNA (siSCR) #3 and #4 were pooled together.
- https://cdn.elifesciences.org/articles/91729/elife-91729-supp1-v1.xlsx
-
MDAR checklist
- https://cdn.elifesciences.org/articles/91729/elife-91729-mdarchecklist1-v1.pdf
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Single-cell ‘omic profiles of human aortic endothelial cells in vitro and human atherosclerotic lesions ex vivo reveal heterogeneity of endothelial subtype and response to activating perturbations
eLife 12:RP91729.
https://doi.org/10.7554/eLife.91729.3
