The Ethics Committee of Guangzhou University of Chinese Medicine approved all the experimental procedures in accordance with the Guide for the Care and Use of Laboratory Animals (Approval No. 20240301051).

Adult male C57BL/6 mice (aged 8–10 weeks, weight between 20 and 25 g, SPF grade) were acquired from the Animal Experiment Center of Guangzhou University of Chinese Medicine. Initially, the mice were randomly allocated into two groups for ChIP-seq: the middle cerebral artery occlusion (MCAO) model group and the electroacupuncture group. The RNA-seq data were derived from our prior sequencing data (Wu et al., 2021), and the mice in this dataset were categorized into three groups: the sham model group, the MCAO model group, and the electroacupuncture group, with three mice in each category. The animals were subsequently randomly distributed into the following three groups (sham group, MCAO model group, and MCAO model + acupuncture group) to validate the core super-enhancer target genes identified from the integrated RNA-seq and ChIP-seq analysis and to preliminarily investigate the epigenetic modification mechanisms of super-enhancers regulated by electroacupuncture in ischemic stroke. All the mice were maintained under regulated temperature and humidity conditions with a 12-h light–dark cycle and unrestricted access to food and water.

The MCAO model was established via the technique previously described (Engel et al., 2011; Bertrand et al., 2017). In brief, the mice were initially anesthetized via the inhalation of isoflurane to minimize discomfort. A midline incision was subsequently made under a stereomicroscope to expose the left common carotid artery, left internal carotid artery, external carotid artery, and vagus nerve. These structures were then separated and ligated as necessary. A suture plug measuring approximately 11 ± 0.5 mm in length was inserted into the internal carotid artery through the incision to prevent blood flow to the middle cerebral artery. One hour after ischemic occlusion, the suture plug was gradually withdrawn to restore the blood supply and induce ischemia–reperfusion injury. The mice in the sham surgery group underwent identical surgical procedures; however, MCAO was not performed. The criteria for evaluating the success of the MCAO model involve assessing the neurological function of animals using the Longa Score Scale (0–5 scales) (Longa et al., 1989). Mice in the MCAO model that receive a score ranging from 1 to 3 points are selected for inclusion in the study for further research. Comprehensive criteria for determining the success of the MCAO model, along with the Longa Score for each group, were detailed in Supplementary Table S1 and Supplementary Figure S1.

The acupoints “Baihui” and “Dazhui,” along with “Zusanli” and “Quchi,” were selected as targeted intervention acupoints on the basis of animal acupoint selection criteria outlined in “Names and Locations of Acupoints Commonly Used in Laboratory Animals Part 3: Mice” (Acupuncture-Moxibustion, 2021). Following standard sterilization protocols using iodine and alcohol, a disposable sterile acupuncture needle (0.16 × 13 mm) was aseptically inserted into each designated acupoint by the operator. An electroacupuncture apparatus was subsequently connected: “Baihui” and “Dazhui” formed one group for electrical stimulation, whereas “Zusanli” and “Quchi” constituted a second group. The electrical wave applied was continuous with a current set at 1 milliampere and a stimulation frequency of 2 Hz (Chang et al., 2018; Xing et al., 2018); slight vibrations were used to assess the response of the tissue surrounding each acupoint during the treatment sessions. Each mouse received acupuncture intervention for 20 min per session over a period of 7 days.

An initial behavioral evaluation of the mice was conducted 1 day prior to the MCAO surgery, with a follow-up assessment performed 7 days post-surgery. The specific methodologies employed for these behavioral assessments are detailed below.

All mice underwent careful anesthesia, followed by cardiac perfusion with chilled phosphate-buffered saline (PBS), and were subsequently fixed in a 4% formaldehyde solution. The brains were promptly extracted and subjected to overnight fixation in the same 4% formaldehyde solution at 4°C. Dehydration was then performed using a sucrose gradient. From each brain, with a primary focus on the striatum, for immunofluorescence analysis to evaluate the infarcted area. Specifically, the sections were incubated overnight at 4°C with a rabbit anti-NeuN polyclonal antibody (Proteintech, United States; 1:500), followed by incubation at room temperature with a goat anti-rabbit secondary antibody (Abcam, United States; 1:1000). Imaging was conducted using a Nikon confocal microscope. The percentage of the infarcted area within the mouse brain was quantified using the following formula: Infarct area % = (Area of neuronal loss in each section)/(Total brain area of each section) × 100%. Additionally, to identify and quantify the surviving neuronal cells, fluorescent images of the tissue sections were captured and analyzed using a confocal microscope. Image J software was employed to determine the number of positive neuronal cells in both the affected cortical and striatal regions.

