# RNA-Seq Expression Analysis of Chronic Asthmatic Mice with Bu-Shen-Yi-Qi Formula Treatment and Prediction of Regulated Gene Targets of Anti-Airway Remodeling

Authors: Jie Cui; Zexi Lv; Fangzhou Teng; La Yi; Weifeng Tang; WenQian Wang; Wuniqiemu Tulake; Jingjing Qin; Xueyi Zhu; Ying Wei; Jingcheng Dong
Journal: Evidence-Based Complementary and Alternative Medicine
 (2021)

---

## Abstract

Airway remodeling is one of the typical pathological characteristics of asthma, while the structural changes of the airways in asthma are complex, which impedes the development of novel asthma targeted therapy. Our previous study had shown that Bu-Shen-Yi-Qi formula (BSYQF) could ameliorate airway remodeling in chronic asthmatic mice by modulating airway inflammation and oxidative stress in the lung. In this study, we analysed the lung transcriptome of control mice and asthmatic mouse model with/without BSYQF treatment. Using RNA-sequencing (RNA-seq) analysis, we found that 264/1746 (15.1%) of transcripts showing abnormal expression in asthmatic mice were reverted back to completely or partially normal levels by BSYQF treatment. Additionally, based on previous results, we identified 21 differential expression genes (DEGs) with fold changes (FC) > (±) 2.0 related to inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPK signal pathways which may play important roles in the mechanism of the anti-remodeling effect of BSYQF treatment. Through inputting 21 DEGs into the IPA database to construct a gene network, we inferred Adipoq, SPP1, and TNC which were located at critical nodes in the network may be key regulators of BSYQF's anti-remodeling effect. In addition, the quantitative real-time polymerase chain reaction (qRT-PCR) result for the selected four DEGs matched those of the RNA-seq analysis. Our results provide a preliminary clue to the molecular mechanism of the anti-remodeling effect of BSYQF in asthma.


---

## Body

## 1. Introduction

Asthma, a chronic inflammatory respiratory disease, is caused by complex factors and affects approximately 300 million people of all ages worldwide [1]. The pathogenesis of asthma is characterized by airway inflammation and airway hyperresponsiveness (AHR), and airway remodeling is the main pathophysiological feature [1, 2]. Airway remodeling is usually observed in asthma and can be driven by pathways that are partly independent of airway inflammation. It includes airway smooth muscle (ASM) mass, peribronchial collagen deposition, goblet cell, and glandular hyperplasia and angiogenesis [3, 4]. The progression of asthma is associated with the functional changes of the airways [5], whereas the complexity of airway function and structure changes in asthma has impeded the development of novel asthma targeted therapy.Asthma-related morbidity results from immune imbalances caused by the release of inflammatory mediators including reactive oxygen species (ROS), cytokines, and growth factors [4]. ROS could induce cell damage, change the physiological function of structural cells, and play a key role in initiation as well as amplification of inflammation in asthma [6, 7]. Meanwhile, accumulating evidence has suggested that ROS is involved in airway remodeling in asthma. ROS can enhance release of transforming growth factor-β1 (TGF-β1) and vascular endothelium growth factor (VEGF), which contribute to subepithelial fibrosis and airway remodeling [6]. Mitochondria are one of the important sources of basal endogenous ROS production in the lung [8]. Furthermore, airway abnormalities, particularly secretion of various cytokines, leads to the activation of intracellular signaling pathways which are involved in airway inflammation as well as remodeling process. The PI3K/AKT pathway has been proved to play an important role in regulating cell proliferation, growth, differentiation, and metabolism [9]. As well, the mitogen-activated protein kinase (MAPK) pathway is known to regulate a variety of biological processes like cell growth and proliferation, chemotaxis, degranulation, and other processes [10]. Therefore, PI3K/AKT and MAPKs signal pathways are considered as therapeutic targets in airway remodeling of asthma.Traditional Chinese medicines (TCMs) were derived from thousands of years of clinical use in China. TCMs containing multiple bioactive ingredients are potential novel resources for asthma treatment drugs. BSYQF is used in the clinical treatment of asthma in China and is composed ofAstragalus membranaceus (Fisch.) Bunge, Rehmannia glutinosa Libosch, and Epimedium brevicornu Maxim. The chemical fingerprint of BSYQF contains at least 16 chemical components including icariin, acteoside, catalpol, leonuride, calycosin, and epimedin A. This indicates that multiple compounds in BSYQF may deliver an integrated anti-asthma effect through multiple targets and their associated molecular pathways.The anti-remodeling effects of BSYQF including the inhibition of ASMC proliferation and peribronchial collagen deposition in chronic asthmatic mice have been demonstrated in our previous study [11]. We also preliminarily explored the anti-remodeling mechanism of BSYQF, which may be related to its anti-inflammatory, antioxidant, and mitochondrial structure restoration. To further investigate the molecular mechanism of BSYQF, we used lung samples obtained from previous experiments and took advantage of high-throughput whole-transcriptome analyses to explore molecular mechanisms targeted by BSYQF. BSYQF is a multicompound formula with wide molecular mechanisms. Given this complexity of multiple targets and their associated molecular pathways, increased understanding of the mechanisms of BSYQF might be best derived from the available experimental results. Hence, according to our previous result, we mainly focused on inflammatory, oxidative stress, and mitochondria in GO enrichment and PI3K/AKT and MAPKs signal pathways which are thought to play an important role in asthma airway remodeling to identify differentially expressed genes (DEGs) which may be therapeutic targets of BSYQF.


