Network based analysis of microarray gene expression profiles in response to electroacupuncture.

Electroacupuncture (EA) has been extensively considered as a tool for treating diseases and relieving various pains. However, understanding the molecular mechanisms underlying its effect is of high importance. In this study, we performed a weighted gene co-expression network analysis (WGCNA) on data collected from a microarray experiment to investigate the relationship underlying EA within three factors, time, frequency and tissue regions (periaqueductal grey (PAG) and spinal dorsal horn (DH)) as well as the biological implication of gene expression changes. Gene expression on rats in PAG-DH regions induced by EA with 2 Hz and 100 Hz at l h and 24 h were measured using microarray technology. The WGCNA was performed to identify distinct network modules related to EA effects. To find the biological function of genes and pathways, the gene ontology (GO) Consortium was applied and the gene-gene interaction network of top genes in important modules was visualized. We identified one network module (466 genes) which is significantly associated with time, another module (402 genes) significantly related to frequency, and three modules each consisting of 1144, 402 and 3148 genes that are significantly associated with tissue regions. Furthermore, meaningful biological pathways were enriched in association with each of the experimental factors during EA stimulation. Our analysis showed the robustness of WGCNA and revealed important genes within specific modules and pathways which might be activated in response to EA analgesia. The findings may help to clarify the underlying mechanisms of EA and provide references for future verification of this study.

Introduction

Electroacupuncture (EA), an electronic developed version of manual acupuncture, has been widely applied as an alternative therapy for treating diseases and relieving various pains. The central nervous system (CNS) regions, PAG and DH are known to play a critical role in EA-induced analgesia.1, 2, 3 The EA stimulation of Chinese acupoints, known as Sanyinjiao (SP6) and Zusanli (ST36), are commonly used for managing pain, a wide range of neurological diseases, diseases of adrenocortical function, inflammation, the gastrointestinal system and the immune system. Gene expression profiling applying microarray technology could explore the unknown molecular signature of EA-induced alteration in the body. Wang and colleagues (2014) have shown the effect of EA in the modulation of a neural-immune network in the CNS by activating a wide range of genes in a specific region, time and frequency. Although their experiment was performed with multiple factors (time, frequency, region), data analysis was conducted in a univariate manner. We assume that EA stimulation should cause expression changes not only for single genes, but also for different cluster of genes. Hence, we used data from the study of Wang et al. and performed a network-based test by applying the weighted gene co-expression network analysis (WGCNA). WGCNA is a systems biology approach that clusters and examines groups of genes whose expression profiles are highly correlated defined as modules. The identified gene clusters or modules are each summarized by a representative eigengene and genes within each module are (1) tested for their enrichment of known biological pathways for functional interpretation and (2) analyzed for connectivity among the genes to define novel gene expression networks characterized by high connectivity hub genes to serve as important molecular targets. In recent literature, WGCNA is performed as a promising systems biology approach for analyzing the correlation pattern of genes in microarray experiments. It is used for identifying a few gene modules or clusters (module eigengenes) and relating them to clinical traits or external factors rather than individually relating thousands of genes to them. This is a profound contribution in alleviating the multiple testing problem in microarray data analysis as it focuses on testing the group of genes with similar expression rather than on single gene expression.

This paper reports the identified meaningful pathways and potential hub genes by WGCNA which may contribute to further elucidating the underlying mechanism of EA-induced alteration in the transcriptome.

Material and methods

Samples

The data used in this paper were collected by Wang et al., 2014, where EA needles were inserted in the hind legs of 39 male Sprague-Dawley (SD) rats, at two locations, one at SP6 and the other at ST36. The Microarray mRNA expression data, comprising 11,444 rat genes or Expressed Sequence Tags (EST) across two CNS, PAG-DH tissue regions (the L5 and L6 spinal cord) induced by EA with 2 Hz and 100 Hz at l hour and 24 h were taken from the Gene Expression Omnibus (GEO) database under accession numbers GSE58803 and GSE21758. Since multiple probes were mapped to the same gene, we summarized the expression for each gene by calculating their mean. This resulted in 10,447 genes for subsequent analysis. In the study, rats were given the frequency of either 2 Hz or 100 Hz for 30 min, nociceptive testing, returned to home cages for 1 h then were killed at 1 h and 24 h. The sample from each CNS (PAG or DH) included: 9 rats who received 2 Hz and were killed at 1 h after the end of EA, 10 rats received 2 Hz and were killed at 24 h after the end of EA, 10 rats received 100 Hz and were killed at 1 h after the end of EA, 10 rats received 100 Hz and were killed at 24 h after the end of EA.

