Multipronged Therapeutic Effects of Chinese Herbal Medicine Qishenyiqi in the Treatment of Acute Myocardial Infarction

Background: Based on global gene expression profile, therapeutic effects of Qishenyiqi (QSYQ) on acute myocardial infarction (AMI) were investigated by integrated analysis at multiple levels including gene expression, pathways involved and functional group. Methods: Sprague-Dawley (SD) rats were randomly divided into 3 groups: Sham-operated, AMI model (left anterior descending coronary artery ligation) and QSYQ-treated group. Cardiac tissues were obtained for analysing digital gene expression. Sequencing and transcriptome analyses were performed collaboratively, including analyses of differential gene expression, gene co-expression network, targeted attack on network and functional grouping. In this study, a new strategy known as keystone gene-based group significance analysis was also developed. Results: Analysis of top keystone QSYQ-regulated genes indicated that QSYQ ameliorated ventricular remodeling (VR), which is an irreversible process in the pathophysiology of AMI. At pathway level, both well-known cardiovascular diseases and cardiac signaling pathways were enriched. The most remarkable finding was the novel therapeutic effects identified from functional group analysis. This included anti-inflammatory effects mediated via suppression of arachidonic acid lipoxygenase (LOX) pathway and elevation of nitric oxide (NO); and amelioration of dyslipidaemia mediated via fatty acid oxidation. The regulatory patterns of QSYQ on key genes were confirmed by western blot, immunohistochemistry analysis and measurement of plasma lipids, which further validated the therapeutic effects of QSYQ proposed in this study. Conclusions: QSYQ exerts multipronged therapeutic effects on AMI, by concurrently alleviating VR progression, attenuating inflammation induced by arachidonic acid LOX pathway and NO production; and ameliorating dyslipidaemia.


INTRODUCTION
Acute myocardial infarction (AMI) is the most frequent cause of heart failure (HF) (Mahmood et al., 2014;Maranhão and Tavares, 2015). Although evidence-based therapies have been established, AMI-induced HF remains a leading cause of global mortality (Gaziano et al., 2006). Herbal formulas of traditional Chinese medicine (TCM) have been effectively used in treatment of AMI with few side-effects (Ferreira and Lopes, 2011). One of the commonly used TCM for AMI is Qishenyiqi (QSYQ). QSYQ originates from the recordings of classic medical textbook, "Treatise on Exogenous Febrile Disease." It has been effectively used in treatment of several diseases for thousands of years. It is prepared from 6 herbs, including 2 star herbs: Astragalus membranaceus (Fisch.) Bunge ("huang-qi" in Chinese) and Salvia Miltiorrhiza Bunge ("danshen" in Chinese), and 4 adjunctive herbs: Lonicera japonica Thunb., Scrophularia aestivalis Griseb., Aconitum fischeri Rchb., and Glycyrrhiza uralensis Fisch.
Previous studies indicated that QSYQ exerts its therapeutic effects on AMI through different target pathways in a synergistic manner (Wang et al., 2012a(Wang et al., ,b, 2013Qiu et al., 2014). However, the major limitation of these studies is that the effects of QSYQ on different pathways were not examined simultaneously. TCM is an ancient medical practice system which emphasizes the integrity of entire human body and it usually exerts therapeutic effects via multiple targets or pathways . However, determination of multiple complex pathways and their interactions remains a challenge for TCM research. In recent years, multiple high-throughput omics technologies, including transcriptomics, proteomics, metabolomics and metagenomics, have been increasingly applied in TCM research (Gu and Chen, 2014). To meet the development and application of omics technologies in TCM research, routine bioinformatics strategies, which focus on differentially expressed genes (DEGs) and enriched pathways, have been evolved into advanced stages, which need to integrate multiple omics data, and more advanced network-based methods to identify hub genes and functional groups regulated by drug intervention (Berg, 2014). In our study, we obtained a well characterized highthroughput RNA sequencing expression profile (Zhu et al., 2011), complemented by Western blot, immunohistochemistry (IHC) and measurement of plasma lipids, to identify global gene expression changes induced by QSYQ treatment on AMI. We also developed a new strategy to analyse the significance of functional group based on the priorities of keystone genes, to better understand the molecular mechanisms of these gene expression changes.

Experimental Animals
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.).

Model Preparation of AMI and Grouping
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 . 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).

Echocardiographic Assessment of Left Ventricular (LV) Function
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 endsystole (LVIDs), ejection fraction (EF) and fractional shortening (FS) were detected by echocardiography (Sharrett et al., 2001;Wang et al., 2015).

Hemodynamic Measurements in Rats
LV performance was measured in rats after they were anesthetized with 2% isoflurane as described earlier . 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).

DGE-Tag Gene Expression Profiling
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 TM 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 wellcharacterized transcriptome profile of total 10,515 mRNAs was established.

Gene Co-expression Network (GCN) Analysis
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 QSYQregulated genes.

