Genomewide analysis of rat periaqueductal gray-dorsal horn reveals time-, region- and frequency-specific mRNA expression changes in response to electroacupuncture stimulation.

Electroacupuncture (EA) has been widely applied for illness prevention, treatment or rehabilitation in the clinic, especially for pain management. However, the molecular events that induce these changes remain largely uncharacterized. The periaqueductal gray (PAG) and the spinal dorsal horn (DH) have been verified as two critical regions in the response to EA stimulation in EA analgesia. In this study, a genetic screen was conducted to delineate the gene expression profile in the PAG-DH regions of rats to explore the molecular events of the analgesic effect induced by low-frequency (2-Hz) and high-frequency (100-Hz) EAs. Microarray analysis at two different time points after EA stimulation revealed time-, region- and frequency-specific gene expression changes. These expression differences suggested that modulation of neural-immune interaction in the central nervous system played an important role during EA analgesia. Furthermore, low-frequency EA could regulate gene expression to a greater degree than high-frequency EA. Altogether, the present study offers, for the first time, a characterized transcriptional response pattern in the PAG-DH regions followed by EA stimulation and, thus, provides a solid experimental framework for future in-depth analysis of the mechanisms underlying EA-induced effects.

Electroacupuncture (EA), evolving from traditional manual acupuncture, has often been used to enhance mechanical stimulation using modern electronics to achieve good effects on pain, relaxation, circulation and muscle1. EA stimulation causes changes in the properties of multiple targeted organs or tissues at the physiological, morphological, and functional levels. In EA, needles are inserted into acupoints to provide precisely pulsed electro-stimulation. The frequency of electro-stimulation is an important parameter known to influence the effects of EA. Many clinical trials and laboratory studies have demonstrated that low- and high-frequency EAs have different therapeutic effects23, suggesting that the frequency-related effects may work through different mechanisms.

To obtain the optimal curative effect in the clinic, the neural mechanisms involved in EA with different frequencies have attracted increasing attention from researchers. First, the central pathways mediating low- and high- frequency EA have been proven to differ. Evidence from neuropharmacological studies has revealed that low-frequency EA activates the arcuate nucleus of the hypothalamus, and high-frequency EA activates the parabrachial nucleus to suppress nociceptive transmission23. However, the periaqueductal gray (PAG) and the spinal dorsal horn (DH) are co-activated by low- and high-frequency EA. Second, different frequency EAs involve different neurochemical mechanisms, even in the same central regions. Low-frequency EA of 2 Hz accelerates the release of enkephalin, β-endorphin and endomorphin, while high-frequency EA of 100 Hz selectively increases dynorphin release in the central nervous system (CNS)3. Meanwhile, Silva et al.4 reported that the analgesia intensity of 2 Hz EA depends on noradrenergic and muscarinic mechanisms, and the effective duration depends on both noradrenergic and serotonergic descending mechanisms, as well as spinal GABAergic modulation. In contrast, the intensity of 100 Hz EA involves spinal muscarinic and GABAB mechanisms, while the effective duration depends on spinal serotonergic, muscarinic, and GABAA mechanisms. These results indicate that different frequency EAs have different effects by acting on multiple targets. However, the molecular mechanisms underlying the EA effects are largely unknown.

In traditional Chinese medicine (TCM), Zusanli (ST36) and Sanyinjiao (SP6) are two commonly acupoints used for the treatment of a wide range of health conditions. It has been shown that EA stimulation of these two acupoints has effects on pain, neurological disorders, inflammation and gastrointestinal function567. In addition, in rats, these acupoints are anatomically similar to those in humans8. Therefore, we stimulated the ST36 and SP6 points of rats with low- and high-frequency EA and examined the gene expression in the PAG-DH region with a cDNA array to elucidate the molecular mechanisms underlying EA effects.

The regulation of gene expression in the CNS detected by microarray offers a powerful tool for exploring functional changes by measuring expression levels for thousands of genes, and the results provide insight into the mechanisms and features of external stimulation9. Gene expression analysis of healthy neural tissues induced by EA stimulation provides a framework to explore the general and specific physiological characteristics of EA with different frequencies on different tissues1011. The aim of this study is to identify the gene expression profiles in the PAG and DH in naïve rats after the end of low- and high-frequency EA stimulations at the 1-hr and 24-hr time points to elucidate the molecular regulatory mechanisms underlying EA effects.