Network construction

WGCNA R software package was applied to identify the co-expression modules. A general framework for this method is explained in Fig. S1. WGCNA calculates the Pearson correlation between every pair of genes and as microarray data can be noisy, the adjacency matrix is converted to a matrix of connection strength (similarity matrix) by raising to a power (β). Equation (1) defines the adjacency () between the expression profile of genes i and j:where x is the gene expression profile and the β ≥ 1 is a soft threshold parameter which is defined in accordance with a scale-free topology criterion and can be considered for emphasizing on strong correlation and punishing weak correlations., Studies have shown that many networks in different domains are approximately scale-free networks., We selected β = 6 as the network followed an approximate scale-free topology criterion by Zhang and Horvath. In addition, β = 6 is suggested in WGCNA as a default value which has shown to satisfy the network to be approximately scale-free.

Module detection

To detect the gene modules (clusters) with similar expression patterns, average linkage hierarchical clustering was performed, and gene modules were illustrated in distinct colors. WGCNA assigned the genes that are not found in any modules in a grey module as they are not co-expressed. Topological overlap (TO) of two genes represents their similarity based on the other genes that are connected to them (shared neighbors)., We constructed the topological overlap matrix (TOM) of dissimilarity using correlation expression values and used as input for hierarchical clustering of genes. Then gene modules were found, using a dynamic tree cutting algorithm. In WGCNA, the module eigengene (ME) is considered as a representative of the entire expression level of a given module which is estimated as the first principal component of the module.

Association of modules with traits

To identify the association of MEs to external factors, we applied the linear mixed-effect module from lme4 package in R as we need to account for the correlation in the samples. The following model was used to test the association of each ME with the experimental factors:

In (2), y is the ME, the random factor SID is simply sample ID. The model estimates marginal effects of acupuncture time, frequency, region, and their interaction effects (frequency-dependent time effect, region-dependent time effect) on ME which is a summary expression of all genes in module or cluster. Note that, (2) estimates and tests the association of each experimental factor with gene expression while adjusting for other factors including the random effect from repeated measurements on same individual (SID).

Hub genes detection and enrichment analysis

For each gene in a module we calculated module membership (MM) by correlating ME with expression profile of gene i in WGCNA is defined to estimate the importance of each gene in a module. The larger the absolute value of , the more similar the gene is to the eigengene of a module, i.e. the more important the gene i is in a module.

To quantify the significance of all the modules for each of the traits we calculated GS as the minus log of the p-value from the linear mixed-effects module for testing the association of each gene expression with the factors. This paves the way to find modules that are highly associated with the factors. In WGCNA, the mean of the gene significance of all the genes in a module is considered as the module significance (MS). The higher the mean of GS in a module, the more important the module is to the factors. Genes with the highest MM were chosen as input for enrichment analysis as they are natural candidates for further investigation., Additionally, we visualized the gene-gene interaction network of the most connected genes in the interesting modules (by computing their TOs) in visANT to discover the hub genes in interesting modules.

Furthermore, genes identified in interesting modules were submitted to GO Consortium to determine whether gene expression changes in time, frequency and tissue regions are significantly enriched with regards to known pathways. The GO Consortium computes the Fisher's exact test and then corrects the p-values by calculating the false discovery rate (FDR).

Results

First, we cleaned the data by removing genes with more than 10% missing values and imputed the rest missing values with the mean of the whole expression level for those genes. This filtering yielded 6811 genes to use for network construction. Then we performed hierarchical clustering on samples to detect outliers. As a result, one sample appeared as an outlier and was removed (Fig. S2). The dendrogram of the samples with their corresponding heatmap showing their gene expression patterns for time, frequency and tissue region is indicated in Fig. S3.

