High Throughput mRNA Sequencing Reveals Potential Therapeutic Targets of Tao-Hong-Si-Wu Decoction in Experimental Middle Cerebral Artery Occlusion

Background: Experimental and clinical studies have shown that Tao-Hong-Si-Wu decoction (THSWD) improved neurological deficits resulting from Middle Cerebral Artery Occlusion (MCAO). However, the mechanisms of action of THSWD in MCAO have not been characterized. In this study, the mRNA transcriptome was used to study various therapeutic targets of THSWD. Methods: RNA-seq was used to identify differentially expressed genes (DEGs). MCAO-induced up-regulated genes (MCAO vs. control) and THSWD-induced down-regulated genes (compared to MCAO) were identified. Intersection genes were defined as up-regulated differentially expression genes (up-DEGs) identified as MCAO-induced gene expression that were reversed by THSWD. Genes down-regulated by MCAO and up-regulated by THSWD were grouped as another series of intersections. Biological functions and signaling pathways were determined by gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses. In addition, several identified genes were validated by RT-qPCR. Results: A total of 339 DEGs were filtered based on the 2 series (MCAO vs. control and MCAO vs. THSWD), and were represented by genes involved in cell cycle (rno04110), ECM–receptor interaction (rno04512), complement and coagulation cascades (rno04610), focal adhesion (rno04510), hematopoietic cell lineage (rno04640), neuroactive ligand–receptor interaction (rno04080), cocaine addiction (rno05030), amphetamine addiction (rno05031), nicotine addiction (rno05033), fat digestion and absorption (rno04975), glycerophospholipid metabolism (rno00564), and others. The protein–protein interaction (PPI) network consisted of 202 nodes and 1,700 connections, and identified two main modules by MOCDE. Conclusion: Cell cycle (rno04110), ECM–receptor interaction (rno04512), complement and coagulation cascades (rno04610), focal adhesion (rno04510), hematopoietic cell lineage (rno04640), and neuroactive ligand–receptor interactions (rno04080) are potential therapeutic targets of THSWD in MCAO. This study provided a theoretical basis for THSWD prevention of neurological deficits resulting from intracerebral hemorrhage.


INTRODUCTION
Traditional Chinese medicinal (TCM) herbs have been widely used in China and other Asian countries to combat various diseases for thousands of years. It is extremely important to elucidate the molecular mechanisms of TCM using modern advanced technologies (Efferth et al., 2007;Terrano et al., 2010;Youns et al., 2010). Fortunately, high throughput sequencingbased transcriptome analysis provides an approach for rapidly distinguishing altered mRNAs following drug treatment and thus contributes to understanding and interpretation of mechanisms of action of TCM (Liu and Guo, 2011;Suo et al., 2016). Combination therapy with individual herbs to form specific formulas (Babhulkar, 2009) has been used in TCM for approximately 2,500 years with the goal of increasing therapeutic efficacy and reducing adverse effects (Wang et al., 2008). Wang et al. (2017) used RNA-SEQ to study the pharmacological mechanism of Orexin-A in the treatment of stroke. RNA-seq technique was used to study the expression of transcriptome in rat stroke model treated with Mongolian medicine Eerdun Wurile by Gaowa et al. (2018). Tao-Hong-Si-Wu decoction (THSWD), a formula of TCM herbs, has been used to prevent and treat cerebrovascular diseases in clinical practice in China. This formula contains six commonly used herbs, including Semen prunus (Taoren in Chinese, TR), Flos carthami (Honghua in Chinese, HH), Rehmanniae Radix Praeparata (Shudi in Chinese, SD), Radix Angelicae sinensis (Danggui in Chinese, DG), Paeonia lactiflora (Baishao in Chinese, BS), and Rhizoma ligustici Chuanxiong (Chuanxiong in Chinese, CX) (Yeh et al., 2007). In previous studies, THSWD protected against cerebral ischemia through PI3K/AKT and Nrf2 signaling pathways (Han et al., 2015). Another study showed that THSWD could regulate angiogenesis and exert neuroprotective activity against cerebral ischemia-reperfusion injury.
In recent decades, the expression of mRNA performed by microarray or high-throughput RNA-seq has been used to reveal molecular mechanisms and explore biomarkers for diagnosis and prediction, especially in complex diseases such as cancer, diabetes and stroke [10][11][12].
To comprehensively elucidate potential molecular mechanisms, our present study applied an RNA-seq strategy to characterize the key molecular mechanisms of THSWD. At first, we performed RNA-seq in rat brain hemispheres of control, Middle Cerebral Artery Occlusion (MCAO), and THSWD groups. EdgeR and venny were applied to identify up-DEGs and down-DEGs. And then, clusterProfiler, GO, and pathway annotation of up-DEGs and down-DEGs was performed. PPI networks were structured based on STRING database, visualization was carried out using Cytoscape, and modules were identified by MCODE. Finally, key DEGs were verified by qPCR.