In brief, RNA was extracted from brain tissue via TRIzol, and quality was assessed with a NanoDrop and Agilent Bioanalyzer. mRNA was purified with oligo(dT) beads and converted to cDNA, followed by library construction per the manufacturer’s standards. The PCR products were denatured and circularized to create the final library, which was amplified to produce DNA nanoballs (DNBs). DNBs were loaded onto a nanoarray for single-end 50-base reads on the BGISeq500 platform. Transcript sequences were normalized to fragments per kilobase of transcript per million mapped reads (FPKM), and differential gene expression was analyzed with DESeq2. Genes with Q ≤ 0.05 or FDR ≤ 0.001 were identified as differentially expressed genes (Cheng et al., 2019). The fold change (FC) value was employed to distinguish between upregulated and downregulated genes.

Brain tissue samples (2–3 g) from each group were collected and placed in a PBS solution. The samples were subsequently ground and filtered through a cell strainer. The samples were treated with 1% formaldehyde to induce cross-linking, followed by the addition of glycine to terminate the cross-linking reaction. The cross-linked DNA–protein complexes were sonicated, after which agarose gel electrophoresis was performed for detection. First, protein G beads were washed before the addition of ChIP dilution buffer and the target protein antibody H3K27ac to form the antibody-target protein–DNA complex. Multiple washes were conducted to remove non-specifically bound chromatin and purify the complex. Crosslinking was then reversed, and the resulting DNA fragments were purified. Finally, the blunt ends of the DNA fragments were repaired, followed by the ligation of sequencing adapters. The PCR products were amplified, and selected fragments were used to construct the sequencing library, which was subsequently sequenced via the Illumina NovaSeq 6,000 PE150 platform.

Cutadapt (v2.5) was used to trim adapters and filter out low-quality sequences. The sequencing reads were aligned to the mouse genome mm10 with Bowtie2 (v2.3.5.1) via default parameters. Regions of H3K27ac enrichment were calculated with MACS2 (v2.1.2) using the parameters m 5 50. The MACS peaks of H3K27ac served as enhancers for identifying super-enhancers.

The enhancers were concatenated, and super-enhancers were discerned via ROSE¹ as previously described (Whyte et al., 2013). In summary, this algorithm merges constituent enhancers if they lie within a specified distance and ranks the enhancers on the basis of their input-subtracted H3K27ac signal. It subsequently distinguishes super-enhancers from typical enhancers by pinpointing an inflection point of the H3K27ac signal in relation to the enhancer rank. ROSE was executed with a stitching distance of 12,500 bp, meaning that we permitted enhancers within 12,500 bp to be amalgamated. Furthermore, peaks within ±2,500 bp of an annotated transcription start site were omitted from stitching. For enhancer annotation, we utilized the homer (v4.10.4) utility annotatePeaks.pl, which links peaks to their adjacent genes.

The GO and KEGG pathway enrichment analyses of the genes linked to super-enhancers were conducted via the clusterProfiler R package (v3.6.0). A p-value ≤0.05 was considered indicative of enriched GO terms and KEGG pathways.

