Introduction

The blood-brain barrier (BBB), which is crucial for drug delivery to the brain, mainly comprises specialized brain endothelial cells with tight junctions, a basement membrane, neurons, pericytes, and astrocytes. Its chemical composition and physical properties create an adaptable interface. Single-cell RNA sequencing (scRNA-seq) can systematically identify cell types and post-intervention molecular changes by analyzing the gene expression profiles of individual cells. This enables researchers to infer cell type-specific interactions and predict the molecular basis and functional significance of these ligand-receptor interactions.

Our previous research revealed that electroacupuncture (EA) employing the parameters of “2/100 Hz, 3 mA, 6 s-6 s, 40 min” at the head stimulation points (Baihui and Shuigou acupoints) yields a notably effective BBB opening effect. Nevertheless, the mechanism remains largely unknown (Zhang et al., 2018; Zhang et al., 2020; Ma et al., 2022; Lin et al., 2023; Ma et al., 2024; Dai et al., 2025). Here, we addressed this challenge by presenting a comprehensive single-cell molecular map of the rat frontal cortex, which covers approximately 98,338 cells from 10 samples, including freshly excised cortical brain tissue after EA intervention and responsive tissue from the control group. We characterized 23 clusters in the rat frontal cortex, including eight annotated cell types.

Methods

Animal model

Male Sprague–Dawley rats (Specific Pathogen Free), 8–10 weeks old and weighing 180–220 g, were obtained from Shanghai Sipul-Bike Experimental Animal Company. All mice were housed in the Zhejiang Chinese Medical University Laboratory Animal Research Center under a 12-h light/dark cycle with food and water given ad libitum. Rats were group-housed and were randomly assigned to experimental groups.

Materials and methods

Grouping and interventions

Rats were randomly assigned to the EA and CON groups (n = 5 per group). The EA group received EA; the CON group received the same non-stimulated binding treatment. Sterile disposable needles from Beijing Zhongyan Taihe Medical Equipment Co., Ltd. (Beijing, China) were used for EA. The GV20 (Baihui) and GV26 (Shuigou) acupuncture points were selected, with needle dimensions of 25 mm x 0.13 mm for GV20 and 16 mm × 0.07 mm for GV26. A custom relay, connected to an acupuncture point stimulator (HANS-200A, Nanjing Jisheng Co., Ltd., Nanjing, China), activated the needles in a 6 s on/off cycle. The current was set at 3 mA at a frequency of 2/100 Hz, and the stimulation lasted for 40 min.

Cerebral cortex issue preparation

At the conclusion of the experiments, rats were anesthetized using intraperitoneally administered pentobarbital (50 mg/kg). Subsequently, normal saline solution was circulated through the left side of the heart. The brain was removed via decapitation, and the frontal cortex tissue was immediately separated and placed in an ice-cold sterile RNase-free culture dish with an appropriate amount of calcium-free and magnesium-free 1 × phosphate-buffered saline (PBS). The tissue was then sliced into 0.5-mm³ pieces and washed with 1 × PBS. Non-target tissues, such as the blood cells and brain white matter, were carefully removed.

Single cells were isolated using a previously described method (Chi et al., 2024). Standard cDNA amplification and library construction were performed. The libraries were sequenced on an Illumina NovaSeq 6000 system (paired-end, 150 bp) by LC-BioTechnology Co., Ltd. (Hangzhou, China) with a minimum depth of 20,000 reads per cell.

Quantification and statistical analysis

Bioinformatics analysis

Sequencing results were demultiplexed and converted to the FASTQ format using Illumina bcl2fastq (version 2.20). The Cell Ranger pipeline (v6.1.1; <https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger>) handles sample demultiplexing, barcode processing, and single-cell 3’ gene counting. From five EA and five CON samples, 98,338 single cells were captured using the 10X Genomics Chromium Single Cell 3’ solution. The output was analyzed using Seurat (version 3.1.1) for dimensional reduction, clustering, and scRNA-seq data analysis. Quality control filtered cells with fewer than 500 expressed genes, UMI counts below 500, and >25% mitochondrial DNA-derived gene expression; as a result, 80,713 cells were included in the analysis.

To enhance data visualization, we used Seurat for the dimensionality reduction of all cells, followed by t-distributed stochastic neighbor embedding (t-SNE) for two-dimensional projection. The process included 1) applying the LogNormalize method from Seurat’s “Normalization” function to determine the gene expression values; 2) conducting principal component analysis (PCA) of normalized expression values, selecting the top 10 principal components (PCs) for clustering and t-SNE analysis; 3) using the weighted shared nearest neighbor (SNN) graph-based clustering method to identify clusters; and 4) selecting marker genes for each cluster using the Wilcoxon rank-sum test (default parameters “bimod” and the likelihood ratio test) via Seurat’s FindAllMarkers function. Marker genes were chosen based on their expression in more than 10% of cells within a cluster, with an average log value (fold change) greater than 0.25 (default parameter: 0.26).

Bioinformatics analysis was performed using the OmicStudio tools, which are accessible at https://www.omicstudio.cn/tool.

Doublet removal

