The siRNA for HSPB1 (siHSPB1) and non-targeting control siRNA (siNC) were purchased from Genepharma (Suzhou, Jiangsu, China). The sequence of siNC was: 5^′-UUCUCCGAACGUGUCACGUTT-3^′ (sense). The sequence of siHSPB1 was: 5^′-CGGACGAGCUGACGGUCAATT-3^′ (sense). We also used the blank cell as background to assess the knockdown level of HSPB1.

The transfection reagent Lipofectamine™ RNAiMAX (13778030, Invitrogen, Carlsbad, CA, USA) was used for siRNA transfection into HUVECs (DFSC-EC-01, Zhongqiaoxinzhou Biotech, Shanghai, China). The cells were cultured at 37 °C with 5% CO₂ in ECSM (ZQ-1304, Zhongqiaoxinzhou Biotech, Shanghai, China) with 10% fetal bovine serum (FBS) (10099-141, Gibco, Shanghai, China). The HUVECs were tested by short tandem repeat (STR) profiling for validation, and tested negative for mycoplasma. Transfected cells were harvested after 48 h for the subsequent experiments.

RT-qPCR was performed to validate the knockdown efficiency of HSPB1. Complete DNA (cDNA) was synthesized according to standard procedures using the Bio-Rad S1000 instrument (Hercules, CA, USA) and Hieff™ qPCR SYBR® Green Master Mix (Low Rox Plus, YEASEN, Shanghai, China). We used the 2^{-Δ⁢Δ⁢CT} method [23] to calculate the concentration of each transcript, with the level of HSPB1 mRNA normalized to that of GAPDH mRNA.

To evaluate the protein level of HSPB1 in HUVEC cells, the experimental procedure for Western blot was conducted as previously described [24]. The membranes were incubated with the primary antibodies to HSPB1 (1:1000, A0240, Abclonal, Wuhan, Hubei, China) and ACTIN (1:1000, 20536-1-AP, Proteintech, Wuhan, Hubei, China), and then with HRP-conjugated secondary antibody. The secondary antibody (anti-rabbit, 1:5000, SA00001-2, Proteintech, China; or anti-mouse, 1:5000, AS003, ABclonal, China) was detected using enhanced chemiluminescence (ECL) reagent (Bio-Rad, 170506, Hercules, CA, USA).

The Cell Counting Kit-8 (CCK-8) assay kit (HY-K0301, MedChemExpress, Monmouth Junction, NJ, USA) was used to evaluate the proliferation of HUVECs. These were seeded into a plate with 96-wells at a density of 10⁴ cells per well (6 replicates) and cultured in an incubator for 0, 24, 48, and 72 h. The 10 µL of CCK-8 solution was added to the wells at each time point and then incubated for a further 3 h. The absorbance of each well at a wave length of 450 nm was then measured using a microplate reader (ELX800, Biotek, Winooski, VT, USA). OD450 values were used to evaluate cell proliferation.

The Annexin V-APC/7-ADD kit (KGA1026, KeyGEN BioTECH, Nanjing, Jiangsu, China) was used to detect apoptosis according to the user instructions. We seeded 10⁶ HUVECs into plates, cultured for 24 h, and then transfected with siHSPB1 or siNC for 48 h. The cells were subsequently mixed with 5 µL of Annexin V-allophycocyanin (V-APC) and incubated at room temperature for 5 min, followed by incubation with 5 µL of 7-aminoactinomycin D (7-AAD) reagent for 5 min. The two cell groups were then analysed by flow cytometry (FACSCanto, BD, Franklin Lake, NJ, USA) to assess apoptosis.

The total RNAs were extracted from the siNC and siHSPB1 HUVEC groups (3 replicates) using the TRIzol method. RNA concentration and purity were measured using Smartspec Plus (A260/A280; BioRad, Hercules, CA, USA). After passing the quality inspection, 1 µg of total RNA was used to construct an RNA-seq library using the VAHTS Universal V8 RNA-seq Library Prep Kit for Illumina (NR605, Vazyme, Nanjing, Jiangsu, China). The ligation product was then purified and the 300–500 bp size fraction was selected. The products were then amplified, purified, and stored at –80 °C. High-quality RNA sequences were obtained using the NovaSeq 6000 sequencing platform (Illumina, San Diego, CA, USA) for paired-end 150 nt sequencing.

“Atherosclerosis” was used as the keyword to retrieve and select the appropriate dataset from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE104140 dataset was downloaded for further analysis, and has not been used in any previous published study.