Key super-enhancer-driven genes were identified via the random forest algorithm, with core genes selected via the R package “randomForest.” To find the optimal error rate and the ideal number of stable trees, we calculated the error rate for each tree from 1 to 500. The random forest method was subsequently employed to screen for key super-enhancer-driven genes, utilizing the minimum Gini (MDG) in the random forest algorithm to assess classification accuracy and dimensional significance values. The top genes that presented the highest MDG scores were designated core super-enhancer genes. Furthermore, the Friends analysis method is extensively utilized for the identification of critical genes. In order to pinpoint critical core super-enhancer genes, the R package GOSemSim (Yu, 2020) was employed to compute the functional correlations among the primary super-enhancer-driven genes. The genes with the highest rankings were designated as key super-enhancer genes. Subsequently, we conducted an intersection analysis of the top three significant genes identified by both the random forest and Friends analysis methods to determine the core super-enhancer-driven genes. Additionally, the signal peaks of H3K27ac associated with key super-enhancer-target genes were visualized and analyzed utilizing the Integrative Genomics Viewer (IGV, version 2.16.2).

The Coltron algorithm we devised was employed to identify essential regulatory circuits, utilizing H3K27ac ChIP-seq read data aligned with the mouse genome, along with the ChIP-seq peaks generated by MACS (Zhang et al., 2008) and an enhancer list generated by ROSE (see Footnote 1) (Whyte et al., 2013). The process of identifying core transcriptional regulatory circuits (CRC) is as follows: (1) determine the nearest active genes using the outcomes from the super-enhancer regions and the alignment results; (2) exclude the segment of the nearest active genes that function as transcription factors and identify the transcription factors governed by super-enhancers; (3) gather transcription factor motifs from the TRANSFAC and JASPAR databases; and (4) utilize FIMO in MEME62 to examine the sequences of super-enhancer elements for motifs. The p-value threshold for FIMO is established at 1*10^⁻⁴. (5) When a transcription factor is capable of binding to the super-enhancer element that regulates it, it is classified as a self-regulating transcription factor. Transcription factors assigned to super-enhancers, whose component set encompasses at least three instances of DNA sequence motifs corresponding to their own protein products, are designated self-regulating transcription factors. From the group of self-regulating transcription factors, those capable of binding to the super-enhancers of other self-regulating transcription factors are predicted via the same criteria outlined above. All the possible fully interconnected self-regulating circuits are subsequently constructed through recursive identification. After multiple configurations of fully interconnected self-regulating circuits of transcription factors are established, the most representative fully interconnected self-regulating circuit is chosen as the CRC model (Saint-Andre et al., 2016; Feng et al., 2023).

Damaged tissues from each group of mice were frozen in liquid nitrogen and then pulverized to prepare protein extracts. A specific amount of primary antibody targeting the gene of interest was added to the lysed tissues for western blotting (WB) protein quantification. Brain sections from the affected areas of each group of mice were subjected to immunofluorescence. A primary antibody from rabbits that targets the gene of interest was added to the brain tissue sections. This was followed by the addition of a corresponding fluorescent secondary antibody for immunofluorescence. The fluorescence of the target gene was examined via laser confocal microscopy, and the expression levels in each group were quantified with Fiji ImageJ².

The data were analyzed via GraphPad Prism 9 (GraphPad Software, San Diego, CA, United States) and R version 4.2.1³. For normally distributed data with homogeneous variances across groups, the results are presented as the means ± standard errors of the means, and one-way ANOVA was conducted, accompanied by Bonferroni correction for post hoc comparisons. In instances where variances were not homogeneous, Dunnett’s T3 test was employed; a p-value of less than 0.05 was considered statistically significant.

Before the commencement of MCAO surgery, there were no significant differences in the time required to remove the adhesive tape among the three experimental groups, suggesting that the groups were comparable at baseline (all p > 0.05). After a period of 7 days, the electroacupuncture group exhibited a significantly shorter duration for tape removal compared to the MCAO model group, indicating improved motor function in the electroacupuncture-treated mice. This observation suggests that electroacupuncture may mitigate neurological impairments and enhance motor performance following ischemic stroke (Figure 1A). The foot slip test is utilized as a primary measure of limb motor coordination in the context of ischemic stroke. Prior to MCAO surgery, statistical analyses indicated no significant differences in foot slip percentages among the three groups of mice, confirming their comparable baseline characteristics (p > 0.05). Seven days post-injury, the MCAO + EA group showed a reduced incidence of foot slips compared to the MCAO control group, with statistically significant differences observed (p < 0.05, Figure 1B). This suggests that the intervention may have a beneficial effect in improving motor coordination in mice subjected to MCAO.

