Studies were performed according to the “Guide for the Care and Use of Laboratory Animals” published by National Institutes of Health (NIH Publications No. 85-23, revised 1996) and with approval of the Animal Care Committee of Beijing University of Chinese Medicine. A total of 45 male Sprague-Dawley (SD) rats with weights of 240 ± 10 g in Specific Pathogen Free (SPF) grade were selected (purchased from Beijing Vital River Laboratory Animal Technology Co. Ltd.).

Total 45 rats were randomly divided into 3 groups (15 rats per group): Sham-operated, AMI model and QSYQ-treated group. Rats in AMI model and QSYQ-treated groups received ligation surgery of left anterior descending (LAD) coronary artery (Wang et al., 2012b). Left thoracotomy between third and fourth intercostal space was performed on rats. After exposing the cardiac tissues, LAD was ligated with a sterile suture (Shuangjian, Shanghai, P.R. China) 1 mm below the left atrium. The thorax was then closed layer by layer. After thoracotomy, rats were warmed on a heated blanket. One day after the surgery, rats in QSYQ-treated group received a dosage of 2.33 g/kg concentrated QSYQ via oral gavage. QSYQ is a combination of six Chinese herbs that includes A. membranaceus (Fisch.) Bunge, S. Miltiorrhiza Bunge, L. japonica Thunb., S. aestivalis Griseb., A. fischeri Rchb. and G. uralensis Fisch. in the ratios of 10:5:3:3:3:2 (Wang et al., 2012b). It was derived by mixing the residue of A. membranaceus (Fisch.) Bunge with all S. Miltiorrhiza Bunge, L. japonica Thunb., S. aestivalis Griseb., A. fischeri Rchb. and G. uralensis Fisch., followed by extraction with hot water (twice, 2 h each). The water extract was then concentrated to form a paste, and the ethanol is added. After 24 h, the filtration is collected to form the final product (Wang et al., 2012b). Quality of QSYQ was evaluated and established by its fingerprints with HPLC-IT-TOF method (shown in Figure S1). The dosage is ideal to exert definite effect for QSYQ on AMI rat model according to our previous study (Li et al., 2014). Rats in AMI model group were only treated with the same volume of water as QSYQ-treated group (Wang et al., 2012a). Rats in Sham-operated group underwent the same procedure like AMI model rats but had no actual ligation of LAD. The mortality rate of surgery was 26.7%. After surgery, there were 2, 6, and 4 rats which died inadvertently in Sham-operated, AMI model and QSYQ-treated groups respectively. Most of these rats died on the day of surgery or the following day, likely due to acute pump failure or lethal arrhythmias. After 28 days, eight rats from each group were assessed by echocardiography following overnight fasting. The cardiac tissues of all rats were harvested for subsequent molecular biology experiment (13, 9, and 11 rats in Sham-operated, AMI model, and QSYQ-treated groups respectively). Among them, 9 rats were selected randomly for digital gene expression sequencing (3 rats per group). Total RNA of cardiac tissues was extracted with TRIzol Reagent (Gibco-BRL, Paisley, UK). The concentration of RNA was measured with Nano Drop 2000 (Thermo Scientific, USA).

The high-frequency, high-resolution digital imaging platform from the Vevo 2100 system (VisualSonics Inc., Canada) was applied to generate 2D M-mode and B-mode echocardiography. After 28 days of QSYQ treatment and before harvesting the cardiac tissues, 8 rats in each group were assessed by echocardiography. The independent operator was from scientific center of Beijing University of Chinese Medicine. LV internal diameter at end-diastole (LVIDd), LV internal diameter at end-systole (LVIDs), ejection fraction (EF) and fractional shortening (FS) were detected by echocardiography (Sharrett et al., 2001; Wang et al., 2015).

LV performance was measured in rats after they were anesthetized with 2% isoflurane as described earlier (Li et al., 2014). A terminal hemodynamic study was performed 28 days after surgery. Nine rats were included in each group. Right carotid artery, LV systolic and end-diastolic, diastolic and mean aortic pressures, Max dP/dt, Min dP/dt, and heart rate were recorded with a system of PowerLab ML880 (AD Instrument, Australia).

Cardiac tissues of ischemic marginal zone were collected after 28 days of QSYQ treatment as described above. After extracting from cardiac tissues, mRNA was purified and reverse transcribed into cDNA with Oligo (dT) magnetic beads. The cDNA was then digested with 2 endonucleases (NlaIII and MmeI) and ligated with adaptors. After linear PCR amplification of cDNA, the library of each sample was sequenced with Illumina HiSeq™ 2000. Dirty tags were filtered out. All clean tags were aligned to rat reference sequences and unambiguous tags were annotated. The number of clean tags belonging to each gene was computed as raw gene expression count, which could be normalized to the total number of transcripts in each sample per million clean tags (TPM normalization). Finally, a well-characterized transcriptome profile of total 10,515 mRNAs was established.