By applying the WGCNA R package, we chose the soft threshold power β = 6 as it led to the lowest power that reached the highest value of R2 = 0.9 threshold for the first time (Fig. S4). As an alternative way, we visualized the log-log plot of connectivity versus the frequency of connectivity which showed that the connectivity distribution followed a power law, thus, the scale-free network criterion was approximately fulfilled (Fig. S5). Next, through the dynamic tree cutting algorithm, 25 distinct co-expression modules containing 58 to 2442 genes were detected and 37 uncorrelated genes were discovered in grey module as they were not assigned to any co-expression modules. In addition, a hierarchical clustering of the MEs was made based on the TOM dissimilarity. The MEs clustering with the heatmap plot of MEs based on their adjacencies suggested us to merge the highly correlated modules (Figure S6. a, b). After merging the highly correlated modules, colored dark red in heatmap plot Figure S6. b, 22 MEs existed, and we made the gene dendrogram obtained from an average linkage hierarchical clustering with the list of assigned modules in different colors. Fig. S7 illustrates the results from both dynamic trees cut and merged dynamic. The list of 22 MEs with the number of genes in each are provided in Table S1 .

To identify modules that are significantly associated with time, frequency, region, time-frequency and time-region interactions, we performed linear mixed-effects module (equation (2)) where we obtained yellow, brown, black and blue as the most significant modules. Fig. S8 illustrates the coefficient and p-value for each, from the summary test statistic. In addition, we compared the MS among the modules for each of the trait factors (Fig. 1). The result showed that the yellow module consisting 466 genes had the highest relevance to time, the black module (402 genes) as the most significant module related to frequency, three modules brown (1144 genes), black (402 genes), and blue (3184 genes) had the most relevance to region, the black module as the most important module identified for interaction of time with frequency. As for the time-region interaction, three significant modules (black, blue and brow) were detected. Furthermore, scatter plot of MM versus the GS of the modules for each of the factors was made to find those with high correlation and significant p-values (Figs. S9–13). Additionally, we summarized the list of genes for each factor along with the number of genes with positive and negative coefficient from regression in Tables S2–6. The scatter plot of MM versus GS in interesting modules for each of the experimental factors are illustrated in Fig. 2. To understand the biological mechanisms associated with the genes in important modules we conducted the enrichment analysis (Table 1 and S7-10). The visualized plot of gene-gene interaction networks for yellow, black, brown and blue modules are illustrated in Fig. 3. Additionally, the top list of genes with high intramodular connectivity in each network are displayed in Tables S11–14.

Fig. 1

Box plot showing the gene significance across models for each of the traits and interactions. The x-axis shows the modules in different colors and the y-axis shows the mean of gene significance for all the genes included in a module. The higher the mean gene significance in a module, the more significantly associated the module to factors.

Fig. 2

Correlation plot based on the module membership (MM) in X-axis and gene significance (GS) in Y-axis for time, frequency, interaction of time and frequency, Region and interaction of time and region in interesting modules. The correlation value and p-value are shown on top of each plot.

Table 1

Top-ranked GO terms enriched by genes in each module and associated experimental factor.

Experimental factor (module color)	GO biological process	Fold Enriched	No. of Genes	FDR P-value
Time (yellow module)	intracellular protein transport (GO:0006886)	6.94	13	5.54e-04
	regulation of binding (GO:0051098)	7.95	8	9.36E-03
	golgi to plasma membrane transport (GO:0006893)	31.21	4	9.65E-03
	glucan metabolic process (GO:0044042)	31.86	4	1.01E-02
	cellular glucan metabolic process (GO:0006073)	31.86	4	1.07E-02
	glycogen metabolic process (GO:0005977)	31.86	4	1.15E-02
	energy reserve metabolic process (GO:0006112)	26.37	4	1.63E-02
	regulation of synaptic plasticity (GO:0048167)	10.57	6	1.773e-02
	regulation of neuron differentiation (GO:0045664)	5.08	10	1.843e-02
	regulation of protein binding (GO:0043393)	10.02	6	2.52e-02
	polysaccharide metabolic process (GO:0005976)	21.85	4	2.52E-02
	golgi vesicle transport (GO:0048193)	9.29	6	2.81E-02
	vesicle-mediated transport to the plasma membrane (GO:0098876)	19.86	4	3.22E-02
