Introduction

Infantile hemangioma (IH), the most common benign vascular neoplasm in infants, affects approximately 5 to 10% of infants by the end of their first year.1 This condition exhibits a higher occurrence in female infants, premature babies, those with low birth weights, and in instances where gestation is marked by placental abnormalities.2 The tumor is characterized by a distinctive life cycle consisting of three stages: the proliferative phase, the involuting phase, and the involuted phase.3 Proliferation occurs during early infancy, specifically the first 4–6 months, while gradual spontaneous involution or regression starts by 1 year of age.4 The involuted phase marks the cessation of the tumor growth, leaving behind scar tissue, telangiectasia, or skin conditions such as redundant or anetodermic skin.5 While the majority of IH cases (approximately 85%) undergo spontaneous regression, a subset necessitates medical intervention based on the tumor’s anatomical location and associated risks.6 Previously, systemic corticosteroids served as the primary therapy for IH. However, the contemporary treatment of choice involves the use of the nonselective beta-adrenergic receptor (ADRB) antagonist, propranolol, albeit its exact mechanism of action is still unclear. The side effects of propranolol include bradycardia, hypotension, hypoglycemia, and disturbed sleep.7 Further, there is a potential for propranolol resistance. Consequently, understanding the molecular mechanisms underlying IH and discovering novel treatments remains an area of significant importance.

The identification of novel therapeutic targets is a critical component of drug development, particularly for tumors. Recent biotechnological innovations offer transformative potential for elucidating disease mechanisms, identifying novel drug targets, and propelling a paradigm shift in pharmaceutical research. The Connectivity Map (CMap) represents one of the first publicly accessible resources that provide a comprehensive record of the transcriptional responses of human cells in response to chemical and genetic perturbations.8 This tool quantifies the connections between diseases, genes, and drugs, thus identifying drugs that can reverse the impact on a large number of genes affected by the disease. In contrast, the Comparative Toxicogenomics Database (CTD) operates as a digital ecosystem that bridges toxicological data on substances, genes, phenotypes, and diseases. Unlike CMap, CTD focuses on the toxicological aspects of various compounds and their impact on genes and diseases.9 Both these robust databases were utilized in our research, contributing significantly to our study’s findings and conclusions.

To elucidate the pathogenic underpinnings of IH and identify putative therapeutic agents, we leveraged the GSE127487 dataset10 from the Gene Expression Omnibus (GEO). We subsequently reanalyzed this dataset using the bioinformatics tools referenced in the following sections. Our analysis involved several steps, beginning with an exploration of differentially expressed genes (DEGs). Next, we conducted an enrichment analysis and implemented a weighted gene co-expression network analysis (WGCNA).11 Construction of a protein-protein interaction (PPI) network revealed high-connectivity hub genes, which were queried against the CMap and CTD databases to explore potential novel therapeutic agents for IH. Furthermore, to validate our computational predictions, we also conducted a series of in vitro experiments. A detailed outline of the study design is presented in Figure 1. This research followed protocols approved by the Review Board of Plastic Surgery Hospital, Chinese Academy of Medical Sciences.

Figure 1 A flowchart depicting the schematics of the workflow pipeline.

Abbreviations: GEO, Gene Expression Omnibus; PCA, Principal Component Analysis; WGCNA, Weighted Gene Coexpression Network Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein‒protein interaction; CMap, Connectivity Map; CTD, Comparative Toxicogenomics Database.

Materials and Methods Reagents and Cell Culture