Keystone Gene-Based Group Significance Analysis
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 QSYQtreated and AMI model groups) were extracted to form a subnetwork. 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 QSYQregulated genes, and estimate significance levels. There are 3 key steps in this strategy.
Step 1: Generation of Priority List for Keystone QSYQ-Regulated Genes PageRank scores in GCN were scaled in the range of 0-1. The minimum of PageRank scores was subtracted from each PageRank score. Each PageRank score was then divided by the range of PageRank scores. A priority list L consists of scaled PageRank scores of keystone QSYQ-regulated genes in subnetwork. Genes in list L were prioritized in decreasing order of scaled PageRank scores.
Step 2: Calculation of Enrichment Score (ES) for Each Functional Group We calculated ES that reflects the degree to which a functional group C is overrepresented at the top of priority list L. Given the list L, ES was calculated by walking down list L. During the walking, when we encounter a gene in group C, we increase a running-sum statistic (hits score "P hit "); when we encounter a gene not in group C, we decrease it (penalty score "P miss "). "P hit " and "P miss " were calculated as formulas (1) and (2), respectively. In the formula (1), s j represents scaled PageRank score of gene g j. In the formula (2), N diff means the count of genes in difference set between group C and list L; N genelist and N group denote the count of genes in list L and group C, respectively.
Finally, ES is the maximum of P hit -P miss . If genes in a group C are randomly distributed across list L, or distributed at the bottom of list L, ES(C) will be relatively low or even negative. On the contrary, if they are concentrated at the top of list L, ES(C) will be high. In our strategy, both count (P miss ) and significance degree (P hit ) of keystone genes were taken into considerations.
Step 3: Estimation of Significance Level of ES by Monte Carlo Simulation For a functional group C, the statistical significance (nominal P value) of ES was estimated over a null distribution, which was generated from 10,000 times of Monte Calro simulations. The main steps are as follows: (a) Randomly sample N C genes from the background network to form a group C random and re-compute ES random . N C is the count of genes in group C, (b) Repeat step (a) for 10,000 times, and create a null distribution of enrichment scores ES NULL , (c) Estimate nominal P value for ES of group C from ES NULL with the positive portion, that is P(C) = #(ES NULL ≥ observed ES)/10,000, which #(ES NULL ≥ observedES) represents the count of ES NULL that is larger than the ES of group C.
The estimated significance levels (P-values) of all the functional groups were adjusted for multiple hypothesis testing with Simes-Hochberg method (Jacquez, 1996). IHC was performed with cell & tissue staining kits (R&D Systems, Inc., USA) and followed by the protocol . 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).

Statistical Analysis
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.

Model Evaluation by Echocardiography and Hemodynamic Measurements
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 Shamoperated 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.

Global Therapeutic Effects of QSYQ on AMI and Identification of Keystone QSYQ-Regulated Genes
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 QSYQregulated 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 Shamoperated 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 FIGURE 1 | The assessments of cardiac function by echocardiography and hemodynamics. (A) Echocardiographic analysis was performed in Sham, AMI, and QSYQ (2.33 g/kg) groups for 28 days after LAD, n = 8 per group. (B) Hemodynamics analysis was carried out in each group. Data shown are the mean ± SD, n = 9 per group. *P < 0.05; **P < 0.01; ***P < 0.001. # DBP and Min dP/dt are negative values.

Significance Analysis of Functional Group Based on Keystone QSYQ-Regulated Genes
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 QSYQresponsive 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.

QSYQ Suppressed Inflammatory Response Induced by Arachidonic Acid LOX Pathway (Group 1)
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.

QSYQ Enhanced NOS3 Expression and Elevated NO Production (Group 6)
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.

QSYQ Ameliorated Dyslipidaemia via Up-Regulating Fatty Acid Oxidation (Group 3)
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).

QSYQ Inhibited Cardiac VR (Top Keystone QSYQ-Regulated Genes and Group 4)
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 lipopolysaccharideinduced 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.

QSYQ Accelerated BCAAs Degradation and Other Energy Metabolism (Analyses of Group 2, 3, and 5)
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).

DISCUSSION
Therapeutic Effects of QSYQ: Anti-Inflammation, Ameliorate Dyslipidemia, and VR Progression