Pearson's correlation coefficient was used as co-expression measure for constructing GCN. For a given coefficient of a gene pair, Student asymptotic p-value was calculated and then adjusted with false discovery rate (FDR). GCN consists of gene pairs with FDR < 0.05.

PageRank centrality (Brin and Page, 2012) is an indicator to identify “hub genes” in GCN. It is the extension of degree centrality, which indicates significance of nodes in a network. PageRank, together with targeted attack on network analysis (Achard et al., 2006), was used to identify keystone QSYQ-regulated genes.

Since the enriched KEGG pathways have close linkages between each other in KEGG map, we merged them to form a background network, using KEGGgraph R package. From this network, all QSYQ-responsive genes (DEGs of comparison between QSYQ-treated and AMI model groups) were extracted to form a sub-network. Functional groups were divided from the sub-network using “leading eigenvector” method (Newman, 2006) with igraph R package. A new strategy was further developed to calculate enrichment degrees of functional groups with keystone QSYQ-regulated genes, and estimate significance levels. There are 3 key steps in this strategy.

Cardiac tissues were lysed and homogenized using RIPA buffer (Applygen Technologies Inc., Beijing, China). Total protein was extracted from the tissue homogenate. Concentration in each sample was detected by BCA assay kit (Applygen Technologies Inc., Beijing, China). Equal amounts of protein were separated with SDS-PAGE and transferred onto NC membranes (Applygen Technologies Inc., Beijing, China). The following primary antibodies were used: MMP-2(1:500 dilution; Abcam, Cambridge, UK), NOS3(1:2,500 dilution; BD Transduction Laboratories, US), ALOX15(1:1,000 dilution; Abcam, Cambridge, UK), CPT2(1:1,000 dilution; Abcam, Cambridge, UK), CD36(1:1000 dilution; Abcam, Cambridge, UK), LPL(1:500 dilution; Abcam, Cambridge, UK), BCKDHA(1:1,000 dilution; Abcam, Cambridge, UK), DBT(1:500 dilution; Abcam, Cambridge, UK), GAPDH(1:10,000 dilution; Abcam, Cambridge, UK). Three rats in each group were detected for 3 biological replicates.

IHC was performed with cell & tissue staining kits (R&D Systems, Inc., USA) and followed by the protocol (Li et al., 2014). 6 rats in each group were incubated with primary antibody (1:200 dilution; Anti-Collagen III antibody, Abcam, Cambridge, UK; 1:200 dilution; Anti-Collagen V antibody, Abcam, Cambridge, UK) for 12 h. 5 visual fields were selected for each sample and analyzed with Image-Pro Plus 6.0 software to determine the integrated optical density (IOD) in the positive stained areas.

Twenty-eight days after surgery, animals were sacrificed after anaesthetization and blood was taken through abdominal aorta. The blood was centrifuged at 8,000 g for 10 min and supernatant was used for detection of plasma indicators. Six rats were included in each group. Plasma levels of total cholesterol (TC), triglyceride (TG), high density lipoprotein (HDL) and low density lipoprotein (LDL) were measured by automatic biochemical analyzer (HITACH17080, Japan). Plasma levels were quantified by comparing with standards of TC, TG, HDL, and LDL from Sekisui chemical Diagnostic reagents (Sekisui chemical co., LTD, Japan).

DEGs were detected by DESeq2 R package, with FDR < 0.10 (for reference of Love et al., 2014) and absolute value of log2 fold change >0.58. KEGG enrichment analysis was performed using “clusterProfiler” R package with P < 0.05. In experimental validation, all data were presented as mean ± standard deviation. Statistical analysis was carried out on one-way analysis of variance (ANOVA) and Dunnetts' test. P < 0.05 were considered statistically significant. Statistical analysis was performed in R 3.1.3 software.

Twenty-eight days after surgery, echocardiography showed down-regulation of EF and FS in AMI model group, with 46.68 and 61.67% decrease respectively as compared with Sham-operated group. This was accompanied by increase in LVIDd and LVIDs (53.14 and 147.49% increments, respectively), which suggested cardiac hypertrophy in AMI development. After treatment with QSYQ, EF, and FS were recovered by 28.53 and 37.20%, respectively. LVIDd and LVIDs of QSYQ-treated group were recovered by 17.55 and 24.85%, respectively (Figure 1A and Figure S2).

As shown in Figure 1B, in AMI model group, 2 indicators reflecting LV systolic function—systolic blood pressure (SBP) and Max dP/dt were reduced by 25.47 and 36.02% respectively, compared to Sham-operated group. In QSYQ-treated group, these indicators were significantly improved (5.8 and 33.61% increments, respectively). Diastolic blood pressure (DBP) and Min dP/dt, which were associated with LV diastolic function, in AMI model group were decreased by 46.47 and 37.16% respectively, compared to Sham-operated group. QSYQ treatment enhanced DBP and Min dP/dt by 63.78 and 30.41% respectively. No difference was observed in heart rate among 3 groups. Echocardiography and hemodynamic measurements demonstrated marked protective effects of QSYQ on improving cardiac function and hemodynamics.

