The microarray datasets GSE160717 (circRNA), GSE31568 (miRNA), GSE61741 (miRNA), GSE24519 (mRNA), and GSE66360 (mRNA) were obtained from the Gene Expression Omnibus (GEO) database.¹ GSE160717 is based on the GPL21825 platform (074301 Arraystar Human CircRNA microarray V2). GSE31568 and GSE61741 are based on the GPL9040 platform (febit Homo Sapiens miRBase 13.0). GSE24519 is based on the GPL2895 platform (GE Healthcare/Amersham Biosciences CodeLink Human Whole Genome Bioarray). GSE66360 is based on the GPL2895 platform ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array). These datasets satisfy the following characteristics: (1) Organism: Homo sapiens; (2) Sample source: Blood; and (3) Grouping: AMI group and healthy controls. Table 1 shows the detailed information of these five microarray datasets. These datasets are widely used for bioinformatics analyses to explore the molecular mechanisms associated with AMI (17–19).

Raw sequencing data of GSE160717 (circRNA), GSE31568 (miRNA), GSE61741 (miRNA), and GSE24519 (mRNA) were downloaded from the abovementioned platform and preliminary differential expression analysis between AMI samples and control was performed using GEO2R.² Benjamini-hochberg method was used to correct adj. P for potential false positive results. FDR < 0.05 and |logFC| > 1.0 was set as the cut-off criteria of DERNAs. CoDEmiRNAs in GSE31568 and GSE61741were identified using Venn diagram web tool.³ Based on relative expression information in raw sequencing data, the heatmap of DEcircRNAs were produced using the package gg plots of R software (version: x64 3.2.1) (20).

Using circBase database,⁴ the circRNA name in GSE160717 containing six digits was converted into circBase ID containing seven digits. circRNAs that do not match in the annotation file were further blatted to complete the ID conversion. The miRDB⁵ database was used to convert the miRNA ID in GSE31568 and GSE61741 from the previous version to the latest version.

CircBase is a circRNA-related database that records circRNA-related annotation information, such as sequence, gene and genomic location (21). DEcircRNAs in GSE160717 were entered into the circBase database to get their annotation information. CSCD⁶ database is the first comprehensive cancer-specific circRNA database which can predict the microRNA response element (MRE) sites, RNA binding protein (RBP) sites, and potential open reading frames (ORF) (22). Structure of these DEcircRNAs were predicted and visualized using CSCD.

Based on DEcircRNAs and CoDEmiRNAs, the circRNA-miRNA interactions were predicted using website tools of circBank⁷ and circInteractome⁸ (23, 24). Based on CoDEmiRNAs and differentially expressed genes (DEGs), the miRNA–mRNA interactions were predicted using website tools of miRDB,⁹ Targetscan (version7.2¹⁰), and miRWalk¹¹ (25–28). According to these regulatory relationships, the up-regulated circRNAs and down-regulated circRNAs were used to construct a circRNA-miRNA-mRNA regulatory network respectively. The network diagram wsa visualized by Cytoscape software (version 3.7.1).

The Database for Annotation, Visualization, and Integrated Discovery (DAVID; version6.8¹²) was used to perform Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs in these two circRNA-miRNA-mRNA regulatory networks, respectively (29). P value < 0.05 was set as the threshold to identify statistically significant pathways.

In order to better understand the molecular mechanism of circRNA’s involvement in transcriptional regulation, a protein–protein interaction network (PPI network) of DEGs in the circRNA-miRNA-mRNA network was constructed. The String database (version10.0¹³) was used to predict the interaction between proteins, and a PPI network was established by Cytoscape software (version 3.7.1) subsequently.

The scores of Maximal Clique Centrality (MCC) algorithm was set as the criteria (30, 31), and the top10 Hub genes with high connectivity in the PPI network was screened out using the Cytohubba plug-in of Cytoscape software. In addition, the MCODE plug-in of Cytoscape software was used to predict meaningful modules with default parameters as thresholds. Finally, the key pathways of these modules were predicted by KOBAS (version 3.0¹⁴) database.

Based on hub genes identified in PPI network, corresponding miRNAs and circRNAs were selected from circRNA-miRNA-mRNA regulatory networks mentioned above. R software ggplot2, ggalluvial, and dplyr package was used to draw sankey diagrams (32). The heat map of circRNAs in this sankey diagram were produced using the package gg plots of R software (version: x64 3.2.1).

Hub gene and hub gene-related miRNAs in the sankey diagram were considered as candidate biomarkers for AMI. Based on their raw expression data in datasets GSE31568 (miRNA), GSE61741 (miRNA), GSE24519 (mRNA), and GSE66360 (mRNA), Graphpad Prism (version:8.0) software was used to plot the receiver operator characteristic curves and calculate the area under the curve (AUC) respectively to evaluate the predictive ability of these biomarkers for AMI. Meanwhile, hub genes and miRNAs with AUC ≥ 7.0 and P < 0.05 were screened as potential biomarkers. Their relative expression levels in these datasets were further evaluated by box plots. Graphpad Prism 9 software was used for analysis and mapping to compare the control and AMI groups. Comparison between the two groups was performed by independent sample T test. P < 0.05 was considered statistically significant.

