Abstract

Background

   Cardiovascular diseases (CVDs) are the leading global cause of death,
   with atherosclerosis as the primary cause. Chronic inflammation,
   endothelial dysfunction, and the role of molecules like nitric oxide
   and reactive oxygen species are crucial in this context. Our previous
   research indicated that cilostazol and ginkgo biloba extract could
   enhance the ability of endothelial cells to dissolve blood clots, but
   the effects of cilostazol on monocytes remain unexplored.

Method

   This study utilized peripheral blood mononuclear cells from 10 healthy
   donors, treated ex vivo with cilostazol. RNA-sequencing,
   over-representation analysis, xCell stromal cell analysis, and Gene Set
   Enrichment Analysis were employed to investigate the gene expression
   changes and biological pathways affected by cilostazol treatment.

Results

   The study identified specific gene sets and pathways that were enriched
   or reduced in response to cilostazol treatment, providing insights into
   its effects on monocytes and potential therapeutic applications in CVD.
   The analysis also revealed the potential impact of cilostazol on the
   stromal cell compartment, further broadening our understanding of its
   multifaceted role.

Conclusion

   The findings offer a nuanced understanding of the advantages and
   mechanisms of cilostazol in CVD, uncovering novel therapeutic targets
   and strategies to enhance the clinical application of cilostazol and
   contributing to the broader implications of this therapy in
   cardiovascular health.

