Integrated analysis of miRNA-mRNA regulatory network and functional verification of miR-338-3p in coronary heart disease

Coronary heart disease is a cardiovascular disease with high morbidity and mortality. Although great progress has been made in treatment, the prognosis is still very poor. Therefore, this project aims to screen potential diagnostic markers and therapeutic targets related to the progression of coronary heart disease. A total of 94 overlapping differentially expressed mRNAs and 70 differentially expressed miRNAs were identified from GSE20681, GSE12288, GSE49823, and GSE105449. Through a series of bioinformatics methods and experiment, we obtained 5 core miRNA-mRNA regulatory pairs, and selected miR-338-3p/RPS23 for functional analysis. Moreover, we found that RPS23 directly targets miR-338-3p by dual luciferase assay, western, and qPCR. And the expression of miR-338-3p and RPS23 is negatively correlated. The AUC value of miR-338-3p is 0.847. Downregulation of miR-338-3p can significantly inhibit the proliferation and migration of HUVEC. On the contrary, overexpression of miR-338-3p promoted the proliferation and migration of HUVEC. In addition, the interference of RPS23 expression can reverse the regulation of miR-338-3p on HUVEC proliferation. In conclusion, miR-338-3p/RPS23 may be involved in the progression of coronary heart disease, and miR-338-3p may be a diagnostic biomarker and therapeutic target for coronary heart disease.


Introduction
Coronary heart disease is one of the most dangerous cardiovascular diseases, which mainly include myocardial infarction, angina pectoris, and sudden death (Case and Waksman 2020;Infante et al. 2020). There are many causes of coronary heart disease, such as obesity, genetics, hyperlipidemia, smoking, and so on (Malakar et al. 2019). At present, the treatment of coronary heart disease is coronary artery bypass surgery. However, the prognosis is very poor, and the mortality rate is still high (Caliskan et al. 2020;Zheng et al. 2020). Therefore, it is necessary to explore the pathogenesis of coronary heart disease and discover new effective biomarkers.
In recent years, with the development of high-throughput technologies, many disease-related databases have been established, such as Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) (Toro-Domínguez et al. 2019;Hutter and Zenklusen 2018). The construction of disease regulation network based on these databases can well screen effective biomarkers for disease diagnosis and prognosis. Zhao et al. (2020) found that through comprehensive bioinformatics analysis of GEO database, six genes related to immune cells in the progression of atherosclerosis have been discovered. The GEO database contains many gene maps in a variety of diseases, which can be used to screen differentially expressed genes in diseases and establish gene regulatory networks (Suzuki et al. 2019). It can also be combined with drug databases to screen effective drug targets and ingredients for diseases. For example, Liang et al. (2020), based on system pharmacology and GEO database mining, revealed the therapeutic mechanism of xuefu zhuyu decoration on atherosclerotic cardiovascular diseases. In addition, the complete function of endothelial cells is an indicator and prerequisite for vascular health, which can inhibit the development of coronary heart disease. The proliferation of vascular endothelial cells is a sign of improvement in coronary heart disease (Yubero-Serrano et al. 2020;Pu et al. 2020;Pelliccia et al. 2020;Rhee et al. 2018). Therefore, human umbilical vein endothelial cell (HUVEC) was selected to construct an in vitro verification model for functional experiments.
In this project, we expect to screen potential biomarkers from GEO database of coronary heart disease RNA expression profile through a series of bioinformatics methods and experiments, and to preliminarily analyze the role of key biomarker in coronary heart disease progression.