Results

Analgesic effects induced by EA

As shown in Fig. 1, two-way repeated measures ANOVA revealed significant main effects for EA stimulation (F(2, 174) = 3.74, p < 0.05), time (F(2, 174) = 11.31, p < 0.0001), and an EA stimulations × time interaction (F(4, 174) = 3.819, p < 0.001). No significant differences were found in the baseline of the tail flick latency (TFL) among the three groups. EAs at 2 Hz and 100 Hz significantly increased the TFL at ten minutes after the end of EA stimulation (Fig. 1). However, there was no significant difference at the 24-hr time point after the end of EA administration among groups with and without EA (Fig. 1).

Figure 1

Analgesic effect induced by 2 Hz and 100 Hz EA in rats.

The analgesic effects of EA on acute thermal pain were quantified using the tail flick latency (TFL) test. Both 2 Hz and 100 Hz EA significantly increased the TFL at ten minutes after the end of EA stimulation. Data are represented as the mean ± SEM. **p < 0.01, ***p < 0.001 in comparison with restraint group. (two-way repeated measures ANOVA followed by Bonferroni's Test).

Transcriptome modulation by EA

Gene expression profiles were investigated across two CNS regions (PAG and DH) at two time points after the end of EA stimulation. The global transcriptomic profiling with all the detected genes stimulated by either 2 Hz or 100 Hz EA at the 1-hr and 24-hr time points after the end of EA stimulation was clustered using a hierarchical clustering strategy. As shown in Fig. 2a, all samples with the largest gene variances in the one-to-one gene set were clustered perfectly, first by time factor, then by tissue region, and finally by frequency of EA. This expression pattern suggested that changes by EA stimulation were most strongly time-dependent, then region-specific, and finally frequency-dependent. The results obtained by hierarchical clustering were also obtained independently with principal component analysis (PCA) (Fig. 2b). These findings indicated that multiple factors underlie EA effects. Therefore, we explored the time, region, and frequency factors of EA effects in our subsequent analysis.

Figure 2

Global gene expression in the PAG-DH regions with either 2 Hz or 100 Hz EA at the 1-hr and 24-hr time points was shown by (a) hierarchical clustering and (b) principal component analysis (PCA).

Time-dependent regulation of gene expression by EA

The first part of the analysis focused on time-dependent regulation of gene expression in both PAG-DH regions by EA stimulation. The gene expression profiles in the PAG-DH of rats at the 1-hr and 24-hr time points after the end of EA stimulation were compared with the corresponding data of restrained controls, and differentially expressed genes (DEGs) with p value < 0.01 and false discovery rate ≤ 0.01 were further analyzed. There were 2756 regulated genes at the 1-hr time point and 2828 regulated genes at the 24-hr time point. The changes in expression level for most DEGs at both time points were at log ratios of −0.6 to 0.6, suggesting that most transcriptional changes induced by EAs were subtle (fold change ≤ 1.5) (Fig. 3). In these DEGs, only 1191 were co-regulated at both time points, and 288 genes had oppositely regulated directions.

Figure 3

Frequency distribution of expression log ratios after EA stimulation relative to control.

Histogram of differentially expressed genes (RankProd analysis, p value < 0.01 and false discovery rate ≤ 0.01) induced by EA stimulation at (a) 1-hr time point (n = 2756 genes) and (b) 24-hr time point (n = 2828 genes). The change in expression level for the most differentially expressed genes at both time points were at log ratios of −0.6 to 0.6 (fold change ≤ 1.5).

As for the function of the differentially expressed genes, DEGs (both at t = 1 hr and t = 24 hr) were subjected to the disease/phenotype web-PAGE to identify enrichment of potential biological processes. Some general cellular metabolic and development processes were co-enriched at both time points. However, DEGs were comprehensively enriched in immunity/stress response-related gene sets (two gene sets Z-score < 0; nine gene sets Z-score > 0) and neural function-related gene sets (six gene sets Z-score < 0; three gene sets Z-score > 0) at the 1-hr time point (Fig. 4a). Only three and five altered gene sets were also enriched at the 24-hr time point (Fig. 4b).

Figure 4

Enrichment of gene sets altered by EA.

Gene sets were classified into three functional categories based on the gene set name derived from the MGI database. The altered gene sets were enriched with PAGE analysis as described in the Methods. (a and b) Time-dependent altered gene sets. (c and d) Region-specific altered gene sets. (e and f) Different-frequency EAs altered gene sets (the differentially expressed genes, Group A, was significantly regulated not only in the 2 Hz vs. 100 Hz comparison but also in the 2 Hz and 100 Hz vs. restraint group comparison; Group B showed significantly different gene expression only when 2 Hz was compared with 100 Hz).

Region-specific regulated genes by EA

Because the main effects induced by EA stimulation persisted only several hours after withdrawal of the needle in the physiological state, gene expressions at the 1-hr time point after the end of EA administration were further used to explore the region- and frequency-EA effects in the subsequent analysis3. At the 1-hr time point, 2088 and 2291 genes in DH and PAG, respectively, were regulated by EA stimulation, while 1044 of them overlapped in both regions.

The web-PAGE analysis showed that more gene sets were enriched in the PAG (49 enriched gene sets) than in the DH (31 enriched gene sets) (Fig. 4c, d). Consistent with time effects by EA stimulation, gene sets enriched in the DH and PAG were mainly related to cellular metabolic and development processes, neural function, and immunity/stress response.

In neural function-related gene sets, 9 enriched gene sets in the DH region could be further divided into three sub-categories, including synaptic transmission (Reduced Long Term Depression, Enhanced Long Term Potentiation, Abnormal Synaptic Vesicle Number, Abnormal Inhibitory Postsynaptic Currents, and Abnormal Excitatory Postsynaptic Currents), behavior (Weakness, Impaired Passive Avoidance Behavior, and Abnormal Grip Strength), and morphology (Decreased Brain Weight) (Fig. 4c). However, in the PAG region, 23 enriched gene sets could be divided into four sub-categories, including synaptic transmission (Decreased Paired Pulse Facilitation, Abnormal Prepulse Inhibition, Enhanced Long Term Potentiation, Abnormal Excitatory Postsynaptic Currents, Abnormal Inhibitory Postsynaptic Currents), behavior (Decreased Eating Behavior, Polydactyly, Abnormal Conditioned Taste Aversion Behavior, Weakness, Abnormal Cued Conditioning Behavior, Decreased Anxiety Related Response, Hyperactivity, and Impaired Passive Avoidance Behavior), morphology (Purkinje Cell Degeneration, Abnormal Sensory Neuron Projections, Abnormal Purkinje Cell Dendrite Morphology, Ectopic Purkinje Cell, Increased Susceptibility to Neuronal Excitotoxicity, Decreased Susceptibility to Neuronal Excitotoxicity, Hippocampal Neuron Degeneration, and Astrocytosis), and nervous system phenotype (Abnormal Hippocampus Function, Increased Thermal Nociceptive Threshold) (Fig. 4d).

Ten gene sets related to the immunity/stress response were enriched in the DH and PAG regions. Although most DEGs in this category differed between the DH and PAG regions, the expression levels of two-thirds of DEGs were up-regulated (Fig. 4c, d).

Different-frequency EAs with different gene expression

Compared to time and region factors, the effect of EA frequency was relatively subtle. To comprehensively elucidate the different effects between 2 Hz EA and 100 Hz EA, we directly compared the microarray data from the 2Hz-1 group (rats that received 2 Hz EA and were killed at 1-hr after the end of EA administration) with data from the 100Hz-1 group (rats that received 100 Hz EA and were killed at 1-hr after the end of EA administration) to find the DEGs between 2 Hz and 100 Hz EA in the DH and PAG regions. There were 1384 DEGs in the DH region for the 2Hz-1 group vs. 100Hz-1 group comparison, while there were 1458 DEGs in the PAG region. These genes could be divided into two classes with different gene expression characters. The first, Group A, was significantly regulated not only in the 2 Hz vs. 100 Hz comparison but also in the 2 Hz and 100 Hz vs. restraint group 1 (R1, which did not receive electrical current stimulation and was killed at the 1-hr time point after the end of restraint) comparison (759 genes in the DH and 804 genes in the PAG), indicating that these genes differed in the degree of regulation after EA stimulations. The other, Group B, showed significantly different gene expression only when 2 Hz was compared with 100 Hz (625 genes in the DH and 654 genes in the PAG).

Of the ten enriched neural function-related gene sets, all except for one gene set (Kinked Neural Tube), presented a down-regulated pattern when 2 Hz EA was compared with 100 Hz EA in the DH region (Z-scores < 0) (Fig. 4e). These neural function-related gene sets in the DH region included behavior categories (Abnormal Spatial Learning and Impaired Passive Avoidance Behavior in Group A; Tremors, Impaired Righting Response, Hypoactivity, and Impaired Coordination in Group B), and neurogenesis (Kinked Neural Tube, Decreased Brain Weight, and Abnormal Microglial Cell Morphology in Group A; Decreased Purkinje Cell Number in Group B). Similarly, in the PAG region, the gene expression pattern of the behavior-related categories, including Decreased Exploration in New Environment, Abnormal Grip Strength, Abnormal Grip Strength, Myoclonus, Tonic Clonic Seizures, and Impaired Passive Avoidance Behavior in Group A and Abnormal Contextual Conditioning in Group B, was down-regulated after when 2 Hz EA compared with 100 Hz EA (Z-scores < 0), except the gene sets of Impaired Balance (Z-scores > 0) (Fig. 4f). In the PAG region, the neurogenesis-related categories were only enriched in Group A (Abnormal Myelin Sheath Morphology (Z-scores > 0) and Absent Corpus Callosum (Z-scores < 0)), and the synaptic transmission-related categories were only enriched in Group B (Abnormal Nerve Conduction (Z-scores > 0) and Abnormal Inhibitory Postsynaptic Currents (Z-scores < 0)) (Fig. 4f).

Of the immunity/stress response–related gene sets, all except for one gene set (Abnormal Tumor Necrosis Factor Physiology) presented a down-regulated pattern after 2 Hz EA compared to 100 Hz EA in the DH region (Z-scores < 0) (Fig. 4e). However, in the PAG region, seven genes sets were up-regulated (Liver Inflammation, Increased Circulating Corticosterone Level, and Abnormal Tumor Necrosis Factor Physiology in Group A; Thymus Atrophy, Abnormal Interleukin Physiology, Increased Eosinophil Cell Number, and Abnormal Chemokine Physiology in Group B) (Z-scores > 0), and four gene sets were down-regulated (Spleen Hyperplasia in Group A, and Increased Susceptibility to Viral Infection, Decreased Susceptibility to Autoimmune Diabetes, and Vasculitis in Group B) (Z-scores < 0) (Fig. 4f).

Based on the enriched gene sets for the neural function and immunity/stress response, we identified the concordance of regulated directions in DEGs of Group A for the 2 Hz vs. 100 Hz comparison and region-regulated DEGs in the above analyzed. According to the index of consistency, we found that 2 Hz EA strengthened the regulation of neural function-regulated genes in the DH region (10/11) and the PAG region (12/18) (Table S1). For the immunity/stress response, different frequency EAs had different regulations in different regions. In the DH, 2- and 100-Hz EA had different up-regulated DEGs (Table S1). However, 2 Hz EA enhanced the regulation on immunity/stress response-related genes in the PAG (Table S1).

Validation of microarray data with qRT-PCR

Because the EA stimulations were effective in regulating the neural and immune function, we analyzed 23 interesting genes implicated in neural function and immunity/stress response with qRT-PCR to validate the data obtained by microarray analysis. In general, the expression of most genes analyzed by qRT-PCR was highly consistent with the results of microarray analysis (Table S2–S4).

Discussion

Previous studies have shown that, in healthy volunteers and normal animals, EA gives rise to biological responses that regulate the physiological functions of the body111213. Furthermore, the performance of surgical operations under acupuncture anesthesia was reported several decades ago1415. Before surgical incision, the use of EA was able to reduce intra-operative narcotic drugs consumption and alleviate postoperative side-effects in patients. However, little attention has been paid to its physiological influences, especially on the changes in mRNA levels following EA treatment under normal conditions. Although a large number of studies have been carried out to examine the effects of EA under different morbid conditions, the regulatory mechanism of EA remains unclear. We hoped to accurately analyze the general effects of EA under normal conditions and clarify the regulatory mechanisms involved to distinguish the effects of EA on specific active substances under morbid conditions. Therefore, naïve rats were used in this study.

Although EA has been extensively investigated for its preventive and therapeutic effects on various disorders, the complexity of the response to EA stimulation in mammals hinders the understanding of the mechanisms of EA's effectiveness using conventional methods. In the current study, microarray mRNA expression profiling was adopted to provide a comprehensive view of the genes induced by 2 Hz and 100 Hz EA in the adult rat PAG-DH regions at different time points. Previous studies revealed that EA was able to regulate gene mRNA expression at 0–3 h after 30 min EA stimulation, and the mRNA expression of some of these genes could reach a peak level at 24 h1617. Based on these findings, we chose 1 hr and 24 hr as the time points for detecting the effects of EA. Both PCA and unsupervised hierarchical clustering analysis showed a strong regular change pattern in gene expression after EA stimulations with a time-, region-, and frequency-dependent effect (Fig. 2). These results suggested that the EA stimulation triggered highly orchestrated biological processes in the body that involved multiple roles in different physiopathological processes.

First, time is a dominant factor in EA stimulation. In gene expression analysis, the profiles revealed that gene expression changes had a long-term effect. At the 1- and 24-hr time points, 2756 and 2828 genes were significantly regulated after the end of EA, respectively (Fig. 3). Furthermore, most transcriptional changes induced by EA were subtle (fold change < 1.5) (Fig. 3), which indicated that EA displayed physiological adjustment on healthy subjects. Compared with the transcriptomic regulation, acupuncture analgesia was a relatively short-lasting acute effect on healthy subjects. In the present study, both 2 Hz and 100 Hz EA produced good analgesic efficacy at 40 min after needle insertion, and the duration of electro-stimulation had no effect at the 24-hr time point after EA stimulation (Fig. 1). This result was consistent with previous studies reporting that acupuncture stimulations could elevate the pain threshold after needle insertion on healthy volunteers or naïve animals, which could persist over 30 min after withdrawal of the needle181920. Notably, in addition to gene expression regulation, we should also pay attention to the change in protein levels by EA stimulation, especially to the acute effects of EA stimulation. For example, the levels of endomorphin-2 and dynorphin in spinal perfusate were significantly increased at the end of EA stimulations and evoked a significant increase in mu-opioid receptor (MOR) binding potential at the 0.5 hr time point at the end of EA stimulations1321.

Using the disease/phenotype web-PAGE analysis, we found changes in the transcriptome. Some of these changes were specifically related to neural functions and immunity/stress response, and others were related to general metabolic processes (Fig. 4a, b). Consistent with the acupuncture analgesia effect, the neural functions and immunity/stress response-related gene sets were more enriched at the 1-hr time point than that of the 24-hr time point. In fact, many genes detected in these gene sets at the 1-hr time point have been linked to the modulation of nociceptive transmission, such as Bbs4, Camk2a, Gabbr1, Gabrd, Gfap, Hrh1, Kcna1, Kif5a, Oprl1, and Slc6a1222324252627282930. The four genes Gabbr1, Gabrd, Kif5a, and Slc6a1 belong to the neurotransmitter gamma-aminobutyric acid (GABA) system, suggesting that the GABA system might play an important role in EA analgesia. Additionally, previous studies have shown that Gfap and Oprl1 are involved in the analgesic effect of EA. For example, EA treatment suppressed expression of spinal protein levels of Gfap and proinflammatory cytokine production to inhibit hyperalgesia in an inflammatory pain rat model31. EA analgesia was enhanced in transgenic Oprl1 knock-out mice32. Furthermore, Gfap and Oprl1, as well as Bbs4, Kif5a, Slc6a1, were still significantly regulated at the 24-hr time point. In addition, acupuncture has been proposed as an important alternative therapy for Alzheimer's disease (AD)33. Strikingly, a group of genes (Apbb1, Abca2, Aplp2, and App) was significantly down-regulated at the 1-hr time point. These genes were closely related to β-Amyloid peptide (Aβ) generation, which is one of the major components of senile plagues in AD343536. Moreover, the gene expressions of Apbb1 and App were still repressed at the 24-hr time point. Apbb1 interacted with App to regulate Aβ generation from App36. Although numerous acupoints are used to treat Alzheimer's disease, ST36-SP6 has long been used in TCM as main acupoints to treat AD in clinical and animal studies3738. This result suggests that EA stimulation results in comprehensive decreases in the expression of genes related to Aβ generation, especially Apbb1/App signaling, to induce significant neuroprotective effects in AD.

Second, different regions exhibit great variation in gene regulation during the EA response. Although the regulation of 1044 genes overlapped in the PAG-DH regions at the 1-hr time point after the end of EA administration, there were still 1044 and 1247 genes regulated in the DH and PAG regions, respectively. The altered expression patterns of these genes in the DH or PAG represented a generalizable molecular response to EA, with many of these genes showing matching directional changes.

The spinal DH and PAG have different physiological functions. The DH participates in somatosensory and EA stimulation information reception and regulation. Therefore, 55.6% of neural function-related enriched gene sets in the DH region were related to synaptic transmission, only 33.3% related to behavior and 11.1% related to morphology (Fig. 4c). The PAG is an anatomic and functional interface between the forebrain and the lower brainstem that has important functions, such as pain and analgesia, fear and anxiety, reproductive behavior, and cardiovascular activity. Corresponding to PAG functions, neural function-related enriched gene sets were enriched in synaptic transmission (21.7%), behavior (34.8%), morphology (34.8%) and nervous system phenotype (8.7%) (Fig. 4d). These results suggest that EA stimulation could regulate the genes of specific gene sets in different CNS regions related to specific neural functions.

Based on the literature, we found that the nociceptive transmission and modulation-related genes regulated by EA stimulation had different patterns in the PAG and DH regions. In the PAG region, the Increased Thermal Nociceptive Threshold gene set (Htr7, Ntsr2, Ppm1f, and Prkar1b) was especially enriched. In addition, within neural function-related enriched gene sets, the Cdk5, Grm4, Hrh1, and Th genes were significantly down-regulated in the PAG region and have been shown to be involved in nociceptive transmission and modulation29394041. However, in the DH region, except for the Mtap6 gene in neural function-related enriched gene sets, other pain modulation-related DEGs were mainly enriched in immunity/stress response-related gene sets, including Ccr5, Cd74, Cxcl10, and Dusp142434445. These results reinforce the idea that orchestrated actions of neural-immune system interactions by EA stimulations play an important role in EA analgesic effects. Furthermore, different CNS regions may differently regulate and coordinate neural-immune system activity by EA stimulations to produce a highly orchestrated biological network and affect the outcomes of pain perception.

Third, despite the similar analgesic effects between 2 Hz and 100 Hz EA, their gene expression profiles remained relative distinct, demonstrating that low- and high-frequency EA stimulation can produce similar effects through distinct mechanisms. The 2 Hz EA and 100 Hz EA showed different effects on neural-immunity-related gene expression. On one hand, many neural function-related genes were down-regulated by 2 Hz EA in both DH and PAG regions compared with 100 Hz EA. Among these neural function-related genes, some of them, being co-regulated by 2- and 100-Hz EA compared to the control group, were suppressed by 2 Hz EA to a greater degree (Table S1). On the other hand, although of different immune-related genes were up- or down-regulated by 2 Hz EA and 100 Hz EA in the DH and PAG regions compared with the control group, the expression patterns of these genes in the PAG region that were co-regulated by 2- and 100-Hz EAs were enhanced by 2 Hz EA to a greater degree (Table S1). Previous studies showed that EA at 2-Hz was more effective than 100-Hz EA in nervous and immune system diseases in clinical and animal experiments. For example, 2 Hz EA stimulation generated a long-term depression in the DH, suppressing cold hypersensitivity for more than 24 h in rats with neuropathic pain, whereas 100 Hz did not have this effect46. Concerning the recovery of motor function after ischemic stroke, 2-Hz EA was much more effective than 120-Hz EA in all surveyed parameters, including the latency, central motor conduction time and amplitude of motor evoked potentials47. As for our results, EA at 2 Hz was more effective than EA at 100 Hz, which may be partly due to the fact that 2-Hz EA produced more effective transcript expression amplitude.

EA could activate different biological responses by transcriptional or post-translational modifications or by protein level mechanisms. Acupuncture effective biomolecules, including genes, proteins, polysaccharides, lipids, nucleic acids, and primary and secondary metabolites, are produced by a living organism in response to acupuncture. As mentioned above, protein has the capacity to react faster than the transcriptomics response to EA stimulations13161721. Therefore, the similar analgesic effect in the acute stage of 2-Hz and 100-Hz EA stimulations in this study may be related to non-transcriptional mechanisms.

In summary, we showed that EA stimulation activated a wide variety of genes in a time-, region- and frequency-specific manner. The present data indicate that the regulatory effect of EA could be achieved through its modulation of a neural-immune network in the CNS. Low-frequency EA was more effective than high-frequency EA, which might be partly due to the fact that low-frequency EA has more powerful gene expression regulation capability than high-frequency EA. Thus, this work establishes a baseline profile of EA stimulation on the DH-PAG regions for future genetic, electrophysiological, physiopathological and behavior studies that will give insight into the mechanisms involve in EA effectiveness.

Methods

Animals

All experiments were performed on male Sprague–Dawley rats, obtained from the Experimental Animal Center, Peking University, weighing 200–220 g at the beginning of the experiment. Animals were housed in a 12 h light/dark cycle with food and water available ad libitum. The room temperature was maintained at 22 ± 1°C and relative humidity at 45–50%. Rats were handled daily during the first three days after arrival. All experimental procedures were approved by the Animal Care and Use Committee of Peking University Health Science Center.

EA stimulation

EA stimulations were performed as described previously48. In brief, stainless steel needles, 0.3 mm in diameter and 3 mm in length, were bilaterally inserted in the hind legs, one at the acupoint ST36 and the other at the acupoint SP6. Constant-current square-wave electrical stimulation generated by a HANS LH 800 programmed pulse generator (manufactured by Astronautics and Aeronautics Aviation, Peking University) was administered via the two needles for a total of 30 min. The frequency of stimulation was set to either 2 or 100 Hz. The intensity of stimulation was increased stepwise from 0.5 to 1.0 and then 1.5 mA for each 10 min step. To exclude the effect of stress caused by animal restraint and needle insertion, the restraint group underwent the same manipulation except for the administration sham EA without electrical current stimulation.

The rats were treated and divided into six groups: 2 Hz group 1 (2Hz-1) received 2 Hz EA and was killed at the 1-hr after the end of EA administration; 2 Hz group 24 (2Hz-24) received 2 Hz EA and was killed at the 24-hr after the end of EA administration; 100 Hz group 1 (100Hz-1) received 100 Hz EA and was killed at the 1-hr after the end of EA administration; 100 Hz group 24 (100Hz-24) received 100 Hz EA and was killed at the 24-hours after the end of EA administration; restraint group 1 (R1) and restraint group 2 (R24) did not receive electrical current stimulation and were killed at the 1-hr and 24-hr time points after the end of restraint, respectively.

Nociceptive testing and statistical analysis

Nociceptive threshold was assessed by recording the TFL test49. Radiant heat from a focused light beam (3 mm diameter) produced by a 12.5 W projector bulb was applied directly to the junction between proximal the 2/3 and distal 1/3 of the tail. The projector bulb was turned off as soon as the rat flicked its tail and the latency was simultaneously recorded using a digital timer with an accuracy of 0.1 s. The voltage of the stimulation was adjusted to 12 V and room temperature was carefully monitored at 22 ± 1°C to minimize the possible influence of ambient temperature on TFLs during the test. We used a cut-off latency of 15 s to avoid possible damage to the superficial tissue of the tail. The average of three successive TFL determinations (pre-EA TFL) before EA stimulation was recorded as basal latency. The TFL ten minutes after the end of EA stimulation was also assessed as EA latency. The results were presented as the mean ± SEM and were analyzed with two-way repeated measures ANOVA followed by Bonferroni's Multiple Comparison Test.

RNA extraction and cDNA microarray hybridization

Each rat was sacrificed by decapitation and the PAG and DH tissues of the fifth and sixth lumbar (L5 and L6) spinal cord were quickly removed and stored immediately in cold RNAlater (Qiagen, Hilden, Germany) at −80°C until later experimentation.

Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and purified with RNeasy column (Qiagen, Valencia, CA, USA). RNA quality was assessed with a Lab-on-chip Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). Homemade cDNA microarray platform containing 11,444 rat genes/ESTs (Expressed Sequence Tags) was used. Microarray manufacture, experiment procedure and data extraction strategy were performed as previously described5051. Equal amounts of RNA from the same CNS regions (PAG or DH) of R1 group (n = 5) or R24 group (n = 5) were pooled and labeled with Cy3 as a reference. Each RNA sample from a different CNS region of groups 2Hz-1 (n = 9), 2Hz-24 (n = 10), 100Hz-1 (n = 10) and 100Hz-24 (n = 10) was individually labeled with Cy5. In the microarray experiment, each Cy5-cRNA for the 2Hz-1 and 100Hz-1 groups was cohybridized with Cy3-labeled cRNA pools of the R1 group, and each Cy5-cRNA for the 2Hz-24 and 100Hz-24 groups was cohybridized with Cy3-labeled cRNA pools of the R24 group.