Frequency, Time-frequency interaction and Region (black module)	purine nucleoside monophosphate metabolic process (GO:0009126)	18.24	8	2.70E-05
	drug metabolic process (GO:0017144)	9.67	11	4.72E-05
	ATP metabolic process (GO:0046034)	21.54	8	6.19E-05
	Oxidative phosphorylation(GO:0006119)	42.92	6	6.55E-05
	Cellular respiration (GO:0045333)	23.04	6	3.40E-04
	respiratory electron transport chain (GO:0022904)	31.84	5	6.03E-04
	aerobic respiration (GO:0009060)	29.84	4	6.50E-03
Region and Time-region interaction (brown module)	organic substance metabolic process (GO:0071704)	1.89	48	6.31E-04
	regulation of cellular localization (GO:0060341)	4.61	13	7.67E-03
	ammonium transport(GO:0015696)	23.23	5	9.15E-03
	response to stress(GO:0006950)	2.40	25	1.68E-02
	organophosphate metabolic process (GO:0019637)	4.24	12	2.55E-02
Region and Time-region interaction (blue module)	metabolic process (GO:0008152)	1.70	180	1.48E-14
	organic substance metabolic process (GO:0071704)	1.72	172	6.24E-14
	response to antibiotic(GO:0046677)	3.98	26	3.95E-06
	regulation of oxidative stress-induced cell death (GO:1903201)	9.92	10	6.80E-05
	regulation of oxidative stress-induced intrinsic apoptotic signaling pathway(GO:1902175)	18.85	7	9.45E-05
	glutathione metabolic process(GO:0006749)	10.15	7	1.53E-03
	negative regulation of response to oxidative stress(GO:1902883)	9.60	7	2.00E-03
	response to nicotine (GO:0035094)	7.73	8	2.09E-03
	hydrogen peroxide metabolic process(GO:0042743)	11.03	6	3.56E-03
	mRNA transcription (GO:0009299)	14.50	5	4.78E-03
	vasodilation(GO:0042311)	10.77	5	1.24E-02
	response to hydroperoxide (GO:0033194)	13.71	4	2.21E-02
	hydrogen peroxide catabolic process(GO:0042744)	11.17	4	3.66E-02
	inner mitochondrial membrane organization (GO:0007007)	11.17	4	3.67E-02
	positive regulation of erythrocyte differentiation (GO:0045648)	11.17	4	3.69E-02
	antibiotic catabolic process (GO:0017001)	7.39	5	4.17E-02
	regulation of neuron apoptotic process(GO:0043523)	3.24	11	4.18E-02
	strartle response(GO:0001964)	10.40	4	4.32E-02

Fig. 3

Gene-Gene interaction network of the highly connected genes in a) yellow, b) black, c) brown and d) blue modules. The dark red nodes in the networks represent the hub genes related to EA.

Discussion

In the present study, we introduced the WGCNA method to gene expression data on an acupuncture experiment to look for clusters of genes with similar co-expression patterns. Among the 22 distinct expression modules, 466 genes clustered in the yellow module showed significant association with time during EA stimulation (Fig. 1, Fig. 2, S8). In fact, the high correlation between MM and GS in this module (cor = 0.62, p-value = 7.8e-51) indicated that genes with high MM might have the potential to be biologically significant in association with time when EA is applied (Fig. 2). The analysis of 181 genes with high MM resulted in potentially important enrichment findings such as regulation of synaptic plasticity, regulation of neuron differentiation, intracellular protein transport and some general metabolic processes (Table 1, S7). These findings are in accordance with the Wang's study. Moreover, we visualized the gene-gene interaction network in the yellow module for the top list of genes with high TO (Fig. 3a, Table S11) and found interesting mRNA hub genes (GenBank ID: NM_019288, Alias: APP), APBB1, APBB3, ARC, TNF, PRSS2, GNAO, NGFR and NAE1 related to Alzheimer's disease (AD) and serotonergic synapse which were up-regulated at two different time points. Two genes APP and APBB1 have already been reported as closely related to β-Amyloid peptide (Aβ) generation, an important component of senile plagues in AD. The hub gene SNCA is related to AD and Parkinson disease (PD) where the effect of EA on PD patients has been reported in several studies., The other hub genes (NM_080586 (GABRG1)), OPRM1 and GABARAP were up-regulated and associated with neuroactive ligand-receptor interaction, GABAergic synapse, morphine addiction and nicotine addiction, which might play a role in EA analgesia. NM_178105 (M6A) was up-regulated and involved in differentiation and migration of neuronal stem cells, positive regulation of filopodium assembly, and synapse assembly. Another up-regulated gene BG672052 (Mapb1) is related to neuron differentiation.

