Transcriptomic Analysis Exploring the Molecular Mechanisms of Hanchuan Zupa Granules in Alleviating Asthma in Rat

Objective To investigate the molecular mechanisms of HCZP treatment of asthma. Materials and Methods Thirty Sprague Dawley (SD) rats were divided into normal, asthma, and HCZP groups (n = 10). The asthma model was sensitized by 1 mg ovalbumin (OVA)/aluminum hydroxide Al(OH)3mixture and then challenged with 1% aerosolized OVA for four weeks. Rats in the HCZP group received 10.08 g/kg/d HCZP for four weeks during OVA challenge. Then, lung tissues of rats in each group were collected for RNA sequencing. Moreover, the expression level of some core genes was detected by using western blotting and immunohistochemistry. Results Inflammatory cell infiltration and pathological damage of the lungs improved in the HCZP group. Compared with the asthma group (0.049 ± 0.002 mm2/mm; 0.036 ± 0.006 mm2/mm; and 0.014 ± 0.001 mm2/mm), total wall thickness (0.042 ± 0.001 mm2/mm), inner wall thickness (0.013 ± 0.001 mm2/mm), and smooth muscle layer thickness (0.012 ± 0.001 mm2/mm) significantly decreased in the HCZP group. Bioinformatics analysis showed that hub genes such as bradykinin receptor B2 (Bdkrb2) and CD4 molecule (Cd4) had different expression patterns between model and HCZP groups. Two transcription factors, forkhead box Q1 (Foxq1) and nuclear factor of activated T cells 2 (Nfatc2), served important regulatory roles in asthma. Compared with the model group, Bdkrb2 protein expression increased and Nfatc2 protein expression decreased in the HCZP group. Discussion and Conclusion. HCZP could alleviate asthma via regulating the expression of several hub genes, which might serve as therapeutic targets for asthma. However, the mechanism of these genes will be studied in the future.


Introduction
Asthma, a common chronic respiratory disease in children and adults, is characterized by reversible airway inflammation obstruction and bronchial hyperresponsiveness [1]. Approximately 315 million people worldwide are affected by asthma, among which about 10% have severe or uncontrolled asthma attacks [2]. Despite the advances in diagnosis and treatment, it remains a serious global health problem [3]. e long-term goal of asthma management is to achieve asthma control and minimize the risk of exacerbation [4]. Currently, pharmacotherapy of asthma mainly includes bronchodilators and anti-inflammatory drugs, such as β 2 -agonists and glucocorticoids [5]. However, previous evidence indicates that β 2 -agonists are associated with an increased risk of myocardial infarction, congestive heart failure, and sudden cardiac death [6]; meanwhile, glucocorticoids possess extensive immunosuppressive activity and potentially serious side effects and may further promote human metapneumovirus infection which is capable of eliciting inflammatory responses [7]. us, it is necessary to explore safe and effective drugs for asthma treatment.
Traditional Chinese medicine (TCM), widely reported to be a complementary and alternative therapy for asthma attacks, can alleviate airway hypersensitivity and inflammatory cell infiltration in patients with asthma [8,9].
Transcriptome sequencing techniques are widely used in molecular biology research [14]. Yick et al. [15] explore the cellular and molecular pathways in asthma using transcriptomic analysis (RNA-seq). ey found that asthma group and normal group had different transcriptomic profiles; in addition, genes such as pendrin and periostin were differentially expressed between asthma and normal, which might be relevant for the pathogenesis and treatment of disease. Nevertheless, genes involved in the HCZP treatment in asthma have not been discovered. Nowadays, the application of RNA-seq in the study of TCM attracts more and more attention of researchers. Due to the lack of genomic data and gene sequence information, the medicinal plants need a large amount of genetic information to analyze the gene function in the whole level [16]. At present, the whole-genome sequencing of most herbs cannot be detected, so it is a quick way to compare the gene sequences and identify the expressed genes by constructing the transcription database [17]. For instance, Jiao et al. [18] revealed the anti-inflammation mechanism of Ma Huang Tang on acute bronchial asthma in mice using RNA-seq. erefore, in the present study, we constructed a rat model of asthma to investigate the mechanism of HCZP treatment in asthma. e lung tissues from rat in the normal, model, and HCZP groups were respectively extracted and used for RNA-Seq. Our study may provide new insights into the molecular mechanisms of HCZP treatment and provide candidate biomarkers for targeted therapy of asthma.