By comparison of AMI model with Sham-operated group, 868 DEGs were identified and regarded as AMI-pathogenic genes, including 418 up-regulated and 450 down-regulated genes (Figures S3A,D). The comparison of QSYQ-treated with AMI model group revealed 2,733 DEGs (Figures S3B,D), consisting of 1,344 up-regulated and 1,389 down-regulated genes. These 2,733 DEGs were considered as QSYQ-responsive genes. Out of these 2 comparisons, 500 genes overlapped and showed reverse differential expression, which were considered as QSYQ-regulated genes. These genes included 221 inversely up-regulated and 279 inversely down-regulated genes under QSYQ treatment (Figures S3C,D). As shown in Figure S3C, QSYQ exerted a significant modulation on global gene expression of AMI toward the normal expression level, which accounting for therapeutic effects of QSYQ on AMI. GCN was constructed with expression profile of 500 QSYQ-regulated genes in 3 groups of samples. PageRank (Brin and Page, 2012) further characterized the centrality of genes in GCN. Targeted network attack analysis (Figure S4) showed that the top 421 genes in PageRank ranking were the hub genes of GCN. Thus, they were considered as the “keystone QSYQ-regulated genes.”

A list of top 20 keystone QSYQ-regulated genes and their biological functions are given in Table S1. It is intriguing to note that most genes are related to extracellular matrix (ECM). ECM abnormality is well known to be associated with VR, which supports well relevance of our method for selection of keystone QSYQ-regulated genes. Expression of these genes was significantly up-regulated in AMI model compared with Sham-operated group. Furthermore, their increased expression was significantly reversed by QSYQ treatment.

Keystone QSYQ-regulated genes were mainly enriched in cardiovascular diseases and cardiac signaling pathways, including dilated cardiomyopathy, hypertrophic cardiomyopathy and β-adrenergic receptor (β-AR)-mediated signaling pathway (Figure S5).

Significance of functional groups was evaluated by our new strategy based on priorities of keystone QSYQ-regulated genes (Materials and Methods). Table S2 demonstrates significance analysis for 15 functional groups. Significant functional groups were identified with adjusted P < 0.005 (Group 1–6 in Table 1 and Figure 2). Some QSYQ-responsive genes (not belong to functional groups), which participate main KEGG pathways of significant functional groups, were also included into analysis to help understand therapeutic effects of QSYQ. Key genes of functional groups and indicators of plasma lipid and lipoprotein levels were validated and showed in Table 2, Table S3, and Figure 3. Results of significant functional group analysis were shown as follows.

Inflammatory response induced by arachidonic acid lipoxygenase (LOX) pathway was activated in AMI and markedly suppressed by QSYQ treatment (Figure S6). Expression of important genes in arachidonic acid LOX pathway were significantly up-regulated in AMI process, including arachidonate 15-lipoxygenase (ALOX15) and matrix metallopeptidase 2 (MMP-2) (Table S3). After QSYQ treatment, their expression was significantly decreased. Western blot of ALOX15 and MMP-2 (Figures 3A,B) further confirmed the regulatory effects of QSYQ on suppression of arachidonic acid LOX pathway in AMI.

Production of nitric oxide (NO) from arginine metabolism catalyzed by nitric oxide synthase 3 (NOS3) was markedly increased after QSYQ treatment compared with AMI model ones (Table S3), which was also observed by Western blot (Figure 3B). Western blot also showed decreased expression of NOS3 in AMI model compared with Sham-operated group.

Blood biochemical analyses were performed to evaluate plasma lipid levels. Compared with Sham-operated, levels of TC, TG, and LDL in AMI model group increased by 72.76, 306.25, and 140.00%, respectively. HDL level decreased by 35.70% in AMI model compared with Sham-operated group. After treatment with QSYQ, levels of TC, TG and LDL decreased by 35.46, 46.15, and 84.37% respectively, as compared with model group. HDL level was up-regulated by 28.49%. These findings suggested that QSYQ ameliorated dyslipidaemia in rats with AMI (Table 2).

As an important pathway to remove lipids from heart, fatty acids oxidation was significantly inhibited in AMI and reversely up-regulated after QSYQ treatment (Figure S7). The affected genes included lipoprotein lipase (LPL), cluster of differentiation 36 (CD36), acyl-CoA synthetase long-chain family member 1 (ACSL1) and carnitine palmitoyltransferase 2 (CPT2) (Table S3). Western blot confirmed significant reverse expression changes of 3 key enzyme genes: LPL, CD36, and CPT2 (Figure 3C).