1. Introduction

   Cardiovascular diseases (CVDs) are the leading cause of death globally,
   significantly contributing to health loss and high healthcare system
   costs [[33]1]. Atherosclerosis, characterized by arterial plaque
   buildup, is the primary cause of CVD. Chronic inflammation, linked with
   endothelial cell dysfunction, appears to govern the disease process
   [[34]2,[35]3], necessitating clinical research into anti-inflammatory
   and immunomodulatory treatments. Inflammation is a key factor in
   atherosclerosis and CVD. The roles of immune-active cells,
   pro-inflammatory cytokines, chemokines, and lipid mediators are crucial
   in this context [[36]4]. Dead cells and the oxidized forms of
   low-density lipoprotein are particularly abundant in atherosclerosis
   [[37]5]. Oxidized low-density lipoprotein, known for its immunogenic
   properties, can activate endothelial cells, leading to monocyte
   adhesion and differentiation into macrophages [[38]6,[39]7]. Within
   atherosclerotic plaques, macrophages tend to accumulate and exhibit
   reduced migratory ability, leading to unresolved inflammation and
   progression to more complex plaques [[40]8,[41]9]. Macrophages
   contribute to inflammation by secreting pro-inflammatory mediators and
   matrix-degrading proteinases. Dying macrophages release lipid
   components and tissue factors, forming a prothrombotic necrotic core, a
   key component of unstable plaques that can rupture.

   The “Response-to-Injury Hypothesis” by Ross and Glomset outlines the
   potential process from endothelial dysfunction to atherosclerosis
   [[42]10]. This hypothesis has been instrumental in understanding the
   underlying mechanisms that lead to atherosclerotic plaque development.
   A key characteristic of endothelial dysfunction is the loss of nitric
   oxide (NO) bioactivity [[43]3], a critical molecule involved in
   vascular homeostasis. NO plays a vital role in maintaining vascular
   tone, inhibiting platelet aggregation, and preventing the adhesion of
   monocytes to the endothelium. The loss of NO bioactivity can lead to a
   cascade of events contributing to atherosclerosis. One of the factors
   closely related to NO production is the generation of reactive oxygen
   species (ROS). ROS are molecules that play a dual role in the vascular
   system. On one hand, elevated ROS levels are associated with reduced NO
   bioavailability [[44]11]. This reduction in NO bioavailability can lead
   to vasoconstriction, increased platelet aggregation, and enhanced
   adhesion of monocytes to the endothelium, which are key events in the
   initiation of atherosclerosis. ROS may also cause an increase in
   monocyte-endothelium adhesion [[45]12], promoting the inflammatory
   response that characterizes atherosclerosis. On the other hand, ROS are
   essential for macrophage differentiation [[46]13], a process central to
   the immune response and atherosclerotic plaque progression. Macrophages
   are key players in the inflammatory response, and their differentiation
   and function are tightly regulated by various signaling pathways,
   including those involving ROS.

   Cilostazol, characterized by its antithrombotic, vasodilatory,
   antimitogenic, and cardiotonic properties, has been broadly
   investigated because of its multifaceted efficacy in CVD management
   [[47]14]. Its application in CVD is substantiated by its demonstrated
   superiority over clopidogrel in reducing the incidence of recurrent
   ischemic strokes in those with non-cardioembolic ischemic stroke
   [[48]15]. Furthermore, the combination of cilostazol and aspirin is
   comparable or superior to that of aspirin and clopidogrel in preventing
   stent restenosis [[49]16]. Regarding cardiovascular event prevention,
   cilostazol has been associated with a notably low rate of serious
   events, including myocardial infarction, observed in only 1% of the
   patients treated with cilostazol [[50]17]. Lin et al. further
   elucidated the drug's efficacy in significantly reducing the risk of
   major coronary events, major adverse cardiovascular and cerebrovascular
   events, and symptoms of angina pectoris [[51]18]. The influence of
   cilostazol on circulating endothelial progenitor cell counts, which are
   pivotal markers in cardiovascular health, also surpasses that of
   aspirin, suggesting an additional therapeutic benefit [[52]19].
   Moreover, combining cilostazol with aspirin may mitigate biological
   resistance to aspirin in ischemic stroke [[53]20]. Crucially,
   cilostazol significantly reduces the risk of cerebrovascular events in
   patients with atherothrombosis, without a corresponding increase in
   bleeding risk, thereby offering a safety advantage over other
   antiplatelet agents [[54]21]. This safety profile is particularly
   relevant for patients with a history of hypertension or
   gastrointestinal bleeding, where cilostazol may be a safer alternative
   than clopidogrel [[55]22].

   Our previous research focusing on endothelial cells discovered that
   both cilostazol and ginkgo biloba extract could enhance the expressions
   of endothelial cell thrombomodulin (TM) and endothelial NO synthase
   (eNOS) by activating the Krüppel-like factor 2 (KLF2) axis. This
   activation subsequently increased the endothelial cells' ability to
   dissolve blood clots [[56][23], [57][24], [58][25]]. These findings
   provide valuable insights into the potential therapeutic applications
   of cilostazol and ginkgo biloba extract in cardiovascular health.
   However, the literature has not explored whether cilostazol has any
   other effects on monocytes. This knowledge gap prompted our current
   study, aiming to shed light on the broader implications of cilostazol
   in CVD and to enhance the precision of its clinical application. In
   this study, we utilized peripheral blood mononuclear cells (PBMC) from
   healthy donors and treated them ex vivo with cilostazol. We
   comprehensively analyzed the gene expression changes induced by
   cilostazol treatment through RNA sequencing. Over-representation
   analysis was employed to identify the biological pathways significantly
   affected by the treatment. Additionally, we utilized xCell immune and
   stromal cell analysis to investigate the potential impact of cilostazol
   on the stromal cell compartment. Gene Set Enrichment Analysis (GSEA)
   further allowed us to identify specific gene sets and pathways that
   were enriched or depleted in response to cilostazol treatment. The
   findings from this study are expected to provide a more nuanced
   understanding of the advantages and mechanisms of cilostazol in the
   context of CVD. By exploring the effects of cilostazol on monocytes, we
   hope to uncover novel therapeutic targets and strategies that could
   enhance the efficacy and precision of cilostazol in clinical settings.

2. Material and methods