Occasionally, two or more cells are encapsulated in the same gem during single-cell sequencing, resulting in abnormal barcodes that could affect the analysis (called doublets or multilets). Therefore, we used Doublet Finder software to identify and remove double and multiple cells to ensure the quality of the results, with the 7.5% doublet formation rate guideline recommended by 10X Genomics for single-cell analysis, aiming to identify the most likely doublets.

Classification of single-cell clusters

The single-cell subset classification analysis involved several steps. First, low-quality cells were removed, and expression values were normalized using the Log Normalize method in Seurat’s “Normalization” function. Next, PCA was performed on the normalized expression values to reduce the dimensionality and variables. The top 10 PCs were selected for clustering and subpopulation analyses. Finally, the Seurat software used a graph-theory-based clustering algorithm to divide the cells. This process included constructing a k-nearest neighbor clustering graph based on the Euclidean distance from significant PCs, optimizing intercellular distance weights using the Jaccard similarity and identifying cell clusters using an SNN module-optimized clustering algorithm.

Analysis of marker genes

Seurat used a bimodal likelihood ratio statistical test to analyze differentially expressed genes (DEGs) across various cell populations to identify genes that were upregulated in different clusters. The screening criteria for upregulated genes were as follows: 1) genes expressed in more than 10% of cells in the target or control subpopulations; 2) P-value ≤ 0.01; and 3) gene expression fold change (logFC) ≥ 0.26, meaning that the fold increase in gene expression was ≥ 2^0.26.

Cell-type annotation

Cell annotation results were defined using SingleR by employing a reference dataset derived from Mouse RNAseq Data. The correlation between the gene expression profiles of each cell and the dataset was calculated by assigning each cell a corresponding identification category. We also tallied the cell types identified within each cluster and their respective percentages, selecting the cell type with the highest correspondence for each cluster as the cell-type identification result.

Differential gene expression analysis

For differential analysis, the bimod algorithm was chosen, and the following parameters for screening DEGs were used: 1) P < 0.01, 2) logFC ≥ 0.26, and 3) the threshold for gene expression in at least one type of cell was 0.1.

Functional enrichment analysis

For the Ḍ significance enrichment analysis, we initially mapped all DEGs that were differentially expressed by various terms. We then calculated the number of genes belonging to each term and used hypergeometric tests to identify Gene Ontology entries significantly enriched in DEGs compared to the entire genome background. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis followed the same strategy. Enrichment scatter plots and bar charts were generated using OmicStudio tools at https://www.omicstudio.cn/tool.

Cell-cell communication analysis

CellPhoneDB was used to analyze cellular communication by predicting ligand-receptor interactions based on receptor expression in one cell type and ligand expression in another. For each gene, we calculated the percentage of cells expressing and the average expression of genes. Ligand-receptor pairs were included only if the proportion of cells expressing both genes exceeded 10%. For the polymers, the subunit with the lowest average expression represented the receptor in the statistical analysis. We compared all cell types in pairs within the filtered dataset by randomly arranging the cells into a new population with a default of 1000 arrangements. We determined the mean expression of ligands in a randomly arranged cell population and the average expression of receptors in interacting cell types and then calculated their mean. This process was repeated several times to obtain a normal mean distribution. Next, we calculated the actual average number of ligand-receptor pairs in the original cell population. The significance of the receptor-ligand pair was determined by the proportion of the calculated average values equal to or exceeding the actual average. We then sorted the interactions between cell types based on the number of significant ligand-receptor pairs for manual screening of biologically relevant interactions. As CellPhoneDB only recognizes human genes or proteins, we used the Ensemble tool to convert genes into their homologous counterparts before conducting further analyses.

Results

ScRNA-seq and clustering analysis

We used 10X Genomics technology to analyze the mRNA levels of individual cells isolated from the frontal cortex of rats (Figure 1a). We analyzed 98,338 individual cells from 10 different biological samples. A clear boundary between cells and non-cells was visible. The environmental background was very clean, with almost no environmentally free RNA (Figure 1b). We eliminated cells with potential double droplets, characterized by an unusually high number of detected genes, and unhealthy cells with typical mitochondrial mRNA loads exceeding 10%.

Figure 1.

Cells were isolated from two groups of 10 tissue samples (n = 5 per group), a blank control group (CON) and a specific mode EA stimulation group (EA). The cells were clustered using the Seurat package, and the t-distributed Stochastic Neighbor Embedding algorithm was used for visualization. Twenty-three clusters (cluster 0–cluster 22) were obtained by cell subgroup classification. The proportion of each cluster in the CON and EA groups is shown in Figure 1c and d.

Identification of cell types in different clusters in the CON and EA groups

The SingleR package was used to identify the predominant cell types. Cells were annotated according to the expression of canonical cell class markers (Figure 2a). We identified eight cell types (Figure 2b). In addition to endothelial cells, we identified astrocytes, microglia, oligodendrocytes, T cells, border-associated macrophages (BAMs), fibroblasts, and granulocytes in this dataset (Figure 2a and b). The combined analysis of cell clustering (Figure 1d) and cell identification showed that endothelial cells were enriched in two subpopulations, astrocytes were enriched in six subpopulations, and microglial marker genes were enriched in 10 subpopulations. Oligodendrocytes, T cells, BAMs, fibroblasts, and granulocytes each exhibited a single subcluster (Figure 2b and c).

