Systems biology analysis of longitudinal functional response of endothelial cells to shear stress

Abstract

Blood flow and vascular shear stress patterns play a significant role in inducing and modulating physiological responses of endothelial cells (ECs). Pulsatile shear (PS) is associated with an atheroprotective endothelial phenotype, while oscillatory shear (OS) is associated with an atheroprone endothelial phenotype. Although mechanisms of endothelial shear response have been extensively studied, most studies focus on characterization of single molecular pathways, mainly at fixed time points after stress application. Here, we carried out a longitudinal time-series study to measure the transcriptome after the application of PS and OS. We performed systems analyses of transcriptional data of cultured human vascular ECs to elucidate the dynamics of endothelial responses in several functional pathways such as cell cycle, oxidative stress, and inflammation. By combining the temporal data on differentially expressed transcription factors and their targets with existing knowledge on relevant functional pathways, we infer the causal relationships between disparate endothelial functions through common transcriptional regulation mechanisms. Our study presents a comprehensive temporally longitudinal experimental study and mechanistic model of shear stress response. By comparing the relative endothelial expressions of genes between OS and PS, we provide insights and an integrated perspective into EC function in response to differential shear. This study has significant implications for the pathogenesis of vascular diseases.

Keywords: RNA-seq; endothelial cells; shear stress; systems biology; time-series analysis.

Fig. S1.

Number of DE genes. (A) Number of DE genes using the ST condition (hour 1) as control. (Upper) The number of up-regulated genes. (Lower) The number of down-regulated genes. In each panel, at each time point, the left and right bars show the number of DE genes for OS vs. ST and for PS vs. ST comparisons, respectively. (B) Number of DE genes in OS vs. PS. In each panel, at each time point, the left and right bars show the number of DE genes that are down-regulated and up-regulated in OS vs. PS, respectively.

Fig. 1.

OS vs. PS log fold-change (LFC) data projected onto the custom pathway consisting of genes and mechanisms exhibiting differential response between OS vs. PS. The heat maps below the gene nodes show the time course of transcriptional changes, representing, from left to right, hours 1, 2, 3, 4, 6, 9, 12, 16, 20, and 24. (A) G1-to-S transition pathway. G1-to-S transition is dependent on the E2F1 activation through RB phosphorylation, which is facilitated in part by CDK2 (bound to cyclin E) and CDK4 (bound to cyclin D). The CDK proteins must be activated by interacting with the CDK-activating kinase. CDKs can also be inhibited by proteins such as CDKN2D. (B) Reconstructed pathway of oxidative stress and superoxide metabolism. Oxygen in the cell is converted to superoxide, hydrogen peroxide, and finally water. Reactive oxygen species are also known to activate HIF1A, a marker of hypoxia. The metallothionein MT1X, which is thought to be hypoxia-responsive, is also shown. (C) AP-1 family of genes and the antioxidant product NQO1. (D) NF-κB and NF-κB target genes.

Fig. S2.

(A) Clustered heatmap of the log₁₀ ratio of the P values in OS vs. ST and in PS vs. ST. Positive values indicate that PS has a smaller P value than OS; thus the pathway is more significantly enriched in PS. Negative values indicate the opposite. (B) When the heatmap is divided into two clusters (via hierarchical clustering), the smaller cluster shows a distinct difference in enrichment of several cell-cycle pathways starting around hour 6. This is particularly pronounced with many pathways relevant to S phase.

Fig. S3.

Endo-MT and related pathways. (A) Raw P values of custom-made Endo-MT pathways through GSEA. See the main text for the composition of these pathways. (B) TGF-β signaling pathway based on pathways in KEGG and Wikipathways (82).

Fig. 2.

(A) The largest contiguous portion of the TF-to-target network for TFs that are down-regulated in OS vs. PS. (B) The largest contiguous portion of the TF-to-target network for TFs that are up-regulated in OS vs. PS. Gene targets are shown as unlabeled nodes. TFs were chosen based on cluster analysis of expression data. Gene targets were chosen based on the presence of an entry in TRANSFAC, as well as whether the target was differentially expressed in at least one time point.