Two genes, BC072495 (RAB5A) and NM_21252 (OM38), with a trend of up-regulation are associated with Motor Neuron Diseases (MNDS). Although the role of EA in the analgesic effect of MNDS is unknown, one study revealed the role of EA for treating quadriplegia and showed that applying the regular body acupuncture could contribute to the movement of hands, fingers, feet, and toes of patients. Hub genes NM_022852 (IDX1) and SLC2A4 are involved in insulin secretion, type 2 diabetes mellitus (T2DM) and maturity onset diabetes of the young which are meaningful considering that some studies have shown that EA might be considered as a new insulin sensitizer, which may control the epidemics of obesity and T2DM. NM_053457 (CLDN11) displaying a trend of up-regulation is linked to tight junction pathways. It has been shown that EA treatment improves neurological function associated with regulation of tight junction proteins in rats. Furthermore, other hub genes in the network might be novel as they have not been reported elsewhere.

The black module consisting of 402 genes was the most interesting module from both association of ME with frequency and mean of GS across all the modules (Fig 1, Figs. S7 and 9). The MM and GS of genes in this module were highly correlated (cor = 0.68, p-value = 7.2e-56) (Fig. 2). Thus, we examined the list of genes with high MM (174 genes) in GO and we obtained ATP metabolic process, oxidative phosphorylation, respiratory electron transport chain, cellular respiration, aerobic respiration namely as meaningful findings that explain the biological interpretation of important genes (Table 1, S8). Moreover, biologically meaningful hub genes from the black network were explored (Fig. 3b). These up-regulated genes include (NM_001025291(MBP)), KCNMA1, CDC42, CALM2, PRKCA and MAPK1 which were enriched in MAPK signaling pathway, cAMP signaling, neurotrophin signaling, T cell receptor signaling, dopaminergic synapse, long-term depression, amphetamine addiction and phosphatidylinositol signaling. The hub gene UBC in the network belongs to the PPAR signaling pathway. One study showed that this pathway was co-regulated by EA in both hypothalamus and epididymal white adipose tissue which contributes not only to lipid metabolism but also to gluconeogenesis and thermogenesisto. Another hub gene (NM_198788 (SDHD)) was up-regulated and involved in citrate cycle, oxidative phosphorylation, non-alcoholic fatty liver diseases, Alzheimer, Parkinson and Huntington diseases. The other hub gene J01435 (Cytochrome C oxidase subunit 3) was also up-regulated and related to oxidative phosphorylation. The rest of the hub genes have not yet been reported to be regulated by acupuncture.

Furthermore, we assessed changes in gene expression on time-frequency interaction during EA. As a result, the black module was the most significant with the highest mean of gene significance (Figs. 1, 2 and S8) as well as the highest correlation between MM and GS (cor = 0.54, p-value = 8.3e-32). Additionally, a negative coefficient of time-frequency interaction on gene expression levels from the regression was identified. This means increases in time and frequency during EA stimulation from 1 h to 24 h and 2 HZ to 100 HZ result in decreases in the expression of genes associated with biological pathways in the black module. Notably, this important phenomenon was not discovered by Wang et al. in their univariate analysis.

Hub genes in the black module in time-frequency interaction showed an opposite regulation compared to the association with frequency. For instance, NM_001025291, NM_198788 and J01435 were down-regulated. This suggests that increases in time and frequency of EA stimulation result in decrease in the expression of genes related to MAPK signaling pathways, citrate cycle, oxidative phosphorylation, non-alcoholic fatty liver diseases, AD, PD and Huntington disease.

