The microarray expression data of peripheral blood samples was obtained from GEO (http://www.ncbi.nlm.nih.gov/geo) and GEO Series accession number GSE48060. This microarray profiling was acquired from human peripheral blood samples of 31 AMI patients and 21 controls. 52 samples were profiled using the chip-based platform GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Suresh et al., 2014).

The use of gene expression data from whole blood instead of heart tissue samples may limit the ability of our approach to detect signatures in the specific diseased tissue. However, peripheral blood is more easily accessed than heart tissue and may provide specific insights into the immune and inflammatory response of AMI exacerbations.

Microarray data analysis was performed using R software and Bioconductor 3.3 (http://www.bioconductor.org/). By quantile normalization, base 2 logarithm conversions, background correction and normalization were performed on original expression data by Robust Multiarray Averaging (RMA) (Gautier et al., 2004); as a result 54,675 mRNAs were carried over for further analysis (Supplement Table 1).

We screened the significant DEGs from 31 AMI patients and 21 controls by Limma (linear models for microarray data) package of R 3.1.3 software (Sanges et al., 2007). The Fold change > 1.2 and P <0.05 was regarded as a threshold.

Cluster and TreeView are software programs for analyzing and visualizing DNA microarray data or other genomic data sets. In this study, Cluster analysis was used to differentially-expressed genes data and then visualize the data through TreeView.

GO analysis was utilized to explain the primary function of the DEGs according to the GO database, which is the crucial functional classification of NCBI (Ashburner et al., 2000; Gene Ontology, 2006). Fisher's exact test was used to calculate the significance level (P-value) of each GO term to screen out the significance GO term of the DEGs enrichment. While P <0.05 was considered to be significant.

Analogously, pathway analysis was performed to discover the significant pathway terms of the DEGs according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000). We used the Fisher's exact test to identify the significant pathway terms; a P <0.05 was considered significant (Kanehisa et al., 2004; Draghici et al., 2007).

The GO hierarchy was a directed acyclic graph (DAG), in which each term had a defined link with one or more other terms. A GO-Tree was constructed based on the GO DAG to supply clear data navigation and visualization. P <0.01 were selected for statistical significance of GO terms in GO analysis, and the GO-Tree was constructed using the up- and down-regulated DEGs in order to outline the functions that influenced AMI (Zhang et al., 2004).

Experimental genes also participated in many significant GO analyses. We organized the mutual regulation and affiliations between all GO analyses into the database based on a hierarchy of GO. By building a functional relationship network, it was possible to summarize the impact of the experimental function groups, as well as internal affiliation significant features. We considered the biological process (BP) terms of GO analysis to build a network of functional regulation (P <0.01).

The KEGG database included signal transduction, metabolism, cell cycle, membrane transport pathways and information of their interactions. The genes we selected may be related to two or more signaling pathways. Due to the same genes in different pathways, overlapping between pathways was inevitable. We selected the genes in enriched biological pathways and used Cytoscape for graphical representations of pathways (Shannon et al., 2003). By constructing the pathway action network, the signal transfer relations between the pathways was explored at a macro level. In order to look for the core pathways and regulatory mechanisms affected by the disease from multiple significant pathway terms, P <0.05 of pathway terms in pathway analysis were selected to construct a pathway action network.

Gene co-expression network analyses (Kim et al., 2001) were carried out to reveal the interrelation between the DEGs, based on their normalized signal intensity in “AMI” and “Control” profiles. The Pearson correlation was calculated for each pair of genes, and the significantly correlated pairs were used to construct the co-expression network (Prieto et al., 2008). To trace the core regulatory genes in the networks, k-core scores were applied to simplify the analysis of graph topology. The k-core of a given gene indicated its core or nodal status linking other genes in the network to “k” (Ravasz et al., 2002; Barabási and Oltvai, 2004). Consequently, we applied k-core scores to identify genes with the highest networking degrees as key regulatory genes.

WGCNA is a common algorithm for constructing a co-expression network. We utilized the similarity of gene co-expression to define a network. When using S_mn to represent the gene co-expression similarity between genes m and n, we can then apply the power adjacency function for correlating adjacency of genes: S_mn = |(1+cor(xm,xn))/2|β. In data processing, the genome-wide gene expression data was preliminarily filtered, followed by measuring the consistency of gene expression profiles by Pearson correlation, and finally using the power adjacent function to Pearson correlation matrix, data was transformed into weighted gene co-expression networks.

Intramodular connectivity (IC) describes the correlation degree of a gene with other genes in a given module, which can be understood as a measure of module membership (MM). Network module (Module) is a cluster of closely interconnected genes. During module detection, the adjacency matrix (which is a measure of topology similarity) is transformed into the topological overlap matrix (TOM), and modules are detected by cluster analysis (Zhao et al., 2014). We then analyzed the importance of genes by the t-test to determine whether the modules were associated with myocardial inflammation. Module eigengene (ME) refers to the first principal component gene of module expression matrix. It is considered the most representative of the module genes, which has important biological significance.

The online tool PANTHER was used to analyze the functional enrichment of network module genes that are associated with myocardial inflammation (P <0.05) (Mi et al., 2013).

Hub genes are described as the genes most closely associated with disease; they are connected to the highest degree of a series of genes in a module. To a certain extent, they are used to determine the character of modules, the hub genes of modules often have more biological significance compared with the hub genes of global networks (Jeong et al., 2001; Goh et al., 2007; Liang and Li, 2007). A gene is considered to be a hub gene if it has a unique character, such as high GS, high MM and high IC in the network (Zhang and Horvath, 2005). Gene significance (GS) showed the different IC of a gene in various networks, and MM described the significance of genes in the network. The IC of a gene suggested the connectivity with the other genes within network. We were therefore also able to identify the hub genes in modules through the GS, MM, and IC.

JASPAR (http://jaspar.genereg.net) is an open-access database storing curated, non-redundant transcription factor (TF) binding profiles representing transcription factor binding preferences as position frequency matrices for multiple species in six taxonomic groups (Mathelier et al., 2016). Cytoscape is a visualization data integration package for biological networks based on the Java language (Smoot et al., 2011). We imported DEGs into the JASPAR database to obtain the interaction between transcription factors and their target genes (TG). Then, we used Cytoscape software to visualize these relationships and finally obtained the transcription factor-target gene regulatory network.

The research object of this study with a total of 20, all are hospitalized patients with AMI in the department of cardiology of first hospital of Jilin university, and each was invited to join the study patients are required to sign a consent form. Our research protocol was approved by the ethical committee of the First Hospital of Jilin University. Among the 20 patients, 12 males, and 8 females, aged 57 ± 9 years. Inclusion criteria: chest pain or distress within 24 h duration > 30 min, and myocardial enzymes CK-MB and cTNT higher than the normal range online. While excluding myocarditis and other diseases caused by chest pain or distress; exclude the patients with history of renal failure; patients with advanced liver disease; patients with malignant tumors and patients with other inflammatory diseases, such as psoriasis, rheumatoid arthritis and so on.

According to the microarray results, the 8 most dysregulated mRNAs were chosen for further validation by RT-qPCR in AMI patients vs. healthy controls. Blood samples were collected and PBMCs were isolated. Total RNAs were isolated from PBMCs using Trizol reagent (Invitrogen), and complementary DNA was synthesized using the TransScript®Frist-Strand cDNA Synthesis SuperMix(Transgen, China)according to the manufacturer's instructions. Primer sets for selected genes were designed by Sangon Biotech (Shanghai, China); their sequences and reaction conditions are available in Table 1. Each sample was run in triplicate in 96-well plates using LightCycler®96 and FastStart Essential DNA Green Master (Roche Diagnostics GmbH, Germany). Quantification cycles (Cq) were calculated using the fit point method (LightCycler®96 Software, Version 1.1 provided by Roche). The expression data were normalized to the reference glyceraldehyde-3-phosphate dehydrogenase (Gapdh). All experiments (sample collection, preparation and storage, primer design, qPCR normalization) were performed according to the MIQE guidelines.

We used R software and corrected batch effect by Bayesian methods, then removed the probes without corresponding annotation information. To determine the expression values of each gene, we used multiple GSE48060 probes corresponding to the median expression value of that gene, and finally obtained the expression profile of 54,675 genes (Supplement Table 1).

By using algorithms provided in the Limma package, we calculated the data and obtained with lists of differentially expressed genes. A total of 551 DEGs were identified between the AMI group and the control group, including 164 upregulated and 387 down-regulated genes (P <0.05, Fold change > 1.2, Figure 1A, Supplement Table 2). Hierarchical clustering analysis was obtained for the 551 DEGs from the 52 samples of the AMI and control groups. The general gene expression patterns were evidently different in the two groups by TreeView (Figure 1B).

To reveal AMI–related biological processes, we conducted a functional enrichment analysis. GO enrichment analysis of 551 DEGs was performed. The most GO terms of biological processes were associated with inflammation, including the inflammatory response, immune response, cellular defense response and chronic inflammatory response in the AMI group (Figure 2A). In the cellular component category, enriched GO terms were mainly associated with lipid particles, the external side of the plasma membrane and the extracellular space. In the molecular function category, GO terms enriched for DEGs in AMI included chemokine activity, GTP binding and carbohydrate binding (Table 2). We constructed gene function regulatory networks (GO-Tree) for the significant GO terms (P <0.01) of the category of biological processes in GO analysis in order to explore the intrinsic link among gene function (Figure 2B). We then found the hierarchical tree relationships between gene function significantly (Supplement Table 3). The analysis showed that the gene function cascade eventually induced inflammation during AMI.

To screen the significant enrichment of DEGs in pathway terms, Fisher's exact test was used to calculate the significance level of the pathway (P <0.05). According to KEGG databases, we performed pathway annotation of DEGs and obtained DEGs involved in all pathway terms (Figure 2C, Table 3). The DEGs of AMI were enriched mainly in inflammatory bowel disease (IBD), the chemokine signaling pathway, natural killer cell mediated cytotoxicity, toll-like receptor signaling pathway and cytokine-cytokine receptor interactions. As there was cross-talk between the pathway terms, in order to search the transitive relation between signaling pathways, we constructed pathway action networks (pathway-act-network) for the significant pathway terms (P <0.05) (Figure 2D, Supplement Table 4). There were interactions in multiple pathways, demonstrated by the pathway-act-network, which ultimately impacted activation of multiple pathway terms including the JAK-STAT signaling pathway, toll-like receptor signaling pathways and antigen processing and presentation.

We performed a gene co-expression network analysis to investigate the phenotypic change of genes associated with AMI (Kumari et al., 2012; Villa-Vialaneix et al., 2013). Twelve genes from the AMI and control groups were chosen as key regulatory genes (|Dif_Kcore| > 30), and it was shown that there was a significant change in the expression pattern of genes in the AMI group (Supplement Table 5), namely CNN2 (calponin 2), CRYZ (quinone oxidoreductase), SULT1A1 (sulfotransferase 1A1), SULT1A2 (sulfotransferase 1A2), PRMT2 (protein arginine methyltransferase 2), ATP1B1 (ATPase, Na⁺/K⁺ transporting, beta 1 polypeptide), CX3CR1 (CX3C chemokine receptor 1), GCH1 (GTP cyclohydrolase 1), INSIG1 (insulin-induced gene protein), CXCL5 (C-X-C motif chemokine 5), GBP3 (guanylate-binding protein 3) and HEG1 (heart development protein with EGF-like domains 1).We hypothesized that the 12 key regulatory genes are likely to be closely related to the occurrence and development of AMI (Figure 3).

However, traditional studies of gene co-expression patterns mainly utilize the Pearson coefficient to describe the correlation between genes, and tend to determine the presence of a co-expression network using the Pearson coefficient and FDR threshold (also known as hard threshold) (Butte and Kohane, 2000; Carter et al., 2004; Prieto et al., 2008). The two genes are considered to be connected when the correlation coefficient of two genes is equal to or greater than this threshold (for example, 0.8). However, there is a significant limitation to this method, i.e., we have no evidence to identify whether the coefficient 0.8 or 0.79 of two genes has a significant difference. The above algorithm cannot avoid this situation. WGCNA is different from the traditional gene-gene correlation coefficient matrix, and introduces a soft-threshold method that avoids this limitation. In WGCNA, the correlation coefficient in genes of a gene co-expression matrix is obtained by weight calculation. The standard of weighting is the connections between genes, which should conform to the distribution of scale-free networks, obeying a power-law distribution. In short, the correlation coefficient of each gene pair was the β-th power of the exponent operation, which was the weight calculation and β was known as the soft threshold.

We selected a suitable weighted parameter of adjacency function, which is the soft-threshold β, before constructing the weighted co-expression network. After the calculation, we selected the correlation coefficient close to 0.8, soft-threshold β = 9 to construct gene modules using the WGCNA package.

After determining the soft threshold, a total of 551 DEGs were used to construct weighted gene co-expression networks. According to the basic idea of WGCNA, we calculated the correlation matrix and adjacency matrix of the gene expression profile of the AMI and control groups, and then transformed them into a topological overlap matrix (TOM), and obtained a system clustering tree of DEGs on the basis of gene-gene non-similarity (dijω = 1− ω_ij). Together with the TOM, we performed the hierarchical average linkage clustering method to identify the gene modules of each gene network (deep split = 2, cut height = 0.99). In both the AMI and control groups, a total of five gene modules were recognized by the dynamic tree cut (Figure 4A). Genes that not belong to any modules were housed in the gray module. The gray gene modules were ignored in this study.

For a preliminarily test to assess whether the network was acceptable, we derived the corresponding genes of each module, and PANTHER was used for the functional enrichment analysis of GO terms, P <0.05 (Mi et al., 2013). For the blue module, genes involved in the inflammatory response (GO: 0006954, P = 1.81E-07), immune response (GO: 0006955, P = 1.54E-05), defense response (GO: 0006952, P = 5.12E-05) and immune system process (GO: 0002376, P = 1.35E-03) were significantly enhanced (Table 4). It is suggested that this module is closely related to the occurrence and development of the immune and inflammation responses in AMI. Meanwhile, by analyzing the brown module, it is evident that the module gene is mainly concentrated in the immune response (GO: 0006955, P = 8.83E-06), regulation of lymphocyte chemotaxis (GO: 1901623, P = 1.05E-05), defense response (GO: 0006952, P = 2.03E-05), response to stimulus (GO: 0050896, P = 9.99E-05) and immune system process (GO: 0002376, P = 1.53E-04) (Table 4). These results indicate that the module may be closely related to the immune response during AMI. The gene modules that we obtained have biological significance, so we named the blue and brown modules as immune modules. The functional enrichment analysis of blue and brown modules indicates that the immune response and inflammation would immediately follow the onset of AMI, which is consistent with previous reports (Fang et al., 2015).

A major goal of our research was to analyze the degree of correlation between genes and disease, and to ascertain the significance of genes in the corresponding modules. We defined GS as the association of a gene with immune and inflammatory responses. Then we calculated the mean GS of all genes in the module, which is module significance (MS), and described subsequently the relationship between modules and the inflammatory response (Langfelder and Horvath, 2008). Figure 4B showed the correlation between the various modules and immune/inflammatory responses, which shows that the most relevant modules to the immune and inflammatory responses were the blue and brown modules. For a particular gene, to evaluate its involvement in the given module, we analyzed the correlation between its expression profile and ME profile in all samples. We calculated the MM of all the genes in the corresponding module to learn their significance within the given modules. We identified modules (blue and brown) that showed a high correlation with immune and inflammatory reactions, and GS and MM of each gene are shown in Figures 4C,D.

The above results suggested that the occurrence and development of AMI were closely related to the immune and inflammatory responses. To screen genes that were most relevant to immune and inflammatory responses, we constructed the gene co-expression network of immune module genes. There was a positive correlation between MM and IC in the modules, but not GS. Therefore, we selected the hub genes according to the IC and MM values of each gene module (Table 5). Eight hub genes were selected, including NCF4, AQP9, NFIL3, DYSF, GZMA, TBX21, PRF1, and PTGDR. Combined with previous reports, we hypothesized that these hub genes were closely relevant to immune and inflammatory responses when AMI occurs. The co-expression network of module genes was then constructed to obtain the degree of connection between genes in the module (Figures 5A,B).

Based on JASPAR database, the transcription regulatory network of DEGs was constructed and included 66 nodes (Figure 6A). There were 30 transcription factors, and the core transcription factors were SATA4, BCL11B, GATA3, and TBX21, which have the regulatory relationship with the most target genes. (TG > 10, Supplement Table 6). Combined with the results of WGCNA analysis, we found that the transcription factor TBX21 was a hub gene. We screened all target genes associated with the transcription factor TBX21 and constructed a transcriptional regulatory network centered around TBX21 (Figure 6B). Significantly, TBX21 was directly related to the transcription factors SATA4, BCL11B, and GATA3.Interestingly, these four transcription factors were also directly related to PRF1, and PRF1 was also a hub gene.

RT-qPCR was used to validate the microarray data. To verify the main conclusion drawn from the microarray results for peripheral blood samples, the expression levels of genes coding for NCF4, AQP9, NFIL3, DYSF, GZMA, TBX21, PRF1, and PTGDR were determined. Overall, the RT-qPCR results were qualitatively consistent with the results of the microarray analysis. However, the RT-qPCR analysis tended to give higher up-regulation levels than those calculated from the microarray data. The highest change was found for TBX21, PRF1, and NFIL3 (Figure 7).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.