## 2. Materials and Methods

### 2.1. Animals

Lung tissue samples were obtained from our recently published study that tested the effects of BSYQF on airway remolding in murine chronic asthma [11]. Female BALB/c mice were sensitized and challenged with ovalbumin for 8 weeks to establish chronic asthmatic model as described previously [11]. BSYQF was prepared as described previously [12]. Briefly, Astragalus membranaceus (Fisch.) Bunge, Rehmannia glutinosa Libosch, and Epimedium brevicornu Maxium were decocted in the ratio of 3 : 2 : 1.5(w/w). From day 14, mice in BSYQF group were oral administrated 20g raw herbs/kg body weight, and mice in control and asthma groups were oral administrated with saline. At the end of experimental period, mice were euthanized and lung tissues used to extract RNA were harvested immediately and frozen in liquid nitrogen.


### 2.2. Tissue Sample RNA Isolation

For RNA-seq, we used 3 lung tissue samples from control mice, 3 from asthmatic mice, and 3 from BSYQF treatment mice. Total RNA was isolated using Trizol reagent as described previously [13]. The quality of RNA was checked by 1% agarose gels electrophoresis, and quantity was determined with NanoPhotometer spectrophotometer (IMPLEN, CA, USA). Optical density values were confirmed an A260: A280 ratio above 1.9. RNA integration number (RIN) was measured using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) to confirm RIN above 7.


### 2.3. Library Preparation and Transcriptome Sequencing

We used NEBNext® UltraTM RNA Library Prep Kit for Illumina (NEB, USA) to construct sequencing libraries according to the manufacturer’s recommendations. Briefly, using poly-T oligo beads, the mRNA was purified from total RNA. In NEBNext First-Strand Synthesis Reaction Buffer, fragmentation was performed using divalent cations. M-MuLV Reverse Transcriptase (RNase H) was used to synthesize the first strand cDNA, and DNA polymerase I and RNase H were used to synthesize the second-strand cDNA. Via exonuclease/polymerase activities, and remaining overhangs were transformed into blunt ends. The NEBNext adaptor with the hairpin loop structure was ligated and prepared for hybridization after the 3 ′-end of the DNA fragment was identified. Using AMPure XP system (Beckman Coulter, Beverly, USA), we purify the library fragments and the cDNA fragment with length of 250–300 bp was selected preferentially. Finally, PCR products were purified using AMPure XP system, and on the Agilent Bioanalyzer 2100 system, the library quality was evaluated. TruSeq PE Cluster Kit V3 -cBot-HS (Illumia, San Diego, CA, USA) was preformed to cluster the index coded samples according to the manufacturer's protocol on the cBot Cluster generation system. After cluster generation, sequenced library was prepared on Illumina Novaseq6000 platform to generate 150 bp paired end reads.


### 2.4. Quality Control, Mapping, and Quantification of RNA-Seq Reads

Firstly, clean reads were obtained by removing reads containing adapter or ploy-N as well as low quality reads from raw reads. Meanwhile, the Q20, Q30, and GC contents of clean data were calculated. All downstream data were analysed on the basis of high quality cleaning readings. STAR was performed to align clean reads to reference genome. In the mapping speed, STAR outperformed other aligners by a factor of >50 and improved alignment sensitivity and precision [14]. The HTSeq V0.6.0 count was mapped to the reading for each gene. The value of the FPKM (number of segments per thousand base exon; mapping per million segments) is calculated based on the length of the gene as well as the read count mapped to the gene.


### 2.5. Differential Expression Analysis

After the significant analysis and FDR analysis, only genes under |log2FC| > 1 andp value < 0.05 criteria were included in the analysis. DESeq2 algorithm was applied to filter the differentially expressed genes.


### 2.6. Functional Analysis

The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to describe the function of each DEGs. We applied the ingenuity pathway analysis (IPA, Redwood City, CA, USA) to construct gene interaction network.


### 2.7. Real-Time PCR

The primers used in this study were purchased from BioTNT (Shanghai, China). Total RNA from treated, asthma, and control lung samples was extracted by the TRIZOL method. cDNA synthesis was performed using cDNA Reverse Transcription Kit (Takara) according to the manufacturer’s instructions. qRT-PCR was performed using TB Green Premix Ex Taq Kit (Takara). We used GAPDH gene as an internal reference to normalize the expression of target genes. Reactions were run on the 7500 Real-time PCR system (ABI, USA), and the data were analysed using its software. Ct (threshold cycle) was used to determine the target gene expression levels.


### 2.8. Statistical Analysis