2.1. IRB approval and PBMC treatment

   This study was conducted with the approval of the institutional review
   boards at both the Tri-Service General Hospital and the Taipei Veteran
   General Hospital, under IRB number 2012-03-001AC. In accordance with
   the ethical guidelines set forth in the 1975 Helsinki Declaration, all
   participants in this study provided their written informed consent.
   Furthermore, all experimental procedures and protocols involving
   animals were not only approved by the Institutional Animal Care and Use
   Committee (IACUC) of National Yang-Ming University but were also in
   full compliance with the ARRIVE guidelines.

   Ten healthy donors were recruited for the study, and blood samples were
   obtained. PBMC were isolated from the blood samples using Ficoll-Paque
   Plus (Sigma-Aldrich, Inc., St. Louis, MO, USA, Cat.: GE17-1440-02) in
   Leucosep™ tubes (Greiner Bio-One, Kremsmünster, Austria), according to
   the manufacturer's instructions. Following isolation, the PBMC were
   cultured in RPMI 1640 medium, supplemented with 10% Fetal Bovine Serum
   (FBS, Thermo-Fisher Scientific Inc., Waltham, MA, USA) and 1%
   penicillin/streptomycin. The cells were treated with 30 μM cilostazol
   and an appropriate vehicle control and were incubated for 24 h at
   standard culture conditions.

2.2. RNA sequencing and differentially expressed genes (DEG) analysis

   Treated PBMCs were subjected to RNA extraction and analysis via
   next-generation sequencing. The extraction process was facilitated by
   the use of TRIzol Reagent (Thermo-Fisher Scientific). Following
   extraction, the RNA samples were quantified and analyzed through the
   NanoDrop device (Thermo-Fisher Scientific), as well as the Agilent 2100
   Bioanalyzer system (Agilent Technologies). The visualization and
   integrity of the RNA were evaluated using 1% agarose gel
   electrophoresis. Subsequently, we purified the mRNA using both Poly(A)
   mRNA Bead Isolation Kits and rRNA Removal Kits. First-strand cDNA
   synthesis was accomplished using ProtoScript II Reverse Transcriptase
   following mRNA fragmentation with First Strand Synthesis Reaction
   Buffer and Random Primer. We then synthesized Second Strand DNA by
   leveraging the Second Strand Synthesis Enzyme Mix. DNA fragments were
   then polished using End Prep Enzyme Mix. Adaptors were attached to both
   ends of the fragments via a process of dA-tailing and T-A ligation. A
   420bp fraction of fragments was subsequently isolated through the
   application of magnetic beads. A 13-cycle PCR was then performed using
   P5 and P7 primers. The resulting library was quantified with a Qubit
   3.0 fluorometer and sequenced on an Illumina HiSeq 4000 platform. The
   resulting sequencing data in fastq format was processed and generated
   by GENEWIZ®. Transcript abundance for each sample was estimated by
   separately running two iterations using both Salmon (v0.91) and
   Kallisto (v0.46.2) in combination with the GRCh38 human reference
   transcript from GENCODE v41 [[59]26]. These tools quasi-mapped the
   sequencing reads using a 31bp k-mer index. Subsequently, a third
   iteration was conducted using Salmon (v0.91) in combination with the
   human reference transcript downloaded from Ensembl
   (Homo_sapiens.GRCh38.cdna.all.fa.gz) [[60]27]. The tximport package
   (v1.24.0) facilitated the importation of transcript abundance estimates
   from the files generated by Salmon and Kallisto for analysis in the R
   (v4.2.2) statistical programming environment. Transcript abundance
   estimates were then aggregated into gene expression counts at the gene
   level. The resultant raw counts were merged for use in downstream
   analyses. DEGs for cilostazol treatment were identified using the EdgeR
   package (v3.40.2). Criteria for selection of differentially expressed
   genes: |log2 fold change| ≥ 1; Benjamini-Hochberg (B–H) adjusted
   P-value <0.05.

