Uncovering emergent phenotypes in endothelial cells by clustering of surrogates of cardiovascular risk factors

Abstract

Endothelial dysfunction (ED) is a hallmark of atherosclerosis and is influenced by well-defined risk factors, including hypoxia, dyslipidemia, inflammation, and oscillatory flow. However, the individual and combined contributions to the molecular underpinnings of ED remain elusive. We used global gene expression in human coronary artery endothelial cells to identify gene pathways and cellular processes in response to chemical hypoxia, oxidized lipids, IL-1β induced inflammation, oscillatory flow, and these combined stimuli. We found that clustering of the surrogate risk factors differed from the sum of the individual insults that gave rise to emergent phenotypes such as cell proliferation. We validated these observations in samples of human coronary artery atherosclerotic plaques analyzed using single-cell RNA sequencing. Our findings suggest a hierarchical interaction between surrogates of CV risk factors and the advent of emergent phenotypes in response to combined stimulation in endothelial cells that may influence ED.

Figure 1

Overall workflow for the identification and validation of emergent phenotypes associated to the combined exposure to surrogates of CV risk factors such as hypoxia (CoCl₂), Oxidized lipids (OxPAPC), inflammation (IL-1β) and OSS on human coronary artery endothelial cells (HCAEC). 1 Experimental design for in vitro exposure of HCAEC to individual or combined stimuli. 2 Global gene expression and gene regulatory network analyses. 3 Orthogonal validation using scRNA-seq data from human coronary atherosclerotics plaques (dataset GSE131778). The cartoons were created using the Mind the Graph platform (www.mindthegraph.com) and data presented in the manuscript to illustrate the overall workflow.

Figure 2

Gene expression profile revealed the hierarchical contribution of the CVD risk factors CoCl₂, OxPAPC, IL-1β and OSS and key TFs relevant to endothelial function. (a) Principal Component Analysis (PCA) of gene expression profile from each sample showed that the PC1 explained 20.4% of the variation while PC2 16.2%. (b) Hierarchical clustering heatmap using Euclidean distance as a measured parameter revealed details of hierarchical contribution by each condition according to their distance to the LSS. (c) Volcano plot showed the differentially expressed genes (DEGs) when comparing the clustered risk factors vs. LSS. The DEGs were considered as adjusted p-value < 0.05 and |Log2foldchange|> 1.3. (d) TF identified among the DEGs and their respective log2 foldchange. Genes in blue are downregulated and in red are upregulated.

Figure 3

Dysregulated DEGs distribution on individual and clustered risk factors followed by GO map of biological process enriched for the 620 DEGs of the clustered risk factors vs. LSS. (a) Venn diagram of DEGs from each condition. The DEGs were considered as adjusted p-value ≤ 0.05 and |Log2foldchange|≥ 1.3; black circle and black arrow regard the unique DEGs from the clustered risk condition; dashed green circle and green arrow regard to the commonly shared genes of the clustered risk factors with each stimulus; dashed black line for unique stimuli-dependent genes, 146 DEGs for OSS, 90 for IL-1β and 37 for OxPAPC. (b) Heatmap of 349 DEGs shared among the DEGs from each stimulus. (c) Heatmap of the similarity matrix of enriched GO terms (biological process) for the clustered risk factors DEGs. Redundancy analysis was applied to identify and keep the most representative term from the redundant terms (cutoff = 0.7 and adjusted p-value ≤ 0.05). (d) Canonical enriched terms in CoCl₂ & IL-1β & OxPAPC & OSS vs. LSS (121 DEGs) and the source of the stimulus individual (represented by colours light blue, yellow and orange), common (grey), more than one stimulus, and if it is emergent, meaning that the gene is only altered if all four stimuli are combined (green). Gene score and enrichment risk score were calculated by − log₁₀(adjusted p-value).

Figure 4

TF-DEGs regulatory network of the combined stimuli condition. (a) Using the TRRUST database (version 2), the transcription factors (TFs) were mapped to their published transcriptional targets. The identified transcriptional network was further filtered based on whether the TF and their targets were a DEG in CoCl₂ & OxPAPC & IL-1β & OSS vs. LSS (adjusted p-value ≤ 0.05 and |Log2foldchange|≥ 1.3). The diamond shape is TF and the square is the DEG target. The edges mean arrow (activation), circle (unknow) and bar (repression). Node size according to the degree. Genes in blue are downregulated, genes in red upregulated. (b) TF-DEGs regulatory network and the condition which they also are differentially expressed. Green is a modulated emergent gene; Orange is an OxPAPC modulated gene; Light blue is an IL-1β modulated gene; Yellow is an OSS modulated gene; Gray is when the gene is modulated in more than one condition.

Figure 5

scRNA-seq of human coronary artery atherosclerotic plaque identified two endothelial cells (EC) clusters which are defined by the clustered risk factors DEGs. (a) UMAP visualization of clustering identified 17 cell populations (n = 10,671 cells). (b) KEGG enrichment analysis of the two EC clusters using the cluster marker differentially expressed genes between a cluster and all remaining cells. DEG were considered for adjusted p-value ≤ 0.05 and |Log2foldchange|≥ 0.25. (c) Heatmap of the 103 clustered risk factors DEGs expression levels from each endothelial cell of the two clusters followed by the condition in which they are differentially expressed. (d) Violin plot of emergent DEGs which are differentially expressed in EC2 vs. EC1. (e) Schematic illustrating the key dysregulated pathway in the clustered risk factors which are recapitulated in human atherosclerotic plaque. (e) Created in the Mind the Graph platform (www.mindthegraph.com).