Figure 2.

Evaluation of DEGs and pathways related to EA intervention on a holistic level

Based on the differences between the EA and CON groups, the cells with large clustering differences between the groups and those closely related to the composition and function of the BBB (endothelial cells, astrocytes, and microglia) were selected for differential analysis to further analyze the mechanism of EA intervention in the BBB.

Endothelial cells

Endothelial cells play a crucial role in BBB function. These cells create a continuous network of tight and adhesive connections along their contacts, forming a selective barrier that regulates the movement of substances inside and outside the brain. The DEGs were analyzed to identify important genes and pathways involved in maintaining this barrier. The expression of a total of 147 genes was significantly increased, and that of 160 genes was markedly downregulated (logFC = 0.26, P < 0.01). The expression fold changes of the top 10 genes most significantly upregulated and downregulated are shown in Figure 3a, and all DEGs are listed in Supplementary File 1.

Figure 3.

To perform enrichment analysis to evaluate the functionality of the 307 genes, we used the Gene Ontology (GO) terms, which encompassed biological processes (BPs), cellular components (CCs), and molecular function (MF), demonstrating the enrichment of 42 GO terms related to the BBB (P < 0.05, count ≥ 3) (Figure 3c). KEGG analysis yielded nine pathways that may be associated with BBB function: tight junctions, endocytosis, focal adhesion, the phosphatidylinositol 3-kinase (PI3K)/Akt signaling pathway, the mitogen-activated protein kinase (MAPK) signaling pathway, the hypoxia-inducible factor 1 (HIF-1) signaling pathway, the Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling pathway, ferroptosis, and cell adhesion molecules (Figure 3b).

Computational gene set enrichment analysis (q < 0.05) revealed eight significantly enriched unique pathways and identified the following BPs as being affected by EA: complement and coagulation cascades, ribosomes, lysosomes, cardiac muscle contraction, the rig I-like receptor signaling pathway, proteasomes, antigen processing and presentation, and cytokine-cytokine receptor interactions (Figure 4).

Figure 4.

DEGs in astrocytes

DEG analysis identified 24 downregulated and three upregulated genes, and the top 10 most upregulated and top three most downregulated genes are shown in Figure 5a. GO was performed. KEGG analyses revealed that DEGs were mainly enriched in calcium-related terms, adenosine triphosphate (ATP) and adenosine diphosphate synthesis and metabolism, and the nicotinamide adenine dinucleotide dehydrogenase complex. These molecules are prominent in certain pathways that may be connected to BBB function and permeability, including positive regulation of vascular permeability, cellular response to the vascular endothelial growth factor stimulus, endothelial cell chemotaxis, and the platelet-derived growth factor receptor signaling pathway (Figure 5b and c). Furthermore, these molecules were also involved in the regulation of the metabolic process of reactive oxygen species, the response to hypoxia, the response to hyperoxia, cell respiration, and the regulation of gaseous respiratory exchange.

Figure 5.

DEGs in the microglia following EA intervention

Among the DEGs observed in the microglia under EA intervention, the top 10 genes that were most significantly upregulated and downregulated are shown in Figure 6a. EA produced an immediate damage-associated molecular pattern response and the co-chaperone of the heat shock protein 70, DNAJ protein family (Dnaja1 and Dnajb1), related to the inflammatory response (Cao et al., 2014; Kovacs et al., 2017). According to the GO analysis, the terms associated with DEGs were mainly classified into five categories: calcium signaling, MAPK-related pathways, chemokine pathways, extracellular function, and angiogenesis (Figure 6c). Figure 6b shows the results of functional enrichment in the KEGG pathway database, which may be related to BBB function: the MAPK signaling pathway, endocytosis, HIF-1 signaling pathway, PI3K-Akt signaling pathway, tumor necrosis factor (TNF) signaling pathway, cell adhesion molecules, Wnt signaling pathway, and tight junctions.

Figure 6.

Endothelial cell subsets following EA intervention

Characteristic gene expression of EC_Clusters

The single-cell analysis identified 23 EC subgroups, of which four clusters (clusters 0, 2, 4, and 5) were dominant between the EA and control groups (Figure 7a). Therefore, further GO and KEGG analysis was performed to evaluate the characteristic gene expression patterns of these four subgroups (Figure 7c–f), and the bubble maps of the top 10 genes expressed in these four subgroups are shown in Figure 7b.

Figure 7.

The upregulated and downregulated genes in these four subgroups were further analyzed (EA vs. CON). The enriched pathways of DEGs in EC_Cluster0 and EC_Cluster2 exhibited a high degree of overlap (Figure 7c and d) and included infectious disease, the immune system, the MAPK signaling pathway, fluid shear stress, and atherosclerosis. In contrast, cluster 0 was enriched in platelet activation, whereas cluster 2 was enriched in the TNF signaling pathway genes. GO and KEGG results showed that the characteristic expression genes of EC_Cluster4 were mainly involved in immune system-related terms (Figure 7e). Similarly, the DEGs of EC_Cluster5 are primarily enriched in pathways related to the immune system, infectious diseases, and the regulation of cell population proliferation. (Figure 7f). In contrast to EC_Cluster4, EC_Cluster5 was also enriched in the interleukin (IL)-17 signaling pathway, fluid shear stress, and the atherosclerosis pathway.