As demonstrated in Figure 1C, the MCAO model exhibited significant neuronal loss compared to the sham surgery group. However, after a seven-day treatment period, electroacupuncture significantly decreased the extent of infraction areas relative to the MCAO model group (p < 0.05). Furthermore, neuronal populations within the brain tissue were evaluated using NEUN staining. The ischemic brain injury model was associated with neuronal loss in the cortical and striatal regions. Nonetheless, electroacupuncture treatment resulted in a significant reduction in neuronal loss in both the cortical and striatal regions compared to the model group (p < 0.05) (Figures 1C–F). Therefore, electroacupuncture treatment appears to promote neuronal survival following cerebral ischemic injury, providing neuroprotective effects.

This study aimed to investigate potentially differentially regulated super-enhancer-driven genes modulated by electroacupuncture in ischemic stroke through ChIP-seq analysis. Considering that super-enhancers play crucial cis-regulatory roles within cells, we utilized the ROSE tool to identify super-enhancers and predict their nearest associated genes on the basis of chromosomal spatial positioning. In total, 77 super-enhancer-associated genes were identified in the electroacupuncture intervention group, whereas 403 genes were detected in the MCAO model group (Figures 2A,B). DeepTools software was employed to analyze and visualize the signal intensity across super-enhancer regions for each experimental group (Figures 2C,D).

The GO and KEGG pathway analyses were performed to examine the differentially super-enhancer target genes between the experimental group (EA) and the model group. The results revealed significantly enriched BP terms, including ‘regulation of synaptic plasticity,’ ‘axon genesis,’ ‘regulation of neurogenesis,’ ‘neuron projection extension,’ and ‘dendrite development’. Additionally, GO cellular component (CC) enrichment analysis revealed significantly enriched GO CC terms such as ‘postsynaptic specialization,’ ‘postsynaptic density,’ ‘asymmetric synapse,’ ‘neuron to neuron synapse,’ and ‘cell leading edge’. The top five enriched GO molecular function (MF) terms included ‘actin binding,’ ‘GTPase regulator activity,’ ‘nucleoside-triphosphatase regulator activity,’ ‘guanyl-nucleotide exchange factor activity,’ and ‘actin filament binding’. Furthermore, KEGG signaling pathways such as Axon guidance, Dopaminergic synapse, Glioma, Regulation of lipolysis in adipocytes, Sphingolipid signaling pathway, Wnt signaling pathway, Endocytosis and HIF-1 signaling pathway were significantly enriched in super-enhancer-related genes. The specific details are shown in Figures 3A,B.

Since super-enhancers are associated with gene transcription, genes that exhibit abnormal expression may represent key pathogenic factors in various diseases. The comparison of super-enhancer-driven genes between the model group and the electroacupuncture group revealed that 41 genes were uniquely regulated by super-enhancers in the electroacupuncture group [namely differential SE-related genes (gained super-enhancers after EA)], whereas 367 genes were regulated solely by super-enhancers in the model group [namely differential SE-related genes (lost super-enhancers after EA)] (Figure 3C). The genes regulated by super-enhancers in the electroacupuncture group, along with those lost super-enhancers by electroacupuncture compared with those in the model group, may serve as potential target genes for electroacupuncture. Super-enhancers facilitate gene expression by enriching transcription factors, which leads to the upregulation of critical genes. We compared the upregulated genes with those associated with super-enhancers. RNA sequencing analysis revealed the expression profiles of the differentially expressed genes within each group, revealing 2,702 upregulated genes in the model group and 53 upregulated genes in the electroacupuncture group. The analysis identified key regulatory genes associated with super-enhancers by examining the overlap between differentially expressed super-enhancer-related genes and upregulated genes in RNA-seq data. Specifically, this involved two comparisons: first, the intersection of genes associated with lost super-enhancers following EA with those upregulated in the model group relative to the EA group; and second, the intersection of genes associated with gained super-enhancers post-EA with those upregulated following the electroacupuncture intervention. Notably, the analysis revealed no overlapping genes between the super-enhancers gained post-EA and those upregulated in the RNA-seq data following the EA intervention. However, a comparison of genes associated with lost super-enhancers post-EA and upregulated genes in the MCAO group from RNA-seq identified 18 overlapping differentially expressed genes, which may serve as super-enhancer target genes, potentially exerting a suppressive role in the regulation of ischemic stroke via electroacupuncture. And these 18 selected super-enhancer-driven differentially expressed genes underwent were further analyzed and validated to confirm their significance (Figures 3D,E and Supplementary Table S2).