Animals.
irty specific-pathogen-free (SPF) male Sprague Dawley (SD) rats weighing 160-200 g were purchased from the Laboratory Animal Center of Hubei Provincial Center for Disease Control and Prevention (Wuhan, China). e rats were kept in the Laboratory Animal Center of Wuhan Hospital of Traditional Chinese Medicine (Wuhan, China) under standard condition (12 h light/dark cycle and 21 ± 1°C) and were acclimated to the conditions for 7 days before the experiment. All animal protocols were approved by the Ethics Committee and Animal Management Committee of Wuhan Hospital of Traditional Chinese Medicine (Approved on 2017-06-02), which conformed to the Guide for the Care and Use of Laboratory Animals published by the US National Institutes of Health [19].

Model Establishment and Drug Treatment.
e rats were randomly divided into three groups: normal, model, and HCZP groups, with 10 rats in each group. HCZP was provided by the Preparation Center of the Base of Traditional Chinese Medicine of Wuhan Hospital of Traditional Chinese Medicine. In this study, the dose of HCZP administered to each rat was selected based on the clinical use of HCZP. Specifically, human doses of HCZP were converted to rats doses according to the human-rat coefficient (0.018) of skin surface area [20]; thus, HCZP was administered intragastrically at a dose of 10.08 g/kg. Rats in model and HCZP groups were sensitized on day 1 and day 8 by intraperitoneal injection with 1 mL of OVA/Al(OH)3 normal saline mixture (containing 10 mg OVA and 100 mg Al(OH) 3 ) as described previously [21], while rats in normal group were intraperitoneally injected with 1 mL of normal saline. Subsequently, the rats in model and HCZP groups were challenged with 1% OVA for 30 min every other day from the 15th day through the 42th day; the rats in normal group were atomized with normal saline from day 15 to day 42. Notably, during the stimulus period (from day 15 to 42), rats in the HCZP were intragastrically given HCZP (10.08 g/ kg) 1 h before atomization, while rats in the model and normal groups were given equal volume of normal saline.
Following a 4-week treatment, rats in each group were given the last challenge after 24 h of withdrawal. Within 2 h after the final challenge, rats were anesthetized with an intraperitoneal injection of 3% pentobarbital sodium. en, the lung tissues were extracted and the residual blood on the surface of the tissues was eliminated using normal saline. e left lung was fixed in neutral formalin and stored at 4°C for hematoxylin-eosin (HE) staining and immunohistochemistry (IHC), while the right lung tissues were stored at −80°C for transcriptome sequencing.

H&E Staining.
After modeling, the pathological changes in the lungs were observed using H&E staining. Lung tissues were fixed overnight in 4% neutral formalin, embedded in paraffin, and cut into 5 μm sections [22]. e sections were dewaxed with xylene and washed with gradient alcohol and distilled water. en, they were stained with hematoxylin and eosin. Finally, sections were sealed with neutral gum and photographed under a microscope (Olympus IX73, Tokyo, Japan). Moreover, H&E images were analyzed by Image-Pro plus software. Briefly, the indicators of bronchial basal perimeter (Pbm), total wall area (WAt), inner wall area (WAi), and smooth muscle area (WAm) were measured. ese measured values were standardized by Pbm, and WAt/Pbm, WAi/Pbm, and WAm/Pbm represented the total wall 2 Evidence-Based Complementary and Alternative Medicine thickness, inner wall thickness, and smooth muscle layer thickness, respectively [23].