Human cardiomyocyte AC16 cell line were purchased from BeNa Culture Collection (Beijing, China). The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Gibco) with 10% fetal bovine serum (FBS, Gibco) and 1% penicillin-streptomycin (Sigma, St. Louis, MO, United States) at 37°C with 5% CO₂.

The AC16 cardiomyocytes were exposed to 24 h of hypoxia (1% O₂, 5% CO₂, and 94% N₂) in DMEM to generate a hypoxia model. AC16 cells under normoxia condition all the time were regarded as the control group.

According to the manufacturer’s protocol, total RNA was isolated using FastPure^® Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Nanjing, China), and the mRNA and circRNA cDNAs were synthesized with Hifair^® II 1st Strand cDNA Synthesis Kit (gDNA digester plus) (Ye Sen, Shanghai, China). The qRT-PCR reactions were carried out in the Bio-Rad CFX96 Real-time PCR Detection System, using Hieff UNICON^® Universal Blue qPCR SYBR Green Master Mix (Ye Sen, Shanghai, China). GAPDH was used as the reference genes and the relative fold change was evaluated by the method of 2^–ΔΔCt. All primers were presented in Supplementary File 3. Comparison between the two groups was performed by independent sample T test. P < 0.05 was considered statistically significant.

The microarray datasets GSE160717 (circRNA), GSE31568 (miRNA), GSE61741 (miRNA), and GSE24519 (mRNA) were obtained from the Gene Expression Omnibus (GEO) database. Based on GEO2R software analysis and screening criteria mentioned above, DERNAs in each dataset were identified. There were 83 DEcircRNAs, including 50 upregulated and 33 downregulated, in GSE160717. There were 249 DEmiRNAs, including 101 upregulated and 148 downregulated, in GSE31568. There were 143 DEmiRNAs, including 88 upregulated and 55 downregulated, in GSE61741. Finally, there were 1204 DEGs, including 766 upregulated and 438 downregulated, in GSE24519. Supplementary File 1 presents detailed results of differential expression analysis. The volcano plots of these DERNAs in each dataset are shown in Figures 2A–D. The online venn tool was used to identify CoDEmiRNAs in GSE31568 and GSE61741. Finally, 58 upregulated and 51 downregulated CoDEmiRNAs were found (Figures 2E,F). The top 5 high-expressed DERNAs and the top 5 low-expressed DERNAs in each dataset are shown in Table 2.

The DEcircRNA’s name in GSE160717 containing six digits was converted into circBase ID containing seven digits (Supplementary File 2). Figure 3A is the Heatmap of these DEcircRNAs. DEcircRNAs were entered into the circBase database to get their annotation information. The annotation information of the top 5 high-expressed DEcircRNAs and the top 5 low-expressed DEcircRNAs are shown in Table 3. Structure of these DEcircRNAs were predicted and visualized using CSCD database. The structural patterns of the top 3 high-expressed DEcircRNAs and the top 3 low-expressed DEcircRNAs are shown in Figures 3B–G.

In order to better understand the molecular mechanism of circRNAs involved in the process of AMI. The up-regulated circRNAs and down-regulated circRNAs were used to construct a circRNA-miRNA-mRNA regulatory network respectively. Figure 4A is an “up-regulated circRNA”-“down-regulated miRNA”-“up-regulated mRNA” network, with 134 nodes and 166 edges, including 19 up-regulated circRNAs, 28 down-regulated miRNA and 87 up-regulated mRNAs. On the contrary, Figure 4B is a “down-regulated circRNA”-“up-regulated miRNA”-“down-regulated mRNA” network, with 72 nodes and 68 edges, including 18 down-regulated circRNAs, 15 up-regulated miRNA and 39 down-regulated mRNAs. These RNA interactions may provide new insight into the mechanism underlying AMI.

Gene Ontology knowledgebase and KEGG database were used to characterize the functional roles of the DEGs in the circRNA-miRNA-mRNA Regulatory Network mentioned above. Figure 5 list the top10 enriched GO terms and KEGG pathways. The biological process (BP) category of the GO analysis results showed that these DEGs are significantly enriched in the term of “positive regulation of transcription, DNA-templated,” “negative regulation of apoptotic process,” and “calcium ion transmembrane transport,” etc. (Figure 5A). For GO cellular component (CC) analysis, the top 4 significantly enriched terms are “cytoplasm,” “neuron projection,” “neuronal cell body,” “actin cytoskeleton,” (Figure 5B). The top four significantly enriched molecular function (MF) terms include “RNA binding,” “enzyme binding,” “microtubule binding,” and “calcium channel activity” (Figure 5C). Furthermore, “FOXO signaling pathway,” “Carbohydrate digestion and absorption,” “T cell receptor signaling pathway,” “MAPK signaling pathway,” “Insulin resistance,” “cAMP signaling pathway,” and “mTOR signaling pathway” are pathways of significant enrichment in KEGG analysis (Figure 5D).