Electroacupuncture altered EC_Clusters of interest

To examine the DEGs and pathways associated with EA exposure from the four subsets of EC, we compared the significantly upregulated or downregulated genes between the EA and CON groups and identified the most representative enriched pathways and processes (Supplementary Files 2 and 3). Neither the upregulated nor downregulated genes of EC_Cluster0 were independent of other clusters; that is, the differential genes of EC_Cluster0 represented the common influence of EA intervention, suggesting that this is the basic response of endothelial cells to the effect of EA.

Thirty-four genes were significantly upregulated, and 11 genes were significantly downregulated by EA in EC_Cluster0 (Supplementary File 2b). Pathway analysis revealed downregulation of cellular responses to extracellular stimuli and inflammatory responses (Supplementary Files 1, 2a and 2d). Genes in pathways that involve leukotriene biosynthetic processes, metabolic processes, and synthase activity were downregulated (Supplementary File 1). We found that EA-related genes related to interferon signaling were upregulated by EA (Supplementary Files 2b and 4), which affects brain endothelial function and barrier integrity (Oshima et al., 2001; Zhang et al., 2005; Mbofung et al., 2017). Interestingly, Stip1 regulated the Hsp90ab1 function (Li et al., 2014), and its mRNA was also regulated by EA in EC_Cluster0 (Supplementary File 4). Notably, Hsp90b1 and Hsp90ab1 were enriched in the fluid shear stress and atherosclerosis pathways (Supplementary File 2c).

The upregulated genes and corresponding functions of EC_Cluster2 were highly correlated with those of EC_Cluster0, with the addition of increased Psap, Mt-co3, Pdia4, and Atrx.1 expression (Supplementary Files 2g, h, j, and 5a) related to the extracellular matrix region and active transmembrane transporter activity (Supplementary File 6). Only 20 genes were downregulated in EC_Cluster2 (Supplementary File 5b). Among them, Nfkbia, Cxcl10, and Cxcl2 had the highest number of enriched pathways, mainly including the NF-κB signaling pathway, IL-17 signaling pathway, TNF signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway (Supplementary Files 2f, g, i, and 7).

After EA intervention, 15 unique upregulated genes and 13 downregulated genes were identified in EC_Cluster4 (Supplementary File 5). Importantly, downregulated gene enrichment analysis identified pathways related to BBB transport and barrier function. We also observed that terms and pathways related to the immune response were enriched (Supplementary Files 3a, 3b, 3d, and 8). The upregulated DEGs were mainly related to the inflammatory response and involved a variety of cytokine-related signaling pathways. Simultaneously, several typical signaling pathways were also observed, including the ERK cascade, MAPK, and PI3K pathways. Finally, the EA intervention in EC_Cluster4 also resulted in changes in vasculogenesis, astrocyte cell migration, pinocytosis, and calcium-mediated signaling (Supplementary Files 3b, 3c, 3e, and 9).

Among the downregulated genes in EC_Cluster5, only six genes overlapped with other subgroups, and the remaining 104 genes were downregulated (Supplementary File 5b). The downregulated genes were primarily related to ribosomes, suggesting active changes in protein synthesis and processing (Supplementary File 10). Among the upregulated genes in EC_Cluster 5, 21 coincided with other subgroups, and the remaining 44 genes (Supplementary File 5a) were enriched in the regulation of the establishment of the BBB, sodium-dependent organic anion transmembrane transporter activity, and sprouting, which are important for BBB formation and material transport, in addition to chemokine activity and immune cell chemotaxis (Supplementary Files 3g, 3h, 3j, and 11).

Nature of the microglial cell subsets and the effects of the EA intervention

Analysis of the microglial cell subsets yielded 15 cluster subsets, and marker genes showed different microglial subsets (Supplementary File 12); therefore, these 15 clusters were considered distinct clusters of microglial cells. We selected six clusters (MG_Cluster0, 1, 2, 3, 5, and 6) with large changes and annotated marker genes that were differentially expressed among these cell populations. GO annotation identified the top 20 enrichment terms with the most significant enrichment, whereas KEGG revealed the top 20 pathways with the largest number of genes and the 20 pathways with the smallest P-values. A Venn diagram displays the overlap between these pathways and entries (Supplementary Files 13 and 14).

Characteristic DEGs of MG_Clusters

