Transcriptomic Insights into Different Stimulation Intensity of Electroacupuncture in Treating COPD in Rat Models

Abstract

Background

Electroacupuncture (EA), with varying stimulation intensities, has demonstrated therapeutic potentials in both animal and clinical studies for the treatment of chronic obstructive pulmonary disease (COPD). However, a comprehensive investigation of the intensity-related effects, particularly 1mA and 3mA of EA, and the underlying mechanisms remains lacking.

Methods

A COPD rat model was established by prolonged exposure to cigarette smoke and intermittent intratracheal instillation of lipopolysaccharide. EA treatment was administered at acupoints BL13 (Feishu) and ST36 (Zusanli), 20 minutes daily for 2 weeks, with intensities of 1mA and 3mA. EA effectiveness was evaluated by pulmonary function, histopathological change, serum level of inflammatory cytokines, and level of oxidative stress markers in serum and lung tissues. Transcriptome profiling and weighted gene co-expression network analysis (WGCNA) were performed to reveal gene expression patterns and identify hub genes. Real-time quantitative PCR (RT-qPCR) and Western blot (WB) were performed to detect the mRNA and protein expression levels, respectively.

Results

EA at both 1mA and 3mA exerted differing therapeutic effects by improving lung function and reducing inflammation and oxidative stress in COPD rats. Transcriptome analysis revealed distinct expression patterns between the two groups, functionally corresponding to shared and intensity-specific (1mA and 3mA) enriched pathways. Eight candidate genes were identified, including Aqp9, Trem1, Mrc1, and Gpnmb that were downregulated by EA and upregulated in COPD. Notably, Msr1 and Slc26a4 exclusively downregulated in EA-1mA, while Pde3a and Bmp6 upregulated solely in EA-3mA. WGCNA constructed 5 key modules and elucidated the module–trait relationship, with the aforementioned 8 genes being highlighted. Additionally, their mRNA and protein levels were validated by RT-qPCR and WB.

Conclusion

Our results demonstrated that 1mA and 3mA intensities induce distinct gene expression patterns at the transcriptional level, associated with shared and 1mA vs 3mA-specific enriched pathways. Genes Mrc1, Gpnmb, Trem1, and Aqp9 emerge as promising targets, and further studies are needed to elucidate their functional consequences in COPD.

Keywords:

• chronic obstructive pulmonary disease
• electroacupuncture
• intensity
• transcriptome profiling
• weighted gene co-expression network analysis

Introduction

Characterized by chronic inflammation, alveolar destruction (emphysema), and bronchiolar obstruction, chronic obstructive pulmonary disease (COPD) is an irreversible yet preventable disease, ranking as the third leading cause of global mortality.¹ Besides genetic factors, COPD is predominantly triggered by chronic exposure to inhaled particulate matter, mainly cigarette smoke (CS) and air pollutants. CS induces airway inflammation and oxidative stress in the lungs. Even after the cessation of CS exposure, inflammation may persist, fostering a microenvironment conducive to chronic tissue damage through elevated pro-inflammatory cytokines, such as interleukin-6 (IL-6), interleukin-1β (IL-1β), and tumor necrosis factor-α (TNF-α),² leading to the infiltration of macrophages, lymphocytes, neutrophils, and other inflammatory cells into small airways. Simultaneously, oxidative stress disrupts redox balance, impairing antioxidant enzymes like superoxide dismutase (SOD) and glutathione peroxidase (GSH-Px), thereby exacerbating the inflammatory process and perpetuating a detrimental loop in COPD pathophysiology.³ Although cytokines and oxidants/antioxidants are recognized as key players in COPD pathogenesis,⁴ a deeper understanding of the molecular mechanisms and identification of therapeutic targets are of utmost importance to address the unmet clinical needs. Recent advances in the integration of transcriptome profiling and Weighted Gene Co-expression Network Analysis (WGCNA) enable robust identification of differentially expressed genes (DEGs) and construction of module eigengenes (MEs) and networks based on the similarity profiles and functions of these DEGs.⁵ This facilitates the identification of novel hub genes and enriched pathways responsible for COPD pathogenesis in patients.^6,7

At present, bronchodilators and anti-inflammatory drugs are the mainstays of pharmacotherapy for COPD; however, long-term use may cause side effects such as tachycardia and tremors.⁸ Clinical studies have demonstrated that electroacupuncture (EA) effectively enhances respiratory function and exercise tolerance in COPD patients by suppressing inflammatory reactions, relieving small airway obstruction, and improving ventilation.^9–12 Nevertheless, the stimulation intensity for EA in the treatment of COPD patients is not standardized, with individual patient tolerance often determining intensity. In contrast, animal studies frequently utilize stimulation intensities of 1mA and 3mA to assess EA effectiveness in COPD rat models, resulting in alleviated airflow limitation and obstructive pulmonary dysfunction.^13–16 Despite the absence of standardized intensity in clinical settings for COPD patients, researchers have extensively explored the distinct regulatory effects of varying EA intensities on experimental pain¹⁷ and cerebral blood flow in the prefrontal cortex of healthy individuals,¹⁸ as well as on clinical symptoms among patients with other diseases.^19,20 In animal studies, significant efforts have been made to further clarify the intricate relationship between EA intensity and its therapeutic effects, along with exploring the underlying presumptive mechanisms.^21,22 For instance, applying low (0.5mA) and high (3mA) intensities at acupoints of ST36/ST25 in pre- and post-lipopolysaccharide (LPS) treated mice revealed opposite inflammatory responses mediated by distinct sympathetic pathways.²² These findings suggested the important role of intensity related to EA efficacy and associated mechanisms.