A large proportion of top keystone QSYQ-regulated genes are related with ECM. Of note, abnormal elevated expression of these ECM-related genes was effectively inhibited by QSYQ treatment. These ECM-related genes include transforming growth factor beta 1 induced transcript 1 (TGFB1I1) and lipopolysaccharide-induced TNF factor (LITAF) (Table S3). ECM abnormality is regarded as a hallmark of VR. Moreover, in functional group 4, it was also revealed inverse expression alterations in ECM genes induced by QSYQ treatment (Figure S8), including collagen type III alpha 1 (COL3A1) and collagen type V alpha 2 (COL5A2) (Table S3). IHC was applied to confirm the inverse alterations of COL3A1 and COL5A2 (Figure 3D), which supported the regulatory effects of QSYQ to ameliorate VR in AMI.

Up-regulation of branched-chain amino acids (BCAAs) degradation by QSYQ treatment was identified in our study (Figure S9). BCAA homeostasis is controlled by mitochondrial branched-chain alpha-keto acid dehydrogenase complex (BCKDC). 4 genes (BCKDHA, BCKDHB, DBT, and DLD) encoding the catalytic subunits of BCKDC, were suppressed in AMI model rats compared with Sham-operated ones (Table S3). After QSYQ treatment, 2 of them were significantly up-regulated: BCKDHA and DBT, which were confirmed by Western blot (Figure 3E). The regulatory effects of QSYQ on BCAAs degradation may mediate cardiovascular contractility and compliance (Details in figure notes of Figure S9). In addition, other energy metabolism including glycolysis (PRKAB1 and ALDOA), TCA cycle (SUCLG1 and SUCLA2) and creatine metabolism (CKM and CrT) were also accelerated by the treatment of QSYQ (Table S3).

Multipronged therapeutic effects of QSYQ were summarized in Figure 4 and Table S4. VR caused by inflammation and dyslipidaemia is considered as a progressive yet irreversible process and the most major pathological manifestation of AMI (Sharrett et al., 2001; Tabas and Glass, 2013). In present study, VR markers, such as COL3A1, COL5A2, and MMP-2, were elevated in AMI progression. QSYQ treatment alleviated VR through counter-acting the aforementioned events. Furthermore, novel therapeutic effects of QSYQ were proposed in our study. QSYQ can inhibit inflammatory response by down-regulating arachidonic acid LOX pathway and elevating production of NO. Meanwhile, QSYQ can ameliorate dyslipidaemia through elevating CD36-CPT2-LPL fatty acid oxidation. Improvement of fatty acid metabolism effectively increases energy supply to cardiac contractility and relaxation in AMI and thus evaluates EF value. In summary, QSYQ concurrently alleviated VR progression, attenuated inflammation induced by arachidonic acid LOX pathway and NO production, and ameliorated dyslipidaemia in the treatment of AMI, thus improved cardiac function and hemodynamics. TCM has upheld the holistic therapeutic philosophy for more than 2000 years. Many Chinese herbal formulae are multi-targeting in the treatment of diseases. Our observations here support a multi-targeting mechanism for QSYQ in the treatment of AMI.

Many traditional strategies for network group significance analysis focus on topological properties. For instance, MCODE algorithm evaluates group significance with a cluster score, to emphasize densely connected clusters (groups). However, density cannot reflect the degree to which the group is relevant with regulatory effects of QSYQ in our study. Our strategy introduces quantitative information (PageRank scores of keystone QSYQ-regulated genes) to character the regulatory effects of QSYQ when analysing the significance of functional group. To some extent, our strategy is a particular enrichment test method which is analogous to Gene Set Enrichment Analysis (Subramanian et al., 2005).

Most enrichment test methods adopt overlap statistics (Berriz et al., 2003; Doniger et al., 2003; Zhong et al., 2004; Subramanian et al., 2005), which only concern the count of overlap between the important genes and genes of a certain pathway or ontology term. However, they do not take advantage of the priorities of important genes and give equal weight to every gene. Compared with these enrichment test methods, our strategy differs in two regards. First, it underscores both the priorities and count of keystone genes in functional group. When genes of a group mainly distribute on the top of priority list of keystone genes, the term P_hit is larger and then ES is larger. Meanwhile, if there are more keystone genes in the group, corresponding P_miss is smaller and then ES is larger. Second, our strategy assesses the significance level of ES by Monte Carlo simulation that sample genes from the background network. This process preserves the enrichment degree of observed functional group with QSYQ-responsive genes, and thus, provides an accurate null distribution. Keystone gene-based group significance analysis has a broad application in evaluation of group significance with priorities of genes.

The reviewer JL declared a shared affiliation, though no other collaboration, with several of the authors YW, CL, LL, and WW to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review. The other 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.