The BP, CC, and MF pathways with the most significant (minimum Q-value) enrichment of the top gene of MG_Cluster1, as well as the top 20 pathways with the highest number of these genes, are displayed in Supplementary Files 13b, 14b, and 15–17. These entries involved the membrane dynamics, G protein-coupled receptor signaling pathway, cytokine-mediated signaling pathway, and extracellular space (Supplementary Files 13, 14a, and 15). KEGG analysis showed cytokine-cytokine receptor interaction, arachidonic acid metabolism, ECM-receptor interaction, signaling pathways regulating stem cell pluripotency, platelet activation, and JAK-STAT signaling (Supplementary Files 13b, 14b, 14c, 16, and 17). The uniquely enriched entries and pathways of MG_Cluster0 were mainly related to the intracellular anatomical structure, responses of cells to various cytokines, and regulation of transcription factors (Supplementary Files 13a, 14a, and 18). The important pathways were the FoxO and NF-κB signaling pathways (Supplementary Files 13a, 14b, 14c, 18, and 19). In contrast to the other clusters, the functions of MG_Cluster2 included the extracellular region, cellular response to organic substances, and lysosome and chemokine signaling pathways (Supplementary Files 13c, 14, 20–23). MG_Cluster3 had important functions that differed from those of the other clusters, summarized as negative regulation of amyloid fibril formation, positive regulation of tau protein kinase activity, and low-density lipoprotein particle receptor binding (Supplementary Files 13d, 14a, and 24). The total number of entries and pathways representing the unique functions of MG_Cluster5 was the largest (Supplementary File 14), and its functions mainly focused on the immune response (Supplementary Files 13e, 14a–c, 25–27). For MG_Cluster6, the top enriched pathways in MG_Cluster6 were not unique to other clusters (Supplementary Files 14b, 14c, 27–29), and GO enrichment analysis showed that they were mainly related to ATP binding, transcription, and protein folding (Supplementary Files 13f, 14b, 14c, and 28).

DEGs following electroacupuncture intervention in the MG_Clusters

We analyzed the coincidence and non-coincidence of the upregulated and downregulated DEGs in the six clusters (Supplementary Files 30a and 31a). Similarly, we performed an intersection analysis of GO entries with more than two DEGs and KEGG pathways with a Q-value ≤ 0.05 (Supplementary Files 30b, 31b, 31l).

Upregulated genes and enrichment analysis of MG_Clusters

Regarding the quantity of upregulated genes, MG_Cluster3 had the highest number (n = 52), followed by MG_Cluster6 (n = 38); the numbers in MG_Cluster5 and MG_Cluster0 were 14 and 13, respectively; the number in MG_Cluster2 was seven, while the number in MG_Cluster1 was only one (Supplementary File 30a). In contrast, GO analysis showed that MG_Cluster5 had the largest number of entries, independent of the other clusters (n = 53), followed by MG_Cluster6 and MG_Cluster0, with 18 and 13 entries, respectively. MG_Cluster2, MG_Cluster1, and MG_Cluster3 had similar numbers of independent enrichment entries; MG_Cluster2 had eight entries, and the remaining two clusters had seven entries (Supplementary File 30b). Regarding the KEGG pathways, for MG_Cluster0, all activated pathways were similar to those of the other clusters. In contrast, the number of upregulated gene enrichment pathways in MG_Cluster1 that do not overlap with other clusters was the largest (13 pathways). The enrichment of upregulated genes in MG_Cluster6 occurred in six pathways, whereas the featured pathways in MG_Cluster5, 3, and 2 were 3, 2, and 1, respectively (Supplementary File 30i).

Next, based on the above results, we analyzed the unique GO and KEGG enrichment of the six clusters of interest to study the functional changes in each cluster after the EA intervention. Among them, MG_Cluster0 was related to vascular function. MG_Cluster0 functioned as a focal adhesion and synapse (Supplementary File 30c). GO enrichment in MG_Cluster1 was mainly related to post-translational modifications of proteins. The KEGG results showed that Hsp90aa1 and Hsp90ab1 were also enriched in the PI3K/Akt signaling pathway (Supplementary File 30d and j). In MG_Cluster2, we identified genes upregulated in the cluster after EA intervention, which involved the basement membrane and extracellular region (Supplementary File 30e and j). Upregulated genes in MG_Cluster3 were mainly enriched in GO terms related to axon and myelin. KEGG analysis revealed that this cluster was related to ribosomes (Supplementary File 30f and j). In MG_Cluster5, the upregulated genes were mainly enriched in the cytosol and plasma membrane. Notably, the two genes upregulated in the cluster (Ccl3 and Ccl4) play key roles in and are mainly enriched in calcium and CCR chemokine receptor-related functions, as well as in the toll-like receptor signaling pathway (Supplementary File 30g and k). The enriched pathways in MG_Cluster6 are mostly related to mitochondria after EA intervention (Supplementary File 30h and k).

Downregulated genes and enrichment analysis of MG_Clusters

We observed that among the six clusters of independently downregulated genes, MG_Cluster3 and MG_Cluster6 involved far more genes than the other clusters, with 82 and 68 genes, respectively. The number of pathways identified in the GO and KEGG enrichments of MG_Cluster6 were also higher than those of the other five clusters (Supplementary File 31a, b, and l), including the cellular response to various stimuli. Another noteworthy pathway was the fluid shear stress and atherosclerosis pathway (Supplementary File 31h and k), because the fluid shear stress generated by directional shear of blood flow directly acts on brain microvascular endothelial cells.