2.3. Over representative analysis and visualization

   Details on differentially expressed genes (DEGs) are provided in the
   supplementary table. To interpret the potential biological implications
   of these DEGs, over-representation analysis (ORA) was carried out to
   uncover significantly enriched biological processes, molecular
   functions, and cellular components. The DEGs were categorized into
   three clusters based on stringent criteria: Cluster “>1” with an
   absolute log2 fold change greater than 1 and FDR less than 0.05,
   resulting in 546 up-regulated and 631 down-regulated genes. Cluster
   “>2” with an absolute log2 fold change greater than 2 and FDR less than
   0.05, resulting in 129 up-regulated and 125 down-regulated genes.
   Cluster “>3” with an absolute log2 fold change greater than 3 and FDR
   less than 0.05, resulting in 54 up-regulated and 32 down-regulated
   genes. Gene Ontology (GO) enrichment analysis for biological processes
   (BP), molecular functions (MF), and cellular components (CC) was
   performed on these three clusters using the compareCluster function of
   the clusterProfiler package (v4.4.4) in R [[61]28]. The fun was set to
   “enrichGO” with ontologies (ont) set as "BP", "MF", and "CC", and the
   showCategory parameter set to 3. To streamline the representation of
   similar gene sets, the “simplify” function was utilized with a cutoff
   of 0.7 by adjusted p-value (by = "p.adjust"). Visualization of the
   enriched GO terms was performed using the “dotplot” function, with the
   dot size determined by the negative log of the FDR (-logFDR) and color
   indicating the gene ratio. For pathway enrichment analysis, a similar
   approach was adopted, but the “fun = enricher” was employed. The
   analysis incorporated KEGG_2021_Human (320 gene sets),
   WikiPathway_2021_Human (622 gene sets), and Reactome_2022 (1818 gene
   sets) databases, which were downloaded from the EnrichR platform
   [[62]29].

   For the visualization of over-representation analysis (ORA), Cytoscape
   (v3.9.1) along with the ClueGo app was utilized, focusing on the 546
   up-regulated and 631 down-regulated genes identified in our study
   [[63]30,[64]31]. The ClueGo app in Cytoscape was set up to reference
   multiple databases, including Gene Ontology biological processes
   (GOBP), Gene Ontology molecular functions (GOMF), Gene Ontology
   cellular components (GOCC), KEGG, Reactome, and WikiPathways with the
   default parameters.

2.4. Gene Expression Omnibus (GEO) database access

   Public datasets were accessed from the National Center for
   Biotechnology Information's Gene Expression Omnibus (NCBI GEO) database
   to validate our findings related to cardiovascular diseases. Dataset
   [65]GSE19151 included 70 adults with prior venous thromboembolism (VTE)
   on warfarin and 63 healthy controls, analyzed using Affymetrix Human
   Genome U133A 2.0 Array (log2 RMA signal) [[66]32]. Dataset [67]GSE28829
   contained samples from atherosclerotic carotid artery segments,
   including 16 advanced and 13 early plaques, analyzed with Affymetrix
   Human Genome U133 Plus 2.0 Array (RMA-calculated signal intensity)
   [[68]33]. Dataset [69]GSE100927 comprised samples from 69
   atherosclerotic and 35 control arteries from deceased organ donors,
   analyzed using Agilent-039494 SurePrint G3 Human GE v2 8 × 60K
   Microarray (normalized signal intensity ln) [[70]34]. All datasets were
   obtained from publicly available GEO database repositories, and the
   corresponding microarray data and patient follow-up information were
   used for analysis, with raw data preprocessed and normalized as
   required for the specific microarray platforms, and appropriate
   statistical methods applied for validation of our findings.

3. Cell subsets enumeration using xCell approach

   For a comprehensive understanding of the cellular composition of the
   samples, the xCell analysis was employed to enumerate various cell
   subsets in the gene expression data [[71]35]. This analysis was carried
   out using the IOBR package (version 0.99.9) in R (v4.2.0), specifically
   employing the deconvo_tme function with the method parameter set to
   “xcell” and arrays set to “TRUE” [[72]36]. The 64 cell types identified
   were categorized into five major classes: Lymphoid, Stem cells,
   Myeloids, Stromal Cells, and Others. This categorization allowed for an
   in-depth analysis of the cell subsets within the samples, enabling a
   nuanced understanding of their respective roles and interactions. A
   heatmap was constructed to visualize the hierarchical clustering of the
   data. Utilizing the Euclidean distance as the metric, the clustering
   was executed using the Morpheus web tool. The raw data were
   standardized using Z-score normalization to ensure comparability across
   the different samples. Statistical analysis was performed using a
   t-test to compare the samples to the control group. The results were
   visualized through violin plots, drawn using the ggplot2 package
   (v3.4.2) in R.