Bioinformatic analysis

Signal intensity normalization within each array was performed by locally weighted scatter plot smoothing (LOWESS) regression to normalize the expression log-ratios for the experiments so that the log-ratios averaged to zero within each array. Further scale normalization between arrays was implemented using the median absolute deviation (MAD) approach. Clustering analysis was carried out using Cluster 3.0 software52. Pearson correlation was used to measure the distance for the data. The data had been submitted to GEO under accession GSE21758 and GSE58803. Clustering results were viewed using software Java TreeView53. PCA was also used to summarize gene expression profiles between groups using R programming language54. The bioconductor package RankProd55 was used to detect DEGs under two experimental conditions (permutation times for the RankProd analysis were set to 1000, and genes with a p value < 0.01 and false discovery rate ≤ 0.01 were defined as DEG). DEGs were applied to Parametric Analysis of Gene Set Enrichment (PAGE)56 using the disease/phenotype web-PAGE (http://dpwebpage.nia.nih.gov)57 to identify key biological functional gene sets derived from the Mouse Genome Informatics (MGI) database. A p value was calculated to test the significance of the cumulative Z-score (positive or negative) of an enriched gene set.

Real-time RT-PCR

The experiment was performed as previously described58. Aliquots of the RNA samples were used in Real-time RT-PCR (qPCR). Two micrograms of each total RNA sample were used for cDNA synthesis using PrimeScript™ RT Master Mix (TaKaRa, Dalian, China). cDNA samples were placed on ice and stored at −20°C until further use. Prior to the analysis, 20 μL of each cDNA sample was diluted with 180 μL of MilliQ water. qPCR reactions were performed with Prism 7900 Sequence Detection System (Applied Biosystems, Foster City, CA). For each reaction, 1 μL of each diluted cDNA sample was added to a mixture containing 10 μL of 2 × SYBR® Premix Ex Taq™ II (TaKaRa, Dalian, China), 1 μL of each primer (5 mM), 0.4 μL ROX Reference Dye, and 7.6 μL of MilliQ water. Cycling conditions were 30 s 95°C, followed by 40 cycles of 5 s at 95°C and 34 s at 60°C. After cycling, a melting protocol was performed with 15 s at 95°C, 1 min at 60°C, and 15 s at 95°C, to control product specificity. The fold change (FC) of target gene cDNA relative to glyceraldehyde-3-phosphate dehydrogenase (Gapdh) endogenous control was determined as follows: FC = 2−ΔΔCt, where ΔΔCt = (CtTarget − CtGapdh)test − (CtTarget − CtGapdh)control. Ct values were defined as the number of the PCR cycles at which the fluorescence signals were detected. The primer sequences are listed in Table S5. Data are presented as the mean ± SEM and analyzed with independent samples t-test or Welch' test if the variances were not equal. The results were considered significant when the two-tailed p value was < 0.05.

Author Contributions

C.C. and G.Z. conceived and designed the experiments. K.W., X.X., J.Q., L.L., R.Z., X.P. and S.X. performed the experiments. K.W., X.X., N.Q. and J.H. analyzed the data. J.H. coordinated the study. K.W., X.X., G.Z. and C.C. co-wrote the paper. All authors discussed the results and commented on the manuscript.