Three significant modules, brown, black and blue were identified to be related to region (Fig. 1). Among them, the brown module had the highest mean of GS and high correlation (cor = 0.78, p-value < 1e-200) of MM versus GS, followed by, black (cor = 0.71, p-value = 6.7e-63) and then blue (cor = 0.63, p-value < 1e-200) (Fig. 2). We assessed genes with high MM in both brown (465 genes) and blue modules (1124 genes) and identified interesting GO terms including those with high fold enrichment and corrected p-value < 0.05: ammonium transport, regulation of cellular localization, organophosphate metabolic process, negative regulation of oxidative stress-induced intrinsic apoptotic signaling pathway, mRNA transcription, response to hydroperoxide, vasodilation, etc (Table 1, S9 and S10). In the network of the brown module (Fig. 3c), the hub gene M21964 (NEFH) was down-regulated which is involved in amyotrophic lateral sclerosis (ALS). Multiple studies have indicated a potential functional relation between EA therapy and anti-neuroinflammation in ALS as well as the effective role of acupuncture for ALS, producing symptomatic relief, respiratory impairment and improving quality of life.20, 21, 22 Two hub genes, CLF1 and ACTB, are linked to the NEFH gene, the former is associated with Fc gamma R-mediated phagocytosis, the latter might be a housekeeping gene. The other hub genes, SLC2A4 (GLUT4), SORBS1 and CBL, are related to insulin signaling pathways with one study showing that “the effect of EA on weight loss may be related to improved insulin resistance caused by changes in the adipocyte size and by reductions in the expressions of neuroprotein Y/agouti-related protein and protein tyrosine phosphatase 1B”. The hub genes OCLN and ITGB1 are connected to the down-regulated gene NM_013081(PTK2) relating to tight junction and cytokine signaling in immune system respectively.

In the network of the blue module (Fig. 3d), we found that the hub gene AY011335 (APP) was down-regulated in association to the region. AY011335 and SNCA are related to AD and serotonergic synapse. Hub genes TGFB1, DACT2 and PRKCA are connected to the down-regulated gene NM_031578 (TESK1). The PRKCA gene is involved in MAPK signaling, GABAergic synapse and dopaminergic synapse and two other hub genes, TGFB1 and DACT2, are related to TGF-beta receptor signaling which regulates in the immune system. The hub gene (NM_017063 (IMPNB)) was up-regulated and enriched in RNA transport pathway. The other hub gene BC62081 (ALDH2) with a trend of up-regulation has links to the glycolysis pathway. Moreover, the list of identified hub genes in the black module showed up-regulation trends in association with region during EA stimulation. Thus, the DH region is characterized by increases of expression levels of those genes in contrast to the PAG region.

In addition, we analyzed the gene expression alteration on time-region interaction in EA stimulation, among all the modules, three significant modules, black (cor = 0.69, p-value = 4.1e-58), blue (cor = 0.54, p-value < 1e-200) and brown (cor = 0.37, p-value = 2e-37) showed the highest correlation between MM and GS and mean of gene significance respectively (Fig. 2, S8 and S13). From the regression of gene expression level on time and region interaction, we achieved negative coefficients in the black and blue modules and positive coefficient in the brown module. This indicates that extended EA stimulation in DH region result in decreased expression of genes associated with detected pathways in the black module. Furthermore, all lists of genes identified from the black, blue and brown networks showed an opposite regulation on interaction between time and region (opposite sign of coefficient) compared to the association with region alone.

Conclusion

In conclusion, our network-based analysis of EA induced transcriptional changes identified important genes and enriched biological pathways regulated by EA and three experimental factors (time, frequency and region) in addition to reported findings. Furthermore, our multivariable model allowed for exploring the interactive effects among the experimental factors revealing novel findings with important biological and clinical implications.

Source of support

This project was supported by the Region of Southern Denmark Research Grant for Alternative TCM acupuncture #13/26014.

Conflicts of interest

The authors declare that they have no conflict of interest.