In this study, we aimed to conduct the first thorough comparative analysis of the therapeutic outcomes and underlying mechanisms associated with EA at intensities of 1mA and 3mA in a COPD rat model. We chose the acupoints BL13 (Feishu) and ST36 (Zusanli) following a systematic review of ancient and modern literature,^23–25 and also based on the principles and practices of traditional Chinese medicine (TCM) for COPD treatment.^26,27 We postulated that the 1mA and 3mA-dependent effects of EA in the treatment of COPD are modulated by substantial transcriptional alterations and their enriched shared and 1mA vs 3mA-specific molecular pathways, as revealed through transcriptome profiling. By constructing MEs through WGCNA algorithm and subsequently validating through existing literature, we attempted to identify potential molecular targets implicated in the EA treatment of COPD.

Materials and Methods

Main Reagents and Instruments

Primary antibodies CD3 (GB12014) and CD68 (GB11067) were purchased from Servicebio (Wuhan, China). β-Actin (ab8226), Aqp9 (ab191056), Mrc1 (ab64693), Msr1 (ab151707) were obtained from Abcam (Cambs, United Kingdom). Trem1 (A16535), Gpnmb (A14270), Pde3a (A17919), and Bmp6 (A4538) were supplied by Abclonal (Wuhan, China). Slc26a4 (DF13409) were obtained from Affinity (Changzhou, China). Elisa and Biochemical kits including IL-10 (E-EL-R0016c), IL-6 (E-EL-R0015c), IL-1β (E-EL-R0012c), TNF-α (E-EL-R2856c), GSH-Px (E-BC-K096-M), SOD (E-BC-K020-M), and malondialdehyde (MDA, E-BC-K025-M) were supplied by Elabscience (Wuhan, China). RNAprep Pure Tissue Kit (DP431), FastFire qPCR PreMix (SYBR Green, FP207), and FastKing gDNA Dispelling RT SuperMix (KR118) were purchased from Tiangen (Beijing, China). Main instruments include whole-body exposure chamber (CSM-100C, Tow-Int Tech, Shanghai, China) and AniRes2005 lung function system (V3.5, Bestlab Company, China). Other reagents and instruments in detail can be found in Supplementary Table 1.

Animal Ethics and Study

All procedures were approved by the Animal Ethics Committee of Chengdu University of TCM in accordance with the NIH Guide for the Care and Use of Laboratory Animals. Male Sprague-Dawley (SD) rats aged 12–14 weeks and weighing 200–220g, were purchased from Dashuo Co., Ltd (SCXK(chuan)2020–0030, Chengdu, China). These rats were housed in the animal facility of Chengdu University of TCM (SYXK(chuan)2020–124) supplied with food and water. The housing environment maintained a consistent temperature of 22 ± 1°C and relative humidity of 50%±5%, adhering to a standard 12-hour light/dark cycle. Twenty-four rats were randomly divided into four groups (n=6 per group): normal control group (NC), COPD group (COPD), COPD+1mA-EA group (EA-1mA), and COPD+3mA-EA group (EA-3mA).

Generation of a Rat Model of COPD

Passive CS with LPS injection is used to induce COPD-like inflammation in rats. The rats were placed in a whole-body exposure chamber (CSM-100C, Tow-Int Tech, Shanghai, China) and exposed to the smoke generated from 15 cigarettes (Jiaozi, Sichuan Tobacco Industry Co, China) containing 11 mg of tar and 1.1 mg of nicotine for 1 h, twice daily (in total 2h, 30 cigarettes per day) for 12 weeks. On days 1, 31 and 61, the rats were intratracheally injected 200 µL of LPS (1 mg/mL, L2880, Sigma, St. Louis, USA) and escaped from the CS challenge.^28–30

EA Treatment

EA treatment started on day 77 (16 days after the last LPS injection) and continued for 2 weeks, 20 mins daily. Prior to EA, rats were anesthetized with 1.5% isoflurane (R510-22-16, RWD Life Science Co., Ltd., Shenzhen, China). We chose acupoints BL13 (Feishu, 3 mm below and lateral to the third thoracic vertebrae on the back) and ST36 (Zusanli, 5 mm below and lateral to the anterior tubercle of the tibia) in treating COPD. Previous study exploring acupuncture point selection for COPD patients revealed that BL13 (Feishu) and ST36 (Zusanli) are among the top five most frequently used acupoints.²⁴ Additionally, both clinical^31,32 and animal studies^16,33,34 have confirmed their effectiveness in treating COPD. Two stainless steel needles (0.25 mm × 15 mm; Suzhou Medical Supplies Co. Ltd., China) were inserted into each acupoint. For BL13, the insertion is at a depth of 6 mm unilaterally, and for ST36, it is 7 mm (left and right sides were treated every other day). The needles were connected to an EA apparatus (HuaTuo SDZ-V type, Suzhou, China). Alternating waves of dense-sparse frequencies were applied to treat the COPD rats. The parameters were set at 4/20 Hz with intensity of 1mA or 3 mA, respectively. The anesthetic procedure was also performed on the rats in the NC and COPD groups.