Sample Preparation for UPLC-TOF-MS
The proportion of THSWD compound was TR: HH: SD: DG: BS: CX = 3: 2: 4: 3: 3: 2 (Yeh et al., 2007). Each component was purchased from Anqing Huashi Chinese Herbal Medicine Beverage Co., Ltd. (Anqing City, Anhui Province, China). First, the decoction was extracted with 10× water (10 times the sum of the weight of six components) for 2 h. After filtration, slag was extracted with 8× water (the amount of water added to the total weight of the six components) for 1.5 h. The extract was filtered twice and combined with ethanol to sink conditions. Ethanol containing THSWD at sink conditions was diluted with ethanol, and the solution was filtered using 0.22 µm membranes. 2 µL of filtrate was used for further analysis.
Eleven control compounds were obtained, including gallic acid, rehmannioside D, protocatechuic acid, P-hydroxybenzoic acid, hydroxysafflor yellow A, caffeic acid, amygdalin, chlorogenic acid, paeoniflorin, ferulaic acid, and benzoic acid ( Table 1). The 11 control compounds were all provided by the National Institutes for Food. The purity of the products was greater than or equal to 98%.

UPLC-QTOF/MS E Analysis
The aqueous extract of THSWD was separated using a Waters Acquity UPLC BEH C18 column (2.1 mm × 100 mm, 1.7 µm, Waters Corporation, Milford, MA, United States). Mobile phases consisted of 0.1% aqueous formic acid (A) and acetonitrile (B). Chromatographic separation was performed at 35 • C. Gradient elution with the flow rate of 0.3 mL/min was performed as follows: 3% B at 0-2 min, 3-8% B from 2 to 8 min, 8 to 25% B from 8 to 12 min, 25% B from 12 to 15 min, 25 to 45% B from 15 to 16 min, 45 to 90% B from 16 to 22 min, 90 to 100% B from 22 to 26 min, and 100% B from 26 to 28 min. A Waters SYNAPTG2-Si MS (Waters Corporation, Manchester, United Kingdom) with an ESI source in negative ion mode was used for MS analysis, and leucine enkephalin was used as the accurate  mass calibration solution. The desolvation gas temperature was 350 • C. The flow rates of the cone and desolvation gases were set at 50 and 600 L/h, respectively. The capillary, cone, and extraction cone voltages were set to -2.5 kV. MS E was used for MS/MS analysis with a low collision energy of 6 V and a high collision energy of 20-80 V. The scan range was m/z 50-1200. MassLynx TM v4.1 (Waters Co.) and UNIFI R Scientific Information System v1.7 (Waters Co.) were used for data acquisition and analysis.

Animal Experiments and Samples Collection
Male Sprague-Dawley (SD) rats (verified SPF, 180-230 g) were obtained from the Experimental Animal Center of Anhui Medical University. Rats were randomly divided into normal group, MCAO group and THSWD group with 10 rats in each group. All rat anesthetized using sodium pentobarbital during all surgeries to minimize suffering, and preparation of the MCAO model was performed following our previous study (Han et al., 2015), which also conformed to the NIH Guide for animal care and welfare . The procedures of THSWD extraction and medicine administration were described in the material and methods section "Sample Preparation for UPLC-TOF-MS" and our previous research (Han et al., 2015). Four days after MCAO, the left hemispheres were collected and immediately frozen in liquid nitrogen. The study protocol was approved by the Committee on the Ethics of Animal Experiments of Anhui University of Chinese Medicine (2012AH-036-03).