RNA Sequencing.
Total RNA from each lung tissue sample was extracted using TRIzol reagent (code no. 9109; TaKaRa, Japan). e 1% agarose gel was used to detect RNA degradation and contamination. RNA purity was monitored using the NanoDrop 2000 ( ermo Scientific). en, RNA concentration and integrity were evaluated using Invitrogen Qubit RNA IQ kit (Q33221, ermo Fisher Scientific in Basingstoke, UK) and Agilent 2100 (Agilent Technologies, Palo Alto, CA, USA), respectively. en, rRNA from the total RNA was depleted using Ribo-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA, USA) according to the manufacturers' instructions. e treated RNAs were interrupted randomly to 200-300 bp by adding ion solution, and then the fragmented RNAs were used as templates to construct the cDNA library. e average insert size for the pairedend libraries was 300-400 bp. After construction of libraries, the concentration and size of library were determined by fluorescence quantification and Agilent 2100 bioanalyzer (Agilent Technologies), respectively. Finally, the libraries were sequenced on the Illumina HiSeq platform [24].

Data Preprocessing.
e sequenced data could contain some low-quality bases and incorrectly sequenced bases. In order to filter out unreliable reads, we used the following steps to control the quality of the raw data: (1) the reads with sequencing joints were removed; (2) when the N content of sequencing read exceeded 10% of the base number of read, their paired reads were removed; and (3) when the low-quality (Q ≤ 5) base number contained in any sequencing read exceeded 50% of the read base number, paired read was excluded [25].