The FASTX-Toolkit (Version 0.0.13, The Hannon Lab, New York, NY, USA) was used to discard low-quality reads from the two RNA-seq datasets. HISAT2 2.2.1 (https://daehwankimlab.github.io/hisat2/) [25] was used to align filtered reads onto the human GRCh38 genome. Differentially expressed genes (DEGs) were identified using DESeq2 https://bioconductor.org/packages/release/bioc/html/DESeq2.html [26], with criteria of (fold-change) FC $\ge$3/2 or $\le$2/3, and a p-value 0.01 used to define a significant difference.

Gene Ontology (GO, http://geneontology.org/) was used to describe the functions of gene sets. GO enrichment analysis was conducted to explore the functions of DEGs. These were first mapped to various terms in the GO database and the number of genes associated with each term was counted. Hypergeometric distribution testing was then used to obtain significantly enriched GO terms (p-value 0.01) associated with DEGs based on the GO annotation of the entire genome.

We used ABLas pipeline (ABLife, Wuhan, Hubei, China), as previously described [27], to compare the differences in alternative splicing events (ASEs) and define and quantify the regulated alternative splicing events (RASEs) between groups. To identify RASEs that were different after silencing the HSPB1 gene, Student’s t-test method was used to calculate changes in alternative splicing of the same splicing type for each gene between siHSPB1 and siNC. The Events with significance at a p-value cutoff of 5% were considered to be RASEs.

The Gene Ontology (GO) enriched pathways were identified using the KOBAS 2.0 server [28] (version 2.0, Peking University, Beijing, China). Experimentally quantitative data were presented by mean $\pm$ standard deviation (SD). Student’s t-test was used to calculate the difference between siHSPB1 and siNC groups (Prism version 8.0, GraphPad, Boston, MA, USA) with p-value 0.05 as significant difference. Sequencing data was analyzed using R software (version 4.3.1, The R Foundation, Vienna, Austria) and plotted using GraphPad Prism software (Version 9.5.0, GraphPad, Boston, MA, USA) and R software.

To examine its role in HUVECs, the HSPB1 gene was silenced using the siRNA transfection method (siHSPB1). The mRNA and protein expression levels of HSPB1 were then evaluated to determine the silencing efficiency. The HSPB1 protein (Fig. 1A, Supplementary Fig. 1) and mRNA (Fig. 1B) in HUVECs showed significantly lower levels in the siHSPB1 group compared to the siNC and blank groups (p 0.0001). Then we used the siNC group as the background in following experiments. We next investigated the changes in phenotype induced by HSPB1 knockdown in HUVECs. Flow cytometry analysis revealed that the apoptosis ratio in the siHSPB1 group was significantly lower than in the siNC group (Fig. 1C, p 0.0001). Moreover, CCK-8 analysis showed that the cell proliferation rate in the siHSPB1 group was significantly higher than that in the siNC group (Fig. 1D, p 0.05). These results indicate that HSPB1 can obviously modulate the cellular phenotype of HUVECs.

Fig. 1.

SiHSPB1 significantly affects the proliferation and apoptosis of HUVECs. (A)Western blot results after HSPB1 knockdown; the right panel was the quantitative results (n = 3). * p 0.05, ** p 0.01. (B) RT-qPCR results after HSPB1 knockdown (n = 3). **** p 0.0001. (C) Flow cytometry and bar plot showing the level of apoptosis after HSPB1 knockdown. The right panel was used for statistical analysis (n = 3). **** p 0.0001. (D) Bar plot showing the rate of cell proliferation after HSPB1 knockdown (n = 6). * p 0.05, ** p 0.01. siNC, non-targeting control siRNA; siHSPB1, The siRNA for HSPB1; HUVECs, human umbilical vein endothelial cells; HSPB1, heat shock protein family B member 1; RT-qPCR, reverse transcription and quantitative polymerase chain reaction; APC-A, allophycocyanin-A.

To clarify the potential mechanism by which HSPB1 regulates HUVECs, RNA-seq and bioinformatics analysis was performed on the siNC and siHSPB1 groups of HUVECs. The sequencing results showed a significant decrease (p-value 0.001) in HSPB1 gene expression in the siHSPB1 group (Fig. 2A), indicating successful knockdown. Then we performed principal component analysis (PCA) on all of the expressed genes from the siNC and siHSPB1 groups. PCA result demonstrated a clear separation between the two groups (Fig. 2B), indicating that HSPB1 knockdown altered the pattern of global transcription. DEGs between these two groups were identified using DESeq2. The fold-change (FC) and p-value were used to determine whether a gene was differentially expressed (assessment criteria: FC $\ge$3/2 or $\le$2/3, and p-value 0.01). A total of 608 DEGs, including 423 upregulated genes and 185 downregulated genes, were finally identified (Fig. 2C,D).

Fig. 2.

HSPB1 regulates gene expression in HUVECs. (A) Bar plot demonstrating the changed expression level of HSPB1 from RNA-seq (n = 3). *** p 0.001. (B) PCA result based on the FPKM values of all expressed genes after HSPB1 knockdown, with confidence ellipse drawn for each group. (C) Volcano plot showing the pattern and number of DEGs between siHSPB1 and siNC groups. (D) Heat map showing the expression pattern of all detected DEGs by hierarchical clustering. (E) Bubble plot showing the top 10 enriched GO biological processes amongst the up-regulated DEGs. The two most enriched pathways are highlighted with red font. (F) The same analysis as in (E), but for the down-regulated DEGs. PCA, principal component analysis; DEGs, differentially expressed genes; GO, Gene Ontology; FPKM, fragments per kilobase per million.

GO analysis of the upregulated genes revealed enrichment of the following pathways: inflammatory response, cytokine-mediated signalling pathway, cytoplasmic translation, signal recognition particle (SRP)-dependent co-translational proteins targeting membranes, positive regulation of gene expression, translation initiation, translation, positive regulation of interleukin-10 production, viral transcription, and intercellular signalling processes (Fig. 2E). This indicates that HSPB1 participates in inflammatory responses and in cytokine-mediated signalling pathways. GO analysis of the 185 downregulated genes revealed enrichment for immune system processes, innate immune response, type I interferon signalling pathway, defence response to viruses, antigen processing, exogenous peptide antigen through major histocompatibility complex (MHC) class I, viral response, transporter associated with antigen processing (TAP)-independent presentation, negative regulation of viral genome replication, antigen processing, peptide antigen presentation through MHC class I, I-$\kappa$$\beta$ kinase/nuclear factor kappa $\beta$ (NF-$\kappa$$\beta$), positive regulation of signalling pathways, and regulation of complement activation in biological pathways (Fig. 2F). In summary, the above results indicate that HSPB1 modulates the inflammatory and immune response of HUVECs, both of which are closely associated with the pathogenesis of AS [29].

By analysing the alternative splicing events regulated by HSPB1 (Fig. 3A), HSPB1 knockdown was found to promote exon hopping and suppress cassette exons. These events may be characteristic of the regulation of AS by HSPB1. To investigate the potential influence of HSPB1-regulated AS events (RASEs) and genes (RASGs), we carried out functional enrichment analysis for RASGs. This revealed that RASGs were mainly enriched in the processes of apoptosis, the I-$\kappa$$\beta$ kinase/NF-$\kappa$$\beta$ biological functional pathway, mRNA processing, RNA splicing, protein transport, mRNA splicing during glycosphingolipid metabolism, type I interferon signalling pathway, and the regulation of cell morphology (Fig. 3B). siHSPB1 was found to induce the splicing pattern of genes that participate in the apoptotic process, consistent with the phenotype alteration observed for HUVECs. A total of 26 overlapping genes were found between RASGs and DEGs (Fig. 3C), of which ACIN1, IFI27, PAK4, UBE2D3, and FIS1 were enriched in the apoptotic pathway. Therefore, we further analysed the RASEs in FIS1 and UBE2D3 by siHSPB1. The mutually exclusive 5^′ untranslated region (5pMXE) of UBE2D3 was included by siHSPB1 (Fig. 3D), whereas the 5pMXE of FIS1 was excluded by siHSPB1 (Fig. 3E).

Fig. 3.

HSPB1 regulates alternative splicing pattern of genes in HUVECs. (A) Bar plot demonstrating the HSPB1-regulated alternative splicing events (RASEs). (B) Bubble plot showing the top 10 enriched biological process pathways for regulated alternative splicing genes (RASGs). The apoptotic process is highlighted in red. (C) Venn diagram demonstrating the overlapped genes between RASGs and DEGs. (D) HSPB1 regulates the 5pMXE splicing event for UBE2D3. The left panel shows the IGV-sashimi plot with the sequencing density on transcripts, while the right panel shows the RNA-seq validation of ASEs. Error bars represent the mean $\pm$ SEM (n = 3). * p 0.05. (E) The same analysis as in (D), but showing the 5pMXE for FIS1. Error bars represent the mean $\pm$ SEM (n = 3). ** p 0.01. SEM, standard error of the mean; ASEs, alternative splicing events; 5pMXE, mutually exclusive 5^′ untranslated region; IGV, integrative genomic viewer.

To further decipher the regulatory targets of HSPB1 in AS, we downloaded the RNA-seq data (GSE104140) of human fibroatheroma patients from the GEO database. This study analysed and compared the RNA-seq transcriptome data from early vascular disease (diffuse intimal thickening, SAMP_DIT) and two late AS states (calcification, SAMP_cal and non-calcified fibrous AS, SAMP_no cal). Data for early vascular disease (SAMP_DIT) and late calcified fibrous AS (SAMP_cal) were used for analysis. The results indicated that HSPB1 expression was significantly lower (p-value 0.001) in the SAMP_cal group than in the SAMP_DIT group (Supplementary Fig. 2A), implying that it may have an important role in AS. Moreover, thousands of DEGs were identified between the SAMP_cal vs. SAMP_DIT groups, including 2740 up-regulated and 987 down-regulated genes (Supplementary Fig. 2B,C). An overlap analysis was performed between 608 HSPB1-regulated DEGs and 3716 genes that were differentially expressed according to the GSE104140 RNA-seq data. This identified 171 overlapping DEGs (Fig. 4A, p = 8.87 $\times$ 10^-28, hypergeometric test). GO pathway enrichment analysis for these overlapped DEGs revealed that inflammatory response showed the strongest enrichment (Fig. 4B). Then we performed a heatmap analysis showing the expression pattern of the 22 DEGs from the inflammatory response pathway, including 6 down-regulated DEGs and 16 up-regulated DEGs (Fig. 4C). Finally, we selected 8 DEGs that showed significant differences in expression and were related to inflammatory response pathways. Their expression pattern according to HSPB1 (HUVEC) KD RNA-seq (Fig. 4D) and GSE104140 RNA-seq (Fig. 4E) data is shown using bar graphs.

Fig. 4.

HSPB1 regulates the expression of genes that were changed in the vascular tissue of AS patients. (A) Venn diagram demonstrating overlapped genes between the GSE104140 dataset and post-HSPB1 knockdown of HUVECs. (B) Scatter plot demonstrating the top 10 enriched GO biological processes for the overlapping DEGs. The inflammatory response is highlighted. (C) Heat map of the overlapped genes for the inflammatory response by hierarchical clustering. (D) Bar plot showing the changed expression pattern of overlapped DEGs between the siNC and siHSPB1 groups. Error bars represent the mean $\pm$ SEM (n = 3). ***p 0.001, ** p 0.01. (E) The same analysis as in (D), but for the GSE104140 RNA-seq dataset. Error bars represent the mean $\pm$ SEM (n = 9). *** p 0.001. AS, atherosclerosis.

To further study the role of HSPB1 in regulating alternative splicing in AS, we overlapped 1096 RASGs identified by RNA-seq of HSPB1 (HUVEC) KD genes with 2096 RASGs identified by RNA-seq of the GSE104140 dataset. This revealed a total of 250 overlapping RASGs (Fig. 5A, p = 4.35 $\times$ 10^-38, hypergeometric test), many of which were enriched in the apoptotic process, as shown by GO pathway enrichment analysis (Fig. 5B). Overlapping RASGs from the apoptotic process were displayed through heatmaps, revealing obvious differences between the siHSPB1 and siNC groups (Fig. 5C). Finally, the distribution of the splicing model in siHSPB1 HUVECs and in AS patients, as well as bar plots showing the dysregulated splicing pattern of regulator of G-protein signaling 3 (RGS3) and NOD-like receptor family CARD domain containing 5 (NLRC5), are presented in Fig. 5D,E. A previous study reported that NLRC5 can regulate vascular remodelling by directly inhibiting smooth muscle cell dysfunction via its interaction with peroxisome proliferators-activated receptors $\gamma$ (PPAR$\gamma$) [30]. In summary, interaction analysis of the splicing pattern suggests that HSPB1 participates in the progression of AS by modulating the splicing of important genes.

Fig. 5.

HSPB1 modulates the splicing pattern of genes in the vascular tissue of AS patients. (A) Venn diagram demonstrating the overlapped RASGs between the GSE104140 dataset and the HSPB1 knockdown RNA-seq of HUVECs. (B) Scatter plot showing the top 10 enriched GO biological processes for the overlapping RASGs. The apoptotic process is highlighted. (C) Hierarchical clustering heat map of overlapping alternative splicing events and genes in the apoptotic process. (D,E) HSPB1 regulates alternative splicing for RGS3 and NLRC5. Schematic diagrams and RNA-seq validation for the RASEs. Error bars represent the mean $\pm$ SEM (n = 3 and n = 9 for D and E, respectively). *** p 0.001, ** p 0.01, * p 0.05