Intergroup differences were analysed using one-way ANOVA. Data in the experiments are expressed as mean ± SD. All statistical analyses were accomplished by using PRISM version 7.0 (GraphPad). Ap value of <0.05 was considered statistically significant.


## 2.1. Animals

Lung tissue samples were obtained from our recently published study that tested the effects of BSYQF on airway remolding in murine chronic asthma [11]. Female BALB/c mice were sensitized and challenged with ovalbumin for 8 weeks to establish chronic asthmatic model as described previously [11]. BSYQF was prepared as described previously [12]. Briefly, Astragalus membranaceus (Fisch.) Bunge, Rehmannia glutinosa Libosch, and Epimedium brevicornu Maxium were decocted in the ratio of 3 : 2 : 1.5(w/w). From day 14, mice in BSYQF group were oral administrated 20g raw herbs/kg body weight, and mice in control and asthma groups were oral administrated with saline. At the end of experimental period, mice were euthanized and lung tissues used to extract RNA were harvested immediately and frozen in liquid nitrogen.


## 2.2. Tissue Sample RNA Isolation

For RNA-seq, we used 3 lung tissue samples from control mice, 3 from asthmatic mice, and 3 from BSYQF treatment mice. Total RNA was isolated using Trizol reagent as described previously [13]. The quality of RNA was checked by 1% agarose gels electrophoresis, and quantity was determined with NanoPhotometer spectrophotometer (IMPLEN, CA, USA). Optical density values were confirmed an A260: A280 ratio above 1.9. RNA integration number (RIN) was measured using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) to confirm RIN above 7.


## 2.3. Library Preparation and Transcriptome Sequencing

We used NEBNext® UltraTM RNA Library Prep Kit for Illumina (NEB, USA) to construct sequencing libraries according to the manufacturer’s recommendations. Briefly, using poly-T oligo beads, the mRNA was purified from total RNA. In NEBNext First-Strand Synthesis Reaction Buffer, fragmentation was performed using divalent cations. M-MuLV Reverse Transcriptase (RNase H) was used to synthesize the first strand cDNA, and DNA polymerase I and RNase H were used to synthesize the second-strand cDNA. Via exonuclease/polymerase activities, and remaining overhangs were transformed into blunt ends. The NEBNext adaptor with the hairpin loop structure was ligated and prepared for hybridization after the 3 ′-end of the DNA fragment was identified. Using AMPure XP system (Beckman Coulter, Beverly, USA), we purify the library fragments and the cDNA fragment with length of 250–300 bp was selected preferentially. Finally, PCR products were purified using AMPure XP system, and on the Agilent Bioanalyzer 2100 system, the library quality was evaluated. TruSeq PE Cluster Kit V3 -cBot-HS (Illumia, San Diego, CA, USA) was preformed to cluster the index coded samples according to the manufacturer's protocol on the cBot Cluster generation system. After cluster generation, sequenced library was prepared on Illumina Novaseq6000 platform to generate 150 bp paired end reads.


## 2.4. Quality Control, Mapping, and Quantification of RNA-Seq Reads

Firstly, clean reads were obtained by removing reads containing adapter or ploy-N as well as low quality reads from raw reads. Meanwhile, the Q20, Q30, and GC contents of clean data were calculated. All downstream data were analysed on the basis of high quality cleaning readings. STAR was performed to align clean reads to reference genome. In the mapping speed, STAR outperformed other aligners by a factor of >50 and improved alignment sensitivity and precision [14]. The HTSeq V0.6.0 count was mapped to the reading for each gene. The value of the FPKM (number of segments per thousand base exon; mapping per million segments) is calculated based on the length of the gene as well as the read count mapped to the gene.


## 2.5. Differential Expression Analysis

After the significant analysis and FDR analysis, only genes under |log2FC| > 1 andp value < 0.05 criteria were included in the analysis. DESeq2 algorithm was applied to filter the differentially expressed genes.


## 2.6. Functional Analysis

The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to describe the function of each DEGs. We applied the ingenuity pathway analysis (IPA, Redwood City, CA, USA) to construct gene interaction network.


## 2.7. Real-Time PCR

The primers used in this study were purchased from BioTNT (Shanghai, China). Total RNA from treated, asthma, and control lung samples was extracted by the TRIZOL method. cDNA synthesis was performed using cDNA Reverse Transcription Kit (Takara) according to the manufacturer’s instructions. qRT-PCR was performed using TB Green Premix Ex Taq Kit (Takara). We used GAPDH gene as an internal reference to normalize the expression of target genes. Reactions were run on the 7500 Real-time PCR system (ABI, USA), and the data were analysed using its software. Ct (threshold cycle) was used to determine the target gene expression levels.


## 2.8. Statistical Analysis

Intergroup differences were analysed using one-way ANOVA. Data in the experiments are expressed as mean ± SD. All statistical analyses were accomplished by using PRISM version 7.0 (GraphPad). Ap value of <0.05 was considered statistically significant.


## 3. Results

### 3.1. Gene Expression Analysis