Screening of Differentially Expressed Genes (DEGs).
Raw data were standardized using the TMM algorithm in the edgeR package (version 3.4, http://www.bioconductor.org/ packages/release/bioc/html/edgeR.html) [26] and converted to logCPM value. e differential expression analyses between model vs. normal and HCZP vs. model were performed. en, the corresponding P values and logFC values of genes were obtained.
e P values were adjusted by the Benjamin-Hochberg method (adj. P value) [27]. In addition, genes with adj. P value < 0.05 and |logFC| > 1 were defined as DEGs. Furthermore, DEGs in model vs. normal groups as well as model vs. HCZP were integrated using Venn analysis, and the overlapping genes were thought to be related to HCZP treatment.

Functional Enrichment Analysis of DEGs.
In order to explore the biological function of genes, the functional enrichment analysis was performed.

Construction of Protein-Protein Interaction (PPI)
Network.
e interactions among genes were predicted using the STRING database (version 10.0, http://www.stringdb.org/) [30]. e input gene set was genes related to HCZP treatment and the species selection was Rattus norvegicus (rats). e DEGs were mapped into STRING to obtain PPIs, and a combined score of >0.7 was set as cut-off threshold.

TF Prediction and TF-Target Network Construction.
All the TFs with motifs in the rat were searched via JASPAR 2018 database (http://jaspar.genereg.net/) [34]. en, these TFs were intersected with the HCZP-related DEGs to obtain the differentially changed TFs and corresponding motifs. Furthermore, the binding site analysis of the TFs' motif information and genes' promoter sequence in PPI network were performed using online tool Find Individual Motif Occurrences (FIMO, version 5.0.5, http://meme-suite.org/ tools/fimo) [35]. TF targets with P value < 0.0001 were selected and then visualized by Cytoscape software.

miRNA-Target Interactions Prediction and miRNA-Target
Network Construction. miRNA targets were predicted using four currently available methods from miRWalk 2.0 database (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) [36], including miRWalk, miRanda, miRDB, and Targetscan. e miRNA-target pairs that existed in four databases were selected. Furthermore, miRNAs that simultaneously regulated at least 10 genes were screened, and the miRNA-mRNA network was constructed using Cytoscape software.

Western Blot (WB)
Analysis. WB analysis was performed to verify the results of bioinformatics analysis. e lung tissues of rats were homogenized with radioimmunoprecipitation assay (RIPA) buffer (code no. P0013B, Beyotime, Shanghai, China), which involved the addition of phenylmethanesulfonyl fluoride (PMSF, no. ST506, Beyotime) and centrifuging at 10000 g for 10 min. e supernatant was collected, and the concentrations of protein were determined using the bicinchoninic acid (BCA) assay [37]. e extracted proteins were separated on the SDS-PAGE gel and transferred to polyvinylidene fluoride membrane. e membranes were blocked in 5% skim milk (0.75 g milk powder + 15 mL PBS) at 37°C for 2 h and then incubated with the Bdkrb2 (cat. no. sc136216, Santa Cruz, CA, USA),

Evidence-Based Complementary and Alternative Medicine
Nfatc2 (cat. no. sc7296), and anti-β-actin (cat. no. ab8226, Abcam, Shanghai, China) antibodies at 4°C overnight, respectively. en, the membranes were incubated with secondary antibody at 37°C for 2 h followed by washing with PBST for six times; immunodetection was performed using ECL system. e images were obtained via TanoImage.

Immunohistochemistry.
e 5 μm paraffin slices were used to perform IHC staining. After the sections were dewaxed, the activity of endogenous catalase was inhibited with 3% H 2 O 2 for 10 min at room temperature. e slices were incubated with PBS containing 10% normal goat serum at 37°C for 30 min to reduce nonspecific adsorption, followed by incubation with anti-Bdkrb2 and anti-Nfatc2 antibodies at 4°C overnight. en, sections were washed with PBS and incubated with second antibody at 37°C for 30 min [38]. ereafter, the sections were incubated with SABC at 37°C for 0.5 h and placed in DAB solution for staining [39]. Finally, the slices were sealed with neutral gum and observed and photographed under a microscope.

Statistical Analyses.
All data are presented as means ± standard deviation (SD). Statistical analyses and plotting were performed using GraphPad Prism 5 (GraphPad Software, San Diego, CA, USA). Statistical significance of the differences between two groups was evaluated by using a two-tailed Student's t-test. e results were considered to be statistically significant at P < 0.05.

H&E Staining Analysis.
After treatment, the histopathological alterations of lung tissues of rats in different groups were observed by H&E staining. Compared with the normal group (Figure 1(a)), H&E staining of the lung tissues from rats in the model group indicated that there was inflammatory cell infiltration around the bronchial wall and blood vessels, which was predominantly composed of eosinophils and lymphocytes (Figure 1(b)). In addition, the epithelial cells were obviously hypertrophic, the bronchial tubes narrowed, and the smooth muscles destroyed in the model group. In contrast, HCZP-treated rats showed alleviated symptoms, indicating pulmonary inflammation, goblet cell proliferation, and mucus secretion were significantly reduced (Figure 1(c)). Moreover, WAt/Pbm, WAi/Pbm, and WAm/Pbm in the model group were significantly higher than those in the normal group (normal vs. model: 0.033 ± 0.001 mm 2 /mm vs. 0.049 ± 0.002 mm 2 /mm; 0.011 ± 0.001 mm 2 /mm vs. 0.036 ± 0.006 mm 2 /mm; and 0.005 ± 0.0003 mm 2 /mm vs. 0.014 ± 0.001 mm 2 /mm). By comparison, WAt/Pbm, WAi/Pbm, and WAm/Pbm in the HCZP group were significantly lower than those of the model group (model vs. HCZP: 0.049 ± 0.002 mm 2 /mm vs. 0.042 ± 0.001 mm 2 / mm; 0.036 ± 0.006 mm 2 /mm vs. 0.013 ± 0.001 mm 2 /mm; and 0.014 ± 0.001 mm 2 /mm vs. 0.012 ± 0.001 mm 2 /mm) (Figure 1(d)).

High-roughput Sequencing Data.
Using highthroughput sequencing, more than 20,000,000 clean reads were generated from each sample. e total clean reads ranged from 6.16 to 7.23 G, the Q30 (%) of these samples was about 90%, and the GC content of samples was less than 50%, suggesting that the quality of the sequencing data was high and could be used for subsequent analysis. According to the sequencing data, a total of 16,818 genes were identified from nine samples, accounting for 78.7-89.9% of mapped reads of all rats genes.

Screening of DEGs.
According to the cut-off criterion of P < 0.05 and |logFC| > 1, there were 1,678 DEGs in normal vs. model, including 1,133 upregulated genes and 545 downregulated genes; in addition, 4,461 DEGs were identified between model and HCZP groups, including 1,912 upregulated genes and 2,549 downregulated genes. Moreover, the DEGs in normal vs. model and model vs. HCZP were intersected to identify the genes that were related to HCZP treatment. Finally, a total of 874 DEGs associated with HCZP treatment were obtained (up-or down-regulation of genes was consistent with that in the model vs. normal group).
e expression pattern of HCZP-related genes was basically consistent with that of the normal group, which was different in the model group ( Figure 2). e top five pathways included cell adhesion molecules (CAMs) (rno04514), ECM-receptor interaction (rno04512), antigen processing and presentation (rno04612), type I diabetes mellitus (rno04940), and autoimmune thyroid disease (rno05320) (Figure 3 and Table 1).

PPI Analysis of DEGs.
e identified DEGs were introduced into the STRING database to obtain the PPI pairs, and these PPIs were then visualized using Cytoscape software. PPI network was composed of 282 nodes and 433 edges (Figure 4). e top 10 nodes with degree scores >10 are listed in Table 2. Among these, Bdkrb2 and CD4 might be regarded as hub genes.

TF-Target Regulatory Network and miRNA-Target Network Construction.
To explore the regulatory relationship between TF and genes, the TFs of the 874 DEGs were predicted using JASPAR database. Two differentially upregulated TFs, Nfatc2 and Foxq1, were obtained. eir motifs e transcription factor binding site of Nfatc2 was "TTTTCCA" and of Foxq1 was " * * * * GTTT * * * ." In addition, based on the promoter sequence of the protein corresponding genes in PPI network, the binding sites of TFs were identified using the online tool FIMO. A total of 132 binding sites of Nfatc2 were found, which corresponded to 78 genes. Meanwhile, 85 binding sites of Foxq1 were obtained, corresponding to 48 genes. Subsequently, the TF-target regulatory network was established ( Figure 5(c)).

WB and IHC Analysis.
In order to verify the accuracy of the transcriptome sequencing results, the protein expression levels of two genes (Bdkrb2 and Nfatc2) were detected using WB and IHC. WB showed that the expression level of Bdkrb2 was significantly lower in the model group than in the normal group (0.139 ± 0.004 vs. 0.999 ± 0.095, P < 0.001), while it was markedly higher in HCZP-treated group than in the model group (0.441 ± 0.038 vs. 0.139 ± 0.004, P < 0.001) ( Figure 6(a)). In addition, the level of Nfatc2 expression was significantly higher in the model group than in the normal group (1.806 ± 0.597 vs. 0.999 ± 0.089, P < 0.001), whereas the Nfatc2 level in the HCZP group was observably downregulated compared to that in the model group (1.169 ± 0.081 vs. 1.806 ± 0.597, P < 0.001) (Figure 6(b)).
Next, IHC staining was performed to evaluate the expression levels of Bdkrb2 and Nfatc2. e expression level of Bdkrb2 was significantly decreased in the model group, compared with that in the normal group (0.0006 ± 0.0002 vs. 0.00155 ± 0.00008, P < 0.001). After HCZP treatment, the expression Bdkrb2 was increased (0.00143 ± 0.00006 vs. 0.0006 ± 0.0002, P < 0.001) (Figure 7(a)). Moreover, the expression level of Nfatc2 in the model group was notably higher than that in the normal group (0.0065 ± 0.0007 vs. 0.0033 ± 00004, P < 0.001); however, it decreased significantly after HCZP treatment (0.0031 ± 0.00035 vs. 0.0065 ± 0.0007, P < 0.001) (Figure 7(b)). ese results were consistent with those of RNA sequencing, further confirming the therapeutic effect of HCZP on asthma.

Discussion
Although asthma can be treated with conventional therapies, such as inhaled corticosteroids and bronchodilators, patients with severe asthma still find the condition intractable. Our study investigated the role of potential genes of HCZP, providing new evidence for molecular therapy and the development of novel drugs. In our analysis, a total of 874 DEGs related to HCZP treatment were screened. Functional enrichment analysis showed that these genes were mainly enriched in GO_BP terms such as immune response and negative regulation of T cell proliferation as well as pathways like CAMs. PPI analysis showed that Cd4 and Bdkrb2 were with higher degrees and might be considered as hub genes. We also screened two TFs (Foxq1 and Nfatc2) and four miRNAs (mo-miR-495, mo-miR-742-3p, mo-miR-126b, and mo-miR-672-3p) with higher degrees in the interaction network.
e protein levels of Bdkrb2 and Nfatc2 were verified using WB and IHC, which were in accordance with the results of RNA-Seq. e CD4 gene encodes the membrane glycoprotein of T lymphocytes and is a receptor for human immunodeficiency virus. is protein functions in initiating or enhancing the early stage of T cell activation and may serve as an important mediator of indirect neuronal injury in immunemediated disease [40]. In this study, functional enrichment analysis showed that CD4 was notably associated with T cell activation and adaptive immune response. Smith and Larché [41] indicated that T cell activation in the body might lead to manifestations of chronic allergic inflammation, including bronchoconstriction and hyperresponsiveness. In addition, adaptive immune response, a biological process in which antigen-specific T/B lymphocytes activate, proliferate, and differentiate into effector cells after receiving antigen stimulation, is involved in innate and adaptive immunity [42]. Previous study revealed that the expression of CD4 might be increased during inflammation and thrombosis, altering the immune cell-mediated response and leading to atherosclerosis [43]. In addition, Noble et al. [44] found that CD4 were  Evidence-Based Complementary and Alternative Medicine increased in the lamina propria of inflamed inflammatory bowel disease tissue. Furthermore, Lee et al. [45] reported the relationship between primary immunodeficiency (PID) and asthma, indicating that PID was a significant risk factor for asthma exacerbation. ese evidences indicated that CD4 was closely involved in immune-mediated disease. However, the relationship between asthma and CD4 gene has not been reported. Taken together, we speculated that CD4 might play a role in asthma by affecting immune-related biological functions. Another gene, Bdkrb2, was also linked to the asthma treatment. Bdkrb2 encodes the receptor for bradykinins. Among these, 9AA bradykinin causes of a number of reactions, including vasodilation, smooth muscle spasms, and painful fibrous stimulation [46]. Sabatini et al. [47] found that the Bdkrb2 was abnormally expressed in bronchial fibroblasts from asthmatics, indicating that bradykinins were actively involved in the proliferation and differentiation of bronchial fibroblasts, as well as airway remodeling in asthma through MAPK pathway and EGF receptor transactivation. Meanwhile, Ricciardolo et al. [48] revealed that bradykinin is a vasoactive proinflammatory peptide that could activate plasma and tissue kinin, thereby mediating acute inflammatory response of asthma, such as excessive mucus secretion, smooth muscle contraction, and sensory nerve stimulation. In this study, Bdkrb2 was significantly upregulated in the asthma rats after HCZP treatment, suggesting that HCZP might exert therapeutic effects on asthma by affecting the expression of bradykinins.
Several studies reported the role of TFs in the pathogenesis of asthma. In the present study, two TFs, including Foxq1 and Nfatc2, were closely associated with HCZP treatment. Nfatc2, a member of the activated T cell family, is transferred to the nucleus after being stimulated by the T cell receptor and becomes an active T cell transcription complex, which serves a crucial role in inducing gene transcription during the immune response [49]. A previous study demonstrated that the number of T regulatory cells was increased in Nfatc2-deficient mice, inducing immunosuppression to control experimental asthma caused by allergens [50]. In this study, we observed that the protein level of Nfatc2 was significantly decreased after HCZP treatment, suggesting that HCZP may play a therapeutic role by suppressing Nfatc2 expression. In addition, Foxq1 is a member of the Fox gene family known to be involved in cell cycle regulation, cell signaling, and embryonic development [51]. e relationship between Foxq1 and immune-related diseases has been explored. Ovsiy et al. [52] indicated that Foxq1 participated in the development of chronic inflammatory disease atopic dermatitis by stimulating monocyte movement and increasing proinflammatory potential. However, no studies have reported the roles of Foxq1 in asthma; further studies with biological experiments are needed to illustrate its exact molecular mechanism.

Conclusion
Based on transcriptome analysis, this study demonstrated that HCZP alleviated the symptoms of asthma through various mechanisms. Results revealed that the immune response and T cell differentiation-related genes contributed to the therapeutic effect of HCZP against asthma. Our findings could help identify novel and effective therapeutic targets and provide evidence for a better understanding of asthma pathology.
Data Availability e raw data supporting the conclusions of this manuscript will be made available by the corresponding author, without undue reservation, to any qualified researcher.

Disclosure
Hailong Yin and Yanbo Fan are the co-first authors.