QSYQ Attenuates Inflammation through Suppressing Arachidonic Acid LOX Pathway and Elevating Production of NO
Inflammation mediated by arachidonic acid (AA) metabolism pathway plays a critical role in AMI progression (Levick et al., 2007;Wang et al., 2014).Our study identified significant reverse gene expression changes in AA metabolism, especially in arachidonic acid LOX pathway ( Figure S6). AA is a kind of polyunsaturated fatty acids, which is metabolized by 3 types of enzymes: lipoxygenase (LOX), cyclooxygenase (COX), and cytochrome P-450 (CYP) (Funk, 2001;Moreno, 2009). As shown in Figure S6, in the LOX-catalyzed sub-pathway, ALOX15 and ALOX5 can form hydroperoxyeicosatetraenoic acids (HPETEs), which are further converted to hydroxyeicosatetraenoic acids (HETEs) by glutathione peroxidases. Previous studies have shown that hypoxia increases the proinflammatory enzyme ALOX15 in human carotid plaques (Magnusson et al., 2012). Furthermore, ALOX15 induces production of the reactive signaling molecule 15-HPETE, which is a peroxidised lipid. Of note, increased levels of peroxidised lipids are tightly connected with complex inflammatory diseases such as AMI (Spiteller, 2006). ALOX15 also catalyzes production of HETE. HETE enhances the adhesion of leukocytes to endothelium by activating chemokine production, which is an early event in development of AMI (Hedrick et al., 1999). Ye et al. (2004) has reported that 5-HETE induced the expression of MMP-2, which supports our results. It is noted that matrix metalloproteinases (MMPs) are involved in the degradation of matrix components and contribute to VR progression in AMI. Thus, our results demonstrated that the enhanced expression levels of ALOX15 and ALOX5, contribute to accumulation of pro-inflammation effects by HPETEs and HETEs, and further exacerbate VR progression induced by MMPs. On the other hand, increased expression of ALOX5AP accelerates the accumulation of LTA4 and thus enhances AMI. Interestingly, QSYQ treatment effectively suppresses inflammation by downregulating the key genes including ALOX15 and ALOX5AP in arachidonic acid LOX pathway, and finally inhibits VR indicators, including MMP-2 and MMP-23. Meanwhile, NO synthesized from L-arginine by NOS3 can reduce inflammation and oxidative stress through inhibiting lipid peroxidation and scavenging superoxide anion (Katakami et al., 2014), and thus alleviates VR in AMI progression.
The initial step of fatty acid metabolism is hydrolysis of TG. LPL is the major enzyme responsible for hydrolysis of TG. Moreover, LPL activity is positively correlated with HDL levels (Blades et al., 1993;Tornvall et al., 1995). As an AMI Risk in Communities Study shown, increased AMI risk has strong associations with TC, LDL-C, and TG; while decreased AMI risk is closely related with HDL-C (Sharrett et al., 2001). It was demonstrated in the study that QSYQ decreases plasma TC, TG, and LDL levels and up-regulate LPL expression and plasma HDL levels, which effectively ameliorated dyslipidaemia and decreased AMI risk.
After hydrolysis of TG into fatty acids, the next step is the transportation of fatty acids from extracellular environment to cytoplasm. CD36 as a fatty acid translocase, binds and transports long chain fatty acids. Previous study (Krzystolik et al., 2015) has suggested protective effect of higher soluble CD36 in coronary artery disease patients. Meanwhile higher soluble CD36 concentration is also associated with lower risk of LV hypertrophy (Krzystolik et al., 2015). Thus, overexpression of CD36 induced by QSYQ is associated with lower risk of LV hypertrophy.
Once transported into cytoplasm, fatty acids participate in β-oxidation in mitochondria. It is noted that ACSL1, CPT2, and ACADM genes, which encode key enzymes of fatty acid β-oxidation, were significantly suppressed in AMI progression and reversely enhanced after QSYQ treatment. ACSL1 is highly expressed in oxidative tissues like brown adipose tissue and heart (Ellis et al., 2011). Activated fatty acid normally provide 60-90% of the substrate used by the heart for energy, but when ACSL1 is absent, fatty acid oxidation decreases more than 90% (Ellis et al., 2011). CPT2 together with carnitine palmitoyltransferase 1 (CPT1), controls the rate of long-chain acyl-CoA transporting into mitochondria from cytoplasm.

QSYQ Alleviates VR Progression
VR markers were elevated in AMI progression. QSYQ treatment suppressed expression of VR markers as our results shown. In atherosclerotic lesions, levels of ECM components, particularly fibrillar collagen such as COL3A1 and COL5A2, are elevated (Tyagi, 1999), which were validated in our study. Furthermore, we identified main causes of VR, inflammation and dyslipidaemia were both ameliorated by QSYQ treatment.

Overview of Multipronged Therapeutic Effects of QSYQ on AMI
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.

New Strategy-Keystone Gene-Based Group Significance Analysis
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 QSYQregulated 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.

CONCLUSIONS
Based on comprehensive transcriptome analyses and experimental validation, we identified therapeutic effects of QSYQ on AMI at multiple biological levels, including gene expression, pathways involved and functional group. We concluded that QSYQ concurrently alleviated VR progression, attenuated inflammation induced by arachidonic acid LOX pathway and NO production, and ameliorated dyslipidaemia in the treatment of AMI. Moreover, our study provided a new strategy to analyse the significance of functional group based on the priorities of keystone genes, which has a broad application in pharmacological studies.

AUTHOR CONTRIBUTIONS
Each author has contributed significantly to the submitted work. WW and RZ conceived and designed the experiments. YW, WL, and CL performed the experiments. YW, WL, CL, SS, GJ, LZ, and LL analyzed the data. YW, WL, CL, SS, GJ, LZ, LL, RZ, and WW wrote the paper. All authors read and approved the final manuscript.