Inflammatory Response in COPD Rats

Assessment of Pulmonary Function

Pulmonary function was tested using AniRes2005 lung function system (Bestlab Company, Beijing, China). Briefly, the rats were anaesthetized using 1% pentobarbital sodium (40 mg/kg, P3761, Sigma, St. Louis, USA) followed by tracheal intubation and connection to a plethysmograph. The pressure changes in the plethysmographic chamber were measured through a pressure transducer. Forced expiratory volume in 0.1 s (FEV_0.1) and in 0.3 s (FEV_0.3), and forced vital capacity (FVC) were recorded. The ratios of FEV_0.1 to FVC (FEV_0.1/FVC) and FEV_0.3 to FVC (FEV_0.3/FVC) were calculated. Each rat was measured five times and sacrificed thereafter. Blood samples and lung tissues were collected for subsequent experiments.

Histopathology Examination

Hematoxylin–eosin (HE) staining was performed to observe lung morphological changes. Tissue specimens of the right lower lobe of the lung were fixed with 10% formalin neutral buffer solution, embedded in paraffin, and cut into 3-μm sections. Whole pathological section images were scanned using a digital panoramic scanner (3DHISTECH, PANNORAMIC DESK/MIDI/250/1000 digital scanner, Hungary). Random selection of 5 non-overlapping regions of interest was performed for each section at 200× magnification (n=6 per group). Inflammatory cells, alveoli, and alveolar septa were manually delineated and counted by two observers. The ratio of the total alveolar area to the septum (A/S) was calculated using Image Pro Plus 6.0 (Media Cybernetics, Rockville, MD, USA).

Immunohistochemistry (IHC) Staining

Paraffin-embedded lung tissue sections (3-μm) were prepared for IHC to assess lung inflammation using CD68 (macrophage marker) and CD3 (T lymphocyte marker).^35,36 Briefly, the slices of lung tissues were incubated with primary antibody of CD68 at 1:200 (GB11067, Servicebio, China) and CD3 at 1:500 (GB12014, Servicebio, China) overnight at 4°C. The slices were then incubated with (horseradish peroxidase (HRP)-labeled) secondary antibodies at 1:300 (BL003A, Biosharp, China) for 50 min at room temperature. Positive cells were visualized with DAB (3, 3-diaminobenzidine) and counterstained with hematoxylin. Image acquisition was performed using a microscope (Leica, DM500, Germany). Initially, an overview of the tissue was obtained at 200× magnification. Subsequently, three random fields per section (n=3 per group) were captured at 400× magnification. Fluorescence intensity and area were measured using the Image Pro Plus 6.0 (Media Cybernetics, Rockville, MD, USA). Mean optical density (integrated optical density (IOD)/area) was calculated for each image. The average IOD/area per sample was calculated by averaging the IOD/area values from three captured images.

Level of Cytokines and Oxidative Stress

Enzyme-linked immunosorbent assay (ELISA) Kit (Elabscience, Wuhan, China) was used to detect the serum levels of inflammatory mediators, including IL-6, IL-1β, IL-10, and TNF-α as previously described.³⁷ Biochemical Kit (Elabscience, Wuhan, China) was used to determine the levels of SOD, MDA, and GSH-Px in serum and lung tissues, which are surrogate markers for oxidative stress.

RNA Sequencing (RNA-Seq)

RNA Isolation

Lung tissues (n=5 per group) were excised, and part of the left upper lobe was used for total RNA isolation using RNAprep Pure Tissue Kit (TIANGEN BIOTECH, Beijing, China) following the manufacturer’s instructions. The purity and concentration of the RNA was assessed by Nanodrop Spectrophotometer (NanoDrop Products, CA, USA), and the RNA integrity was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA).

RT-qPCR Analysis

After RNA isolation and quality control, RNA (1μg) was reverse-transcribed into cDNA (TIANGEN BIOTECH, Beijing, China). The PCR reaction mix included 10 μL of cDNA, 10 μL of SyberGreen Mix, and 0.5 μL of forward and reverse primers, for a total volume of 20 μL. Real-time PCR was performed on a BIOER Detection System (Bioer Technology Co. Ltd., Hangzhou, China). Data were calculated using the 2^−ΔΔCt method³⁸ and relative expression levels of the genes were normalized to the desired group. Each sample was run three times in triplicate. All the primers used are listed in Supplementary Table 2.

Library Construction and Sequencing