3.1. Machine learning algorithms

   Supervised machine learning algorithms were employed to predict the
   relative scores of PBMCs being identified as diseased post cilostazol
   treatment. To predict the relative likelihood of PBMCs being classified
   as diseased following cilostazol treatment, we employed a series of
   supervised machine learning algorithms. Machine learning analysis was
   conducted using Orange (v3.35) [[73]37]. The xCell signatures of each
   database were imported into Orange via the widget “CSV File Import”
   feature. For cross-validation, we used the widget “Data Sampler”. We
   created 10 subsets and reserved 1 as an unused subset. We configured
   the widget “Neural Network” with 3,6,3 neurons in the hidden layers and
   employed ReLu as the activation function. The Adam method was selected
   as the solver and regularization was set to α = 0.0001. The maximum
   number of iterations was set to 200, and replicable training was
   enabled. The “Random Forest” widget was configured with 500 trees. The
   option in Growth Control "Do not split subsets smaller than" was set to
   5. For logistic regression, we utilized the widget “Logistic
   Regression”. The regularization type was set to Ridge (L2), and
   strength was set at C = 1. Widget “Test and Score” was used for
   cross-validation with a number of folds set to 10. The "stratified"
   option was selected. ROC Analysis was carried out using the widget “ROC
   Analysis”. The target was set to venous thromboembolism (VTE), advanced
   atherosclerotic plaque, and atherosclerotic artery, respectively. The
   widget “Prediction” was used to calculate the cardiovascular disease
   score for samples of PBMCs both with and without cilostazol treatment.
   The predicted scores were saved using the widget “Save Data” feature.
   The resulting data were then exported to R for visualization and
   statistical analysis using the ggplot2 package (v3.4.2). Specifically,
   bar plots were created to illustrate the predictive scores of
   cardiovascular diseases for the different treatment groups.

3.2. Gene Set Enrichment Analysis (GSEA)

   GSEA was executed to discern gene sets that showed statistically
   significant, concordant differences between cilostazol-treated and
   untreated samples. The analysis was carried out using the gseaplot2
   function within the enrichplot package (v1.16.2). The gene sets were
   derived from various sources, including the vascular endothelial cell
   activation pathway related to growth factor, NO, and ROS in triggering
   vascular inflammation from the Elsevier Pathway Collection; fibrin
   clot-associated gene sets from the Reactome pathway; regulation of ROS
   biosynthetic process and metabolic process from the Gene Ontology
   Biological Process (GOBP); macrophage behavior from the MGI Mammalian
   Phenotype Level 4 2021, filtered by “macrophage”; and macrophage
   biomarker gene sets from the CellMarker database, filtered by
   “macrophage” and “peripheral blood.” All these gene sets are available
   for download from enrichR [[74]29]. A Gene-Concept Network was
   constructed using the cnetplot function in the enrichplot package
   (v1.16.2), and a heatmap-like plot for functional classification was
   created with the heatplot function in the same package.

3.3. Cell viability assay and reactive oxygen species production assay

   The production of ROS in Human Umbilical Vein Endothelial Cells (HUVEC)
   was assessed to gauge cellular oxidative stress status. Initially, cell
   viability was evaluated using the MTT assay (Sigma-Aldrich, Cat.:
   M5655), treating cells with vehicle (DMSO, Thermo-Fisher Scientific;
   Cat.: 036480.K2) and cilostazol at concentrations of 10 μM, 30 μM, and
   100 μM. After a 24-h incubation, optical density (OD) was measured to
   calculate relative cell viability. Subsequently, ROS production was
   quantified using the H2DCFDA (Thermo-Fisher Scientific; Cat.: D399)
   method. This fluorescent dye is oxidized by ROS into a highly
   fluorescent compound, allowing ROS detection within the cells using a
   BD FACSCanto™ II flow cytometry system (Becton, Dickinson and Company,
   Franklin Lakes, NJ, USA). Cells were seeded and allowed to adhere
   before treatment with cilostazol (0.1 μM, 1 μM, 10 μM, and 30 μM), then
   treated with H2DCFDA and incubated under specified conditions for ROS
   production. Following incubation, ROS levels were detected and
   quantified through flow cytometry.