DEGseq analysis results applied on RNA-seq FPKM. The RNA expression result in lung tissue showed significant differences among control (C, control), asthma (M, model), and BSYQF treatment in the mouse model of asthma (T, treatment). As shown in Figure1(a), clustering analysis was used to detect overexpressed genes and underexpressed genes.Figure 1
(a) The clustering heat maps of gene expressions. Red represents the upregulation expressions; blue represents the downregulation expressions. The color brightness is associated with differences in multiples. (b) The Venn diagram shows RNA-seq data for all genes that were regulated in asthma compared to control (M vs C) and by BSYQF treated compared to non-treated asthma controls (T vs M). FC: fold changes. Blue: M vs C, positive FC (FC > 1). Yellow: M vs C, negative FC (FC < 1). Green: M vs T, positive FC (FC > 1). Pink: M vs T, negative FC (FC < 1). (c). Volcano plot (M vs C). (d) Volcano plot (M vs T). Genes of interest are labeled as examples to show changes in expression levels. Horizontal dashed line represents expression valuep=0.05. Genes outside the midline have FC > 1.
(a)(b)(c)(d)


### 3.2. Gene Expression Regulated by Asthma and BSYQF

Gene expression differences between the asthmatic compared to control (A vs C) and asthma compared to BSYQF treatment (A vs T) groups were shown in the Venn diagram in Figure1(b) and Table 1. As the Venn diagram shown, changes in expression levels were observed in 1,899 transcripts in the lung. 271 (180 + 7+84 + 0) of these were regulated by both asthma and treatment. As shown in Table 1, 1,746 transcripts were regulated by asthma with log2 fold changes (FC) ranging from −8.6 to 11.7, p<0.05. Of these, 1,025 were upregulated and 721 were downregulated by asthma. BSYQF treatment alters the expression of 424 transcripts with FC between −5.5 and 7.8 (p<0.05). Among these, 325 were upregulated and 99 were downregulated. By BSYQF treatment, a total of 264/1746 (15.1%) transcripts including 180 upregulated and 84 downregulated transcripts, and their expression levels were reversed to normal levels.Table 1
Differently expressed genes showing the highest fold changes (in either direction) in lung samples of the mice in the control (C), asthmatic model (M), and BSYQF treatment (T) groups.
A. M vs CB. T vs MCGeneFold changeGeneFold changeGeneM vs CT vs MIncreasedIncreasedTimp16.22−3.01Chil411.74A530016L24Rik7.79Reg16.20−3.30Sprr2a311.18Tmem1797.17Dmp15.65−4.15Clca19.86Serpinb117.01Glod55.55−5.52Mcpt18.71Ebf26.85Rnf175.35−5.33Rnase2a8.43Htr3a6.59Dspp5.14−5.11Capn97.48Lep6.52Lep−8.646.52Saa37.42Glycam16.43A530016L24Rik−8.367.79Saa17.23Syt46.39Tmem179−7.577.17Saa47.18Lrtm15.95Lctl−7.365.20Retnlb7.09Nalcn5.91Ebf2−6.776.85A6.93Slc5a85.89Nalcn−6.425.91Cd209e6.93Dsc35.84Adipoq−5.975.80Mmp126.84Adipoq5.80Htr3a−5.586.59Timd26.78Lsmem15.76Myh1−5.574.94Retnla6.64Ugt2b345.57Car3−5.505.10Cbln16.47Cartpt5.54Gm10382−5.404.72Mrgprg6.32Dnmt3c5.51Lrtm1−5.395.95Prss326.27Nefm5.46Slc5a8−5.295.89Timp16.22Slc38a35.42Hist1h4k−5.275.29Reg16.20Pnoc5.39Kcnq3−5.265.22Fgl16.18Hist1h4k5.29Pnoc−5.185.39Btbd166.12Krt6a5.26Dnmt3c−5.165.51Nuggc6.11Serpina1e5.24Retn−5.094.96DecreasedAcp75.24Cfd−5.065.01Lep−8.64Krt55.23Lsmem1−5.035.76A530016L24Rik−8.36Kcnq35.22——Tmem179−7.57BC1170905.21Lctl−7.36Lctl5.20Ebf2−6.77BC0165795.19Gys2−6.69Car35.10Ebf3−6.58Gm54165.09Nalcn−6.42Dync1i15.03——Cfd5.01DecreasedGlod5−5.52Rnf17−5.33Dspp−5.11Gm49320−5.03Azgp1−5.03A, all genes with increased or decreased expression levels by FC > (±)6 (p<0.05) in the asthma group (M) relative to the control group (C); B, all genes upregulated or downregulated by the BSYQF treatment group (T) with FC > =(±)5 (p<0.05) relative to the asthma group (M). C, all genes with FC > (±)5 (p<0.05) that were reverted in their expression towards normal levels by BSYQF treatment. Genes are ranked orderly by FC (positive or negative values).


### 3.3. Volcano Results