The significance ranking of 18 differentially expressed super-enhancer-driven genes was determined via random forests, with the top-ranked super-enhancer-driven genes identified as “HDAC7,” followed by “ARAP3” and “IRF2BP2,” as illustrated in Figure 4A. The Friends analysis revealed that HDAC7 exhibited a robust correlation with other super enhancer-driven genes and may therefore play a pivotal role in ischemic stroke (Figure 4B). Given that super-enhancers are intricately linked to histone acetylation, and the intersection of random forests and Friends analysis indicated that HDAC7 is a crucial super-enhancer-driven gene (Figure 4C). The histone deacetylase HDAC7, upregulated by super-enhancers, is highly expressed in MCAO model mice according to the RNA-seq data, leading to the preliminary hypothesis that HDAC7 could be a pivotal pathogenic factor in ischemic stroke. Moreover, electroacupuncture may address ischemic stroke by modulating the expression of HDAC7, which is influenced by super-enhancers (Figure 4).

Through WB and immunofluorescence, the expression of HDAC7 regulated by super enhancers in both the model group and the electroacupuncture group was confirmed. The findings revealed that the protein and fluorescence expression levels of HDAC7 were significantly greater in the model group than in the electroacupuncture group and the sham surgery group (p < 0.05). Initial validation suggested that HDAC7 was aberrantly overexpressed in the ischemic stroke group, whereas electroacupuncture intervention had a suppressive effect on the fluorescence intensity and protein level of HDAC7 (Figures 5A–D). Furthermore, the IGV plot and H3K27ac peak value analysis revealed that H3K27ac (a marker for super-enhancers) levels were significantly more enriched in the model group than in the electroacupuncture group (Figures 5E,F).

SEs and open chromatin regions are frequently enriched with specific transcription factors (TFs) that collaboratively regulate cellular transcription. To this end, we developed a super-enhancer-mediated transcriptional regulatory model for ischemic stroke by integrating H3K27ac data (Figure 6A). Using ChIP-seq experimental data to predict core transcription factors, we calculated the total degree (including both in-degree and out-degree) of genes encoding CRCs through machine learning. Genes with a relatively high total degree are considered potential core transcription factors. The primary core transcription factors included SOX8, KLF13, FOXK1, FOXG1, RARA, NR2F1, SMAD3, ZBTB16, and KLF7 (Figure 6B). The following sets were among the top-ranked CRCs identified by machine learning: [‘RARA,’ ‘KLF13,’ ‘TCF4,’ ‘FOXG1,’ ‘ZBTB16,’ ‘NR2F1,’ ‘FOXK1,’ ‘MEF2D,’ ‘HIVEP2’], [‘RARA,’ ‘KLF13,’ ‘SMAD3,’ ‘FOXG1,’ ‘ZBTB16,’ ‘NR2F1,’ ‘FOXK1,’ ‘MEF2D,’ ‘HIVEP2’], [‘SOX8,’ ‘KLF13,’ ‘TCF4,’ ‘FOXG1,’ ‘ZBTB16,’ ‘NR2F1,’ ‘FOXK1,’ ‘MEF2D,’ ‘HIVEP2’], [‘SOX8,’ ‘KLF13,’ ‘SMAD3,’ ‘FOXG1,’ ‘ZBTB16,’ ‘NR2F1,’ ‘FOXK1,’ ‘MEF2D,’ ‘HIVEP2’] (Figure 6C). These factors may represent key regulatory circuits of core transcription factors influencing pathogenic genes in the MCAO model (Figure 6D). Detailed information is presented in Figure 6 and Supplementary Table S3.