The drugs entinostat (#HY-12163), sorafenib (#HY-10201), dasatinib (#HY-10181) and sirolimus (#HY-10219) were obtained from MedChem Express (Monmouth Junction, NJ, USA). Hemangioma Endothelial Cells (HemECs) were purchased from Shanghai Honsun Biological Technology Co., Ltd. (Shanghai, China). Cells were cultured in DMEM (#SH30022.01B, HyClone, USA), fortified with 10% fetal bovine serum (FBS; #16000-044, Gibco, USA) and 1% penicillin-streptomycin (#15140-122, Gibco, USA) at 37 °C in a humidified atmosphere containing 5% CO2. The HemECs were digested with 0.25% trypsin and 0.02% EDTA solution. Subsequent propagation occurred when the cell density achieved an 80% confluence. For our varied experimental needs, cells from passages 3 through 6 were selected.

RNA Extraction and qRT-PCR

To verify that the cells obtained were HemECs, we assessed the expression levels of Factor VIII and CD31 in the cells via quantitative Real-Time PCR (qRT-PCR). We extracted the total RNA from HemECs using the Trizol reagent (Invitrogen, CA, USA). Total RNA (1μg) was used to synthesize cDNA using Reverse Transcription kit (#CW2569, CWBIO, China). We performed the qPCR using Power SYBR Green PCR Master Mix (Vazyme, China) on the ABI StepOne Plus real-time PCR system (Applied Biosystems, USA). The relative expression levels were determined using the 2−ΔΔCt method. The expression levels were referenced against Human Umbilical Vein Endothelial Cells (HUVEC) and were normalized considering the expression of GAPDH.

Cell Proliferation Assay

Cell proliferation was assessed using a Cell Counting Kit-8 (CCK-8; Beyotime Biotechnology, China), with the procedure being strictly based on the manufacturer’s instructions. Cells (3×103 cells/well) were plated in 96-well plates and subjected to various concentrations of entinostat, sorafenib, dasatinib, and sirolimus at multiple time points. Subsequently, 10μL of kit reagent, WST-8 [2-(2-methoxy-4-nitrophenyl)-3- (4-nitrophenyl)-5-(2,4-disulfonyl)-2H-tetrazolium], was added to each well, and the plate was incubated for 3 h at 37 °C. The absorbance at a wavelength of 450 nm was measured post-incubation, utilizing a Microplate Reader (Thermo Scientific, USA). Cell viability was determined as a percentage, comparing the number of viable cells in the compound-treated groups against those in the untreated control group. To ensure the reliability of the data, each experiment was independently executed five times.

Flow Cytometry

To assess the induction of apoptosis in HemECs by entinostat, sorafenib, dasatinib, and sirolimus, we employed an Annexin V-FITC/PI double-staining analysis. Following a 72-hour treatment with specified concentrations of entinostat (10μM), sorafenib (10μM), dasatinib (0.1μM), and sirolimus (10μM), we collected and washed the cells with preheated PBS. In the control group, 0.1% DMSO was included as a vehicle control. We evaluated the presence of phosphatidylserine on the surface of apoptotic cells using the Annexin V-PI Apoptosis Detection Kit (#C1065M, Beyotime Biotechnology, China) and a flow cytometer (BD Accuri C6 Plus, US). Data analysis was performed using FlowJo software. Each cell experiment was conducted with five replicates to ensure the reliability and reproducibility of our results.

Transwell Migration Assay

We assessed the migratory potential of HemECs employing a Transwell system (6.5 mm, 0.4 µm pore polycarbonate membrane; #3413, Corning Cat, USA). HemECs (2.5 x 105) were plated in the upper chamber using DMEM supplemented with 0.5% BSA (#10735078001, Roche, Swiss). We filled the lower chamber with the standard culture medium, which comprised DMEM and 10% FBS. Subsequently, we added either 10μM of entinostat, 10μM of sorafenib, 0.1μM of dasatinib, or 10μM of sirolimus to both the upper and lower chambers, and proceeded with incubation for 48 hours. Thereafter, cells that had not migrated or invaded were fixed on the upper side of the filter with methanol for 20 minutes and stained with a 0.1% crystal violet solution for 20 minutes. We then imaged and quantified the remaining cells on the upper chamber using an inverted microscope (Nikon Eclipse Ci-S) with 200 × objective lens. Five replicates were used for each cell experiment.

Tube Formation Assay

Matrigel (Corning Biocoat, Cat #354243) was prepared and stored at 4°C overnight to ensure proper gelation. On the day of the experiment, the Matrigel was thawed, and the 96-well plates, along with pipette tips, were pre-cooled at 4°C for 30 minutes to maintain the Matrigel in a liquid state during handling. A total of 50 μL of Matrigel (10 mg/mL) was added to each well of the pre-cooled 96-well plates and incubated at 37°C for 30 minutes to allow the Matrigel to solidify. Following the incubation, the tube formation assay was performed. The formed tubular structures were captured using a microscope, and the extent of tube formation was quantified using ImageJ software. Each experimental condition was tested in triplicate to ensure reproducibility and accuracy of the results.

Data Acquisition and Preprocessing

To uncover the transcriptomic differences in IH patients, the GSE127487 (Stone A., 2019)10 dataset in the GEO database (https://www.ncbi.nlm.nih.gov/geo/) was obtained, and the sample platform was GPL10558 (Illumina HumanHT-12 V4.0). Our study involved a total of 28 samples, which were divided across five groups with the respective sample sizes of 5, 6, 6, 6, and 5. We proceeded to rectify the background and normalized the data using the normalizeBetweenArrays algorithm within the limma package to construct a mRNA expression matrix.12 A correlation diagram of mRNA expression levels between samples was applied to affirm the reliability of the experiment and the accuracy of sample selection. We employed principal component analysis (PCA) to assess intergroup variances and intragroup sample duplications. Based on the PCA results, we identified and excluded abnormal samples from further analysis.

Differential Gene Analysis and Enrichment Analysis

The R package limma was used to identify DEGs with a cutoff value of adjusted p value < 0.05 and |log2−fold change| > 1.12 The visualization of DEGs was accomplished through volcano and heatmap plots by the ggplot213 and pheatmap14 packages in R language, respectively.

We further explored the biological functions of the shared DEGs using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.15,16 GO provides a comprehensive gene function classification system, dividing genes into three categories: biological process (BP), cellular component (CC), and molecular function (MF). KEGG, on the other hand, facilitates a systematic visual analysis of gene functions and metabolic pathways. Both GO and KEGG pathway analyses were performed using the clusterProfiler package.17 Significant GO and KEGG terms were identified using thresholds of q-value < 0.05 and p-value < 0.05. We also generated networks within the Cytoscape software environment using the Enrichment Map plugin, with an uncorrected p-value threshold of 0.005, a false discovery rate (FDR) cutoff of 0.1, and an overlap coefficient threshold of 0.1.18,19 We employed Gene-Set Enrichment Analysis (GSEA) to evaluate gene expression matrix files, enabling us to retain more information than conventional GO and KEGG enrichment analysis and to avoid potential loss of relevant genes.20 To validate potential cellular functions and signaling pathways, the clusterProfiler R package was employed to perform GSEA, with c2.all.v7.5.entrez.gmt set as background (https://www.gsea-msigdb.org/gsea/msigdb).

Weighted Gene Coexpression Network Analysis (WGCNA) Analysis and Enrichment Analysis

We conducted Weighted Gene Coexpression Network Analysis using the WGCNA package in R.11 WGCNA is a widely used technique for identifying correlations between genes and pinpointing characteristic and intra-modular hub genes. This analysis also assists in estimating measurement values for topological properties and module membership. The soft threshold was selected using the function pick soft threshold, a criterion based on the approximate scale-free network that allows the constructed network to be more consistent with the power-law distribution. For this study, we defined modules using the following module-cutting parameters: height = 0.25, deepSplit = 2, and minModuleSize = 30. Similar gene modules were merged with a threshold of 0.25. We conducted GO and KEGG analyses of the tissue type-related module using the clusterProfiler package and Cytoscape software, as detailed earlier.

Enrichment Analysis and Protein‒Protein Interaction Network Analysis of the Overlapping mRNAs of WGCNA and Differential Gene Analysis

We selected the overlapping mRNAs from WGCNA and differential gene analyses. Enrichment analysis was performed using the clusterProfiler, ggplot2, and GOplot packages.21

For the identification of key genes in IH development, we performed a protein-protein interaction (PPI) analysis for these intersecting genes. This was facilitated using an online database known as Search Tool for the Retrieval of Interacting Genes (STRING), which examines the interconnectivity between the genes.22 DEGs with a minimum required confidence score ≥0.7 were selected to build a comprehensive network model, which was subsequently visualized in Cytoscape software.18

Module Analysis and Hub Gene Identification

As outlined in Figure 1, to further dissect the underlying molecular mechanisms of IH, we employed the MCODE plugin, a molecular complex detection tool integrated into Cytoscape, to identify key modules within our dataset.23 The operational parameters were set as follows: Degree Cutoff = 2, Node Score Cutoff = 0.2, K-Core = 2, Max.Depth = 100. The most significant and largest module was defined by an MCODE score > 4. This module analysis not only provided validation for our bioinformatics approach but also highlighted additional biological insights that guided the subsequent steps of our study.

With the aid of the cytoHubba plugin in Cytoscape, we used the MCC prediction algorithm to identify the top 100 hub genes.24 Extended Polydimensional Immunome Characterization (EPIC),25 the Microenvironment Cell Population-Counter (MCP-counter),26 Tumor Immune Dysfunction and Exclusion (TIDE),27 and gene signature enrichment–based xCell algorithm28 were used to investigate the correlation between gene expression and cancer-associated fibroblast infiltration of the top 3 genes. Furthermore, we explored the gene expression of the top 3 hub genes in 33 different cancer types by using the gene module in TIMER2.0 (http://timer.cistrome.org/).

Drug Identification

To identify potential small-molecule therapeutic candidates, we uploaded the Infantile Hemangioma (IH) gene signature into the Connectivity Map database (http://www.broadinstitute.org/cmap/). Enrichment scores were calculated within a range from −1 to +1, symbolizing the similarity between the gene profile of IH tissue and various drugs. A positive raw connectivity score (raw_cs) value, closer to +1, suggested that a small molecule could induce similar mRNA expression in IH cells, whereas a negative raw_cs value, closer to −1, suggested that the small molecule could inhibit the progression of IH. Additionally, we searched the CTD database for the top 100 hub genes and searched for chemicals associated with infantile hemangioma. Finally, we visualized the intersections between the IH targets and the predicted top drugs utilizing the R sankeywheel package.

Results Data Acquisition and Preprocessing

The GSE127487 dataset was normalized to obtain a mRNA expression matrix comprising five controls (Nor group), six 6-month infant IH samples (IH6 group), six 12-month infant IH samples (IH12 group), six 24-month infant IH samples (IH24 group) and five IH samples from patients treated with propranolol (Prop group) (Figure 2A and B). Systematic and dataset-specific bias was markedly diminished after preprocessing. Principal component analysis was employed to visualize the disparity between the IH and control groups before (Figure 2C) and after data normalization (Figure 2D). Upon normalization, samples from distinct groups demonstrated notable differences as depicted in Figure 2D, while samples within the same group (except for the GSM3635438 samples) displayed a high degree of uniformity. According to Figure 2D, the mRNA signature of the IH24 group most closely resembled that of healthy tissues, followed by the IH12 group, Prop group, and IH6 group. It is recognized that IH progresses through three unique stages. However, determining a clear boundary between these stages can be complex due to various influencing factors such as the lesion’s shape, location, size, depth, and genetic factors. Considering the likelihood of GSM3635438 being an involuted lesion, we elected to exclude it from our study.

Figure 2 Data preprocessing. (A and B) Standardization of gene expression data before (A) and after standardization (B). *Indicates the outlier detection by boxplots. (C and D) Principal component analysis of the compendium of IH tissue transcriptomics before (C) and after normalization (D).

Differential Gene Analysis and Enrichment Analysis

Figure 2D indicates that the genetic profile of the IH6 group is the most distinct from that of the healthy tissues, suggesting a more accurate differential gene analysis between the IH6 group and the Nor group. A total of 1203 DEGs were extracted, including 802 upregulated genes and 401 downregulated genes. A heatmap was constructed using the pheatmap R package, providing a visual representation of the top 100 DEGs between the IH6 and Nor samples (Figure 3A). In addition, a volcano plot was employed to display the cluster of DEGs (Figure 3B). Comprehensive information regarding the 100 most significant DEGs can be found in Table S1 and Figure S1.

Figure 3 Differential gene analysis and enrichment analysis. (A) Heatmaps showing gene expression values for the 100 most significantly differentially expressed genes (DEGs) between the IH6 and Nor groups. (B) Volcano plot showing the distribution of log2−fold change and − log10−P value of all quantified transcripts between the IH6 group and the normal group. Blue circles: adj.p.value <0.05, log2−fold change < −1; red circles: adj.p.value <0.05, log2−fold change >1. (C and D) GO (C) and KEGG (D) enrichment analyses of the top 20 most significantly enriched molecular functions of DEGs. (E) Enrichment analysis identifying over 500 significantly enriched gene functions that were clustered using EnrichmentMap and AutoAnnotate in Cytoscape to identify the key gene functions. Nodes represent individual GO terms, with size relating to the log2−fold change value in each term and the color indicating the functional category. Edges represent connections between nodes that share genes. (F−G) GSEA identifying the PI3K/AKT/MTOR signaling pathway (F) and CGMP/PKG signaling pathway (G) as significant. Black vertical lines denote the positions of each gene from the PI3K/AKT/MTOR or CGMP/PKG pathway gene sets within the ranked list of genes. The green curve represents the Enrichment Score (ES) plot, illustrating the cumulative ES value across the differential gene ranking. The heatmap depicts genes highly expressed in the IH6 group in red, and those highly expressed in the control group in blue. The grey area plot indicates the signal-to-noise ratio for each gene, with genes above the zero cross positively correlated with the pathway, and those below negatively correlated.

Enrichment analysis identified over 500 significantly enriched gene functions (Table S2). The most significant GO enrichment consequences of DEGs were focused on epithelial cell proliferation, actin filament organization, cell-substrate adhesion, and extracellular matrix organization (Figure 3C). Further, these genes were enriched in the PI3K/AKT/MTOR signaling pathway, RAP1 signaling pathway, and CGMP/PKG signaling pathway (Figure 3D). To facilitate a more straightforward interpretation of functional enrichment, a network of GO terms was created using EnrichmentMap. DEGs were significantly enriched in ‘actin filament’, ‘platelet aggregation’, and ‘cell leukocyte adhesion’ (Figure 3E). In addition, GSEA was performed to further validate that the PI3K/AKT/MTOR signaling pathway (Figure 3F) and CGMP/PKG signaling pathway (Figure 3G) were differentially enriched among the DEGs.

Weighted Gene Coexpression Network Analysis and Enrichment Analysis

WGCNA was employed to construct a co-expression network and discern co-expression modules. The soft threshold was validated sufficiently, converging to a scale-free topology with a value of 7 (Figure S2A). The Dynamic Tree Cut algorithm provided the initial set of modules, after which related modules were merged. Consequently, a total of 22 modules were identified (Figure 4A). Further analysis was conducted to examine the correlation between each module and the sample type (tissue and treatment) (Figure 4B). The blue module, comprising 1780 mRNAs, exhibited the strongest correlation with infantile hemangioma (r = 0.81) (Figure S2B).

Figure 4 WGCNA and enrichment analysis. (A) Hierarchical cluster analysis detecting coexpression clusters with corresponding color assignments. Each color represents a module in the constructed gene coexpression network by WGCNA. (B) Module–trait relationships between 22 modules and clinical phenotypes. Each row represents a color module, and every column represents a clinical trait. Each cell contains the corresponding correlation and p value. (C and D) GO (C) and KEGG enrichment (D) analyses of the top 20 most significantly enriched molecular functions of the blue module. (E) Enrichment analysis identifying over 500 significantly enriched gene functions that were clustered using EnrichmentMap and AutoAnnotate in Cytoscape to identify the key gene functions. Nodes represent individual GO terms, with size relating to the log2−fold change value in each term and the color indicating the functional category. Edges represent connections between nodes that share genes. (F−I) GSEA identifying the PI3K/AKT/MTOR signaling pathway (F), CGMP/PKG signaling pathway (G), RAP1 signaling pathway (H) and RAS signaling pathway (I) as significant. Black vertical lines denote the positions of each gene from the pathway gene sets within the ranked list of genes. The green curve represents the Enrichment Score (ES) plot, illustrating the cumulative ES value across the differential gene ranking. The heatmap depicts genes highly expressed in the IH6 group in red, and those highly expressed in the control group in blue. The grey area plot indicates the signal-to-noise ratio for each gene, with genes above the zero cross positively correlated with the pathway, and those below negatively correlated.

We used these 1780 mRNAs for further gene enrichment analysis. It demonstrated that the blue module was concentrated on mesenchyme-related processes, including cell-substrate adhesion, regulation of angiogenesis, regulation of vasculature development, and mesenchyme development (Figure 4C and Table S3). KEGG pathway analysis showed that these genes were mainly enriched in the PI3K/AKT/MTOR signaling pathway, RAS signaling pathway, RAP1 signaling pathway, and CGMP/PKG signaling pathway (Figure 4D). To facilitate understanding, a detailed network of GO terms was generated using EnrichmentMap (Figure 4E). GSEA further highlighted the importance of PI3K/AKT/mTOR (Figure 4F), CGMP/PKG (Figure 4G), RAP1 (Figure 4H), and the RAS signaling pathway (Figure 4I).

Enrichment Analysis and Protein–Protein Interaction Network Analysis of the Overlapping Genes of WGCNA and Differential Gene Analysis

A total of 674 overlapping genes between WGCNA and differential gene analysis were identified. GO annotation analysis indicated that these genes participated in various processes, including angiogenesis and cell−substrate adhesion (Figure 5A). Additionally, they were involved in actin binding and functioning as extracellular matrix structural constituents. The cnetplot of the five most statistically significant biological processes (Figure 5D), cellular components (Figure 5E), and molecular function GO terms (Figure 5F) illustrated the correlation between differentially expressed genes (DEGs) and GO terms. KEGG enrichment analysis suggested that the overlapping genes were significantly enriched in several KEGG pathway terms, including the RAP1 and PI3K/AKT/MTOR, RAS/MAPK, and CGMP/PKG signaling pathways (Figure 5B). The PI3K/AKT/MTOR signaling pathway was further confirmed by GSEA (Figure 5C). In the PPI network analysis, interactions between proteins expressed from the overlapping genes were observed, comprising 273 genes as illustrated in Figure 6.

Figure 5 Enrichment analysis of the overlapping genes of WGCNA and differential gene analysis. (A) GOBubble plot of GO enrichment annotation containing biological process (red), cellular component (blue), and molecular function (green) of the gene ontology terms. The y-axis is the log−p value, and the x-axis is the Z score. The bubble size in each GO term is considered proportional to the DEG number from the dataset. (B) KEGG enrichment analyses of the top 10 most significant pathways; the inner ring is a bar plot where the bar height indicates the significance of the term (−log10 P value), and the color indicates the z score. The outer ring displays scatterplots of the expression levels (log2−fold change) for the genes in each term. (C) GSEA identifying the PI3K/AKT/MTOR signaling pathway as significant. (D−F) Cnetplot of the five most statistically significant biological process (D), cellular component (E) and molecular function (F) GO terms showed the correlation among genes and GO terms.

Figure 6 Protein–protein interaction network analysis of the overlapping genes of WGCNA and differential gene analysis. The interaction network between proteins coded by DEGs comprised 273 nodes and 379 edges in PPI analysis of the overlapping genes of DEG and WGCNA. The red and blue circle nodes represent upregulated and downregulated DEGs, respectively.

Module Analysis and Hub Gene Identification

Through MCODE analysis of the protein-protein interaction network, a total of four modules with MCODE scores exceeding 4 were obtained (Figure 7A). Module 1, with the highest score of 5, consisted of five genes, including THSD1, ANAMTS9, SBSPON, SEMA5A, and SEMA5B. Module 2, containing GNB4, GNA12, LPAR4, LPAR6, and GNG11, scored 4.5. Module 3, with a score of 4, included PDE1A, PDE1B, GUCY1A2, and GUCY1B. Module 4, also with a score of 4, comprised BST2, IFIT2, IFIT3, and ISG15.

Figure 7 Module analysis and hub gene identification. (A) Four cluster modules extracted by MCODE. (B) Correlation between CD34, FN1, and PECAM1 expression and cancer-associated fibroblast immune infiltration. (C) The expression status of CD34 in different tumor types was visualized by TIMER2 (*P < 0.05; **P < 0.01; ***P < 0.001). Red represents the tumor group, and blue represents the control group. (D) Gene expression of CD4 in different IH tissues. (E) The expression status of FN1 in different tumor types was visualized by TIMER2 (*P < 0.05; **P < 0.01; ***P < 0.001). (F) Gene expression of FN1 in different IH tissues. (G) The expression status of PECAM1 in different tumor types was visualized by TIMER2 (*P < 0.05; **P < 0.01; ***P < 0.001). (H) Gene expression of PECAM1 in different IH tissues.

Using the MCC algorithm in cytoHubba plugin, we identified a total of 100 genes as hub genes (Table S4). The top three significant genes were CD34, FN1, and PECAM1. Previous research has discovered that cancer-associated fibroblasts in the stroma play a role in the regulation of various tumor-infiltrating immune cells.29 To evaluate the relationship between cancer-associated fibroblast infiltration and CD34, FN1, and PECAM1 expression in different tumors, we used the EPIC, MCP-counter, Xcell and TIDE algorithms. CD34 and FN1 expression were positively correlated with cancer-associated fibroblast infiltration in different malignancies, while PECAM1 expression had a negative correlation (Figure 7B). The expression pattern of the top 3 genes in tumor tissues was then investigated. CD34 mRNA expression was increased in various tumor tissues compared to corresponding normal tissues (Figure 7C). Tumor tissues of cholangiocarcinoma (CHOL), glioblastoma multiforme (GBM), kidney renal clear cell carcinoma (KIRC), and liver hepatocellular carcinoma (LIHC) had significantly higher CD34 expression than corresponding normal tissues (Figure 7C). Meanwhile, significantly decreased CD34 expression was observed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), colon adenocarcinoma (COAD), kidney chromophobe (KICH), kidney renal papillary cell carcinoma (KIRP), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ) and uterine corpus endometrial carcinoma (UCEC) tumor tissues (Figure 7C). Furthermore, when compared to that in primary SKCM tumor tissues, CD34 expression was also significantly elevated in metastatic SKCM tissues (p < 0.001) (Figure 7C). Figure 7D illustrates the mRNA expression in IH tissues and normal tissues. Figure 7E–H demonstrates the mRNA expression of FN1 and PECAM1 in 33 different types of tumors and IH tissues. To further validate these findings specifically in the context of infantile hemangioma, we performed qPCR analysis to assess the expression levels of CD34, FN1, and PECAM1 in HemECs. The results confirmed that these genes are significantly expressed in HemECs, as shown in Figure S3.

Identification of Drugs

Chemicals with opposing mRNA signatures were identified when compared to compounds from the CMap database, which contains over 1 million gene expression signatures sourced from the L1000 high-throughput assay. Propranolol stands as the gold standard treatment for IH. Thus, we utilized its raw_cs value (−0.3827) as the significant compound cut-off. In total, 4231 drugs were extracted according to the CMap database. The 100 hub genes were uploaded into the Comparative Toxicogenomics Database, and the specific medications that could reduce or increase the expression levels of these genes were chosen. A total of 285 target drugs targeting at least five genes were predicted. Ultimately, we identified 48 intersecting drugs (aflatoxin B1, amiodarone, belinostat, calcitriol, capsaicin, catechin, chir 99021, curcumin, cytarabine, dasatinib, dexamethasone, diclofenac, disulfiram, doxorubicin, entinostat, etoposide, geldanamycin, genistein, hydrocortisone, irinotecan, ivermectin, methotrexate, paclitaxel, panobinostat, pioglitazone, piroxicam, progesterone, pyrazolanthrone, resveratrol, rotenone, simvastatin, sirolimus, sorafenib, sulforaphane, sunitinib, tamibarotene, tamoxifen, temozolomide, testosterone, thalidomide, topotecan, troglitazone, U 0126, valproic acid, vincristine, vorinostat, wortmannin, and zidovudine) via the CMap and CTD databases, suggesting their potential to inhibit the expression of the IH signature by targeting the top 100 hub genes. Upon meticulous examination, specific compounds were ruled out based on their inherent properties that could compromise the safety of potential treatments. For example, aflatoxin B1 was excluded for its strong carcinogenic property. Similarly, rotenone, a known inhibitor of mitochondrial complex I, was removed.30 Geldanamycin (the first-discovered HSP90 inhibitor) and troglitazone (an oral antihyperglycemic agent) were also ruled out due to their high hepatotoxicity.31,32 Ultimately, 44 drugs were identified through the CMap and CTD databases (Figure 8B). Remarkably, most of these candidate compounds have been utilized in cancer treatments (Table 1). Among them, corticosteroids and sirolimus have already been used for treating complicated IH in patients intolerant to propranolol, further validating the credibility of these drug candidates for IH.33 The Sankey diagram revealed the correlation between the top five drugs and their corresponding targets (Figure 8A).

Table 1 Forty-Four Significant Chemicals for IH Based on Genes with Enriched Behavior Features

Figure 8 Identification of drugs. (A) Sankey diagram revealing the correlation between the drugs and targets. The nodes represent drug treatments or genes. Each node signifies a distinct entity within the analyzed system. The flows in the diagram correspond to specific transitions or relationships between these entities. The width of each flow quantitatively indicates the magnitude of the respective transition. (B) The molecular structure of 44 small component drugs targeting hub genes in IH by the CMAP and CTD databases.

The Inhibitory Effect of Entinostat, Sorafenib, Dasatinib, and Sirolimus on the Proliferation, Apoptosis, Migration and Angiogenesis of HemECs

In our investigation, we identified a total of twelve molecular targeting agents which fall into four primary modes of action: PI3K/AKT/MTOR pathway inhibitors, represented by key drugs like wortmannin and sirolimus; MAPK pathway inhibitors, predominantly featuring U0126 and sorafenib; multi-target inhibitors such as dasatinib, sorafenib, and sunitinib; and histone deacetylase (HDAC) inhibitors, with primary agents including belinostat, entinostat, panobinostat, and vorinostat. To validate the effectiveness of these classes, we performed experiments utilizing representative drugs from each category, specifically entinostat (HDAC inhibitor), sorafenib (MAPK pathway inhibitor and multi-target inhibitor), dasatinib (multi-target inhibitor) and sirolimus (PI3K/AKT/MTOR pathway inhibitor).

We have acquired mature HemECs cell lines, which were utilized for culturing from commercial sources, ensuring high validity of our experiments. The identity and quality of these HemECs were confirmed through qPCR analysis and flow cytometry with factor VIII (Figure 9A) and CD31 (Figure 9B and Figure S4) markers. Subsequently, we assessed the effects of the four drugs on the proliferation, migration, and apoptosis of HemECs.

Figure 9 Effects of entinostat, sorafenib, dasatinib, and sirolimus on the apoptosis and migration of HemEC cells. (A and B) Expression of CD31 (A) and factor VIII (B) in HUVEC and HemEC cells examined using qRT-PCR. (C) Flow cytometry showed that the employment of entinostat, sorafenib, and sirolimus increased the rate of cell apoptosis in HemECs. (D) The invasion and migration abilities of HemEC cells treated with entinostat, sorafenib, dasatinib and sirolimus were detected using Transwell assay. (E) Statistical analysis of the number of cells that had not migrated or invaded based on the Transwell invasion assay. (F) The tube formation of HemECs treated with entinostat, sorafenib, dasatinib and sirolimus. (G) Statistical analysis of the total lengths of tubes using ImageJ software. *P < 0.05, **P< 0.01, ***P < 0.001, ****P < 0.0001.

To determine if entinostat-induced cytotoxicity involved apoptosis, HemECs subjected to entinostat were stained using Annexin-V/FITC and analyzed via flow cytometry (Figure 9C). We observed an increase in apoptosis to 8.01%. Similarly, sorafenib and sirolimus exposure led to increased apoptosis rates of 12.60% and 3.81%, respectively, after 72 hours. However, upon treatment of HemECs with dasatinib, there was no significant change in the proportion of apoptotic cells compared to the control group, indicating its potential non-apoptotic mode of action. We also investigated the influence of these drugs on cellular migration and invasion using the Transwell system. Our findings revealed that fewer cells remained in the upper chamber for the entinostat, sorafenib, dasatinib, and sirolimus treatment groups compared to the control group (Figure 9D). The utilization of entinostat, sorafenib, dasatinib, and sirolimus significantly inhibited both the migration and invasion capacities of the HemECs (P < 0.001) (Figure 9E). To evaluate the tube-forming effects of entinostat, sorafenib, dasatinib, and sirolimus, we conducted Matrigel tube formation assays by adding these drugs to the HemEC model. Comparative analysis revealed a significant reduction in tubule formation following treatment with entinostat, sorafenib, dasatinib, and sirolimus, as illustrated in Figure 9F and G.

The cells were treated with varying concentrations of these drugs for durations of 0, 24, 48, or 72 hours, and cell viability was assessed using the CCK-8 assay. The initial concentrations were selected based on previously published literature, and the CCK-8 assay was then used to further validate these concentrations. Our data suggested that the proliferative capacity of the cells in the entinostat (Figure 10A), sorafenib (Figure 10B), dasatinib (Figure 10C), and sirolimus (Figure 10D) treatment groups was reduced compared to the control group. The CCK-8 assay further revealed that entinostat (Figure 10E), sorafenib (Figure 10F), dasatinib (Figure 10G), and sirolimus (Figure 10H) inhibited cell proliferation in both dose- and time-dependent manners.

Figure 10 Effects of entinostat, sorafenib, dasatinib, and sirolimus on the proliferation of HemEC cells. (A−D) The proliferative capacity of HemECs treated with entinostat (A), sorafenib (B), dasatinib (C), and sirolimus (D) determined by CCK-8 assay; OD values represent the cell viability; higher OD value reflects stronger cell proliferation ability (*P < 0.05, ***P < 0.001, ****P < 0.0001). (E−H) The cell viability of HemEC was observed using the CCK-8 assay with entinostat (E), sorafenib (F), dasatinib (G), and sirolimus (H).

Discussion

Infantile hemangiomas, which have a prevalence of 5 to 10%,1 are the most prevalent tumors in infancy.2 The study of the transcriptomic etiology and pathogenesis of IH has advanced with the development of high-throughput transcriptomic technology. Unfortunately, research on the identification of target drugs remains limited. Therefore, it is crucial to fully utilize bioinformatics analytic methods while performing in-depth data mining on the IH dataset to obtain valuable transcriptomic information. We searched the GEO database to obtain the published human IH dataset (GSE127487) and performed differential expression analysis, weighted gene coexpression network analysis, and protein-protein network analysis to screen the main genes which most likely to be responsible for the occurrence of IH. A total of 1203 DEGs were found, of which 674 genes were screened through WGCNA. We performed further PPI analysis on the overlapping genes, and the top 100 hub genes of IH were identified. Pan cancer analysis revealed that these genes were closely associated with malignancies and immune infiltration, indicating that these identified molecules have the potential to be used as disease biomarkers for patient-specific treatment and to improve the diagnosis, treatment, and prevention of IH. In this study, CTD and CMap database were employed to identify a total of 44 potential drugs for IH based on the 100 hub genes. Out of the twelve predicted molecular-targeting drugs, four (entinostat, sorafenib, dasatinib, and sirolimus) inhibited the proliferation and migration and induced apoptosis in HemEC cells.

The findings of our investigation were consistent with the clinical course and first-line therapy of IH. The majority of IH growth occurs before 12 weeks of age; after that, it slows down and typically comes to a halt between 4 and 6 months.33 Between the second and sixth years of childhood, 85% to 90% of IHs experience spontaneous involution.34 Our PCA analysis of IH tissues from various time points revealed that 6-month IH tissues exhibited a mRNA signature that was considerably different from normal tissues, corresponding to the clinical course of IH. Furthermore, it has been illustrated that hypoxia-induced proliferation plays an important role in the pathogenesis of IH.35 GO enrichment showed that these genes were significantly enriched in mesenchyme-related functions, such as a response to decreased oxygen levels, mesenchyme development, and regulation of angiogenesis (Figure 4C). Vascular anomalies are believed to originate and progress as a result of excessive cross-activation of the PI3K/AKT/MTOR and RAS/MAPK pathways.36 KEGG enrichment and GSEA analysis of our research further emphasized the core effect of the PI3K/AKT/MTOR and RAS/MAPK signaling pathways (Figure 4D). Additionally, accumulating evidence has suggested that the main risk factors for the development of IH are female sex, low birth weight, and prematurity.2 Girls are affected 2.3–2.9 times more often than boys.2 In line with the gender-specific characteristic of IH, three of the 44 possible medications work by modulating sexual hormones (progesterone, progesterone receptor agonist; tamoxifen, estrogen receptor antagonist; testosterone, androgen receptor agonist). In addition, another three drugs (dexamethasone, glucocorticoid receptor agonist; hydrocortisone, glucocorticoid receptor agonist; sirolimus, mTOR inhibitor) were included as potentially useful systemic treatment options for IH in the recommendations of the European expert group.33 All these findings emphasized the accuracy of our study.

In the targeted drug identification process, we screened medicines based on the raw_cs value of propranolol, which is currently the first choice for infantile hemangioma. Forty-four drugs were identified after drug toxicity screening, including 9 frontline cytotoxic drugs (cytarabine, doxorubicin, etoposide, irinotecan, methotrexate, paclitaxel, temozolomide, topotecan, and vincristine), 12 molecular targeting agents (belinostat, chir 99021, dasatinib, entinostat, panobinostat, sirolimus, sorafenib, sunitinib, thalidomide, U 0126, vorinostat, and wortmannin), 6 natural extracts (capsaicin, catechin, curcumin, genistein, resveratrol, and sulforaphane), 4 glucocorticoids or nonsteroidal anti-inflammatory medications (dexamethasone, diclofenac, hydrocortisone, and piroxicam), 3 sex hormone-related drugs (progesterone, tamoxifen, and testosterone), and 10 other drugs (amiodarone, calcitriol, disulfiram, ivermectin, pioglitazone, pyrazolanthrone, simvastatin, tamibarotene, valproic acid, and zidovudine). Of 44 predicted drugs, glucocorticoids and sirolimus were used as second-line treatments for individuals with complex IH who did not react to propranolol or who exhibited primary contraindications or developed side effects. Despite being second to propranolol in clinical practice, glucocorticoids and sirolimus displayed a more opposing transcriptome signature alteration than propranolol when added to cell models. Systemic corticosteroid therapy, however, was not recommended as a first-line treatment due to its side effects, which included cushingoid facies (71%), stomach irritability (21%), fungal infection (6%), shorter stature (35%) and weight increase (42%).37 Propranolol, the first-line therapy for problematic infantile hemangiomas, is also associated with various adverse events. However, a recent study demonstrated that atenolol, when compared with propranolol, has a significantly lower incidence of adverse events, making it a promising alternative treatment for IH patients requiring systemic therapy.38 It is worth mentioning that six drugs possessed potent anti-inflammatory effects (dexamethasone, diclofenac, hydrocortisone, piroxicam, U0126, and pioglitazone), indicating that inflammation might play a crucial role in the pathogenesis of IH. An increasing number of cytotoxic drugs are being used to treat vascular anomalies, but they are only applied topically due to their systemic side effects. Bleomycin is widely used in the treatment of hemangioma and vascular malformation due to its relatively low toxicity (no myelosuppression, thrombocytopenia, or leukopenia).39 However, pulmonary fibrosis was reported in some oncologic patients receiving a high cumulative intravenous dose of bleomycin.40 Other documented common adverse events included edema and ulceration.41 Thus, the 9 predicted cytotoxic drugs are not the most advantageous prospective treatments for IH due to systemic and topical side effects.

We identified a total of twelve molecular targeting agents, which are probably the most promising medications since they cause fewer adverse drug reactions. Wortmannin and sirolimus targeted the PI3K/AKT/MTOR pathway via separate targets. Wortmannin, a microbial product, is a noncompetitive, irreversible inhibitor of phosphatidylinositol-3-kinase (PI3K),42 while sirolimus inhibits the mammalian target of rapamycin (MTOR). Previous studies have reported that sirolimus is a potential therapeutic target for infantile hemangioma and is currently a second-line option in the treatment of complicated IH. The MTOR pathway, which controls protein synthesis, cell growth, metabolism, and survival, is crucial in the development of cancer.43 Most human malignancies have been found to have an overactive MTOR signaling pathway.44 It has been reported that excessive activation of the MTOR pathway and abnormal cell metabolism jointly lead to the occurrence of vascular anomalies.45 According to our KEGG enrichment analysis, the PI3K/AKT/MTOR pathway communicates with a wide range of upstream and downstream signal transduction pathways, primarily the RAS/MAPK signaling pathway, thus playing a core role in the progression of IH (Figure 11). The identified drugs U0126 and sorafenib are inhibitors of the MAPK cascade via mitogen-activated protein kinase kinase (MEK/MKK)46 and the serine kinase RAF,47 respectively, corresponding to the KEGG enrichment results. Furthermore, dasatinib, sorafenib, and sunitinib are all multitarget inhibitors. Among them, sorafenib is an oral multikinase inhibitor with actions against RAF, FLT3, KIT, PDGFR, RET, and VEGFRs, while dasatinib targets KIT, BCR-ABL, EPHRIN, PDGFR, SRC, and tyrosine kinase. Sunitinib is an oral multitarget tyrosine kinase inhibitor with low molecular weight that inhibits PDGF-Rs, VEGFRs, C-KIT, FLT3, and RET. There has been widespread consensus that molecules interfering simultaneously with multiple targets might be more effective than single-target agents. In the future, we can try multitarget agents to treat IH. CGMP/PKG is another signaling pathway we found to play a core role in the development of IH. It has been reported that CGMP/PKG is essential for the maintenance of vascular homeostasis,48 but their function in the progression of IH remains unreported. We also screened four FDA-approved histone deacetylase (HDAC) inhibitors, including belinostat, entinostat, panobinostat, and vorinostat. The majority of clinical indications for this class of medications are hematologic neoplasms, while they are also being used to treat HIV infection, muscular dystrophies, inflammatory conditions, neurological diseases, and Friedreich’s ataxia.49 Further research is needed to elucidate the biological mechanisms underlying the associations between HDAC inhibitors and IH.

Figure 11 Map of the PI3K/AKT/MTOR signaling pathway. Red boxes indicate the differentially expressed genes.

Our investigation comes with several inherent constraints that should be noted. First and foremost, our reliance on public databases for data acquisition might lead to incomplete clinical-pathological details, introducing possible discrepancies and biases. Additionally, the identification of hub genes and drugs was fundamentally based on bioinformatics tools and repositories. Even though these offer valuable preliminary insights, their true accuracy demands rigorous experimental validation in vivo. Although we undertook cellular experiments, our experiments primarily focused on targeted therapies, leaving other potential therapeutics, such as hormone-related drugs and vitamin D, unexplored at this stage. This limitation highlights the need for broader experimental approaches in future studies. Thorough validations involving animal models and clinical human trials are imperative for conclusive evidence. Furthermore, infantile hemangioma is far from a monolithic or single-entity lesion. It emerges from a convoluted microenvironment, characterized by various cell types, signaling agents, and intercellular interactions that jointly orchestrate the IH’s growth, maturation, and eventual regression. This milieu is not limited to endothelial cells; it also encompasses smooth muscle cells and fibroblasts. Therefore, future research endeavors should delve into understanding the drug impacts on these varied cellular components, enriching our understanding of IH’s intricate nature.

Conclusion

In conclusion, we identified the top 100 hub genes using differential gene analysis, weighted gene coexpression network analysis, and PPI analysis. Bioinformatics analyses showed that hub genes were associated with IH pathogenesis via the PI3K/AKT/MTOR pathway, RAS/MAPK pathway, and CGMP/PKG pathway. A total of 44 drugs were predicted, and 12 of them (belinostat, chir 99021, dasatinib, entinostat, panobinostat, sirolimus, sorafenib, sunitinib, thalidomide, U 0126, vorinostat, and wortmannin) might be candidates for the treatment of IH. Furthermore, our in vitro experiments substantiated the therapeutic promise of entinostat, sorafenib, dasatinib, and sirolimus.

Data Sharing Statement

Datasets analyzed in the current study are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA databases.

Funding

This study was financially supported by the National clinical key specialty construction project (23003), the Plastic Medicine Research Fund of Chinese Academy of Medical Sciences (2024-ZX-1-01), and the National Major Disease Multidisciplinary Diagnosis and Treatment Cooperation Project (No.1112320139). The funding body played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Disclosure

The authors report no conflicts of interest in this work.

References

1. Léauté-Labrèze C, Harper JI, Hoeger PH. Infantile haemangioma. Lancet. 2017;390(10089):85–94. doi:10.1016/s0140-6736(16)00645-0

2. Dickison P, Christou E, Wargon O. A prospective study of infantile hemangiomas with a