In Figure1 C/D, the volcano plot shows transcriptome changes in asthmatic compared to control (A vs C) and asthma compared to BSYQF treatment (A vs T) groups in the lung tissues of mice. The scatter plot showed that the distribution of genes in significance (y axis, –log10 (p value)) vs fold changes (x axis, log2 (normalized fold change (FC)). As shown in Figure 1 C/D, genes are outside the midline of the absolute normalized FC > =(±)1, and –log10 (p value)> = 1.3 are colored red for overexpressed genes and blue for underexpressed genes. The underlying key genes in asthma pathogenesis and those regulated by treatment including Adipoq, HMOX1, SPP1, Cyp2e1, TNC, MB, MPO, Col9a1, Ckmt2, and Erbb4 were taken as examples to illustrate the positions and relationships with other points and midline in the figure. There were small changes in the expression levels of many genes in both directions, but only the FC > =(±)1 and p<0.05 genes that met our screening criteria were meaningful in this study.


### 3.4. BSYQF Treatment Normalized Gene Expression Altered in Asthma

From this analysis (Table1 A/B), we found a small group of genes that showed significant (above 50%) or complete reversal of changes in asthma during BSYQF therapy. These genes are listed in Table 1 C and rank-ordered by FC> =(±)5 (p<0.05). These DEGs were regulated by both asthma and BSYQF and may be closely related to the mechanism of BSYQF’s anti-asthma.


### 3.5. Key Gene Interaction Network Identification

Based on the gene expression analysis and our previous results, 21 genes were selected from DEGs between asthma group and treatment group according to their relevant functions with inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPKs signal pathways which are thought to be related to the action mechanism of BSYQF (Table2). Then, we applied IPA software to identify gene networks about 21 DEGs (Figure 2). 11 of 21 input genes could be located in this network including Adipoq, HMOX1, SPP1, Cyp2e1, TNC, MB, MPO, Col9a1, Ckmt2, Erbb4, and MAPK10. In the network, Adipoq, HMOX1, SPP1, and TNC appeared to be key components which directly or indirectly interact with multiple genes and participate in multiple biological functions.Table 2
Differently expressed gene with FC > (±) 2.0 in asthma group (M) relative to treatment group (T) related to inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPK signal pathways.
Gene nameDescriptionlog2 fold change (A vs T)AdipoqAdiponectin C1Q and collagen domain containing5.80NefmNeurofilament medium polypeptide5.46Car3Carbonic anhydrase 35.10MbMyoglobin3.94Krt15Keratin 153.99Cyp2e1Cytochrome P450 family 2 subfamily e polypeptide 13.71Gdap1Ganglioside induced differentiation associated protein 13.66Cyp4a12bCytochrome P450 family 4 subfamily a polypeptide 12B3.10Cacna2d2Calcium channel voltage dependent alpha 2 delta subunit 22.77MpoMyeloperoxidase2.77Ckmt2Creatine kinase mitochondrial 22.75Erbb4Erb-B2 receptor tyrosine kinase 42.63Mapk10Mitogen-activated protein kinase 102.47PlnPhospholamban2.42ObscnObscurin cytoskeletal calmodulin and titin interacting RhoGEF2.24Cox8bCytochrome c oxidase subunit VIIIb2.23Col9a1Collagen type IX alpha 1−4.69Mt-Nd4lMitochondrially encoded NADH dehydrogenase 4L−2.94Hmox1Heme oxygenase 1−2.48TncTenascin C−2.16Spp1Secreted phosphoprotein 1−2.02Figure 2
The gene networks by IPA. The intensity of green or red coloured shapes indicates the degree of down- or upregulation, respectively.


### 3.6. Verification of Selected DEGs

According to IPA analysis result, Adipoq, HMOX1, SPP1, and TNC were selected for further validation. We performed qRT-PCR analysis on the four genes to confirm the results of DEGs in the lung. As shown in Figure3, a significant reduction of adiponectin (Adipoq) gene expressions in asthma group compared with control group were observed. In addition, the expressions of heme oxygenase 1 (HMOX1), tenascin C (TNC), and secreted phosphoprotein 1 (SPP1) in asthma group compared with control group genes were significantly elevated, while BSYQF treatment completely or partially reversed asthma-induced expression changes of Adipoq, HMOX1, SPP1, and TNC. The results showed that the expression of the four genes from qRT-PCR was consistent with the RNA-seq analysis pattern.Figure 3
qRT-PCR validation of Adipoq, HMOX1, SPP1, and TNC results.p∗<0.05, p∗∗<0.01 vs control; p#<0.05, p##<0.01 vs asthma.


## 3.1. Gene Expression Analysis

DEGseq analysis results applied on RNA-seq FPKM. The RNA expression result in lung tissue showed significant differences among control (C, control), asthma (M, model), and BSYQF treatment in the mouse model of asthma (T, treatment). As shown in Figure1(a), clustering analysis was used to detect overexpressed genes and underexpressed genes.Figure 1
(a) The clustering heat maps of gene expressions. Red represents the upregulation expressions; blue represents the downregulation expressions. The color brightness is associated with differences in multiples. (b) The Venn diagram shows RNA-seq data for all genes that were regulated in asthma compared to control (M vs C) and by BSYQF treated compared to non-treated asthma controls (T vs M). FC: fold changes. Blue: M vs C, positive FC (FC > 1). Yellow: M vs C, negative FC (FC < 1). Green: M vs T, positive FC (FC > 1). Pink: M vs T, negative FC (FC < 1). (c). Volcano plot (M vs C). (d) Volcano plot (M vs T). Genes of interest are labeled as examples to show changes in expression levels. Horizontal dashed line represents expression valuep=0.05. Genes outside the midline have FC > 1.
(a)(b)(c)(d)


## 3.2. Gene Expression Regulated by Asthma and BSYQF

Gene expression differences between the asthmatic compared to control (A vs C) and asthma compared to BSYQF treatment (A vs T) groups were shown in the Venn diagram in Figure1(b) and Table 1. As the Venn diagram shown, changes in expression levels were observed in 1,899 transcripts in the lung. 271 (180 + 7+84 + 0) of these were regulated by both asthma and treatment. As shown in Table 1, 1,746 transcripts were regulated by asthma with log2 fold changes (FC) ranging from −8.6 to 11.7, p<0.05. Of these, 1,025 were upregulated and 721 were downregulated by asthma. BSYQF treatment alters the expression of 424 transcripts with FC between −5.5 and 7.8 (p<0.05). Among these, 325 were upregulated and 99 were downregulated. By BSYQF treatment, a total of 264/1746 (15.1%) transcripts including 180 upregulated and 84 downregulated transcripts, and their expression levels were reversed to normal levels.Table 1
Differently expressed genes showing the highest fold changes (in either direction) in lung samples of the mice in the control (C), asthmatic model (M), and BSYQF treatment (T) groups.
A. M vs CB. T vs MCGeneFold changeGeneFold changeGeneM vs CT vs MIncreasedIncreasedTimp16.22−3.01Chil411.74A530016L24Rik7.79Reg16.20−3.30Sprr2a311.18Tmem1797.17Dmp15.65−4.15Clca19.86Serpinb117.01Glod55.55−5.52Mcpt18.71Ebf26.85Rnf175.35−5.33Rnase2a8.43Htr3a6.59Dspp5.14−5.11Capn97.48Lep6.52Lep−8.646.52Saa37.42Glycam16.43A530016L24Rik−8.367.79Saa17.23Syt46.39Tmem179−7.577.17Saa47.18Lrtm15.95Lctl−7.365.20Retnlb7.09Nalcn5.91Ebf2−6.776.85A6.93Slc5a85.89Nalcn−6.425.91Cd209e6.93Dsc35.84Adipoq−5.975.80Mmp126.84Adipoq5.80Htr3a−5.586.59Timd26.78Lsmem15.76Myh1−5.574.94Retnla6.64Ugt2b345.57Car3−5.505.10Cbln16.47Cartpt5.54Gm10382−5.404.72Mrgprg6.32Dnmt3c5.51Lrtm1−5.395.95Prss326.27Nefm5.46Slc5a8−5.295.89Timp16.22Slc38a35.42Hist1h4k−5.275.29Reg16.20Pnoc5.39Kcnq3−5.265.22Fgl16.18Hist1h4k5.29Pnoc−5.185.39Btbd166.12Krt6a5.26Dnmt3c−5.165.51Nuggc6.11Serpina1e5.24Retn−5.094.96DecreasedAcp75.24Cfd−5.065.01Lep−8.64Krt55.23Lsmem1−5.035.76A530016L24Rik−8.36Kcnq35.22——Tmem179−7.57BC1170905.21Lctl−7.36Lctl5.20Ebf2−6.77BC0165795.19Gys2−6.69Car35.10Ebf3−6.58Gm54165.09Nalcn−6.42Dync1i15.03——Cfd5.01DecreasedGlod5−5.52Rnf17−5.33Dspp−5.11Gm49320−5.03Azgp1−5.03A, all genes with increased or decreased expression levels by FC > (±)6 (p<0.05) in the asthma group (M) relative to the control group (C); B, all genes upregulated or downregulated by the BSYQF treatment group (T) with FC > =(±)5 (p<0.05) relative to the asthma group (M). C, all genes with FC > (±)5 (p<0.05) that were reverted in their expression towards normal levels by BSYQF treatment. Genes are ranked orderly by FC (positive or negative values).


## 3.3. Volcano Results

In Figure1 C/D, the volcano plot shows transcriptome changes in asthmatic compared to control (A vs C) and asthma compared to BSYQF treatment (A vs T) groups in the lung tissues of mice. The scatter plot showed that the distribution of genes in significance (y axis, –log10 (p value)) vs fold changes (x axis, log2 (normalized fold change (FC)). As shown in Figure 1 C/D, genes are outside the midline of the absolute normalized FC > =(±)1, and –log10 (p value)> = 1.3 are colored red for overexpressed genes and blue for underexpressed genes. The underlying key genes in asthma pathogenesis and those regulated by treatment including Adipoq, HMOX1, SPP1, Cyp2e1, TNC, MB, MPO, Col9a1, Ckmt2, and Erbb4 were taken as examples to illustrate the positions and relationships with other points and midline in the figure. There were small changes in the expression levels of many genes in both directions, but only the FC > =(±)1 and p<0.05 genes that met our screening criteria were meaningful in this study.


## 3.4. BSYQF Treatment Normalized Gene Expression Altered in Asthma

From this analysis (Table1 A/B), we found a small group of genes that showed significant (above 50%) or complete reversal of changes in asthma during BSYQF therapy. These genes are listed in Table 1 C and rank-ordered by FC> =(±)5 (p<0.05). These DEGs were regulated by both asthma and BSYQF and may be closely related to the mechanism of BSYQF’s anti-asthma.


## 3.5. Key Gene Interaction Network Identification

Based on the gene expression analysis and our previous results, 21 genes were selected from DEGs between asthma group and treatment group according to their relevant functions with inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPKs signal pathways which are thought to be related to the action mechanism of BSYQF (Table2). Then, we applied IPA software to identify gene networks about 21 DEGs (Figure 2). 11 of 21 input genes could be located in this network including Adipoq, HMOX1, SPP1, Cyp2e1, TNC, MB, MPO, Col9a1, Ckmt2, Erbb4, and MAPK10. In the network, Adipoq, HMOX1, SPP1, and TNC appeared to be key components which directly or indirectly interact with multiple genes and participate in multiple biological functions.Table 2
Differently expressed gene with FC > (±) 2.0 in asthma group (M) relative to treatment group (T) related to inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPK signal pathways.
Gene nameDescriptionlog2 fold change (A vs T)AdipoqAdiponectin C1Q and collagen domain containing5.80NefmNeurofilament medium polypeptide5.46Car3Carbonic anhydrase 35.10MbMyoglobin3.94Krt15Keratin 153.99Cyp2e1Cytochrome P450 family 2 subfamily e polypeptide 13.71Gdap1Ganglioside induced differentiation associated protein 13.66Cyp4a12bCytochrome P450 family 4 subfamily a polypeptide 12B3.10Cacna2d2Calcium channel voltage dependent alpha 2 delta subunit 22.77MpoMyeloperoxidase2.77Ckmt2Creatine kinase mitochondrial 22.75Erbb4Erb-B2 receptor tyrosine kinase 42.63Mapk10Mitogen-activated protein kinase 102.47PlnPhospholamban2.42ObscnObscurin cytoskeletal calmodulin and titin interacting RhoGEF2.24Cox8bCytochrome c oxidase subunit VIIIb2.23Col9a1Collagen type IX alpha 1−4.69Mt-Nd4lMitochondrially encoded NADH dehydrogenase 4L−2.94Hmox1Heme oxygenase 1−2.48TncTenascin C−2.16Spp1Secreted phosphoprotein 1−2.02Figure 2
The gene networks by IPA. The intensity of green or red coloured shapes indicates the degree of down- or upregulation, respectively.


## 3.6. Verification of Selected DEGs

According to IPA analysis result, Adipoq, HMOX1, SPP1, and TNC were selected for further validation. We performed qRT-PCR analysis on the four genes to confirm the results of DEGs in the lung. As shown in Figure3, a significant reduction of adiponectin (Adipoq) gene expressions in asthma group compared with control group were observed. In addition, the expressions of heme oxygenase 1 (HMOX1), tenascin C (TNC), and secreted phosphoprotein 1 (SPP1) in asthma group compared with control group genes were significantly elevated, while BSYQF treatment completely or partially reversed asthma-induced expression changes of Adipoq, HMOX1, SPP1, and TNC. The results showed that the expression of the four genes from qRT-PCR was consistent with the RNA-seq analysis pattern.Figure 3
qRT-PCR validation of Adipoq, HMOX1, SPP1, and TNC results.p∗<0.05, p∗∗<0.01 vs control; p#<0.05, p##<0.01 vs asthma.


## 4. Discussion

In this research, we identified a set of transcriptome changes in the lung of chronic asthmatic mice with and without BSYQF treatment. Lung tissues were taken from our recently published study in which we have demonstrated that BSYQF treatment could ameliorate airway remodeling in asthmatic mice by reducing airway inflammation and oxidative stress in the lung [11]. The results from this study provided us insight into the molecular mechanism of BSYQF’s anti-airway remodeling in chronic asthmatic mice.We identified genes whose expression in asthmatic mice were impacted by BSYQF treatment. 264/1746 (15.1%) of genes showing abnormal expression in asthma group were reverted back to completely or partially to normal levels by BSYQF treatment. These included Adipoq, also known as Adiponectin, which in this study was found to be one of the most prominently downregulated genes in asthma and completely reversed back by BSYQF treatment. It has been shown that Adiponectin is an anti-inflammatory adipokine [15]. Since BSYQF has been shown to effectively regulate many pathological processes in asthma, we proposed that this core set of genes may be related to the pathology of asthma and represent key BSYQF targets.Based on previous result, we identified a set of genes related to inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPKs signal pathways which may play important roles in the mechanism of BSYQF’s anti-remodeling. We screened 21 DEGs with FC > (±) 2.0 that were differentially expressed between asthma and treatment. Additionally, 21 DEGs were input into IPA database to obtain gene networks related to regulatory functions. Gene network analyses by IPA are based on current known pathways. We used IPA because we aimed to elucidate the interaction of the gene expression found in specific functions and pathways which are thought to be related to BSYQF's anti-remodeling. Furthermore, we identified the possible key targets of BSYQF through IPA. In the gene network analyses by IPA, we found that Adipoq, HMOX1, SPP1, and TNC are genes at critical nodes. We also verified the four DEGs selected through qRT-PCR and confirmed that their variation trend were consistent with the RNA-seq analysis.Adipoq, also known as diponectin, is an important adipokine with anti-inflammatory and antioxidative effects [16]. Adiponectin plays a physiological role by activating receptors such as AdipoR1, AdipoR2, and T-cadherin which are expressed on airway epithelial and endothelial pulmonary cells playing functional roles on lung physiology [17, 18]. AdipoRs mediate pleiotropic adiponectin actions through signaling mechanisms including AMPK, AKT, ERK1/2, and P38 [19]. Previous study reported that adiponectin was reduced in the asthmatic murine model [20]. Additionally, adiponectin attenuated the airway inflammation and AHR in asthmatic mice [20]. Another study has also shown that adiponectin can reduce inflammation and suppress oxidative stress in vivo [21]. Adiponectin not only downregulates proinflammatory cytokines including TNF-α, IL-6, and NF-κB but also upregulates anti-inflammatory cytokines such as IL-10 [22]. Benjamin D et al. found that adiponectin deficiency increased allergic airway inflammation and pulmonary vascular remodeling [23]. In our RNA-seq result, Adipoq was significantly downregulated more than five-fold in asthma group and reversed completely after BSYQF treatment which drew our attention. The qRT-PCR results further confirmed the result of RNA-seq. Furthermore, in the gene network analysed by IPA, Adipoq interacted with multiple genes and participated in multiple cellular functions. Hence, the mechanism related to Adipoq expression might provide a new understanding for illustrating the therapeutic target of BSYQF’s anti-remodeling.HOMX1, known as HO-1, is a rate-limiting enzyme that catalyzes the degradation of heme into biliverdin, ferrous ion, and carbon monoxide [24]. Studies have shown that HO-1 has anti-inflammatory, excessive cell proliferation and interstitial fibrosis effects [25]. Previous study has shown that HO-1 could ameliorate airway inflammation by downregulating the tumor necrosis factor receptor (TNFR)1 dependent oxidative stress [26]. In asthma, increased HO-1 expression and activity are considered protective. The expression level of HO-1 was associated with its anti-asthma effect. Our results turned out that HO-1 gene expression was increased in asthmatic mice which was consistent with that in the previous study [27]. However, there are conflicting results about the role of BSYQF in regulating the HO-1 expression. We observed BSYQF treatment reduced HO-1 gene expression by using RNA-seq analysis and qRT-PCR. This implies that it is possible that the effect of BSYQF on HO-1 expression differs between gene and protein levels and further studies are needed.SPP1 (secreted phosphoprotein 1, also known as osteopontin) is a phosphorylated glycoprotein secreted by activated macrophages, leueocytes, and activated T lymphocytes [28]. It is recognized as a key cytokine involved in Th1 cytokine expression and immune cell recruitment at sites of inflammation [28, 29]. It is also a mediator of tissue repair and remodeling [30, 31]. According to RNA-seq and qRT-PCR results, we found the SPP1 expression increased in asthma and decreased after BSYQF treatment which suggested that SPP1 might play important effect on the mechanism of BSYQF’s anti-remodeling.TNC is considered an important gene in asthma pathogenesis [32] and known as an extracellular matrix (ECM) protein whose expression is increased in asthma [32]. TNC plays important roles in cell communication and signal transduction and regulates cell adhesion, migration, proliferation, and differentiation [32, 33]. Study in TNC knockout (KO) mice showed that lung disease phenotype was largely attenuated in the absence of TNC in OVA-induced asthma [34]. Our RNA-seq and qRT-PCR results confirmed that TNC up-regulation in asthma. Additionally, we found that TNC mRNA expression decreased after BSYQF treatment. These results provided that BSYQF may have a beneficial effect on asthma by reducing the TNC expression.


## 5. Conclusion

In this study, we used RNA-seq analysis to obtain the differentially expressed genes. We found expressions of a group of coding genes were altered in asthmatic mice. Treatment with BSYQF reversed the expressions of a specific subset of these genes in asthmatic mice. As BSYQF was previously shown to ameliorate airway remodeling by reducing airway inflammation and oxidative stress in asthmatic mice, we enriched a group of genes related to inflammatory, oxidative stress, mitochondria, PI3K/AKT, and MAPKs signal pathways and used IPA to construct gene network. Based on gene network, we inferred that Adipoq, SPP1, and TNC may be key regulators of BSYQF's anti-remodeling effect. While RNA-seq analysis could not show protein levels and the precise role of these genes needs further verification, they provided initial clues to explore the mechanisms of BSYQF's anti-remodeling in asthma.



---

*Source: 3524571-2021-01-18.xml*