3.4. Flow cytometry detecting macrophage differentiation markers

   Flow cytometry was employed to detect specific markers of macrophage
   differentiation, providing validation of the impact of cilostazol on
   monocyte differentiation into macrophages. The THP-1 acute monocytic
   leukemia cell line was cultured in RPMI1640 growth medium with 10%
   Fetal Bovine Serum (FBS, Thermo-Fisher Scientific) and treated with
   DMSO and cilostazol at a concentration of 30 μM for 24 h. The specific
   markers CD80-APC (Elabscience Biotechnology Co., Ltd., Wuhan, Hubei,
   China, Cat.: E-AB-F1232E), CD163-BV421 antibody (BioLegend, Inc., San
   Diego, CA, USA, Cat.: 333612), and CD206- PE/Dazzle™ 594 antibody
   (Biolegend, Cat.: 321130) were targeted to evaluate macrophage
   differentiation. After the incubation period, the cells were analyzed
   by flow cytometry using the Attune™ NxT Flow Cytometry (Thermo Fisher
   Scientific). to assess the fluorescence intensity of the markers. The
   variation in fluorescence intensity was utilized to quantify the
   expression levels of the specific differentiation markers, thereby
   evaluating the effect of cilostazol on the differentiation process.

4. Results

4.1. Exploring the impact of cilostazol on peripheral blood mononuclear cells

   The present study investigated the effects of cilostazol on PBMCs using
   RNA sequencing ([75]Fig. 1A). PBMCs from 10 healthy donors were treated
   ex vivo with vehicle control (DMSO) or cilostazol. The transcriptomic
   profile of the treated PBMCs was examined in detail by RNA sequencing.
   This high-dimensional transcriptomic data was subsequently analyzed
   using the edgeR package, specifically designed for differential
   expression analysis. This provided a list of genes that showed
   significant changes in expression upon cilostazol treatment. Gene
   Ontology (GO) analysis was used to interpret the differentially
   expressed genes, which covered three domains: biological processes,
   molecular functions, and cellular components, providing a comprehensive
   understanding of the biological context of the genes. Pathway analysis
   was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG),
   Reactome, and WikiPathway databases to gain insight into the complex
   network of biological pathways involving these genes. The Cytoscape
   ClueGO plugin, which provides interactive visual exploration of
   functionally grouped GO and pathway annotation networks, was used to
   present a clear picture of the relationships between enriched terms. In
   addition, xCell, a digital cytometry tool, was used to derive
   cell-type-specific signatures from the gene expression data, providing
   a comprehensive landscape of the cellular composition of the samples
   and helping to identify potential shifts in cell populations upon
   cilostazol treatment. Machine learning approaches were used to
   construct predictive models of treatment response based on gene
   expression profiles. Finally, GSEA focusing on endothelial cell
   function and macrophage differentiation was performed to determine
   whether predefined gene sets showed significant differences between the
   two biological states, thus improving our understanding of the effects
   of cilostazol on PBMCs.

Fig. 1.

   [76]Fig. 1
   [77]Open in a new tab

   Analysis of cilostazol's effects on PBMCs. A) Schematic diagram
   illustrating the study design. PBMCs from healthy donors were treated
   ex vivo with either vehicle (DMSO) or cilostazol, followed by RNA
   sequencing. Subsequent analyses involved differential gene expression
   analysis, Gene Ontology (GO) interpretation, pathway analysis, digital
   cytometry, machine learning approaches, and Gene Set Enrichment
   Analysis (GSEA). B) Volcano plot visualizing the differential gene
   expression between cilostazol and DMSO-treated PBMCs. The plot
   represents the -log10 transformed adjusted p-values (y-axis) against
   the log2 fold change (x-axis). Upregulated genes (548) are depicted in
   red and downregulated genes (632) in blue. C) Heatmap illustrating the
   expression levels of the top 100 genes with the smallest false
   discovery rate (FDR) upon cilostazol treatment. Each row corresponds to
   a gene, and each column represents a sample. The genes (rows) and
   samples (columns) were hierarchically clustered based on their
   expression profiles. The clustering process involved calculating the
   Euclidean distance between each pair of genes or samples and then
   grouping them using the average linkage method. This resulted in a
   dendrogram where closely related genes or samples are placed near each
   other. The color intensity in each cell reflects the expression level
   of the gene in that sample, with red indicating upregulation and green
   indicating downregulation. (For interpretation of the references to