The DEGs obtained from the circRNA-miRNA-mRNA Regulatory Network were introduced into the online database String. After removing the isolated genes without interaction, a PPI network with 105 nodes and 203 edges was established (Figure 6A). The 74 orange red nodes in the network represent up-regulated genes and the 31 blue nodes in the network represent down-regulated genes. Next, based on MCC method, the top 10 hub genes were filtered out using the plug-in Cytohubba in Cytoscape. They were ATP2B2, KCNA1, GRIN2A, SCN2B, GPM6A, CACNA1E, HDAC2, SRSF1, ANK2, and HNRNPA2B1. The PPI network of these 10 hub genes is shown in Figure 6B. On the basis of MCODE, we identified 3 modules in the whole network (Figures 6C–E).

The KOBAS online tool was used for pathway enrichment analysis. The results indicated that “the DEGs in module 1” and “hub genes” are significantly enriched in the pathway of Calcium signaling pathway. The DEGs of module 2 are enriched in the pathway of RNA polymerase. In addition, the DEGs of module 3 are significantly enriched in the pathway of Peroxisome. Table 4 shows the top three significant modules and hub genes in the PPI network.

Figure 7A is the sankey diagram of hub genes, hub gene related miRNAs and circRNAs selected from the circRNA-miRNA-mRNA regulatory network mentioned above. This sankey diagram contains 13 circRNAs,10 miRNAs and 10 hub genes, as well as 24 competing endogenous regulatory relationships. For example, hsa_circ_0023461 might act as a hsa-miR-589-5p sponge and play a regulatory role in acute myocardial infarction by influencing the expression of ATP2B2. In order to further analyze the relative expression of circRNA in this sankey diagram, heat map was used to further demonstrate (Figure 7B).

Receiver operating characteristic curve (ROC) analysis was performed on 10 miRNAs and 10 hub genes in the Sankey diagram and AUC values were calculated. Meanwhile, molecules with AUC ≥ 7 and P < 0.05 in two datasets were considered as potential biomarkers. Finally, hsa-miR-31-3p, hsa-miR-1291, hsa-miR-646, SRSF1 and HNRNPA2B1 were identified as potential markers of AMI. Figures 8A–E shows the detailed results of ROC analysis. GSE31568 (miRNA) and GSE66360 (mRNA) were used to further evaluate the relative expression of these molecules. As predicted by the previous analysis, hsa-miR-31-3p, hsa-miR-1291, hsa-miR-646, and SRSF1 were decreased, while HNRNPA2B1 was increased (Figures 8F–J).

Hub genes (ATP2B2, KCNA1, GRIN2A, SCN2B, GPM6A, CACNA1E, HDAC2, SRSF1, ANK2, and HNRNPA2B1) and Hub gene related circRNAs (hsa_circ_0023461, hsa_circ_0032623, hsa_circ_0030368, hsa_circ_0004561, hsa_circ_0047796, hsa_circ_0092367, hsa_circ_0025522, hsa_circ_0001147, hsa_circ_0004771, hsa_circ_0061276, hsa_circ_0013599, hsa_circ_0045519, hsa_circ_0075503)were further verified in AC16 hypoxia modle. The results showed that 4 Hub genes and 8 Hub gene related circRNAs were differentially expressed. They were hsa_circ_0023461(up), hsa_circ_0030368 (down), hsa_circ_0004561 (up), hsa_circ_0001147 (down), hsa_circ_0004771(down), hsa_circ_0061276 (down), hsa_circ_0045519(up), hsa_circ_0075503(down), ATP2B2(down), GPM6A (down), SRSF1(down), and ANK2(up) (Figures 9A–L). The expression of hsa_circ_0023461, hsa_circ_0004561, hsa_circ_0001147, hsa_circ_0004771, hsa_circ_0061276, hsa_circ_0045519, GPM6A, SRSF1, and ANK2 werre consistent with the predicted result. Therefore, there may be hsa_circ_0001147/hsa-miR-204-5p/GPM6A axis, hsa_circ_0004771/hsa-miR-629-3p/SRSF1 axis, hsa_circ_0061276/hsa-miR-629-3p/SRSF1 axis, hsa_circ_0045519/hsa-miR-298/ANK2 axis, and hsa_circ_0004561 hsa-miR-298/ANK2 axis regulatory relationship involved in myocardial infarction.