Screening of differentially expressed mRNAs and differentially expressed miRNAs
According to the screening criteria of |log2-fold change (FC) |> 1, adjusted p-value < 0.05, and false discovery rate (FDR) < 0.05 (Xie et al. 2022), differentially expressed mRNAs and miRNAs were screened from GSE20681, GSE12288, and GSE49823 by R limma package (Ritchie et al. 2015), respectively. Then we draw the volcano and heat map of the difference mRNA and miRNA through the ggplot2 package and pheatmap package. Finally, the Venn diagram was used to screen the overlapping differentially expressed mRNAs and miRNAs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of overlapping differentially expressed mRNAs
We use the David online tool (http:// david. abcc. ncifc rf. gov/) (Sherman et al. 2022) and Metascape (http:// metas cape. org/ gp/ index. html)  to analyze the GO and KEGG pathways of overlapping differentially expressed mRNAs. Functional analysis is mainly divided into three categories: biological process (BP), cell component (CC), and molecular function (MF). Then we use sangerbox 3.0 (http:// sange rbox. com/ Index) to visualize the results. P < 0.05 was considered statistically significant.

Protein-protein interaction (PPI) network and hub mRNAs identification
The online tool STRING 11.5 (https:// string-db. org/) was used to construct the PPI network, and the confidence of the interaction ≥ 0.4 was considered significant (Szklarczyk et al. 2021). Then, we use Cytoscape v3.9.0 for visualization. Subsequently, molecular complex detection (MCODE) plug-in was used to screen the core genes from the PPI network. The MCODE parameters are set to: degree cutoff = 2, hair-cuff, node score cut-off = 0.2, max depth = 100, k score = 2.

Weighted gene co-expression network analysis (WGCNA)
We preprocess the GSE105449 to remove samples without clinical information, and then perform WGCNA in sangerbox 3.0 (http:// sange rbox. com/ Index). We calculated the median absolute deviation (MAD) of each miRNA, excluded the top 50% of the minimum MAD miRNAs, and eliminated the outlier miRNAs and samples using the R software package WGCNA. Furthermore, WGCNA is used to construct a scale-free co-expression network with the following parameters: height = 0, net type = unsigned, block size = 20,000. Then the adjacency matrix is transformed into a topological overlap matrix (TOM). Then we perform hierarchical clustering to identify the modules, and set the parameters as follows: lock size = 7000, minimodules = 20, deep split = 2, merge cut height = 0.25, hub cut = 0.9. Subsequently, the correlation between the module and the clinical data is analyzed to determine the meaningful clinical module.

mRNA-miRNA interaction network construction
The miRNA-mRNA regulatory network was constructed using the key mRNAs in the PPI regulatory network and miRNAs between differentially expressed miRNAs and WGCNA hub module. The relationship between miRNA and mRNA was analyzed by the online RNA prediction tool starbase 2.0 (http:// starb ase. sysu. edu. cn/) (Li et al. 2014). Finally, we use Cytoscape v3.9.0 to visualize the regulatory network.

RNA extraction and real-time quantitative PCR (RT-qPCR)
A total of 15 peripheral blood of patients with coronary heart disease and peripheral blood of 15 healthy volunteers were collected (supplementary material 1). This study was approved by the hospital ethics committee and was carried out in accordance with the Declaration of Helsinki. In addition, each blood donor provided a written informed consent. Samples were stored at 4 ℃ within 24 h after collection, and total RNA was extracted within 48 h.

Cell proliferation
HUVECs were seeded in 96-well plates at 1.0 × 10 3 cells/ well and cultured for 24 h. Then the corresponding plasmid or small RNA was transfected and incubated for 48 h and 72 h. Then we add 10 μl CCK8 solution (Beyotime, China) to each well and incubate at 37 °C for 1 h. Then we use a multifunctional microplate reader (Thermo, USA) to detect the absorbance value at 450 nm. Process was repeated three times for each hole.

Wound healing assays
HUVECs were cultured in 6-well plates and transfected with proper amounts of plasmid and small RNA. When the cells reached 70-75% fusion (about 24 h), the fused cell layers were cut with the tip of 200-μl sterile pipetting tube, then washed with sterile phosphate buffer solution (PBS, Gibco, USA) three times, and added into serum-free DMEM. After incubation for 48 h, the scratch changes were recorded using a microscope, and the scratch closure rate was analyzed by ImageJ software.

Luciferase reporter assays
The wild-type (WT)-RPS23 3′UTR or mutant (Mut)-RPS23 3′UTR plasmid was co-transfected with miR-338-3p mimic using Lipofectamine 8000 (Invitrogen, USA). After 48 h of incubation, the dual luciferase reporter gene detection kit (Promega, USA) was used to detect luciferase activity, and renal cell luciferase activity was used as an endogenous control.

Western blot
After 48 h of transfection, the transfected cells were collected, and then the protein lysis buffer (Beyotime, China) was used to lysed cells to collect the total protein. The protein was quantified by the BCA kit (Thermo, USA). After 12% SDS-PAGE (Beyotime, China), it was transferred to PVDF membrane (Beyotime, China). Incubate with 5% skimmed milk powder (Beyotime, China) for 3 h, primary antibody (RPS23: ab268130, Abcam, USA; β-actin: AF5003, Beyotime, China) was incubate at 4 °C for 12 h, and secondary antibody (horseradish peroxidase-labeled goat anti-mouse IgG (H + L: A0216, Beyotime, China) was incubated at room temperature for 1 h. Finally, it was developed with western color developing solution (Beyotime, China; P0211) and scanned in an image.

Statistical analysis
GraphPad Prism software was used for all analyses. t test was used for the difference between the two groups, and one-way analysis of variance was used for more than two groups. The data are expressed as mean ± SD. p < 0.05 indicates statistical significance.

Identification of differentially expressed mRNAs and differentially expressed miRNAs
The overall process of this study is shown in supplementary Fig. 1. We use the limma package to screen the differentially expressed mRNAs and miRNAs between the coronary heart disease samples and healthy patient samples. A total of 863 differentially expressed mRNAs (464 upregulated and 399 downregulated) were identified in GSE12288. And 659 differentially expressed mRNAs (501 upregulated and 158 downregulated) were identified in GSE20681. Besides, 70 differentially expressed miRNAs were acquired with 68 downregulated and 5 upregulated in GSE49823. Next, the different mRNA and miRNA expression between coronary heart disease samples and healthy samples were visualized by heatmap and volcano plots ( Fig. 1A and B). Furthermore, 94 overlapping differentially expressed mRNAs were identified by Venn diagram between GSE12288 and GSE20681 (Fig. 1C).

GO and KEGG pathway enrichment analysis of overlapping differentially expressed mRNAs
To determine the biological characteristics of differential expression, the GO and KEGG pathways of overlapping differentially expressed mRNAs were analyzed. As shown in Fig. 2, biological processes (BP) are mainly concentrated in muscle contraction, epithelial cell transformation and proliferation, and immune stress. Cellular component (CC) analysis showed that overlapping genes were mainly concentrated in ribosomal assembly, basement membrane, endoplasmic reticulum, and endocytosis vesicles. Molecular function (MF) is mainly concentrated in protein binding, peptide antigen binding, and signal receptor activity recognition. In addition, KEGG pathway showed that overlapping and differentially expressed mRNAs were mainly enriched in cell cycle, DNA repair, Ras signal pathway, and cAMP signal pathway.

Protein-protein interaction (PPI) network construction and key mRNA analysis
To further understand the interaction between overlapping genes, we constructed a PPI network (Fig. 3A), which included 42 nodes and 230 edges. Then, a hub subnetwork was obtained from PPI network by MCODE plug-in, including 11 nodes and 98 edges (Fig. 3B). The results of GO and KEGG analysis (Fig. 3C) showed that 11 node genes were mainly enriched in immune regulation, immune stress, apoptosis signal, and inflammatory stress process. In addition, KEGG pathway analysis showed that these genes were mainly enriched in cell cycle, DNA repair, WNT signal pathway, cell senescence, interferon signal, and transcriptional regulation pathway.

Weighted gene co-expression network analysis (WGCNA)
To further identify the potential miRNA in coronary heart disease, we performed WGCNA on GSE105449. To ensure that the co-expression network is unsigned, β = 4 was chosen to be soft threshold ( Fig. 4A and B). Then, the average linkage hierarchical clustering method was used to cluster the genes (Fig. 4C). At the same time, we carried out the cluster analysis to the module and obtained 6 modules. Then, we found the turquoise module was the greatest correlation module (Fig. 4D), which included 21 miRNAs. Subsequently, we obtained 2 overlapping hub miRNAs (miR-1296-5p and miR-338-3p, Fig. 4E) from 21 hub miRNAs in turquoise module and 70 differentially expressed miRNAs.
Then, the target genes of 2 miRNA (miR-1296-5p and miR-338-3p) were predicted by starbase, and a total of 549 overlapping targeted mRNAs were screened. Then, we compare the 11 hub mRNAs in the PPI with the 549 overlapping targeted mRNAs. Finally, 3 overlapping abnormal regulation mRNAs (IFITM2, IFITM1, and RPS23, Fig. 4F) were obtained, and the miRNA-mRNA network is visualized by Cytoscape v3.9.0 (Fig. 4G). RPS23 has the highest connectivity in the PPI subnetwork, so RPS23 is selected for follow-up research.

Validation of miR-1296-5p, miR-338-3p, and RPS23
To study the role of key miRNAs and mRNAs in coronary heart disease, we first verified their expression changes in peripheral blood of patients with coronary heart disease. Peripheral blood samples from 15 patients with coronary heart disease and 15 healthy controls were collected for qPCR detection (supplementary material 1). qPCR analysis showed that the expression of miR-1296-5p and miR-338-3p in peripheral blood of patients with coronary heart disease was significantly lower than that of healthy controls (Fig. 4H). On the contrary, the expression of RPS23 was high in the peripheral blood of patients with coronary heart disease ( Fig. 4H), which were consistent with the expression trend of bioinformatics analysis. Subsequently, we found that RPS23 was negatively correlated with the expression of 2 miRNAs (miR-1296-5p and miR-338-3p, Fig. 4I and J). The expression of miR-338-3p was higher in clinical samples and GEO, so miR-338-3p/RPS23 was selected for follow-up analysis.

RPS23 is a direct target of miR-338-3p
According to the online analysis of starbase, miR-338-3p has a binding site (Fig. 5A) with the 3′UTR sequence of RPS23. Dual luciferase experiment results showed that the fluorescence intensity of the miR-338-3p mimic and WT-RPS23 3′UTR co-transfection group was significantly lower than that of the control group. However, the fluorescence intensity of the miR-338-3p mimic and Mut-RPS23 3′UTR co-transfection group did not change significantly (Fig. 5B).
In addition, we also tested the effect of interference miR-338-3p on the expression of RPS23 by qPCR and western. qPCR results showed that miR-338-3p mimic can inhibit the expression of RPS23. In contrast, miR-338-3p inhibitor significantly upregulated the expression level of RPS23 (Fig. 5C). The trend of western results was consistent with qPCR. Overexpression of miR-338-3p reduced the level of RPS23 protein, and inhibition of miR-338-3p increased the amount of RPS23 protein (Fig. 5D and E). These results indicate that RPS23 is a direct regulatory target of miR-338-3p and is negatively regulated.
In addition, ROC analysis was used to evaluate the expression of miR-338-3p between patients with coronary heart disease and healthy samples. The area under the curve of miR-338-3p was 0.7117, indicating that the experimental and analytical results are reliable (Fig. 5F). However, the role of miR-338-3p and miR-338-3p/RPS23 in the progression of coronary heart disease is still unclear.

miR-338-3p contribute to HUVEC cell proliferation and migration
To study the role of miR-338-3p in coronary heart disease, we transfected miR-338-3p mimic and miR-338-3p inhibitor into HUVECs. The qPCR results showed that the expression level of miR-338-3p was significantly upregulated with miR-338-3p mimic, and the expression of miR-338-3p was downregulated in cells transfected with miR-338-3p inhibitor (Fig. 6A). Subsequently, CCK8 results showed that inhibiting the expression of miR-338-3p significantly reduced the proliferation activity of HUVECs. At the same time, overexpression of miR-338-3p promoted the proliferation of HUVECs (Fig. 6B). Similarly, wound healing showed that the miR-338-3p overexpression group had a significant increase in cell migration rate compared with the control group ( Fig. 6C and

Interference with RPS23 can reverse the effect of miR-338-3p on HUVECs
To further explore the biological functions of the miR-338-3p/ RPS23 axis, we interfered with miR-338-3p and RPS23 at the same time, and analyzed the effect of the interference on HUVECs. In correlation experiments and qPCR experiments, we found that miR-338-3p expression and RPS23 expression were negatively correlated (Fig. 4I). And interference with the expression of miR-338-3p can correspondingly affect the expression of RPS23. CCK8 experiment showed that overexpression of RPS23 inhibited HUVEC cell proliferation. On the contrary, knocking out RPS23 can significantly increase the proliferation activity of HUVECs (Fig. 6E). In addition, miR-338-3p mimic-treated cells and transfected with pCDH-RPS23 can significantly reverse the effect of miR-338-3p mimic on the proliferation activity of HUVECs. Similarly, miR-338-3p mimic-treated cells and transfected with si-RPS23 can significantly enhance the effect of miR-338-3p mimic on the proliferation activity of HUVECs (Fig. 6F). These results indicate that overexpression of miR-338-3p mimic promotes the growth of HUVECs by downregulating RPS23. Overexpression of RPS23 can reverse the proliferation regulation of miR-338-3p mimic.

Discussion
In this study, we downloaded mRNA and miRNA expression profile from GEO for analysis. Based on bioinformatics analysis such as PPI, GO, and KEGG, we have obtained 11 key differentially expressed mRNAs, which are mainly involved in immune regulation, apoptosis signal and inflammatory stress process, cell cycle, DNA repair, and WNT signal pathway. Subsequently, the hub miRNA-mRNA pair: miR-338-3p/RPS23 was screened through WGCNA, miRNA-mRNA network, qPCR, and correlation analysis. We determined that RPS23 is the direct target of miR-338-3p through online prediction software and luciferase reporter gene detection. Therefore, we believe that miR-338-3p can target RPS23 and participate in the regulation of a series of cell activities in coronary heart disease. Generally, miRNAs directly bind to mRNA 3′UTR to regulate the gene expression by degrading or inhibiting translation (Ali Syeda et al. 2020). Various miRNAs are involved in the pathogenesis of cardiovascular disease (coronary heart disease) by binding to various target genes. For example, Zhu et al. found that MALAT1/miR-15b-5p/MAPK1 mediates endothelial progenitor cell autophagy and affects coronary atherosclerotic heart disease through the mTOR signaling pathway (Zhu et al. 2019). Lin et al. found that miRNA-451b participates in the progression of coronary The effects of miR-338-3p mimic and miR-338-3p inhibitor on the migration of HUVECs measured using wound healing assay. (E) The effects of si-RPS23 and pCDH-RPS23 on the proliferation activity of HUVECs measured using CCK-8. (F) The effects of co-transfection (si-RPS23 and miR-338-3p mimic; pCDH-RPS23 and miR-338-3p mimic) on the proliferation activity of HUVECs measured using CCK-8. The horizontal line indicates the comparison between groups at both ends, and the non-horizontal line indicates the comparison with the control group. NC, negative control. ***p < 0.001; **p < 0.01; *p < 0.05 heart disease by targeting VEGFA (Lin et al. 2019). Li et al. found that miR-145-5p alleviates cardiac microvascular endothelial cell injury induced by hypoxia/reoxygenation in coronary heart disease by inhibiting the expression of Smad4 . Wang et al. also found that serum miR-3129-5p mediates mTOR to promote the progression of coronary heart disease, suggesting that miR-3129-5p may become one of the biological targets for the treatment of coronary heart disease (Wang et al. 2021). Therefore, we speculate that miR-338-3p may be involved in the pathological process of coronary heart disease by directly binding its target genes. There are many causes of coronary heart disease, and the specific pathogenesis is not very clear. However, a large amount of data indicates that miRNA can promote the occurrence of coronary heart disease by regulating related genes Fazmin et al. 2020). Previous studies have shown that miR-338-3p participates in tumor angiogenesis and promotes angiogenesis (Zhang et al. 2016). In addition, miR-338-3p can regulate the proliferation and apoptosis of vascular endothelial cells in atherosclerotic mice, which indicates that miR-338-3p has a regulatory role in cardiovascular diseases (Yu et al. 2020). Yin et al. demonstrated through miRNA sequencing data that miR-338-3p can induce endothelial cell injury by targeting BAMBI and activating TGF-beta/Smad pathway (Yin et al. 2019). Based on these data, we attempted to analyze the effect of miR-338-3p on HUVECs. Therefore, this study confirmed that miR-338-3p regulates the proliferation and migration of HUVECs by targeting RPS23. Ribosomal protein S23 (RPS23) is a component of the 40S small ribosomal subunit encoded by the RPS23 gene (Wang et al. 2013). Studies have shown that RPS23 plays an important role in cell growth, protein folding, apoptosis, and autophagy. At the same time, RPS23 was found to be located near the mRNA entrance channel and played an important role in identifying the interaction of specific human eukaryotic translation initiation factors (Wang et al. 2013;Katz et al. 2014). But at present, there is no relevant research data of RPS23 in cardiovascular disease. Our research results indicate that miR-338-3p/RPS23 may be a potential biomarker or therapeutic target for coronary heart disease, which may help to explore new therapies for coronary heart disease.
Previous studies reported that miR-338-3p showed tumor suppressor effects in different types of cancer. As an anti-tumor molecule, miR-338-3p is downregulated in hypopharyngeal carcinoma and inhibits the migration and invasion of human hypopharyngeal carcinoma cells by targeting ADAM17 (Hong et al. 2020). miR-338-3p acts as a tumor suppressor in gastric cancer by targeting PTP1B (Sun et al. 2018). In this study, the selected miR-338-3p showed low expression in the patients with coronary heart disease. In vitro, we proved that the low expression of miR-338-3p inhibits the proliferation and migration of HUVECs through gain-and loss-of-function experiments. More importantly, we confirmed that miR-338-3p regulates the migration and proliferation of HUVEC by targeting RPS23. Overexpression of miR-338-3p can reduce the mRNA and protein levels of RPS23 in HUVECs. At the same time, overexpression of RPS23 can reverse the proliferation activity of miR-338-3p mimic on HUVECs. However, this study has some limitations. In this study, peripheral blood samples of patients were collected, and there is a lack of tissue sample verification. In addition, the specific pathway of miR-338-3p/RPS23 axis activation is unknown, and the in vivo experimental data is worthy of further study.

Conclusion
In summary, through bioinformatics analysis, a total of 11 hub mRNAs and 2 hub miRNAs were obtained, and the key pair: miR-338-3p/RPS23 was screened. In vitro functional experiments confirmed that miR-338-3p/RPS23 axis affected the proliferation and migration of HUVECs, indicating that miR-338-3p/RPS23 may be involved in the progression of coronary heart disease. In addition, miR-338-3p may be a therapeutic target for the treatment of coronary heart disease.