RNA Extraction and RNA-Sequencing
Total RNA was extracted using mirVana TM miRNA Isolation kit (Cat# AM1561, Ambion) following the manufacturer's instructions. Quantity and quality of RNA samples were determined using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, United States) and Agilent 2100 Bioanalyzer

Functional Enrichment Analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched using R package to better understand the functions of DEGs (Altermann and Klaenhammer, 2005). In our study, clusterProfiler was applied to analysis of GO terms and KEGG pathways, and top30 GOs and pathway were presented (Ashburner et al., 2000).

PPI Network Construction and Module Analysis
STRING is a database that provides comprehensive information about interactions between proteins, including prediction and experimental interaction data (Szklarczyk et al., 2011). In our study, the STRING tool was used to generate PPIs among the DEGs, based on interactions with combined scores ≥ 0.4. Then, Cytoscape was used to visualize the network (Shannon et al., 2003). The PPI network was used to filter modules based on the Molecular Complex Detection plug-in (MCOD) in Cytoscape (Bader and Hogue, 2003) with the following conditions: degree cut-off = 2, k-core = 2, node score cut-off = 0.2, and max depth = 100. MCODE score ≥ 4 and node ≥ 10 were considered for functional enrichment analysis of the modules. Moreover, GO and KEGG enrichment for DEGs in the four modules was performed using clusterProfiler.

Validation of Differentially Expressed mRNAs From the Sequencing Profile by qRT-PCR
Real-time PCR was performed to verify the RNA-seq results, and GAPDH was used as the internal control. Relative expression of mRNAs was determined using the 2 − CT method. The following six genes were analyzed: Cdk1, Ccna2, Cdc20, Ndc80, Th, and Calb2. Brain tissues from control, MCAO, and THSWD groups were analyzed (15 samples). Frontiers in Pharmacology | www.frontiersin.org

Statistical Analysis
Raw data are shown as mean ± standard deviation (SD). Statistical calculations were carried out using GraphPad Prism 5 (GraphPad Software, La Jolla, CA, United States). A p-value less than 0.05 was considered to be statistically significant.  (Figure 1 and Table 1). These components were sourced from HH, BS, CX, and TR ( Table 1).

Identification of DEGs
To further understand the multifaceted mechanisms of THSWD in MCAO, we performed RNA-seq to obtain the transcriptomes of control, MCAO, and MCAO+THSWD samples. Using bioinformatics analysis, 339 DEGs were filtered (MCAO vs. control and MCAO vs. THSWD, |log2(foldchange)| > 1 and p < 0.05). Figure 2 and Supplementary Table S1 show that 63 down-DEGs were screened from the intersection of 222 down-regulated mRNAs (MCAO vs. Control) and 125 down-regulated mRNAs (MCAO vs. THSWD) (Figure 1 and Supplementary Table S1). In addition, 276 up-DEGs were identified from the intersection of 362 upregulated mRNAs (MCAO vs. Control) and 785 up-regulated FIGURE 5 | PPI network of the DEGs resulting from THSWD treatment of MCAO. Node sizes correlate with node importance; Pink nodes denote up-regulated genes, while green nodes denote down-regulated genes; PPI, protein-protein interaction; DEGs, differentially expressed genes. mRNAs (MCAO vs. THSWD) (Figure 2 and Supplementary  Table S1).

KEGG Analysis of DEGS Altered by THSWD
Pathway enrichment analysis could provide further information regarding functions of genes and their interactions. ClusterProfiler was applied to GO enrichment analysis with 63 down-DEGs and 276 up-DEGs, which were associated with THSWD treatment of MCAO.  KEGG enrichment analysis showed that a total of 45 pathway terms were enriched with 63 down-DEGs, 15 of which were significant with p < 0.05. The top10 enriched pathways were neuroactive ligand-receptor interaction (rno04080), cocaine addiction (rno05030), amphetamine addiction (rno05031), nicotine addiction (rno05033), fat digestion and absorption (rno04975), glycerophospholipid metabolism (rno00564), pancreatic secretion (rno04972), cholinergic synapse (rno04725), and amyotrophic lateral sclerosis (ALS) (rno05014). The top30 pathway-terms with highest enrichment factors are shown in Figure 4A.
A total of 10 modules were obtained using default criteria by MCODE. Modules were listed in descending order by MCODE score. Two modules with MCODE score ≥ 3 and nodes ≥ 10 were named as module 1 and module 2. Two modules were selected for module network visualization (Figure 6). Most DEGs belonging to these modules were up-regulated in the MCAO model group in both modules.
The Results of qRT-PCR Verification of the DEGs As shown in Figure 7, real-time PCR results showed that Cdk1, Ccna2, Cdc20, and Ndc80 were over-expressed in the MCAO group. Meanwhile, Th and Calb2 showed low expression in the MCAO group as determined by real-time PCR and RNA-seq ( Figure 7A). THSWD reversed expression of these differentially expressed genes (DEGs) to different degrees in MCAO ( Figure 7B). These results verified that the RNAseq results were consistent with real-time PCR results, and the RNA-seq results were reliable. Primers are listed in Table 3.

DISCUSSION
We analyzed the protein-coding mRNA expression profile in cerebral hemispheres from experimental and control rats by high throughput RNA-seq. 1,007 DEGs were identified, and GO, KEGG pathway, PPI, and PPI module analyses were performed, providing a better understanding of the pathological mechanisms of MCAO.

Complement and Coagulation Cascades
Using pathway enrichment analysis, 10 DEGs were enriched into the complement and coagulation cascades, such as C3ar1, Plau, C1qb, C1qc, Vwf, Kng1, Cfd, Vsig4, and C5ar1 (Figure 9). Previous studies reported that complement components C3 and C5 can promote MCAO injury by complement cascade and TLR2/NFκB activation (Conductier et al., 2010;Semple et al., 2010;Xu et al., 2017). Activation of complement and coagulation cascades result in an inflammatory response, which includes degranulation, chemotaxis, and phagocytosis. We found that there were 78 abnormally expression genes in MCAO involved in the inflammatory response. THSWD treatment reversed abnormal expression of 18 of the 78 genes. Our study confirmed these findings, and THSWD inhibited C3ar1, C1qb, C1qc, and C5ar1 expression changes resulting from MCAO. Therefore, we speculated that THSWD can reverse MCAO injury by inhibiting complement and coagulation cascades.

Neuroactive Ligand-Receptor Interaction
Using pathway enrichment analysis, we found 12 DEGs involved in neuroactive ligand-receptor interactions: Brs3, Mc3r, Calcr, Chrna6, Glra1, Galr1, Ntsr1, Cckar, Gabre, C3ar1, Htr2b, and C5ar1 (Figure 10). Lan et al. (2014) found that 5-HT and 5-HT1A receptor expression were suppressed in rats with MCAO, and treadmill exercise could blunt this effect. We found that Htr2a and htr2b were up-regulated in the MCAO group, and Htr2b was suppressed by THSWD. We hypothesized that MCAO could induce htr2a and htr2b overexpression, leading to excitability during the early stage of MCAO, and that THSWD could reverse this dysregulation. Previous studies have shown that GABA-A and glycine receptor expression were suppressed in rats with MCAO, and receptor agonists could reverse this effect (Wang et al., 2007;Costa et al., 2016). We found that some subunits of GABA-B, such as Gabrr3, Gabrq, Gabra6, Gabre, and Gabrr1 were dysregulated during MCAO, and Gabre was modulated by THSWD.

AUTHOR CONTRIBUTIONS
XD and DP conceived and designed the study. XD and LH performed the UPLC-QTOF/MS E analysis. XD, LX, and QB performed the animal experiments. XD, WC, and CP analyzed the data and pre-viewed the manuscript. All authors read and approved the final manuscript.