RNA Sequencing Analysis of Gene Expression by Electroacupuncture in Guinea Pig Gallstone Models

Abstract

Background: Clinical studies have shown that electroacupuncture (EA) promotes gallbladder motility and alleviates gallstone. However, the mechanism underlying the effects of EA on gallstone is poorly understood. In this study, the mRNA transcriptome analysis was used to study the possible therapeutic targets of EA.

Methods: Hartley SPF guinea pigs were employed for the gallstone models. Illumina NovaSeq 6000 platform was used for the RNA sequencing of guinea pig gallbladders in the normal group (Normal), gallstone model group (Model), and EA-treated group (EA). Differently expressed genes (DEGs) were examined separately in Model vs. Normal and EA vs. Model. DEGs reversed by EA were selected by comparing the DEGs of Model vs. Normal and EA vs. Model. Biological functions were enriched by gene ontology (GO) analysis. The protein-protein interaction (PPI) network was analyzed.

Results: After 2 weeks of EA, 257 DEGs in Model vs. Normal and 1704 DEGs in EA vs. Model were identified. 94 DEGs reversed by EA were identified among these DEGs, including 28 reversed upregulated DEGs and 66 reversed downregulated DEGs. By PPI network analysis, 10 hub genes were found by Cytohubba plugin of Cytoscape. Quantitative real-time PCR (qRT-PCR) verified the changes.

Conclusion: We identified a few GOs and genes that might play key roles in the treatment of gallstone. This study may help understand the therapeutic mechanism of EA for gallstone.

Figure 1

Representation of the flow chart (a) and the ultrasound examinations of guinea pig gallbladders from the Normal group (b), the Model group (c), and the EA group (d) separately after 6 weeks of lithogenic diet feeding. The arrow in (c) and (d) shows gallstones in gallbladder. EA: electroacupuncture.

Figure 2

DEGs identified: volcano plots (a) and heatmap (c) of 257 DEGs in Model vs. Normal, with 176 upregulated DEGs and 81 downregulated DEGs. Volcano plots (b) and heatmap (d) of 1704 DEGs in EA vs. Model, with 270 upregulated DEGs and 1434 downregulated DEGs. A heatmap of 94 reversed DEGs (e). Venn diagram of EA reversed DEGs (f). Normal: the normal group, model: the model group, EA: the electroacupuncture group, DEGs: differently expressed genes, reversed DOWN DEGs: electroacupuncture-reversed downregulated differently expressed genes, and reversed UP DEGs: electroacupuncture-reversed upregulated differently expressed genes.

Figure 3

GO analysis of 94 EA reversed DEGs. Bar plot (a) of EA-reversed DEGs showed an enrichment score (-log(P value)) of the significant enrichment terms involving all 2 BPs, 3 CCs, and 10 MFs. Each bubble in the bubble plot (b) represents different -log(P value). GO: gene ontology; DEGs: differently expressed genes; BP: biological process; CC: cellular component; MF: molecular function.

Figure 4

PPI network of 94 EA reversed DEGs (a) with a score >0.150, and the top 10 hub genes of EA-reversed DEGs by the cytohubba plugin (b). Each node stands for a gene or a protein, and the edges represent the interactions between the nodes. PPI: protein-protein interaction; DEGs: differently expressed genes.

Figure 5

E2F1 mRNA expression using qRT-PCR. Statistical analysis was performed with one-way analysis of variance (ANOVA) and Tukey's multiple comparisons tests by GraphPad Prism 8.4.2. Compared with the Model group, ^∗∗P < 0.01.