Furthermore, the downregulated genes in MG_Cluster0 were mainly associated with focal adhesion, synapses, the response to angiotensin, and positive regulation of angiogenesis; however, no unique KEGG enrichment pathway was identified (Supplementary File 31c). The G protein-coupled purine nucleotide receptor and platelet activation were the key terms of MG_Cluster1 enriched (Supplementary File 31d and j). In the GO enrichment of MG_Cluster2, nine important genes were downregulated (Supplementary File 31g and l). In the unique GO of MG_Cluster3, we focused on the regulation of T cell and B cell migration. KEGG analysis yielded three pathways, including inositol phosphate metabolism, the C-type lectin receptor signaling pathway, and the phosphatidylinositol signaling system (Supplementary File 31e and l). Finally, the MG_Cluster5 downregulated genes enriched GO terms and pathways were mainly related to immune function (Supplementary File 31f and l).

Coordination of cell communication in the BBB following EA intervention

Because the BBB is orchestrated by the coordinated actions and interactions of several different cell types, cell communication analysis is necessary. First, we compared cell communication between the EA and CON groups (Supplementary File 32a–f) and further investigated the communication network between the two groups of cells and other cells (Supplementary File 32g and h). The number of protein interactions between the different cell types (P < 0.05) was calculated to generate a heat map (Supplementary File 32a and d), and the heat map values were log (natural logarithm) log-transformed (Supplementary File 32e and f). Since endothelial cells are the most direct cell type constituting the physical characteristics of the BBB, we focused on the degree of communication between various types of cells and endothelial cells. Upon comparing the EA and CON groups, we found that BAMs, microglia (MG), T cells, granulocytes, and endothelial cells demonstrated obvious differences in their interactions with endothelial cells (Supplementary File 32a and d). The results suggested that the cells with the most active relationship with endothelial cells were astrocytes, fibroblasts, and endothelial cells, whereas the largest changes between the two groups were granulocytes (GC) (Supplementary File 32b, e, g, and h).

To further explore the interaction between cell types, we analyzed the ligand-receptor pairs of other cells interacting with endothelial cells (Supplementary File 33a). Among the top 20 ligand-receptor pairs in the mean value, there were four pairs of interactions with endothelial cells that were different between the groups, namely CD74_APP, PLXNB2_PTN, LAMP1_FAM3C, and CD74_COPA. TNFRSF21_APP was present only in the CON group, whereas CSF1R_CSF1 was specifically present in the EA group (Supplementary File 33a). In the CON group, the ligand-receptor pair CD74_APP was not observed, suggesting that the EA intervention activated the ligand-receptor interaction of CD74_APP. In the communication between the AC and the EC, the ligand-receptor pair activated by the EA also contained TNFSF12_TNFRSF25, while GDF11_TGFR_AVR2B was blocked. Additionally, the mean interactions of TYRO3_GAS6, TNFSF13_TNFRSF1A, NOTCH1_JAG2, FGFR2_EPHA4, and FGFR2_PTPRR were consistently decreased in the EA group (Supplementary File 33b). The FB-EC communication attracted our attention because there were more ligand-receptor pairs than others. In the EA group, the communication between TGFB3 in the FB and its corresponding receptors in the endothelial cells was activated. Furthermore, apart from the higher mean PLXNB2_PTN in the EA group than that in the CON group, we observed that all other ligand-receptor pairs with mean changes were higher in the CON group (Supplementary File 33c).

Except for the TNFSF12_TNFRSF25 exclusively in the EA group and PLXNB2_SEMA4C in the CON group, in the ligand-receptor interaction between microglia and EC, we observed nine pairs of ligand-receptor pairs promoted by EA (Supplementary File 33d). In terms of interaction between granulocytes and endothelial cells (GC-EC), many ligand-receptor pairs were present only in the EA or CON group (Supplementary File 33e). We then focused on oligodendrocyte and endothelial cell (OC-EC) communication and found that the P-values of most ligand-receptor pairs decreased in the EA group, suggesting that EA rendered these interactions more pronounced. In addition, there were also communications that appeared exclusively in one group (Supplementary File 33f). In terms of the TC and EC communication, TGFB1_TGFbeta_receptor1 and TGFB1_TGFbeta_receptor2 appeared only in the EA group, whereas AR_Testosterone_byHSD17B12, COL11A2_integrin_a1b1_complex, IFNG_Type_II_IFNR, and SEMA4A_PLXND1 were only present in the CON group. The P-value of TGFB1_TGFBR3 was smaller in the EA group, while the P-value of CXCR6_CXCL16 was smaller in the CON group, and three pairs had a larger mean in the EA group (Supplementary File 33g). Finally, EA activated the communication of JAG1_NOTCH4 and TNFSF12_TNFRSF25 and blocked the communication of IGFBP3_TMEM219 and MERTK_GAS6 in BAM-endothelial cells (Supplementary File 33h).

Discussion

