Single-Cell RNA Sequencing Reveals That Adaptation of Human Aortic Endothelial Cells to Antiproliferative Therapies Is Modulated by Flow-Induced Shear Stress.
Free full text in Europe PMC
Abstract
Background
Endothelial cells (ECs) are capable of quickly responding in a coordinated manner to a wide array of stresses to maintain vascular homeostasis. Loss of EC cellular adaptation may be a potential marker for cardiovascular disease and a predictor of poor response to endovascular pharmacological interventions such as drug-eluting stents. Here, we report single-cell transcriptional profiling of ECs exposed to multiple stimulus classes to evaluate EC adaptation.Methods
Human aortic ECs were costimulated with both pathophysiological flows mimicking shear stress levels found in the human aorta (laminar and turbulent, ranging from 2.5 to 30 dynes/cm2) and clinically relevant antiproliferative drugs, namely paclitaxel and rapamycin. EC state in response to these stimuli was defined using single-cell RNA sequencing.Results
We identified differentially expressed genes and inferred the TF (transcription factor) landscape modulated by flow shear stress using single-cell RNA sequencing. These flow-sensitive markers differentiated previously identified spatially distinct subpopulations of ECs in the murine aorta. Moreover, distinct transcriptional modules defined flow- and drug-responsive EC adaptation singly and in combination. Flow shear stress was the dominant driver of EC state, altering their response to pharmacological therapies.Conclusions
We showed that flow shear stress modulates the cellular capacity of ECs to respond to paclitaxel and rapamycin administration, suggesting that while responding to different flow patterns, ECs experience an impairment in their transcriptional adaptation to other stimuli.Free full text
Single-Cell RNA Sequencing Reveals That Adaptation of Human Aortic Endothelial Cells to Antiproliferative Therapies Is Modulated by Flow-Induced Shear Stress
Abstract
BACKGROUND:
Endothelial cells (ECs) are capable of quickly responding in a coordinated manner to a wide array of stresses to maintain vascular homeostasis. Loss of EC cellular adaptation may be a potential marker for cardiovascular disease and a predictor of poor response to endovascular pharmacological interventions such as drug-eluting stents. Here, we report single-cell transcriptional profiling of ECs exposed to multiple stimulus classes to evaluate EC adaptation.
METHODS:
Human aortic ECs were costimulated with both pathophysiological flows mimicking shear stress levels found in the human aorta (laminar and turbulent, ranging from 2.5 to 30 dynes/cm2) and clinically relevant antiproliferative drugs, namely paclitaxel and rapamycin. EC state in response to these stimuli was defined using single-cell RNA sequencing.
RESULTS:
We identified differentially expressed genes and inferred the TF (transcription factor) landscape modulated by flow shear stress using single-cell RNA sequencing. These flow-sensitive markers differentiated previously identified spatially distinct subpopulations of ECs in the murine aorta. Moreover, distinct transcriptional modules defined flow- and drug-responsive EC adaptation singly and in combination. Flow shear stress was the dominant driver of EC state, altering their response to pharmacological therapies.
CONCLUSIONS:
We showed that flow shear stress modulates the cellular capacity of ECs to respond to paclitaxel and rapamycin administration, suggesting that while responding to different flow patterns, ECs experience an impairment in their transcriptional adaptation to other stimuli.
Maintaining vascular homeostasis is fundamental to the regulation of blood flow, hemostasis, solute transport, and immunologic response. Endothelial cells (ECs) are the cornerstone of vascular regulation, comprising the endothelium, the inner layer of blood vessels in contact with the bloodstream. ECs are the main regulators of vascular homeostasis, by (1) sensing flow fluctuations and modifying vascular tone, (2) regulating hemostasis, (3) acting as permeability gatekeepers, and (4) enabling the immune response.1
See accompanying editorial on page 2282
ECs are capable of quickly responding in a coordinated manner to the entire array of stresses present in the vasculature. Biochemical (vascular signaling mediators such as TNF-α [tumor necrosis factor alpha], IL [interleukin]-1β, or Ox-LDL [oxidized low-density lipoprotein]), biomechanical (hemodynamic forces or mechanical stretches such as shear stress, hemodynamic pressure, or vasodilatation/vasoconstriction), and cell-derived environmental stimuli or alterations induced during endovascular treatments (eg, endothelial rupture, drug stimulation) are sensed by ECs and drive ECs to commit to specific reactive programs.2 Because ECs perform so many complementary functions, each potential stimulus may lead to a wide array of overlapping fluctuations in the transcriptome and secretome.3
EC dysfunction in response to pathology disturbs the nonmigrating, nonproliferating quiescent state; thus, EC state can be used as a mirror of vascular integrity and overall health. Studying EC states in different pathophysiological contexts is crucial to (1) add to our understanding of how and why vascular response to pharmaceuticals differs across target lesions, organs, and people and (2) define novel markers of vascular disease, the leading cause of death globally.4 Here, we hypothesize that loss of EC cellular adaptation may be an early indicator of vascular disease and a predictor of poor response to clinical intervention.
Addressing this hypothesis first requires identifying paradigmatic stimuli from distinct classes that are known to induce state changes in ECs. We chose 2 orthogonal classes of stress: mechanically mediated (flow-induced shear stress) and chemically medicated (antiproliferative drugs released from indwelling endovascular implants such as drug-eluting stents). Shear stress is a well-known mechanical stimulus driving EC state.5–7 Variations in shear stress levels induce broad phenotypic switching in ECs.2,6 Alterations in shear exposure are implicated in almost every biological function of ECs, including gene expression,8 endoplasmic reticulum stress,9 mitochondrial function,10 and epigenetic regulation.11 Additionally, ECs are the first cells exposed to pharmacological agents, and as such, antiproliferative drugs like rapamycin and paclitaxel are also known drivers of the EC state. Rapamycin and its analogs are compounds with complex biology,12 most recently assuming prominence in vascular biology by virtue of their use in drug-eluting stents. Here, locally delivered rapamycin influences EC proliferation and migration with uncertain cytotoxicity.13 Rapamycin and related analogs impact the EC state via action on mTOR (mammalian target of rapamycin) and by altering EC calcium dynamics.14 Paclitaxel—a member of the taxane class of drugs—is a compound used to avoid restenosis because it is a potent inhibitor of EC proliferation and angiogenesis by modulating the activity of α/β-tubulin heterodimers (microtubules) that conforms the dynamic cell cytoskeleton.15
Shear stress and rapamycin/paclitaxel are not necessarily truly independent stimuli. For example, rapamycin also influences ECs by direct activation of vasoactive receptor families16 including eNOS and by impact on key endothelial TFs (transcription factors) including KLF217; the eNOS and KLF2 pathways are themselves shear responsive. Activation of mTOR has also been implicated in the atheroprone effect of altered shear stress via its effect on autophagy.18 Thus, shear stress and rapamycin/paclitaxel represent distinct stimuli with broad effects on disease-relevant EC states. Understanding these effects in combination at single-cell resolution will elucidate how complex microenvironments may modulate EC dynamism.
Previously, single-cell transcriptomics has been used to describe EC state in the context of biological processes such as angiogenesis19,20 and environmentally driven processes such as diet alteration21,22 and disturbed flow.23 Some authors have also explored the role of shear stress in EC with single-cell RNA sequencing (scRNA-seq). In this study, we used Seq-Well24—a picowell-based scRNA-seq platform—to build upon previously explored markers of shear stress in ECs25 and define the components of EC cellular adaptation to a combination of pathophysiological shear stresses and chemical stimuli under controlled conditions. This allows us to understand how EC cellular adaptation differs when ECs are committed to different flows in a highly controlled, stimulus-specific analysis.
Using our scRNA-seq profile of aortic ECs under well-controlled pathophysiological flows, we detected differentially expressed genes (DEGs) and defined the TF landscape influenced by flow-induced shear stress. These markers were sufficient to differentiate previously defined distinct subpopulations of ECs in scRNA-seq data from murine aorta, exhibiting varying flow sensitivities in spatially distinct domains.
Next, we investigated the influence of dual-stimulus classes on ECs and defined the modular components of the EC response to rapamycin and paclitaxel under varying flow-induced shear stress. We demonstrated that the application of mechanical stimulation via shear stress alters the EC response to antiproliferative agents. This suggests that ECs change responsiveness to stimulation once committed to a stimulus such as different levels of flow-induced shear stress and may serve to explain the diverse drug responses in ECs observed across different target lesions, organs, and people in a clinical setting.
METHODS
Data Availability
Data presented in this study are available at Gene Expression Omnibus, reference number GSE212388 (token: ovapmmqmntcbpwz). All software required to reproduce figures is open source. Code to reproduce presented figures is available on Github. Please see the Major Resources Table in the Supplemental Material for further information.
EC Culture
A combination of primary human aortic ECs (HAECs; PromoCell) isolated from both the ascending and descending aortas at passage P5 was used for all experiments. The media used for cell culture (Lonza EGM-2 or PromoCell Growth Medium 2) and in the flow-controlled bioreactor was supplemented to a total serum concentration of 7% using fetal bovine serum. Media was additionally supplemented with 100 U penicillin/mL, 100 µg streptomycin/mL, and 1 µg amphotericin/mL.
Multimodal Culture of ECs in a Flow-Controlled Bioreactor
A previously described controllable flow loop bioreactor system was utilized for culture of HAECs under laminar high/physiological and low/altered-flow shear stress conditions.26 Briefly, the flow loop was driven by a peristaltic pump (ISMATEC MCP, product No. 78002-00) with an attached 8-roller, 12-channel pump head (ISMATEC Type IBM733A, product No. 78002-36). ISMATEC PharMed BPT pump tubing with an internal diameter of 2.79 mm (Cole Parmer) was utilized in the flow loop. Silastic laboratory tubing with an internal diameter of 0.188 in/4.78 mm and an outer diameter of 0.312 in/7.92 mm (Dow Corning) connected the pump tubing and the cell culture chamber to a custom glass media reservoir, and 3/16″×3/32″ tubing connectors (Cole Parmer) connected the pump tubing and the Silastic tubing. The non–cell chamber components of the flow loop were preassembled and autoclaved before flow loop assembly.
For the laminar high- and low-flow shear stress experiments, HAECs were seeded in Ibidi µ-Slide I Luer 0.8 flow chambers with ibiTreat proprietary tissue culture coating for all experiments. These chambers have a height of 800 µm, a channel volume of 200 µL, a channel length of 50 mm, a channel width of 5 mm, a reservoir volume of 60 µL, and a growth area of 2.5 cm2 per channel. Ibidis were seeded with 200 µL of HAECs at a 0.6 to 1.2 million cells/mL concentration at least 12 hours before connecting to flow; cell attachment and confluence were confirmed by microscopy before initiating shear conditioning. Pump flow rates were set to achieve shear stresses of 30 and 2.5 dynes/cm2 for high and low shear experiments, respectively.
In this study, turbulent flow and numerically quantified wall shear stress on HAECs was also simulated. HAECs were seeded in a custom-made chamber with a narrow channel (1-mm wide, 2-mm high, and 5-mm long), which transitioned into a wider section (3.6-mm wide, 2-mm high, and 15-mm long). An oscillatory volumetric flow with an amplitude of 160 mL/min and a frequency of 1 Hz was applied at the inlet, leading to shear stresses ranging from 0 to 65 dynes/cm2. The simulation of the chamber was conducted using Fluent 2022 R2 (ANSYS, Inc, Canonsburg, PA), utilizing the viscous k-omega model due to its superior performance in resolving internal and separated flows, as well as local recirculation. Chamber’s walls were modeled as stationary, no-slip boundaries. To account for the entrance region (21.7 mm), both tubes at the inlet and outlet were included in the experiment. The computational domain was meshed using ≈4.5 million polyhedral cells after mesh independence for the solution was confirmed. The cell media was modeled as a Newtonian fluid with a density of 997 kg/m3 and a viscosity of 6.91×10−4 kg·m−1·s−1. At the maximum volumetric flow rate (160 mL/min), the Reynolds number was calculated to be 2564 in the narrower channel, justifying our decision to model the flow as turbulent.
To fabricate the chamber, a mold made of stacked layers of polymethyl methacrylate was laser-cut using an Epilog Laser Fusion 120 W, and polydimethylsiloxane (Sylgard) was poured at a concentration of 1:10 and cured overnight at 70 °C. Plasma bonding between a glass coverslip and the polydimethylsiloxane chamber was achieved using a plasma cleaner (Harrick Plasma; PDC-001) with a pressure controller (Harrick Plasma; PDC-FMG). Two holes were punched into both the inlet and outlet, and tubes were connected and sealed with polydimethylsiloxane before being left overnight at 70 °C. Chambers were autoclaved using a 20/20 cycle and then coated with 1.5 µg/cm2 of fibronectin solution in PBS for an hour. After 3× 10-minute washes in PBS, HAECs were seeded, and the experiment proceeded similarly to the laminar flow experiments.
Fifteen milliliters of media were used in each flow reservoir to ensure a physiological pH for the total duration of the experiments (72 hours). The oxygenation of medium was accomplished by using oxygen-permeable flow loop tubing.
Unimodal Flow Stimulus Experiment
HAECs were cultured using the flow loop bioreactor system under constant high flow–induced (30 dynes/cm2), low flow–induced (2.5 dynes/cm2), or turbulent flow–induced shear for 72 hours. Following shear culturing, cells were isolated from flow chambers and used as input to Seq-Well for scRNA-seq.
Bimodal Flow and Chemical Stimulus Experiment
HAECs were cultured using the flow loop bioreactor system under constant high (30 dynes/cm2), low (2.5 dynes/cm2) or turbulent shear stress for a total of 3 days. Twenty-four hours after the initiation of flow, once ECs have been committed to high, low, or turbulent flow, rapamycin (InSolution, Calbiochem/Sigma-Aldrich) or paclitaxel (Invitrogen, ThermoFisher) dissolved in media was added to two-thirds of the flow loops to achieve a final concentration of 1 and 10 µM, respectively. These concentrations have been previously demonstrated to result in alterations in RNA expression in ECs.17,27 An equal volume of normal media was added to the nondrug flow loops. Following shear culturing (24 hours without drugs, 48 hours with drugs: a total of 72 hours), cells were isolated from Ibidis and used as input to Seq-Well for scRNA-seq.
Single-Cell Transcriptional Profiling of Cultured ECs via Seq-Well
Following flow conditioning, cells were harvested from Ibidi flow chambers, and scRNA-seq was performed via Seq-Well, an scRNA-seq method developed to be appropriate for low-input (≈104 cells) samples.28 Seq-Well S324 using second-strand synthesis was utilized for all samples. Per condition, 10 000 cells were loaded as the input to Seq-Well.
After amplification, sample libraries were prepared for sequencing using the Nextera XT DNA tagmentation/amplification. Libraries were purified with SPRI (solid-phase reversible immobilization) ratios of 0.6× and 0.8× applied to the amplified and tagmented product, respectively. Final libraries had an average fragment size distribution between 400 and 600 bp. Libraries were loaded onto NextSeq 500 at a final loading concentration of 2.2 pM with a read structure as follows: 20 cycles covering the 12-bp cell barcode and 8-bp unique molecular identifiers (read 1); 50 cycles on the mRNA strand (read 2); 8 cycles covering sample-specific barcode sequences (index 1). Each sample was allotted 200 million reads for a target sequencing depth of 10000 to 15
000 reads per cell.
Illumina’s bcl2fastq software (2.20.0.422) was used to demultiplex samples according to sample-specific indices introduced in the post-tagmentation polymerase chain reaction. The Drop-seq Pipeline maintained by the Broad Institute was used to align raw FASTQ files to the human reference genome and to generate gene count matrices using the default parameters. Mapped reads filtered for quality were used to generate digital gene expression matrices for each sample. Unique molecular identifier–collapsed data were used to generate the final counts table used for analysis; digital gene expression matrices were combined from separate array samples by retaining the genes detected in all samples to generate count tables for further analysis.
Computational Analysis and Cell Quality Control
Raw count tables were used as the input for analysis in the Seurat package for R (version 4.1.0).29,30 Cells with at least 500 genes and 10 000 molecules detected were retained for analysis; cells with >5\% mitochondrial reads were excluded from analysis. Following quality control, 18 466 cells were analyzed from the in vitro experiments with an average of ≈9194 reads per cell and ≈3017 genes per cell. Normalized expression values were computed for expression in each cell by scaling each cell to 10 000 total unique molecular identifiers and log-normalizing these values with an added pseudocount of 1 using the Seurat NormalizedData function.29,30 Variable genes were identified using the Seurat FindVariableFeatures function with the VST (variance stabilizing transformation) selection method and 2000 as the number of genes to select as top variable features.29,30
High-depth normal mouse aorta transcriptional profiles were downloaded from the Broad Institute Single Cell Portal. ECs were defined based on the expression of the canonical marker Cdh5. A total of 825 ECs were included in the analysis, with an average of 2077 reads per cell and 5900 genes per cell. Data normalization, feature selection, and scaling were performed as previously described above.
Dimensional Reduction of Flow-Conditioned and Dual-Stimulus EC Data
Dimensional reduction of the shear-conditioned EC scRNA-seq profiles was performed on scaled normalized count data using 30 principal components (PCs) based on examination of the PC significance plot. The Seurat implementation of the uniform manifold approximation and projection for dimensional reduction algorithm was used to project 30 PCs into a 2-dimensional space for demonstration of flow state separation.29,30 Dimensional reduction of dual-stimulus data was accomplished similarly also using 30 PCs.
Differential Expression Analysis
Identification of differentially expressed transcripts between conditions (eg, high/physiological versus low/altered shear states or control versus rapamycin/paclitaxel) was performed using the FindMarkers function in Seurat with a minimum log-fold change threshold of 0.25 and with P values computed using a Wilcoxon rank-sum test. Subpopulation-specific markers following cluster analysis were identified using the FindAllMarkers function in Seurat with a minimum log-fold change threshold of 0.25 and with P values computed using a Wilcoxon rank-sum test. Gene set enrichment results were robust to ≈5% variations in choice of log-fold change threshold. P values were adjusted for multiple hypothesis correction in both cases using a Bonferroni correction based on the total number of genes in the analyzed data set; adjusted P<0.05 defined a significant marker. For gene set enrichment analysis (GSEA) plots, differential expression analysis was conducted using R library DESeq2 from Bioconductor (1.32.0).31 To validate the results and contrast several cluster algorithms to evaluate functional heterogeneity in HAECs, we also computed the results with the scGEAToolbox in MATLAB.32
Gene Set Enrichment for Defined Marker Lists
Gene set enrichment on gene lists describing differentially expressed transcripts between 2 conditions or whole-genome correlation network analysis (WGCNA)–derived modules was performed using the R library BiocManager (1.18.0)33 with a P value cutoff of 0.05. Annotated gene sets were drawn from a custom GMT (gene matrix transposed) file derived from the canonical pathways and gene ontology gene sets in v6.2 of the Molecular Signatures Database.34 These include gene sets from the Biocarta, Kyoto Encyclopedia of Genes and Genomes, Reactome, and Gene Ontology databases.
TF Activity Inference
The known targets of each TF in humans were identified using DoRothEA (1.4.2).35 mRNA expression levels were used to compute TF activity, and only TFs with confidence levels A to C were used for this analysis. DoRothEA regulons were coupled with the statistical method VIPER (1.26.0).36 The resulting TF activity matrix was then analyzed similarly to gene analysis as explained above.
WGCNA-Inspired Unsupervised Gene Module Identification
DEGs between rapamycin or paclitaxel and control conditions were identified in both high and low shear contexts as described above using a corrected significance threshold but no fold change threshold. Gene modules within these drug response gene sets that were coexpressed were identified using functions from the R WGCNA library.37 WGCNA was applied with a soft power of 6 used to create the adjacency matrix, transformation of adjacency matrix into a topological overlap matrix, and clustering of the topological overlap matrix using the cutreeDynamic function from WGCNA with method tree using a minimum module size of 10 as described previously.38 Module significance was estimated using a permutation test as described previously.38 Briefly, genes from the total list of genes used for module discovery were binned by expression, and 10 000 random modules with gene expression distributions like that of the tested module were generated. For each random module, a 1-sided Mann-Whitney U test was performed to compare the distribution of dissimilarity values between the tested module and the random module. The P values of these tests are corrected for multiple hypothesis correction with the Benjamini-Hochberg procedure, and the number of tests with false discovery rate >0.05 (indicating that the random module and tested module have similar dissimilarity distributions) was determined. If fewer than 5% of tests (P<0.05) fail to reject the null hypothesis of the Mann-Whitney U test, the tested module was deemed to have significant correlation structure relative to background.
Correlation and Scoring of Unsupervised Gene Modules Defining Antiproliferative Drug Response Across Shear Conditions
To determine whether significant correspondence existed between unsupervised gene modules identified from paclitaxel/rapamycin responses of high/low flow–conditioned cells, a Fisher exact test for significance of gene overlap was implemented. The significance of overlap size between each pair of modules was computed using the set of all differentially altered genes across both shear conditions as the gene universe. P values were adjusted using the Benjamini-Hochberg procedure. A score was computed for each WGCNA-based module for each cell as the total number of postnormalization transcripts (of the 10 000 normalized total) from genes in the module. These scores were scaled to Z scores based on the mean and SD across all cells.
Pseudotime Analysis
Single-cell trajectory analysis was conducted on RNA sequencing data using the Monocle 3 R package (version 3.3.9).39 First, data were filtered based on drug conditions, and a standard preprocessing workflow was applied, including data normalization, selection of highly variable features, data scaling, and linear dimensional reduction, following previously described methods. To identify distinct cellular groups for each drug/flow condition, the Seurat package’s FindNeighbors and FindClusters functions were used at different resolutions, ensuring at least 2 groups were discerned per condition. Uniform manifold approximation and projection was used as the nonlinear dimensional reduction technique for each drug condition.
To define the trajectory path, the learn_graph function from the Monocle package was utilized. Additionally, pseudotime computation and cell ordering were performed using the plot_cells and order_cells functions from Monocle, respectively.
Immunofluorescence Imaging
For immunofluorescence staining, flow was interrupted after 72 hours, and flow chambers were washed with PBS for 5 minutes. Cells were permeabilized with 0.1% Triton (Invitrogen). After 2× 10-minute washes, cells were blocked for 1 hour with 5% goat serum (Sigma-Aldrich) in PBS. The primary antibody was incubated overnight at 4 °C and washed 3× 10 minutes each before adding the secondary antibody for 2 hours at 37 °C. Cell Tracker (Invitrogen) and DAPI (Invitrogen) solutions were added after 3× 10-minute washes and incubated for 30 minutes. Following 3× 10-minute washes, cells were imaged for ki67 (Cell Signaling) expression in rapamycin-treated cells and α-tubulin (Cell Signaling) expression in paclitaxel-treated cells, and the morphology was evaluated based on the expression of DAPI and Cell Tracker with Cell Profiler (4.2.5).
RESULTS
Shear Stress Commits ECs to an Identifiable State
To understand EC dynamism defined by shear forces and understand cellular adaptation of ECs committed to well-differentiated controlled laminar flows, we used Seq-Well S3 to characterize ≈5994 HAECs subjected to either high/physiological (30 dynes/cm2) or low/altered (2.5 dynes/cm2) laminar shear conditions in a controlled bioreactor system (Figure (Figure1A;1A; Figure S1). Following transcriptional profiling, dimensionality reduction was performed using PC analysis, and results were visualized with uniform manifold approximation and projection for dimension reduction, demonstrating that scRNA-seq profiles clearly resolved high and low shear-conditioned EC populations from each other (Figure (Figure1B).1B). The 3 shear replicates from each condition clustered together, demonstrating reproducibility of shear conditioning and transcriptional profiling independent of cell cycle phase (Figure S2A and S2B). DEGs of each shear condition were identified using the Wilcoxon rank-sum test as implemented in Seurat (Figure (Figure1C).1C). DEGs induced by high physiological shear include a number of genes that have previously been associated with high shear stress in HAECs and other human EC populations including PI16, a peptidase inhibitor41; ADAMTS1, a metalloproteinase42; IGFBP5, a mechanosensor in ECs11 (Figure (Figure1D).1D). LYVE1—typically considered a marker of lymphatic ECs—was also a DEG associated with high shear in this HAEC population. DEGs associated with low shear included THBS1, a thrombospondin43; VWF, a glycoprotein associated to hemostasis44; ESM1, EC secreted dermatan sulfate proteoglycan and related to inflammatory response45; and TFPI2, a serine protease inhibitor46 (Figure (Figure11D).
Shear forces determine the commitment of endothelial cells to distinct states. A, Workflow for the stimulation of human aortic endothelial cells (HAECs) using a custom bioreactor system to control shear stress. Single-cell RNA sequencing (scRNA-seq) was used to assess transcriptional changes in HAECs subjected to high (30 dynes/cm2) or low laminar flow–induced shear stress (2.5 dynes/cm2). B, Uniform manifold approximation and projection (UMAP) visualization of 5994 HAECs from all samples (n=3). Cell types are color coded according to their respective shear states. C, Heatmap displaying the top 20 differentially expressed genes (DEGs) for each flow condition. D, Gene expression levels of 4 relevant markers for each flow condition, depicted in the UMAP. E, Visualization of gene set enrichment analysis results for different flow conditions. The color intensity of the bars represents the adjusted P value. F, Average expression of the top 30 transcription factors that discriminate between flow conditions.
We then conducted shear-dependent GSEA using signature gene sets from the Canonical Pathways and Gene Ontology databases.34 High laminar shear stress significantly upregulated (1) sterol biosynthesis, key in modulating biosynthesis of cholesterol and fatty acids47; (2) the Rho pathway, crucial in the regulation of cell morphology, proliferation, and adhesion48; and (3) regulation of actin filament depolarization, enabling turnover of actin filaments and cytoskeletal restructuring49 (Figure (Figure1E).1E). These pathways are correlated with morphological changes observed in HAECs under immunofluorescence microcopy (Figure S3A and S3C). Top upregulated pathways in low laminar flow included extracellular matrix constituents, proteoglycan-binding pathways, and modulators of adherens junctions (Figure (Figure1E).1E). This GSEA confirms that the selected levels of shear stress in this study correspond to the 2 different regions in the set point theory that determines EC remodeling.50
Shear stress–induced TF activity was inferred using DoRothEA regulons35 to partition shear-induced EC state (Figure (Figure1F),1F), allowing modular definition of shear-dependent states. TFs with activity upregulated in high laminar flow were associated with a wide range of EC functions, including, but not limited to, NO processing (KLF451), redox response (NFE2L26), vasculogenic and angiogenic processes (FOSL152 and HNF4α53), and blood pressure control (VDR54). TFs associated with low flow included NFKB2 (nuclear factor kappa B subunit 2), a mediator of inflammation with multiple links to thrombotic processes in ECs55; TP53, a regulator of EC dysfunction in hyperglycemia and ischemia56; and RUNX2, a mediator of EC migration and invasion during tumor angiogenesis.57
Single-cell profiling in combination with TF activity characterization enabled identification of biological associations among shear-responsive TFs in ECs (Figure (Figure1F);1F); for example, the HNF4α-FOSL1 activity was correlated, a relationship previously demonstrated in human colorectal cancer cells.58 A relationship was also identified between retinoic acid and CDX family members (RARA-CDX2), a link previously found in vertebral patterning role through Hox gene expression.59 The RUNX2-TP53 link identified in the study was found to mediate EC migration and invasion during tumor angiogenesis,57 and a NFE2L2-SREBF2 interaction was previously found to regulate lipid metabolism.60
Taken together, the identification of shear-associated DEGs, shear-driven TF modules, and associations among these TF modules via single-cell profiling enabled a multitiered approach to characterization of shear-induced EC state in this model. Most importantly, analysis of DEGs and TF modules validated the physiological relevance of our in vitro model by recapitulating previously identified EC responses to shear stress, for example, the upregulation of KLF family members such as KLF2, KLF4, and KLF6 (Figure (Figure1F;1F; Figure S4D) in response to high physiological flow.
We additionally investigated the effects of turbulent flow in our bioreactor model. Within the vascular tree, atheroprone regions are known to be exposed to turbulent flow patterns accompanied by recirculation. To simulate this phenomenon, an oscillatory volumetric flow was deliberately induced in our turbulent flow model, resulting in the initiation of flow recirculation as demonstrated by computational fluid modeling (Figure S4A).
To comprehensively assess the impact of various flow conditions, we integrated transcriptional profiles of ECs subjected to both laminar and turbulent flow (Figure S4B and S4C). This endeavor revealed a significant overlap between turbulent flow conditions and the previously scrutinized high/physiological and low/altered-flow states. The coexistence of these conditions in turbulent areas of the vascular tree gives rise to EC heterogeneity and varying EC sensitivities to flow (Figure S4E). Moreover, these findings align consistently with the wall shear stress profiles derived from our computational fluid dynamic analysis (Figure S4A).
Notably, our investigation of the top DEGs influenced by shear stress (Figure (Figure1C)1C) further corroborates the convergence of turbulent flow conditions with the high and low shear stress phenotypes, with an overall average response lying between the high and low shear stress conditions (Figure S4D and S4E). Consequently, these genes hold promise as potential indicators to evaluate ECs’ responsiveness to flow (Figure S4E), thereby validating their relevance in this context. Also, we demonstrated that the high and low laminar shear stress levels previously selected provide a distinct response to evaluate flow-sensitive transcriptional change in ECs and that turbulent flows with recirculation present a combination of both states.
Flow-Sensitive DEGs Discriminate Spatially Distributed EC Populations in Murine Aorta
Validation of DEGs and inferred TF profiles identified in this in vitro model was accomplished using previously described scRNA-seq profiles of ECs isolated from the murine aorta21 (Figure (Figure2A),2A), which identified 3 spatially and functionally distinct EC subpopulations within the aortic vascular tree (Figure (Figure2B).2B). Subpopulations were previously differentiated by markers including VCAM1 and CD36 (Figure (Figure2D)2D) and thought to differ in aortic location; the VCAM1-positive EC1 group was localized to the aortic lesser curvature and the CD36-positive EC2 group to the greater curvature and descending aorta.
Shear-dependent genes reveal spatially distinct subpopulations of endothelial cells (ECs) in the murine aorta with distinct flow sensitivities. A, Experimental methodology used by Kalluri et al21 to isolate cell populations from the mouse aorta. Fluorescence microscopy demonstrated high VCAM1 and low CD36 fluorescence expression in the aortic root EC populations, while low VCAM1 and high CD36 fluorescence expression was observed in the descending aorta EC populations. B, Uniform manifold approximation and projection (UMAP) visualization identifying the 3 main EC clusters present in the mouse aorta. C, UMAP colored according to the expression of relevant shear-dependent differentially expressed genes (DEGs) identified in vitro using human aortic ECs (HAECs). D, Transcriptional expression levels of VCAM1 and CD36 in each EC subpopulation. E, Heatmap demonstrating how relevant shear-dependent DEGs identified in vitro using HAECs can characterize EC subpopulations in the mouse aorta. F, Correlation between shear-dependent TF (transcription factor) response in the mouse aorta and the TF response observed in HAECs.
The top DEGs influenced by flow shear stress in the in vitro model (Figure (Figure1C)1C) were sufficient to discriminate these 3 in vivo subpopulations previously identified using the entire transcriptome (Figure (Figure2B2B and and2E).2E). EC1 and EC3 subpopulations exhibited higher overall sensitivity to flow, with an upregulation of shear stress–sensitive genes compared with EC2 (Figure (Figure2E).2E). Moreover, EC1 and EC3 demonstrated less expressions of high shear–related genes, such as RGS5 (regulation of G protein signaling 5), NDRG1, SLC9A3R2, and KLF4, but greater expression of low shear–related genes, such as THBS1, VVF, BGN, and MGP (Figure (Figure2C).2C). This finding aligns with the idea that the low shear stress level defined in our in vitro experiment correlates with regions of recirculation or aberrant flow, where atheroprone regions are found within the vascular tree. Furthermore, analysis of the inferred TF modules defined using the flow loop bioreactor in the in vivo data revealed an upregulation of low/altered shear–related TFs in the EC1 subpopulation, such as NFKB2, ARNT, MYC, or TEAD4, while EC2 predominantly overexpressed the high/physiological shear-related TFs identified in our in vitro data, such as KLF4, STAT5, STAT6, or WT1 (Figures (Figures1F1F and and22F).
The inferred TF module analysis on murine data revealed that some previously defined in vitro biological associations correlated with in vivo findings (Figures (Figures1F1F and and2F).2F). For instance, the NFE2L2-SREBF2 interaction may be involved in the regulation of oxidative stress and lipid metabolism, both of which were affected by the response of ECs to flow shear stress in our GSEA (Figure (Figure1E).1E). Another relationship found in both data sets was WT1 and STAT6, which may play a role in EC development, inflammation, or immune response and share common or parallel pathways in the regulation of NO production and inflammation mediated by the mechanotransduction response of ECs to shear stress. STAT5B, ESRRA, SP3, and NFE2L2 were identified as low shear–related TFs in both in vitro and in vivo analyses, and they may share signaling pathways with the response of ECs to shear stress through their role in the regulation of the immune system, energy metabolism, and oxidative stress response. Finally, MYC, TEAD4, and ARNT were associated with high/physiological shear stress–related TFs, and their roles in regulating EC proliferation, development, and environmental signaling may be linked to the EC signaling pathways activated by ECs in response to physiological shear stress to maintain vascular homeostasis.
Antiproliferative Therapies Elicit Unique EC Fingerprints When Costimulated With Distinct Flow Conditions
To build on our results in EC populations responding to a single stimulus, we used scRNA-seq to profile ≈6252 HAECs cultured in a bimodal stress experiment under high/physiological or low/altered laminar flow (as defined previously) with or without administration of antiproliferative therapies (Figure S1). This experiment resembles endovascular treatment scenarios where ECs are first committed to high or low levels of shear stress at the endothelium wall by local anatomy and pathology. In the course of clinical treatment, antiproliferative drugs are released by diffusion from a drug-eluting stent/balloon. Given the broad known effects of rapamycin and paclitaxel on EC proliferation and migration,13,61 calcium dynamics,14,62 and autophagy,18,63 as well as impact on flow-responsive elements like KLF217 (rapamycin), their effects in this system were defined both independently of and in combination with shear alteration. For shear-independent analysis, drug effects were evaluated irrespective of shear stress level. Drug concentrations selected in this work have been previously demonstrated to result in alterations in RNA expression in ECs.17,27 Dimensional reduction with PC analysis and visualization with uniform manifold approximation and projection demonstrated clear resolution of cell populations by shear condition and clear resolution of the rapamycin/paclitaxel-treated population under high shear stress conditions but not in low shear stress conditions for rapamycin-treated cells (Figure (Figure3A3A and and33D).
Antiproliferative therapies induce a distinctive signature in endothelial cells (ECs) when costimulated with different laminar flows. A, Uniform manifold approximation and projection (UMAP) visualization of 1625 human aortic ECs (HAECs) stimulated with shear forces and paclitaxel. Cell types are labeled with colors corresponding to the shear flow and drug states. B, Volcano plot illustrating the top differentially expressed genes (DEGs) upregulated or downregulated by paclitaxel when ECs are exposed to any flow condition. C, Dot plot summarizing the top 10 genes upregulated or downregulated by the drug condition (paclitaxel). D, UMAP visualization of 4627 HAECs stimulated with shear forces and rapamycin. E, Volcano plot illustrating the top DEGs upregulated or downregulated by rapamycin when ECs are exposed to any flow condition. F, Dot plot summarizing the top 10 genes upregulated or downregulated by the drug condition (rapamycin). G, Visualization of the gene set enrichment analysis for paclitaxel independent of flow condition. The color intensity of the bars represents the adjusted P value obtained from gene set enrichment analysis (GSEA). H, Visualization of the gene set enrichment analysis for rapamycin independent of flow condition.
First, the response elicited by antiproliferative therapies was analyzed independently of shear-induced state. DEGs upregulated by paclitaxel were involved mainly in encoding heterodimers of α/β-tubulins and microtubules (TUBB2A, TUBA1A, TUBA1B, TUBB, and MAP2), modulating cholesterol biosynthesis (MSMO1), and regulating autophagy (ATF4; Figure Figure3B3B and and3C).3C). DEGs downregulated by paclitaxel include those previously found in high shear-dependent genes such as ADAMTS1, IGFBP5, THBD, RGS5, and HYAL2 (Figure (Figure3C).3C). GSEA confirmed primary pathways activated by paclitaxel were related to sterol biosynthesis and tubulin formation (Figure (Figure3G).3G). Primary pathways negatively enriched by paclitaxel included regulation of EC plasma membrane, such as ankyrin and integrin binding (Figure (Figure3G).3G). Overall, paclitaxel impacts genes related to cell function and vascular integrity (EMCN, by stabilizing the EC barrier and modulating vascular permeability64; SCD, through lipid metabolism65), cell division, cellular intracellular transport, and maintaining cell shape (TUBA1B, TUBB, TUBA1A, and TUBB2A, major members of microtubules of the cytoskeleton66; MAP2, stabilizing microtubules67), apoptosis (ATF4, by modulating amino acid metabolism and antioxidant responses68), and cell proliferation (MALAT1, a long noncoding RNA regulating the expression of other genes69). We substantiated the influence of paclitaxel at the protein level through an analysis of its effects on α-tubulin. Our observations revealed a significant alteration in cell morphology induced by paclitaxel, regardless of whether the cells were subjected to physiological or altered-flow conditions (Figure S3C, S3D, and S3F).
The EC response to rapamycin independent of shear state revealed upregulated DEGs involved in the regulation of kinase A (AKAP12), cell differentiation and proliferation (HMGA1, NDRG1, and IGFBP5), and calcium dynamics (S100A6; Figure Figure3E3E and and3F).3F). Downregulated DEGs by rapamycin included POSTN, which regulates tissue regeneration and wound healing through cell attachment; GJA4 and GJA5, which encode gap junction proteins; and HSPA5 and HSP90B1, which are key in protein folding and stabilization (Figure (Figure3F).3F). GSEA revealed that the main effect of rapamycin in ECs is to repress multiple biological pathways, including regulation of proteasomes, protein folding, angiogenesis, and calcium ion channels (Figure (Figure3H).3H). Overall, rapamycin impacts genes related to cell migration and adhesion processes (LMO7, via regulation of scaffolding/cytoskeletal organization70), cell proliferation (S100A6, by regulating calcium-binding proteins71; NDRG1 and SERPINB1, well-known tumor suppressors72,73; HMGA1, a chromatin-binding protein modulating transcription of various genes74), and endothelial permeability (CLDN5, through tight junction expression75). The impact of rapamycin on HAECs was additionally examined at the protein level, assessing its influence on proliferation by observing the ki67 expression (Figure S3B and S3E). Our findings demonstrate that rapamycin effectively regulates HAEC proliferation at both the transcriptomic and protein levels, leading to a notable reduction in ki67 expression, regardless of the prevailing flow conditions.
EC Cellular Adaptation to Antiproliferative Therapies Is Altered by Flow-Induced Shear Stress
The combined response of shear forces and antiproliferative therapies was analyzed by first computing a hierarchical clustering of single-cell profiles based on a Euclidean distance matrix. A single heatmap was used to represent the response induced by flow state and either paclitaxel (Figure (Figure4A)4A) or rapamycin (Figure (Figure4C).4C). Identifiable shear states are denoted in both therapies, with drug stimulation clearly displaying a unique fingerprint in ECs when costimulated with high shear stress conditions (Figure (Figure4A4A and and4C).4C). However, the pharmacological response is more limited in low shear stress states, in which the single-cell transcriptional profile cannot differentiate cells stimulated with antiproliferative drugs from control cells (Figure (Figure4A4A and and44C).
Shear stress influences the cellular response of endothelial cells to antiproliferative drugs. A, Hierarchical heatmap displaying the Euclidean distances between differentially expressed genes (DEGs) and cells costimulated with paclitaxel. B, Gene expression levels of 5 relevant markers for each flow/paclitaxel condition. C, Hierarchical heatmap showing the Euclidean distances between DEGs and cells costimulated with rapamycin. D, Gene expression levels of 5 relevant markers for each flow/rapamycin condition. E through H, Weighted gene correlation network analysis (WGCNA) applied to DEGs in high and low shear stress conditions with paclitaxel or rapamycin. Candidate gene modules of minimum size are identified with colored bars, and gene modules with significant intramodule correlation strength are labeled. I and J, Identification of corresponding modules in high and low shear stress conditions using overlap significance testing. K and L, Mean WGCNA-based module scores for each experimental condition, scaled and centered across all cells. Ctrl, control condition (high or low shear stress state without drug); High, high shear stress (30 dynes/cm2); Low, low shear stress (2.5 dynes/cm2); Pacli, paclitaxel stimulation; and Rapa, rapamycin stimulation.
Limitation of EC response to antiproliferative therapy when exposed to altered shear stress was corroborated using trajectory-based differential expression analysis. Subclusters defining EC heterogeneity within flow/drug groups (Figure S8A through S8E and S8H through S8L) were subjected to trajectory analysis, with the dominant pseudotime coordinate primarily driven by the sensitivity of ECs to flow-induced shear stress both in rapamycin and paclitaxel-exposed cells (Figure S8F and S8M). Furthermore, we found that ECs have a higher sensitivity to antiproliferative drugs when exposed to high/physiological flows rather than low/altered flows, showing a more significant difference in the average pseudotime response in the presence of high/physiological shear stress (Figure S8F and S8M).
Relevant DEGs were identified for each flow/drug condition (Figure (Figure4B4B and and4D).4D). DEGs upregulated by high shear stress or downregulated by paclitaxel are involved in cleavage of proteoglycans (ADAMTS9), RGS5, and modulating cellular adhesion (PCDH12). Cholesterol biosynthesis (MHGCS1, INSIG1, and ACAT2) and autocrine/paracrine signaling modulation through action on calcium and phosphate ions were related to paclitaxel response in high laminar flow shear stress state (STC2). Few DEGs offered a unique response to identify paclitaxel stimulation in the low shear stress state. To further validate the relevance of these DEGs, pseudotime trajectory analysis revealed significant changes in pseudotime response for ACAT2, ADAMTS9, INSIG1, PCDH12, RGS5, and STC2 in paclitaxel-treated HAECs in the presence of flow (Figure S8G).
Pseudotime trajectory analysis further validates the fluctuations experienced by DEGs in response to both flow shear stress and rapamycin, showing a higher relative expression change in the regions where the pseudotime corresponds to high/physiological flow and rapamycin-treated HAECs (Figure S8N). Altogether, these findings reveal that EC cellular adaptation depends primarily on shear stress levels. More specifically, in low-flow shear stress states, the EC cellular adaptation to chemical stimuli has been diminished, limiting the response to antiproliferative drugs.
To leverage the ability of scRNA-seq to further dissect the components of EC dynamism in response to antiproliferative therapies in an unbiased fashion, a previously described unsupervised WGCNA-based method coupled with a significance test38 was used to identify coexpressed gene modules consisting of at least 10 genes in the response of flow-conditioned ECs to paclitaxel/rapamycin. The paclitaxel/rapamycin response under high shear stress was defined as the set of DEGs between the high shear stress cells with drug and the high shear stress cells without drug, which consisted of 244 genes in paclitaxel-stimulated cells and 324 genes in rapamycin-stimulated cells (Figure S5C). The paclitaxel/rapamycin response under low shear stress was defined as the set of DEGs between the low shear stress cells with and without drugs, consisting of 66 genes in paclitaxel-stimulated cells and 137 genes in rapamycin-stimulated cells (Figure S5C). Of note, in both cases, the high flow-conditioned drug response consisted of more total genes than the low flow–conditioned drug response. The WGCNA-based module identification and significance testing procedure was applied separately to both sets. This process identified 8 distinct coexpressed modules of at least 10 genes in the high shear stress condition and only 3 modules in the low shear stress condition in cells stimulated with paclitaxel (Figure (Figure4E4E and and4F).4F). For rapamycin-stimulated cells, 10 modules were found in high flow and 3 modules under low flow. (Figure (Figure4G4G and and4H).4H). Overall, fewer DEGs and WGCNA-derived modules were found in low flow conditions in either paclitaxel or rapamycin-stimulated cells, showing a reduced complexity of antiproliferative drug response in low shear–conditioned ECs.
The coexpressed transcriptional modules comprising the response of high shear–cultured HAECs to paclitaxel included a tubulin/cytoskeleton assembly and rearrangement set including TUB family genes (TUBA1B, TUBB, TUBA1A, TUBB2A, TUBB6, and TUBA1C), which clearly correlates with changes in the α-tubulin protein expression, modifications in the cell cytoskeleton, and overall morphology when HAECs are exposed to paclitaxel, both under high and low shear stress conditions, which were not observed when cells were treated only with rapamycin (Figure S3C, S3D, and S3F); a blood vessel morphogenesis set including ADAMTS1 and THBD; a cell adhesion and migration/locomotion set, including platelet and EC adhesion molecules (PECAM1) and integrins (ITGA6 and ITGB4); a supramolecular complex and tubulin/cytoskeleton complex set, including microtubule-associated proteins (MAP2); a sterol metabolic/biosynthetic process set including several enzyme coding proteins such as reductases (HMGCR, DHCR7, and DHCR24) or isomerases (IDI1 and EBP); a ribosome/cytoplasmic translation set, including ribosomal proteins (eg, RPL41 and RPL27); a tube development/morphogenesis set and one unannotated gene set without known functional gene enrichment including KLF6, SERPINE2, and ANKRD1 (Table S1). Coexpressed transcriptional modules comprising the simpler response of EC culture under low shear stress and paclitaxel included a tubulin/cytoskeleton assembly a rearrangement set, an extracellular matrix/cell surface set, and a cell adhesion set (Table S2). Twenty-three genes were present in both the high shear response modules and the low shear response modules. Of these, 7 were members of the respective tubulin/cytoskeleton assembly and arrangement modules (High.M1 and Low.M1). Notably, the remainder were distributed among both the other low shear–response modules, but 2 high shear–response modules (High.M3 and High.M7) included no common genes.
Similar analysis was performed on rapamycin data, with the coexpressed transcriptional modules comprising the response of high shear–cultured ECs to rapamycin including an extracellular matrix/cell surface molecule set with GJA4 and GJA5 molecules; a cell adhesion/migration set including TMEM100 and BMPER; an endoplasmic reticulum/antigen processing set including B2M, CALR, and PDIA3; a vasculature development set including TUBB; a peptide/ribosomal synthesis gene set; 2 additional segregated sets also related to cell adhesion and endoplasmic reticulum processing; and a mixed set enriched mainly for regulation of protein metabolic process but including mechanotransduction genes such as PIEZO2. Two unannotated gene sets without known functional gene enrichments included an IL-33/COL1A2/SNX3 module and a PLXNA2/TUBA1C/SLC26A2 module (Table S3). All these modules affecting cell adhesion, extracellular matrix development, and cell migration are indirectly correlated with the well-known effect of rapamycin on cell proliferation, as it affects the cell cycle. These effects are relevant in both high/physiological and low/altered shear stress conditions, as evidenced by the diminished expression of ki67 when cells are treated with rapamycin (Figure S3B and S3E). Coexpressed transcriptional modules comprising the simpler response of ECs cultured under low shear stress and rapamycin included a peptide/ribosomal synthesis gene set, an endoplasmic reticulum processing gene set including CALR, PDIA6, and MMP14, and cell anchoring/junction organization set including GJA4 and GJA5. Thirty genes were present in both the high shear response modules and the low shear response modules (Table S4). Of these, 13 were members of the low shear endoplasmic reticulum model modules (Low.M1) and equally divided in the respective high shear modules (High.M5 and High.M7); and 12 were members of the low shear cell anchoring/junction organization module (Low.M2) and unevenly divided in the respective high shear module (10 genes in High.M3). The remainder were distributed among both the other low shear response modules but 4 high shear response modules included no common genes (High.M1, High.M4, High.M9, and High.M10).
To use this method of decomposing the EC cellular adaptation to paclitaxel/rapamycin to understand the effect of shear on drug response, we identified significant overlaps between the modules from each shear condition with the goal of functionally linking modules under both conditions (Figure (Figure4I4I and and4J).4J). No significant linkages were identified in the paclitaxel data set. The tubulin/cytoskeleton assembly and rearrangement modules, main paclitaxel targets, showed a substantial but nonsignificant overlap in each flow state. Notably, the remaining modules of paclitaxel-stimulated cells (Low.M2 and Low.M3) present in the low shear set had no significant overlaps with multiple modules in the high shear paclitaxel response, including the blood vessel morphogenesis module, the cell adhesion and migration/locomotion module, and the tube development/morphogenesis module. Regarding rapamycin-stimulated cells, a significant overlap was present between the endoplasmic reticulum processing module of the 2 flow responses, as well as with the cell anchoring/junction organization set. Another group with strong correspondence in both flow conditions was peptide/ribosomal synthesis gene sets. A weaker overlap was also present between the cell anchoring/junction organization module in the low shear state and the cell adhesion/migration module in the high shear state.
The differences between the high and low shear responses to paclitaxel/rapamycin were further clarified by scoring expression of each module across all cells and examining differences in mean scores among shear conditions (Figure (Figure4K4K and and4L).4L). In general, rapamycin-induced changes in the modules identified under high shear were blunted under low shear stress, while changes in the modules identified under low shear had similar magnitude in both conditions (Figure (Figure4L).4L). Paclitaxel-induced changes under high shear were also blunted under low shear stress, but changes identified under low shear had either a similar magnitude in both conditions or were significantly decreased in high conditions (Figure (Figure4K).4K). In addition to the mean scores, score distributions were examined across conditions (Figures S6 and S7). Unsurprisingly, the addition of paclitaxel/rapamycin resulted in a significant change in expression of all modules in both shear conditions. However, most of the modules first identified via the paclitaxel/rapamycin response were also significantly altered by shear alone, in keeping with the interrelation of the 2 stimuli.
Finally, we contrasted the combined response of paclitaxel and rapamycin under both high and low shear stress conditions (Figure S5A and S5B). As expected, there were only a few common DEGs in the different flow conditions treated with these antiproliferative mechanisms (Figure S5C). Most of the relevant DEGs showed different expression levels for each drug/flow condition (Figure S5E) due to their distinct mechanisms of action: rapamycin inhibits the mTORC1 complex, leading to cell cycle arrest at the G1 phase, while paclitaxel binds and stabilizes microtubules, resulting in cell arrest in the mitotic phase.
Unsurprisingly, many of the regulated programs revealed with gene enrichment were mainly driven by TUBB family genes (Figure S5D), which encode β-tubulin proteins crucial for microtubules, the target of paclitaxel, and indirectly affected by rapamycin due to the mTOR pathway’s influence on microtubule dynamics as a downstream effector when regulating cell metabolism, proliferation, and growth.
DISCUSSION
Understanding EC cellular adaptation is fundamental to understanding the regulation of blood vessel homeostasis. ECs are exposed to multiple categories of stressors that alter the homeostatic quiescent state of the healthy, intact endothelium. Therefore, elucidating EC cellular adaptation in response to discrete stressors is a requirement for further advances in vascular biomedicine, from finding novel biomarkers for cardiovascular disease prognosis to explaining the high variability in endovascular interventional outcomes for cardiovascular disease. Here, we use scRNA-seq to profile HAEC populations subjected to differing shear stresses found in the human aorta as a paradigmatic mechanical stimulus, both alone and in combination with antiproliferative agents as paradigmatic chemical stimuli. Our results enable definition of EC capacity to respond to these environmental factors with single-cell resolution, which in turn provides new insights into the transcriptional modules comprising the EC state, as well as the coupling among these modules that gives rise to functional cellular adaptation in the context of environmental stimuli singly or in combination. Our results demonstrate the following advances in the study of EC functional dynamism: (1) shear- and drug-induced states can be defined with high resolution via the single-cell EC transcriptome; (2) shear-induced EC states correspond to previously identified in vivo EC heterogeneity; (3) shear stress modulates the EC adaptation to other stimuli; and (4) the low/altered shear-induced EC state limits the complexity of the response to widely used pharmacological regulation. ECs committed to respond to low shear, as found in the atheroprone regions of the aorta, have a reduced ability to alter their state in response to pharmacology, here, to the chief classes of drugs used in vasoactive implants (ie, paclitaxel and rapamycin).
The understanding of changes induced by shear stress has evolved over decades as the linchpin of investigations into the EC state.7 At a population level, our transcriptional profile demonstrates a generalized increase in the expression of ribosomal and mitochondrial transcripts under low shear stress (Figure S2C and S2D). Low shear stress induces endoplasmic reticulum stress and the related unfolded protein response pathway.9 Typically, activation of the unfolded protein response is associated with translational repression, making the association of low shear with ribosomal transcript upregulation in our data set surprising; however, a distinct feedback arm in the unfolded protein response mediated by ATF4 and GADD34 has the effect of restoring protein synthesis and is linked to ribosomal biogenesis via NMP4.76 This is particularly relevant in the setting of long-term endoplasmic reticulum stress, as might be expected during the multiday bioreactor model of shear conditioning used in this study.77 The increase in mitochondrial transcripts may be partially explained by the previous suggestion that mitochondrial biogenesis is protective against endoplasmic reticulum stress.78
A value of extending the characterization of shear-induced EC state to the single-cell level is in placing previous results in context by identifying the top TF-based programs that are most important for defining shear state and to enable investigation of how they are coupled to each other across a population of ECs. Given the complex and multifactorial drivers of the EC state, we accomplished this by conceptualizing shear-driven EC state as the combination of multiple TFs for each flow condition (Figure (Figure1F).1F). Our results demonstrate that a relatively small number of TFs based on their shear-responsive activity can capture shear-induced state sufficiently to differentiate shear-induced states. Although these top TFs do not include the full complexity of the original scRNA-seq profiles, they serve as an interpretable bridge between the complete transcriptome and prior knowledge on shear-responsive elements in ECs. This analysis enables direct comparison of TF modules to determine which factors are most important for definition of shear stress–induced state. FOSL1, HNF4A, VDR, CDX2, and RARA were most important for defining high shear–induced state, while ETV1, NFKB2, RUNX2, and TP53 were the greatest contributors to low shear–induced identity. Understanding endothelial dynamism from the perspective of discrete TFs with single-cell resolution enables a middle ground between the traditional model of dysfunctional or activated endothelium and the molecular biology of individual pathways. This same paradigm can be used to characterize EC adaptation in response to a variety of metabolic stressors in terms of shifts in the network of coupled TF-based activity.
Our single-cell transcriptional profiling also demonstrated that the components of shear state are not all expressed in the same cells but are instead upregulated differentially among different shear-cultured HAECs, as observed previously in human umbilical vein ECs.25 This suggests that the endothelial response to shear stress is not monolithic but instead consists of distinct but correlated transcriptional programs, which may be linked to distinct phases or degrees of shear response in combination with other stimuli, as evaluated in our pseudotime analysis. This heterogeneity has been further investigated as a potential origin of spatial EC heterogeneity observed in vivo,21 revealing well-localized EC subpopulations in the aorta expressing different flow sensitivities. The top flow-sensitive DEGs and TFs found here in vitro enabled the identification of EC subclusters previously defined in vivo.
The components of the endothelial response to antiproliferative drugs are much less well described than those comprising the response to shear stress, although they have been linked to multiple effects on endothelial processes including migration and proliferation.14 From a population perspective, our data on the response of high shear–cultured ECs to paclitaxel and rapamycin provide new insight into their effects on EC functional phenotype, including changes in the expression of transcripts implicated in sterol biosynthesis, tubulin formation and vasculature development in paclitaxel-stimulated cells and endoplasmic reticulum handling, cell adhesion/anchoring, and angiogenesis in rapamycin-stimulated cells. Intriguingly, exposure to low shear stress limits the antiangiogenic potential of both paclitaxel and rapamycin on ECs as identified here both by gene set enrichment on high shear–specific paclitaxel/rapamycin response genes and by identifying discordant genes in corresponding paclitaxel/rapamycin response modules under differing shear conditions. Functional analyses of this type are necessary to understand high-resolution transcriptional data in terms of implications for vascular pathology and therapy. Thus, transcriptional profiling of this type enables both new understanding of the functional effects of a complex multifactorial stimulus and how these effects are modulated by additional stressors.
As with shear stress alone, the characterization of drug (paclitaxel/rapamycin)-exposed populations with single-cell resolution enabled identification of discrete transcriptional modules as components of EC adaptation. ECs subjected to high shear had a more complex response to antiproliferative drugs, consisting not only of more DEGs but also with more distinct transcriptional modules with significant correlation among component genes. Many of the modules identified in the response of high shear–conditioned ECs corresponded to the major EC functions identified in the bulk response, while others had no clearly enriched function. By contrast, the response of low shear–conditioned ECs to antiproliferative drugs consisted of fewer genes with fewer distinct modules with less well-defined functions, with a notably blunted response in genes associated with angiogenesis, vasculature development, or cell migration. In rapamycin-exposed cells, the endoplasmic reticulum–based module in the high shear–exposed cells notably included certain homeostatic functions such as antigen processing, which were lost in the low shear modules, possibly due to engagement of endoplasmic reticulum machinery by low shear–induced unfolded protein response. This result has an important implication for our understanding of stimulus-driven EC functional heterogeneity—namely, that the response to one external stimulus can limit the metabolic response to a second stimulus (rapamycin). Dysfunctional EC processes including endoplasmic reticulum stress and the unfolded protein response may represent forcing states, which dampen shifts driven by other stimuli. This is a vital consideration for the development and understanding of other vasoactive therapeutics, which as a rule are deployed in the context of endothelium already made dysfunctional with vascular disease. This method for module-based understanding of EC adaptation can also be utilized for other combinatorial stimuli—whether pharmacological or pathological—with complex effects on ECs that are as yet poorly understood.
Our study has shed light on the effects of single doses of antiproliferative agents on transcriptional changes in ECs, as well as the impact on proliferation and cytoskeletal protein function. However, it is essential to acknowledge the limitations that may influence the interpretation of our findings. First, the use of a single dose for each agent, while validated for inducing significant changes, may warrant further investigation with multiple doses to strengthen the robustness of our conclusions. Second, though 3 replicates were included to enhance reliability, the lack of access to the genetic backgrounds of donors may limit the generalizability of our results, especially when studying populations with diverse genetic backgrounds. Third, many of the transcriptional programs identified in this work as responsive to flow and antiproliferative drugs include substantial posttranscriptional components, which must be investigated in future work applying proteomic and functional assays. Despite these limitations, our research provides valuable insights and highlights the EC sensibility to flow and the altered capacity of EC to respond to antiproliferative treatments when cells are committed to distinct pathophysiological flows.
In conclusion, we present here the use of scRNA-seq to profile EC populations subjected to controlled external stimuli singly and in combination to elucidate the distinct components of EC cellular adaptation. We demonstrate how shear stress limits EC response to other stimuli, including antiproliferative drugs delivered during endovascular treatments. ECs are above all environmentally sensitive cells, especially driven by flow and flow-induced shear stress, and serve to discern distinct EC subpopulations is the aorta. As such, it seems that once exposed to flow alterations, ECs commit to responding to this stress state, altering their capacity to remain responsive to subsequent stresses, for example, chemical exposures. Continued understanding of the function and dysfunction of distinct endothelial subpopulations and their potential response to vasoactive therapeutics will depend on continued utilization of these techniques for probing ECs, especially at the protein level.
ARTICLE INFORMATION
Sources of Funding
A.G. Salazar-Martín received the support of a fellowship from La Caixa Foundation (ID 100010434; LCF/BQ/AA1 9/11720044) and the MathWorks Engineering Fellowship. These fellowships were used to support his PhD studies at the Harvard/Massachusetts Institute of Technology (MIT) Medical Engineering and Medical Physics Program. A.S. Kalluri was supported by the National Institutes of Health (NIH) National Institute of General Medical Sciences via grants to the Harvard/MIT Medical Scientists Training Program and the Harvard Biophysics Program (T32GM007753 and T32GM008313), as well as the American Heart Association Predoctoral Fellowship (18PRE34060009). A.K. Shalek was supported by the Searle Scholars Program, the Beckman Young Investigator Program, a Sloan Fellowship in Chemistry, and the NIH (2RM1HG006193 and 5U24AI118672). E.R. Edelman was supported, in part, by the NIH (R01GM04903924 and R01 HL161069).
Disclosures
A.S. Kalluri reports compensation for consulting from nference, inc. A.K. Shalek reports compensation for consulting and scientific advisory board membership from Merck, Honeycomb Biotechnologies, Cellarity, Repertoire Immune Medicines, Ochre Bio, and Dahlia Biosciences. The other authors report no conflicts.
Supplemental Material
Major Resources Table
Figures S1–S8
Tables S1–S4
Publication and Licensing Rights for the Graphic Abstract, Figures 1A and 2A and Figures S1 and S4A
Nonstandard Abbreviations and Acronyms
- DEG
- differentially expressed gene
- EC
- endothelial cell
- GSEA
- gene set enrichment analysis
- HAEC
- human aortic endothelial cell
- IL
- interleukin
- mTOR
- mammalian target of rapamycin
- PC
- principal component
- RGS5
- regulation of G protein signaling 5
- scRNA-seq
- single-cell RNA sequencing
- TF
- transcription factor
- TNF-α
- tumor necrosis factor alpha
- WGCNA
- whole-genome correlation network analysis
*A.G. Salazar-Martín and A.S. Kalluri contributed equally.
For Sources of Funding and Disclosures, see page 2279–2280.
Supplemental Material is available at https://www.ahajournals.org/doi/suppl/10.1161/ATVBAHA.123.319283.
REFERENCES
Full text links
Read article at publisher's site: https://doi.org/10.1161/atvbaha.123.319283
Read article for free, from open access legal sources, via Unpaywall:
https://www.ahajournals.org/doi/pdf/10.1161/ATVBAHA.123.319283
Citations & impact
This article has not been cited yet.
Impact metrics
Alternative metrics
Discover the attention surrounding your research
https://www.altmetric.com/details/154700539
Data
Data behind the article
This data has been text mined from the article, or deposited into data resources.
GEO - Gene Expression Omnibus
- (1 citation) GEO - GSE212388
Similar Articles
To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.
Smooth muscle cells orchestrate the endothelial cell response to flow and injury.
Circulation, 121(20):2192-2199, 10 May 2010
Cited by: 40 articles | PMID: 20458015 | PMCID: PMC2887340
Flow pattern-dependent endothelial cell responses through transcriptional regulation.
Cell Cycle, 16(20):1893-1901, 18 Aug 2017
Cited by: 42 articles | PMID: 28820314 | PMCID: PMC5638382
Review Free full text in Europe PMC
A comparative study of the effects of vitamin C, sirolimus, and paclitaxel on the growth of endothelial and smooth muscle cells for cardiovascular medical device applications.
Drug Des Devel Ther, 7:529-544, 24 Jun 2013
Cited by: 12 articles | PMID: 23836963 | PMCID: PMC3699137
Differential effects of orbital and laminar shear stress on endothelial cells.
J Vasc Surg, 41(5):869-880, 01 May 2005
Cited by: 174 articles | PMID: 15886673
Funding
Funders who supported this work.
NHGRI NIH HHS (1)
Grant ID: RM1 HG006193
NHLBI NIH HHS (1)
Grant ID: R01 HL161069
NIAID NIH HHS (1)
Grant ID: U24 AI118672
NIGMS NIH HHS (3)
Grant ID: T32 GM007753
Grant ID: T32 GM008313
Grant ID: R01 GM049039