Fig. S4.

(A) Higher-resolution version of Fig. 2A with gene target names visible. Largest contiguous portion of TF-to-target network for TFs that are down-regulated in OS vs. PS. (B) Higher-resolution version of Fig. 2B with gene target names visible. Largest contiguous portion of TF-to-target network for TFs that are up-regulated in OS vs. PS. (C) Longest paths in TF-to-gene networks. Pathways were found by observing the longest contiguous TF-to-TF pathways in the networks in A and B. Because RXRA, RARG, and the RXRA:RARG complex all independently appear in the network in A, and all target RARB, it was inferred that an RARG:RXRA complex is a part of this regulatory pathway. RARG:RXRA inhibition of EGR1 was inferred from manual literature curation (50).

Fig. S5.

(A) Dendrogram of 500 hierarchically clustered TFs that were selected based on their level of expression. The dendrogram was split into three clusters. Displayed are the mean log fold changes per cluster, which is calculated from the log fold change for each TF from each time point. (B) Cluster analysis of the dendrogram. The selection of three clusters was determined by the observed change in the magnitude of the height cutoff between three and four clusters. (C and D) Dendrogram of the largest of the three clusters with mean log fold change 0.047 (C) along with further cluster analysis of this cluster (D). The selection of five clusters was determined by the observed stability of mean log fold change from five clusters to six, in addition to the drop in difference between cutoffs.

Fig. S6.

Proposed mechanosensitive mechanism of regulation of endothelial morphology. DHH, an integral component of the hedgehog signaling pathway, is up-regulated in PS vs. ST while being down-regulated in both OS vs. ST and OS vs. PS. DHH signaling activates the Gli family of TFs. Among the Gli targets is STMN3, a stathmin protein that is involved in microtubule formation and function. The expression profile of STMN3 exhibits strong similarity with the expression profile of DHH, indicating a potential signal transduction from shear to transcriptional regulation.

Fig. 3.

Relative differences in gene expression in ECs exposed to OS or PS over time. Pathways and genes in red are up-regulated in OS vs. PS, and those in blue are down-regulated in OS vs. PS. The TGF-β signaling pathway is labeled “mixed” because its genes’ expression directions were mixed. Divergence in cell-cycle activity between OS and PS begins by hour 6 with E2F1 up-regulation in OS. This divergent cell-cycle activity may occur through KLF2- and KLF4-mediated repression of E2F1 expression in PS. E2F1 may up-regulate v-ATPase genes under OS, as observed in the intermediate hours postshear. This promotes lysosomal trafficking to the cell periphery, thus activating mTOR and inhibiting autophagy by hour 20. The activation of mTOR by E2F1 may also contribute to divergence in cell-cycle activity through S6K activation (20). Through CEBPB, KLF4 activates JUNB at hour 4, which may contribute to KLF4-induced inhibition of cell-cycle activity and to antioxidative stress activity in PS. Genes that contribute to ROS production are up-regulated in OS vs. PS by hour 2, and other oxidative stress-related genes exhibit changes up to hour 24. This increased ROS production activates NF-κB and promotes several early-response genes pertaining to inflammation such as MCP-1 and VCAM-1. Inflammation-related genes show changes during hours 2–9. KLF2 transcriptionally activates RARG, which forms a heterodimer with RXRA and may repress EGR1 activity. EGR1, along with ROS production, transcriptionally activates HIF1A, which are all up-regulated in OS vs. PS. HIF1A is observed to be up-regulated beginning in hour 4 and may contribute to angiogenesis and to Endo-MT. Endo-MT may occur in OS beginning in hour 12. Oxidative stress and autophagy repression, both of which occur in OS, also contribute to Endo-MT.