RNA-seq library was constructed using DNBSEQ DNB Rapid Make Reagent Kit (V2.0) and quantitatively controlled by Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MD, USA). Then, the library was sequenced using DNBSEQ-T7 (MGI, Shenzhen, China). For analysis, sequencing adapters and low-quality reads were first removed from raw data. Clean reads were aligned to the rat genome sequences using HISAT2 followed by transcript assembly and quantification using StringTie.³⁹ TPM (transcripts per million reads) method was used for normalization, where median TPM values (n=5 per group) of less than 1 were defined as low expression and excluded from the data. Using edgeR, significant DEGs were selected using the criteria: |log₂ FC| > 1, p value < 0.05, coefficient variation < 0.7 and false discovery rate (FDR) < 0.05. The RNA-seq data is deposited in the Bioproject database in NCBI, and the BioProject ID is PRJNA942343.

Functional Analysis of DEGs

Functional analysis of DEGs was performed through DAVID (Database for Annotation, Visualization and Integrated Discovery) online platform (https://david.ncifcrf.gov) for KEGG pathway analysis (Kyoto Encyclopedia of Genes and Genomes) and Gene Ontology (GO) analysis. The WGCNA analysis was carried out using R package, including module identification, network construction, and topological overlap matrix (TOM) computation to find clusters (modules) of highly correlated genes.⁵

Western Blot

Left upper lobe of lung tissues was homogenized in Radio-Immunoprecipitation Assay (RIPA) lysis buffer containing proteinase, phosphatase inhibitors, and phenylmethylsulfonyl fluoride (CR2202124, Servicebio, China). Protein concentrations were then determined using a BCA kit (BL521A, Servicebio, China). Equal amounts of protein (30 μg) were loaded onto 7.5–12.5% SDS-PAGE gel for electrophoresis and subsequently transferred to PVDF membranes. The membrane was blocked by non-fat milk and incubated overnight at 4°C with the following primary antibodies: β-Actin (1:10,000, Abcam, ab8226), Aqp9 (1:10,000, Abcam, ab191056), Trem1 (1:1000, Abclonal, A16535), Mrc1 (1:1000, Abcam, ab64693), Gpnmb (1:2000, Abclonal, A14270), Msr1 (1:5000, Abcam, ab151707), Slc26a4 (1:1000, affinity, DF13409), Pde3a (1:1000, Abclonal, A17919) and Bmp6 (1:1000, Abclonal, A4538). Subsequently, HRP-conjugated secondary antibodies were applied for 1 hour at room temperature. Luminescence detection was achieved by applying a freshly prepared ECL mixture. Finally, the optical density values of the target bands were precisely assessed using the ChemiScope Analysis image processing software.

Statistical Analysis

Data were presented as mean ± standard deviation (SD) and analyzed using Prism 9.0 (San Diego, California, US). If the data were normally distributed and had equal variance among groups, a one-way analysis of variance (ANOVA) was used, and least-significant difference t-test was used for a pairwise comparison; otherwise, Kruskal–Wallis test was used. p < 0.05 was considered statistically significant.

Results

Effects of EA-1mA Vs EA-3mA in Rats with COPD

Pulmonary Function

We established a COPD rat model to test the therapeutic effect of EA on pulmonary function. The schematic diagram of the experimental timeline is illustrated in Figure 1A. The pulmonary function was assessed by FEV_0.1/FVC and FEV_0.3/FVC ratios, the indexes for airway obstruction in clinics and animal models.^40,41 The ratios were significantly reduced in COPD vs NC group, indicating severe airflow limitation in the lung (Figure 1B, p<0.01). Both EA-1mA and EA-3mA intensities were able to improve lung capacity by significantly increasing the ratios (Figure 1B, p<0.01); nevertheless, no marked difference was observed between them.

Figure 1 The study timeline and effects of EA on lung function and morphology. (A) Schematic showing establishment of a COPD rat model and EA treatment schedule. (B) Effects of EA on lung function assessed by FEV_0.1/FVC and FEV_0.3/FVC ratios, A/S ratio, and quantitative analysis of inflammatory cells in NC, COPD, EA-1 and EA-3 groups (n=6 per group). (C) Representative images of HE stained lung sections from each group (n=6 per group). Scale bar: 50μm; *p<0.05, **p<0.01.

Abbreviation: ns, not significant.

HE staining showed thickening of the alveolar septum in the COPD group, along with emphysema formation and inflammatory cell infiltration. Quantitative analysis of the staining found a significantly decreased A/S ratio and increased number of inflammation cells in the COPD vs NC group (Figure 1B, p<0.01). Compared to COPD, EA-3mA intensity significantly increased the A/S ratio (p<0.05) and EA-1mA intensity did not reach the significance. Both intensities inhibited the local accumulation of inflammation cells (p<0.01). Representative photomicrographs are displayed in Figure 1C. Furthermore, IHC staining revealed remarkable infiltration of inflammation cells, such as macrophages (CD68+) and T lymphocyte (CD3+), in the lung tissue of rats in the COPD group compared to the NC group (P<0.001, Supplementary Figure 1). Quantitative analysis demonstrated that both EA-1 and EA-3 treatments significantly reduced the number of macrophages and lymphocytes, as indicated by decreased IOD/area values compared with COPD group (P<0.05).

Proinflammatory Cytokines and Oxidative Stress

As release of inflammatory cytokines and increased oxidative stress are key drivers of COPD pathophysiology, we measured their level changes upon EA treatment. The cytokine results showed that compared to COPD, EA-3mA significantly reduced IL-6, IL-1β, and TNF-α and elevated IL-10 in serum, which were similar to the baseline level of the NC group. No significant difference in the levels of cytokines, except for IL-6, was found between the EA-1mA group and the COPD group (Figure 2A).

Figure 2 Effects of EA on cytokines and oxidative stress markers. (A) Effects of EA on serum cytokines (IL-6, IL-1β, IL-10, and TNF-α) measured by ELISA assay (n=6 per group). (B) Effects of EA on oxidative stress markers (GSH-PX, SOD and MDA) examined in serum and (C) in lung tissues (n=6 per group). *p<0.05, **p<0.01.

Abbreviation: ns, not significant.

Oxidative stress is determined by antioxidant enzymes SOD and GSH-Px, and lipid peroxidation production MDA. As shown in Figure 2B and C, EA-3mA treatment significantly increased SOD and GSH-Px activities, and decreased MDA level both in serum and lung tissues compared to the COPD group (p<0.05). The antioxidant effect of EA-1mA varied, showing significance in serum but not in lung tissues.

Collectively, both EA intensities showed improvements in lung function, anti-inflammatory and antioxidant capacity in treating COPD rats.

Transcriptome Profiling Revealed Gene Expression Changes

Transcriptome Analysis Identified Significant DEGs

Next, we performed transcriptome profiling to characterize gene expression changes regulated by EA treatment. After removing sequencing adapters and low-quality reads, approximately 43.1 million clean reads were generated per sample (n=5 per group). Average read counts per transcript were TPM normalized and log2-ratio transformed. About 94.2% of them were mapped to the rat genome corresponding to 12,679 genes, of which 12,043, 11,885, 11,974, and 12,039 were distributed in the NC, COPD, EA-1mA, and EA-3mA groups, respectively (Supplementary Figure 2). To identify the DEGs, we applied a stringent cutoff criterion (TPM > 1, |log2FC| > 1, p value < 0.05, and FDR < 0.05), resulting in a total of 1101 DEGs displayed in the heatmap after hierarchical clustering (Figure 3A). As predicted, the COPD and NC groups showed different expression patterns, with the genes being discriminated into separate clusters. Strikingly, the EA-3mA group shared similar expression features with the NC group, suggesting that treatment with 3mA could modulate the dysregulated genes in COPD rats and restore their function. However, this strong concordance was not found between the EA-1mA and NC groups.

Figure 3 Results of transcriptome profiling. (A) Hierarchical cluster heatmap showing the expression patterns of differentially expressed transcripts across groups. (B) Venn diagram analysis depicting shared and unique genes among group comparison: downregulated DEGs in COPD vs NC group, upregulated DEGs in EA-1 vs COPD group, and upregulated DEGs in EA-3 vs COPD group. Likewise, shared and unique genes in (C) upregulated DEGs in COPD vs NC, downregulated DEGs in EA-1 group vs COPD group, and downregulated DEGs in EA-3 group vs COPD group. Arrows indicate selected target genes. (D) KEGG pathway enrichment analysis between COPD and NC groups, (E) between EA-1 and COPD groups, and (F) between EA-3 and COPD groups. Red color marks the shared pathways among the three comparisons. (G) The selected GO terms enriched between COPD and NC groups, (H) between EA-1 and COPD groups, and (I) between EA-3 and COPD groups; Biological process (red), cellular component (blue), and molecular function (green).

The genes that inversely changed their expressions in the COPD vs EA group were the target genes associated with EA modulatory function. Therefore, we performed a Venn diagram analysis to determine the overlapping and exclusive presence of genes by pairwise comparison (Figure 3B and C). Compared to the COPD group, EA treatment upregulated 94 genes in the EA-1mA group and 255 genes in the EA-3mA group (Figure 3B), resulting in 32 overlapping genes. Twenty-seven of the 32 genes were downregulated in the COPD vs NC group, suggesting that EA treatment may restore their expression (Supplementary Table 3). Likewise, a consensus of 63 genes (Supplementary Table 4) were identified to be upregulated in the COPD vs NC group, and their expression was inversely downregulated after EA treatment (Figure 3C). Of note, several genes were exclusively presented in each EA group, eg, downregulation of Msr1 and Slc26a4 only in EA-1mA, and upregulation of Pde3a and Bmp6 only in EA-3mA, suggesting 1mA vs 3mA-specific regulation of gene expression.

Function Annotation of DEGs

To comprehensively understand the functions and pathways of the DEGs involved in the pathogenesis of COPD and EA efficacy, GO function and KEGG pathway enrichment analyses were performed by pairwise group comparisons. Top 20 enriched pathways with the number of genes annotated in each pathway are summarized in Figure 3D–F. Selected GO functions, including three GO categories of biological process (red), cellular component (blue), and molecular function (green), are listed in Figure 3G–I and Supplementary Table 5.

As shown in Figure 3D, the DEGs in the COPD vs NC group were mainly enriched in the pathways of, eg, complement and coagulation cascades (rno04610), phagosome (rno04145), circadian rhythm (rno04710), cGMP-PKG signaling pathway (rno04022), platelet activation (rno04611), etc, which have been implicated in the pathophysiology of COPD.^42–45 These mentioned pathways were also enriched in the EA-3mA vs COPD group, indicating their strong association with COPD. Pathways such as cholinergic synapse (rno04725), glutamatergic synapse (rno04724), and axon guidance (rno04360) were specifically enriched in EA-3mA group (Figure 3F).

In contrast to the EA-3mA group, the DEGs in the EA-1mA vs COPD group were more likely in favor of anti-inflammation or inducing immune response, mostly functionally enriched in neutrophil extracellular trap formation (rno04613), phagosome (rno04145), HIF-1 signaling pathway (rno04066), arachidonic acid metabolism (rno00590), and natural killer cell mediated cytotoxicity (rno04650), etc. (Figure 3E)

Validation of COPD-Associated Genes with Published Datasets

To support our findings, we validated the identified COPD-associated genes (1101 DEGs) with published datasets, including 10 animal studies examining gene expression in rat or mouse lung tissue (Supplementary Table 6) and 9 clinical studies evaluating gene expression in lung tissue of COPD patients (Supplementary Table 7). While transcriptional data directly studying EA treatment in COPD patients or animals was unavailable, we examined the regulatory effect of EA by cross-comparing genes validated in the COPD vs NC group.

In the EA-3mA group, 13 genes (5 down-regulated and 8 up-regulated) were confirmed in animal studies, and 32 genes (15 up-regulated and 17 down-regulated) were matched to clinical data. EA-3mA treatment could normalize their expression levels. Similarly, EA-1mA treatment reversely regulated 15 genes (3 down-regulated and 12 up-regulated) that overlapped with the animal datasets, as well as 23 up-regulated and 6 down-regulated genes validated in human studies. Notably, the upregulation of Mrc1 and Gpnmb in the COPD vs NC group was validated in both animal^46,47 and human studies.^48,49 Both EA-1mA and EA-3mA were able to restore their expression levels, suggesting that Mrc1 and Gpnmb may serve as promising therapeutic targets for EA treatment.

WGCNA Analysis Underpinned COPD Phenotype-Associated Genes

We used the WGCNA method to generate modules and identify the module-correlated hub genes. Cluster dendrogram was constructed from the 1101 DEGs identified from 20 samples. To merge the modules with high similarity of feature genes, dynamic mixed shearing methods were used, and accordingly obtained 5 colored MEs, including MEturquoise, MEbrown, MEblue, MEyellow, and MEgrey (Figure 4A). Gene network was visualized using a heatmap plot (Figure 4B). Correlations between the 5 MEs and traits of interest were evaluated using Pearson's correlation coefficient (Figure 4C).

Figure 4 WGCNA analysis, RT-qPCR and WB results. (A) Clustering dendrogram of DEGs revealing 5 module eigengenes (MEs), including MEturquoise, MEbrown, MEblue, MEyellow and MEgrey. (B) Heatmap plot of gene network displaying topological overlap, with light colors denoting low overlap and progressively darker colors higher overlap. Blocks of darker colors along the diagonal are the identified modules. (C) Module-trait relationship by Pearson correlation analysis. Each row corresponds to a ME and each column to a trait. Correlation coefficients (r value, scaled by colors) and corresponding P values are displayed in each rectangle. Positive correlation (Red); Negative correlation (Green); No correlation (White). (D) The qPCR results of Mrc1, Gpnmb, Aqp9, Trem1, Msr1, Slc26a4, Pde3a and Bmp6. Abbreviations: RNE, RNE = 2^−∆∆Ct. (E) The WB results of Mrc1, Gpnmb, Aqp9, Trem1, Msr1, Slc26a4, Pde3a and Bmp6 (n=4). *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001.

Abbreviation: ns, not significant.

The MEblue module was negatively correlated with FEV_0.3/FVC (r=−0.7, p=6e-04) and serum level of SOD (r=−0.77, p=7e-05), whereas it was positively correlated to serum level of MDA (r=0.76, p=9e-05). Genes Aqp9, Trem1, Mrc1, and Gpnmb were found to be co-expressed in this module and downregulated by both intensities. The MEblue and MEbrown modules were negatively correlated with the serum level of SOD, in which Msr1 and Slc26a4 were clustered and their expression could be specifically restored by EA-1mA. The MEturquoise module was negatively correlated with the serum level of IL-1β (r=−0.66, p=0.001), which contained Pde3a and Bmp6 that were specifically restored by EA-3mA.

Subsequently, we quantitatively assessed the mRNA and protein levels of Mrc1, Gpnmb, Aqp9, Trem1, Msr1, Slc26a4, Pde3a, and Bmp6 using RT-qPCR and WB, respectively. The differential mRNA expression profiles of these 8 genes are shown in Table 1. The RT-qPCR and WB results illustrated the same trend as the sequencing results, ie, Aqp9, Trem1, Mrc1, Gpnmb, Msr1, and Slc26a4 were upregulated, whereas Pde3a and Bmp6 were downregulated, in the COPD group relative to the NC group. Notably, EA treatment downregulated Aqp9, Trem1, Mrc1, and Gpnmb. Furthermore, Msr1 and Slc26a4 were exclusively downregulated by EA-1mA, whereas Pde3a and Bmp6 were uniquely upregulated by EA-3mA (Figure 4D and E).

Table 1 The Sequencing Results of Mrc1, Gpnmb, Aqp9, Trem1, Msr1, Slc26a4, Pde3a, and Bmp6

Gene Symbol	Log2FC(COPD vs NC)	P value	Log2FC(EA–1mA vs COPD)	P value	Log2FC(EA–3mA vs COPD)	P value
Mrc1	1.5535	1.5396E–07	−1.7458	4.51E–09	−1.3187	7.39453E–06
Gpnmb	2.3276	1.79E–14	−1.3347	5.74E–06	−1.1828	5.49581E–05
Aqp9	1.6667	9.71E–08	−1.1366	1.71E–04	−1.2811	2.93171E–05
Trem1	1.9775	5.93E–11	−1.6523	2.90E–08	−1.4452	1.1016E–06
Msr1	2.1001	3.59737E–12	−1.1764	6.13181E–05	–	–
Slc26a4	3.2948	2.08673E–25	−1.3314	6.15691E–06	–	–
Pde3a	−1.4851	5.65839E–07	–	–	1.2555	2.1645E–05
Bmp6	−1.1757	6.13532E–05	–	–	1.1369	0.000105384

Discussion

Previous research has convincingly demonstrated the EA efficacy in the management of COPD, yet the influence of different intensities of EA on COPD remains inadequately explored. In this study, we sought to address this knowledge gap by conducting a comprehensive investigation into the intensity–effect relationship of EA, specifically focusing on 1mA and 3mA, and elucidating its associated molecular mechanisms in a COPD rat model. To the best of our knowledge, this study represents the first systematic attempt to unveil the 1mA- and 3mA-EA effects on the transcriptional level, providing novel insights into EA modulation in COPD pathophysiology.

First, we successfully established the rat model of COPD by CS inhalation and LPS injection. Parameter evaluations, including FEV₁/FVC and FEV₃/FVC ratios, cytokine levels of IL-6, IL-1β, IL-10 and TNF-α, as well as antioxidant enzymes, eg, SOD and GSH-Px, collectively indicated severe inflammation and oxidative stress in rat lungs, affirming the fidelity of the COPD model to clinical manifestations. Pathological evaluation using HE and IHC staining further confirmed the morphologic change and inflammatory cell infiltration, providing additional insights into the inflammatory status of the lung tissue in COPD rats. Both 1mA and 3mA intensities improved lung airflow capacity, reduced proinflammatory cytokine levels, and decreased inflammatory cell infiltration, while also elevating antioxidant activity in COPD rats. However, these effects varied slightly between the two intensities, suggesting differential responses in addressing these indicators.

Next, we revealed the plausible mechanisms underlying the observed 1mA and 3mA-effect differences using transcriptome profiling. By hierarchical clustering, EA-1mA and EA-3mA groups displayed individual gene sets, specifically regulating 121 and 158 genes that were reversely expressed in the COPD vs NC group. Venn diagram analysis further identified a small number of genes uniquely presented in each intensity group, such as Msr1 and Slc26a4 in the EA-1mA vs COPD group, and Pde3a and Bmp6 in the EA-3mA vs COPD group. The WGCNA results categorized the DEGs into 5 modules, each establishing a module–trait relationship that was either negatively or positively correlated with lung function (FEV_0.3/FVC), inflammation (IL-1β), and oxidative stress (SOD, MDA). Additionally, the WGCNA results underpinned the potential targets we selected from the analyzed transcriptome data. This integrative approach strengthens the credibility of our identified targets, supporting their relevance in the intricate molecular landscape of EA effects on COPD.

Through functional annotation of overlapping and discrepant genes, KEGG and GO analyses unveiled enriched pathways and biological processes associated with the intensity differences. Common pathways such as vascular smooth muscle contraction, complement and coagulation cascades, and osteoclast differentiation were enriched for both 1mA and 3mA stimulation intensities. Notably, 1mA and 3mA-specific pathways emerged. EA-1mA likely initiated an anti-inflammatory and immune response via the regulation of cytokine production, neutrophil chemotaxis, phagosome, and natural killer cell-mediated cytotoxicity, among others. In the context of inflammation, both innate and adaptive immune systems are activated, leading to accumulation of circulating macrophages, neutrophils, and T lymphocytes.⁵⁰ EA-1mA treatment remarkably reduced inflammatory cell infiltration (including macrophages and T lymphocytes) in lung tissues and cytokine release in the serum of COPD rats, consistent with observations in other studies.^13,23,51 In contrast, EA-3mA was more involved in processes related to intracellular communication and signal transduction via cell adhesion, proteolysis, regulation of angiogenesis, neuroactive ligand–receptor interaction, cholinergic synapse, and cGMP-PKG signaling pathway. We infer that low-intensity stimulation influences on genes responsive to CS- and LPS-induced inflammatory signals, while high-intensity stimulation may not only enhance local anti-inflammatory effect but also achieve a remote effect of suppressing systemic inflammation through the somatosensory-autonomic reflex. Indeed, several studies have reported that electrical stimulation may suppress systemic inflammation partly through mechanisms such as the vagal-splenic sympathetic reflex^52,53 or the vagal-adrenal reflex.^22,54 This is due to the molecular and functional heterogeneities of somatosensory neurons. These neurons, when activated, drive specific autonomic pathways that are influenced by factors such as acupoint selectivity, disease state, and the intensity of stimulation.⁴⁴

Based on the analysis results and literature review, we chose Mrc1, Gpnmb, Trem1, and Aqp9 as target genes. MRC1 (macrophage mannose receptor 1) is a type I transmembrane protein that plays a crucial role in regulating inflammatory responses⁵⁵ and tissue remodeling.^56,57 In experimental colitis, Mrc1-deficient mice had attenuated inflammation compared with wild-type mice.⁵⁸ GPNMB (glycoprotein nmb) is an endogenous glycoprotein that functions as a feedback regulator of proinflammatory responses.^59,60 TREM-1 (triggering receptor expressed on myeloid cells 1), belonging to immunoglobulin superfamily, triggers and amplifies inflammation,⁶¹ especially through synergism with Toll-like receptor signaling.⁶² Inhibition of Trem1 effectively improved the injury in lung tissues of COPD mouse, and reduced the infiltration of macrophages.⁶³ AQP9 (aquaporins 9), an aquaglyceroporin membrane channel, is abundantly expressed in hepatocytes and immune cells, impacting glycerol metabolism and redox balance, immune response, inflammation, and oxidative stress.^64–66 It was significantly upregulated or highly expressed in lung tissues of COPD mice or patients,^46,48,49,63 or in blood samples of COPD patients.⁶⁷ Consistent with them, our results showed that the 4 genes were significantly upregulated after COPD modeling, and EA treatment was able to downregulate their expression to potentially prevent the enhanced inflammation, systemic oxidative stress, and impaired lung function. The RT-qPCR and WB results corroborated their mRNA expression and protein levels, aligning with the sequencing data.

We acknowledge several limitations in our study. Firstly, the relatively small sample size may limit the generalizability of our findings, and a larger sample size is crucial for endorsing our results for eventual clinical translation. Secondly, our study focuses solely on comparing the therapeutic effects and mechanistic differences between 1mA- and 3mA-EA in treating COPD rats. Future investigation should include more other stimulation intensities with different EA frequencies, to rigorously and thoroughly elucidate the dose–effect relationship and underlying mechanisms of EA in the treatment of COPD in rats. Additionally, our investigation lacks consideration of the frequency of EA, which is an important parameter that could impact the observed effects. Thirdly, while the COPD rat model offers insights, it does not fully recapitulate the gene profiles of humans. Therefore, conducting functional studies on the candidate genes is imperative for a more comprehensive understanding of the coordinated networks and pathways they influence.

Conclusions

Our results demonstrated that 1mA and 3mA intensities induce distinct gene expression patterns at the transcriptional level, associated with shared and 1mA vs 3mA-specific enriched pathways. Genes Mrc1, Gpnmb, Trem1, and Aqp9 emerge as promising targets, and further studies are needed to elucidate their functional consequences in COPD.

Abbreviations

COPD, Chronic obstructive pulmonary disease; EA, electroacupuncture acupuncture; WGCNA, weighted gene correlation network analysis; FEV_0.1/FVC, the ratio of forced expiratory volume in 0.1s to forced vital capacity; FEV_0.3/FVC, the ratio of forced expiratory volume in 0.3s to forced vital capacity; A/S, the ratio of the total alveolar area to the septum; IL-6, interleukin-6; IL-1β, interleukin-1β; IL-10, interleukin-10; TNF-α, tumor necrosis factor-α; SOD, superoxide dismutase; MDA, malondialdehyde; GSH-Px, glutathione peroxidase; IOD, integrated optical density; TPM, transcripts per million reads; DEGs, differentially expressed genes; log₂FC, log₂ fold change; FDR, false discovery rate; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; ANOVA, analysis of variance; MEs, module eigengenes; TOM, topological overlap matrix; NC, normal control group; EA-1, electroacupuncture with the current of 1 mA; EA-3, electroacupuncture with the current of 3 mA; RT-qPCR, Real-time quantitative PCR; CS, cigarette smoke; LPS, lipopolysaccharide; TCM, Traditional Chinese Medicine; HE, Hematoxylin-eosin; IHC, Immunohistochemistry; ELISA, Enzyme-linked immunosorbent assay; HRP, horseradish peroxidase.

Data Sharing Statement

The sequencing data in our manuscript have been deposited in the Bioproject database in NCBI, and the BioProject ID is PRJNA942343.

Ethics Approval and Consent to Participate

All procedures involving animals were approved by Animal Ethics Committee of Chengdu University of Traditional Chinese Medicine (2021-10) (Sichuan, China).

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

The authors declare that they have no competing interests.

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (Grant No. 82274664, 82305008) and the and Innovation Team and Talents Cultivation Program of National Administration of Traditional Chinese Medicine (No: ZYYCXTD-D-202003).