Increased expression of inflammatory markers such as C-C motif chemokine ligand (CCL-4 and CCL-3) and TMEM119 and decreased expression of cellular communication network family proteins results in increased BBB permeability (Lee et al., 2017; You et al., 2019; Davidson et al., 2022; Liu et al., 2022; Liu et al., 2024). Previous studies have shown that the AAV-PHP.eB delivery is facilitated by the receptor protein lymphocyte antigen-6E (LY6E) in humans (Fu et al., 2021). EA pretreatment reduces BBB permeability and upregulates the expression of Irf7 in healthy rats (Wu et al., 2019). Moreover, functional analysis in our study identified numerous terms and pathways associated with the BBB. Unexpectedly, the large clustering differences in microglia with and without EA intervention caught our attention. Indeed, microglia may serve as crucial priming cells for EA-induced opening of the BBB. Consequently, future research could focus on this aspect, and our laboratory is currently investigating this phenomenon, even though current research on microglia and BBB is mostly focused on the pathology and neuroinflammation of neurodegenerative diseases (Xu et al., 2015).

Higher basal Hmox1 expression in microglia is accompanied by increased oxidative stress, which directly affects the permeability of the BBB. Furthermore, Hmox1 regulates microglial polarization to the M2 phenotype (Groh et al., 2023). Of the top 10 genes that were downregulated after EA treatment, Mx1, Mx2, Ifi27, and Ly6e are immune-related genes, while Iba1 (encoded by Aif1) is a canonical microglial marker (Gandal et al., 2018). In addition, Pycard and Ltc4s function as mediators in apoptosis and inflammation (McConnell and Vertino, 2000; Agostini et al., 2004; Faustin et al., 2007; Fernandes-Alnemri et al., 2007; Bryan et al., 2009; Hasegawa et al., 2009; Hornung et al., 2009; Taxman et al., 2011; Lu et al., 2014; Guan et al., 2015; Shen et al., 2019; Shen et al., 2021).

The relationship between FK506 binding protein (FKBP) and the immunosuppressive drug FK506 has been examined; FKBP5 and FKBP4 are associated with neurological diseases, including Parkinson’s disease and AD (McConnell and Vertino, 2000; Bryan et al., 2009; Jiang et al., 2023). Notably, the basic neural basis synapses of memory and cognitive function are highly susceptible to synaptic phagocytosis by activated microglia and may improve cognitive function by reducing excessive activation of microglia (Reverte et al., 2024). We also observed enrichment of the pathway associated with cocaine addiction, which has been associated with microglia for synaptic adaptations in the nucleus accumbens synapses during cocaine withdrawal (Victor et al., 2022). Another pathway, thyroid hormone synthesis, is important for the microglial phenotype (Supplementary Files 13d, 14b, 14c, 29, and 34). The trigger receptor expressed in myeloid cell 2 (TREM2) is a cell surface receptor found in macrophages and microglia that responds to disease-related signals to regulate the phenotype of these innate immune cells. TREM2 is regulated by thyroid hormones, and its expression in macrophages and microglia is stimulated by thyroid hormones and synthetic thyroid hormone agonists (thyroid inhibitors).

In addition to considering the overall clustering of endothelial cells and microglia, which are cells with large differences in clustering, subsequent experiments should also focus on changes in the characteristics of the clustering of these cell subpopulations to explore whether there are certain unknown potential key subpopulations that play a critical role in BBB permeability. Similarly, cell communication analysis should not be confined to astrocytes and pericytes, which are typically regarded as cell types that are functionally associated with the BBB or comprised in the BBB. Differences in clustering and gene expression were also present in fibroblasts, BAMs, and T cells, although not so pronounced as to be focused on in this article, but it is possible that they are also involved in the regulation of BBB permeability.

This study had some limitations. First, as this investigation represents the initial analysis of single-cell transcriptome sequencing in the frontal cortex following EA intervention, there exists a limited understanding of the signaling pathways involved in EA-mediated BBB function. Second, a rat model was selected for this study due to considerations of EA operational precision. Third, additional protein analyses and functional validation are necessary to confirm the putative roles of each cellular phenotype. Finally, the possibility cannot be excluded that certain cells or cell types may be lost during tissue dissociation and cell enrichment; indeed, pericytes and neurons, which also play significant roles in the brain, were not enriched in this study. Nevertheless, this single-cell transcriptome sequencing in a rat model to investigate the target of action of EA-opening BBB could facilitate the study of BBB-opening-related mechanisms, which is crucial to numerous current investigations of exogenous drugs crossing the BBB into the brain.

Conclusions

The study of the molecular and cellular mechanisms that govern the function of the BBB is complex and requires a thorough understanding of the classification and molecular characteristics of cells, as well as detailed characterization of the spatial organization and interaction of molecularly defined cell types. This is because the spatial relationship between cells is the primary determinant of intercellular interactions and communication, which are mediated by juxtacrine and paracrine signaling. Our study revealed a unique frontal cortex-specific transcriptome signature and elucidated the cell types and clusters that mediate EA-induced transcriptional changes, including cellular communication. Our study highlights the characteristic changes in brain endothelial cells under EA intervention, as well as the cell types that synergistically participate in the structural and functional changes of the BBB, such as microglia and astrocytes, providing a framework for the mechanism of EA opening the BBB. The development of high-resolution, spatially resolved whole-brain cell atlases will provide invaluable resources for investigating BBB function during EA interventions. Future research in this area is necessary to elucidate the intricate details of intercellular communication and gain a deeper understanding of the mechanisms underlying BBB function.

Availability of data and materials

This study did not generate new unique reagents. The scRNA-seq data files are publicly available in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE272895) under accession number GSE272895.

Acknowledgements

We appreciate the support of the Key Laboratory of Acupuncture and Neurology of Zhejiang Province and the technical support from the Medical Research Center, Academy of Chinese Medical Sciences, Zhejiang Chinese Medical University. We thank Mr. Xushen Chen of the School of Public Health, Zhejiang Chinese Medical University, for his help with the language. We would also like to thank Editage (www.editage.cn) for English language editing.

Additional information

Ethics approval and consent to participate

The Animal Ethics Committee of Zhejiang Chinese Medical University approved the experimental procedures (Approval No.: IACUC-20220507-01). All experiments were performed in accordance with the ARRIVE guidelines. Consent to participate was not applicable.

Authors’ contributions

C.M.: Conceptualization, Data curation, Funding acquisition, Writing – original draft, and Visualization. Z.J., L.G., K.Q., and T.J.: Validation, Software, and Resources. J.Y., Z.P., and Q.C.: Project administration and Formal analysis, and M.D. and Q.Y.: Investigation and Methodology. X.L.: Funding acquisition and Writing – review & editing. All authors read and approved the final manuscript.

Funding

National Natural Science Foundation of China (82474626)

• Xianming Lin

National Natural Science Foundation of China (82174502)

• Xianming Lin

Science and Technology Department of Zhejiang Province (LQN25H270008)

• Congcong Ma

Zhejiang Provincial Health and Family Planning Commission (2022ZX009)

• Xianming Lin

Zhejiang Chinese Medical University (2022FSYYZZ10)

• Xianming Lin

Zhejiang Chinese Medical University (2023FSYYZQ11)

• Congcong Ma

Additional files

Supplemental files. Supplementary File 1: Pathway analysis for genes downregulated by SMES only in EC_cluster0 Supplementary File 2: Differential gene and enrichment analysis results of EC_Cluster0 and EC_Cluster2, related to Figure 7. Supplementary File 3: Differential gene and enrichment analysis results of EC_Cluster4 and EC_Cluster5, related to Figure 7. Supplementary File 4: Pathway analysis for genes upregulated by SMES only in EC_cluster0 Supplementary File 5: Venn diagram showing the differential gene intersection and uniqueness of the four clusters of interest in endothelial cells, related to Figure 7. Supplementary File 6: Pathway analysis for genes upregulated only in EC_cluster2 Supplementary File 7: Pathway analysis for genes downregulated only in EC_cluster2 Supplementary File 8: Pathway analysis for genes downregulated only in EC_cluster4 Supplementary File 9: Pathway analysis for genes upregulated only in EC_cluster4 Supplementary File 10: Pathway analysis for genes downregulated only in EC_cluster5 Supplementary File 11: Pathway analysis for genes upregulated only in EC_cluster5 Supplementary File 12: Subgroup clustering and functional enrichment of microglia, related to Figure 6. Supplementary File 13: The GO and KEGG results of the characteristic expression genes of the six subgroups of microglia, related to Figure 6. Supplementary File 14: Upset map of crossover and uniqueness of the enrichment function of the characteristic genes of 6 MG_Clusters of interest, related to Figure 6. Supplementary File 15: GO analysis for MG_cluster1 top genes only (S≥2) Supplementary File 16: KEGG analysis for MG_cluster1 top genes only (counts top 20) Supplementary File 17: KEGG analysis for MG_cluster1 top genes only (20 smallest P values) Supplementary File 18: GO analysis for MG_cluster0 top genes only (S≥2) Supplementary File 19: KEGG analysis for MG_cluster0 top genes only (counts top 20) Supplementary File 20: KEGG analysis for MG_cluster0 top genes only (20 smallest P values) Supplementary File 21: GO analysis for MG_cluster2 top genes only (S≥2) Supplementary File 22: KEGG analysis for MG_cluster2 top genes only (counts top 20) Supplementary File 23: KEGG analysis for MG_cluster2 top genes only (20 smallest P values) Supplementary File 24: GO analysis for MG_cluster3 top genes only (S≥2) Supplementary File 25: GO analysis for MG_cluster5 top genes only (S≥2) Supplementary File 26: KEGG analysis for MG_cluster5 top genes only (counts top 20) Supplementary File 27: KEGG analysis for MG_cluster5 top genes only (20 smallest P values) Supplementary File 28: GO analysis for MG_cluster6 top genes only (S≥2) Supplementary File 29: KEGG analysis for MG_cluster3 top genes only (counts top 20) Supplementary File 30: Functional enrichment analysis of unique upregulated genes in MG_clusters after electroacupuncture, related to Figure 6. Supplementary File 31: Functional enrichment analysis of unique downregulated genes in MG_clusters after EA, related to Figure 6. Supplementary File 32: The interactions between cells, related to the section “Coordination of cell communication in the BBB following EA intervention.” Supplementary File 33: Cell communication with endothelial cells, related to the section “Coordination of cell communication in the BBB following EA intervention.” Supplementary File 34: KEGG analysis for MG_cluster3 top genes only (20 smallest P values)