The experimental animal protocol was approved by the Animal Experiment Ethics Committee of Zhujiang Hospital of Southern Medical University (LAEC-2020-235) and was carried out following the regulations of experimental animal management of Zhujiang Hospital of Southern Medical University. Specific pathogen-free (SPF) healthy male Sprague Dawley (SD) rats (weighing 280–300 g and 6–8 weeks old) were provided by the Laboratory Animal Center of Southern Medical University (license number: SCXK (Guangdong) 2016-0041). The rats were fed and drank water ad libitum. They were raised in an environment with a temperature of 24 $\pm$ 1 °C, 40 $\pm$ 5% humidity, and a 12 h light-dark cycle. Six rats were randomly divided into the Sham and middle cerebral artery occlusion (MCAO) groups, with three in each.

The rats were weighed and intraperitoneally injected with Avertin (0.2 mL/10 g, reagent center of Southern Medical University, Sigma-Aldrich, City of Saint Louis, MO, USA) [22]. After anesthesia, the rats were subjected to MCAO according to the method described by Longa [23]. Briefly, after exposing the common carotid artery (CCA), external carotid artery (ECA), and internal carotid artery (ICA) on the right side of the rat, a silicone-coated nylon monofilament (diameter, 0.3 mm) was inserted into ECA and passed through ICA until the middle cerebral artery was blocked. After ischemia for 1.5 h, the suture was slowly removed, and the wound was disinfected and sutured. The rats in the sham group were operated on according to the above procedure, but no sutures were inserted.

Rats were randomly selected from each group 24 h after reperfusion for euthanasia [24]. After removing the brains, five 2-mm thick slices were serially transected in the brain matrix device. Brain slices were placed in a 2% 2, 3, 5-Triphenyl-2H-Tetrazolium Chloride (TTC) (T8877, Sigma-Aldrich, Riverdale, NJ, USA) solution and incubated at 37 °C for 15 min in the dark. Finally, sections were photographed, with red representing normal brain tissue and white representing infarcted areas.

Rats in Sham and MCAO groups were euthanized, and their brains were removed to obtain cerebral cortical cells for scRNA-seq. One brain tissue specimen meeting the study requirements was selected from each group. The isolated cerebral cortex was prepared as a single-cell suspension. After gently mixing 10 µL of the cell suspension and 10 µL of 0.4% Trypan Blue Stain (T10282, Thermo Fisher Scientific, Waltham, MA, USA), 10 µL of the mixture was quickly added to Countess® II Cell Counting. Cells were counted using a Countess® II Automated Cell Counter (C10228, Thermo Fisher Scientific, Waltham, MA, USA). The cell concentration in sham group was 1.65 $\times$ 10${}^6$ cells/µL, and 86% were living cells, while the cell concentration in MCAO group was 1.43 $\times$ 10${}^6$ cells/µL and 88% were living cells. If the proportion of viable cells was greater than 80%, it was considered a qualified sample, and subsequent experiments were performed after adjusting the cell concentration to 1000 cells/µL.

Single-cell samples of rat cerebral cortex were sequenced using BD Rhapsody Single Cell Analysis System (Becton, Dickinson and Company, Lake Franklin, NJ, USA) [25, 26]. The prepared single-cell suspension was added to a 20W+ microwell honeycomb plate. Beads with a unique molecular identifier (UMI) and cell barcode in the microplate were combined with the cells to capture and identify single cells. Following cell capture, total RNA was extracted from each well and reverse-transcribed into cDNA. The constructed cDNA library was sequenced on the Illumina NovaSeq 6000 sequencing platform (Illumina Inc., San Diego, CA, USA) to obtain the sequencing data. A count of $\ge$1 indicated that the gene was expressed. After excluding the genes expressed in less than three cells in the sequencing results, unqualified cells were excluded according to the following requirements: (1) Cells with more than 6000 or less than 200 genes were removed. (2) Cells with a mitochondrial gene ratio greater than 10% were removed. (3) The cells were removed with a hemoglobin gene ratio of $>$0.1%.

After processing the sequencing results, Wilcox algorithm was used to compare cell subpopulation differences among the samples. Differentially expressed genes (DEGs) with fold change (FC) $>$2 and p-value 0.05 were obtained in cell subpopulations. The top 40 DEGs were selected for gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis to clarify the biological processes and functions of cell subpopulations, according to the order of p-value from small to large (Fig. 1).

Fig. 1.

Technology road map.

DESeq2 package was used to determine the read counts of differentially expressed genes. The Student’s two-tailed t-test was used to compare the two groups. All statistical analyses were performed using SPSS software (version 20.0, IBM Corp., Armonk, NY, USA), and figures were drawn using GraphPad Prism (version 5.0, GraphPad Software Company, San Diego, CA, USA) and Photoshop CC (version 2019, Adobe Systems Incorporated, San Jose, CA, USA). A p-value 0.05 was defined as statistically significant.

MCAO rat models were constructed and evaluated by TTC staining using scRNA-seq to detect changes in cell subpopulations in the cerebral cortex of acute ischemic rats. Fig. 2A displays that the cerebral infarction area of MCAO rat was white, and the normal tissue was red, indicating that MCAO rat model was successfully constructed. We then isolated the cortex of rats in the sham and MCAO groups and prepared cell suspensions for scRNA-seq (Fig. 2B). After sequencing, the results were preprocessed and normalized using Seurat R package. Unqualified cells were removed (Fig. 2B–D). The results demonstrated that 4,871 cells in the Sham group and 6,836 cells in MCAO group were identified by scRNA-seq (Fig. 2B–D). The filtered cells were subjected to cell cycle analysis. The cell cycle score of each cell was calculated according to the expression of the cell cycle genes to determine the cell cycle of each cell. Visual principal component analysis (PCA) data were used to evaluate cell cycle effects. If the cell cycle effect was excessive, a linear regression was used. The results indicated that the cell cycle effect was insignificant in this experiment, proving that the cell cycle effect had no significant effect on subsequent analysis (Fig. 2E–H).

Fig. 2.

Experimental scheme and data preprocessing. (A) 2, 3, 5-Triphenyl-2H-Tetrazolium Chloride (TTC) staining results after 24 h of middle cerebral artery occlusion (MCAO). The cerebral infarct area was stained white, and the normal tissue was red, indicating that the model was successfully established. (B) Schematic diagram of scRNA-seq workflow. (C) Violin diagram of the number of genes, the proportion of mitochondrial genes, and the proportion of hemoglobin genes in the samples of sham and (D) MCAO groups. (E,F) Principal component analysis (PCA) of cell cycle-related genes and all genes in the Sham group. (G,H) PCA of cell cycle-related genes and all genes in the MCAO group. In the PCA diagram, red dots represent cells in G1 phase, green dots represent cells in G2 or M phase, and blue dots represent cells in S phase.

In this study, the brain cortex cells of sham and MCAO groups were identified using scRNA-seq based on BD Rhapsody system. After unbiased clustering analysis, the distribution of cells derived from the cerebral cortex of rats in sham and MCAO groups in each cell cluster was identified and visualized using t-distributed stochastic neighbor embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction algorithms, respectively (Fig. 3A,C). Subsequently, cell types were annotated using SingleR, and 21-cell clusters with unique gene expression patterns were defined (Fig. 3B,D). Fig. 3E depicts that although the number of cells and gene expression levels were different in each cell cluster in the two groups, the cell clusters were similar in the two groups. Notably, these 21-cell clusters were divided into 11 major cell lineages (Fig. 3F): microglia (clusters 0, 1, and 16), astrocytes (cluster 2), oligodendrocytes (clusters 3, 9, and 10), endothelial cells (clusters 4, 6, 11, 12, and 17), macrophages (cluster 5), ependymal cells (clusters 7 and 20), T cells (cluster 8), fibroblasts (clusters 13 and 19), monocytes (cluster 14), granulocytes/monocytes (cluster 15), and fibroblast activated (cluster 18). In conclusion, scRNA-seq data identified 21-cell clusters with high intercellular heterogeneity in the rat cortex.

Fig. 3.

Twenty-one cell clusters were identified in the sham and MCAO groups in the rat cortex. (A) t-distributed stochastic neighbor embedding (t-SNE) plot of cells in the sham and MCAO groups, where red dots represent cells derived from the sham group, and green dots display cells derived from MCAO group (n = 1). (B) t-SNE plot of 4871 cells in the sham group and the t-SNE plot of 6836 cells in MCAO group. The labels on the right were 21-cell clusters with corresponding colors, of which 0 to 20 represent microglia, microglia, astrocytes, oligodendrocytes, endothelial cells, macrophages, and endothelial cells, ependymal cells, T cells, oligodendrocytes, oligodendrocytes, endothelial cells, endothelial cells, fibroblasts, monocytes, granulocytes/monocytes, microglia, endothelial cells, activated fibroblasts, fibroblasts, and ependymal cells, respectively. (C) Uniform Manifold Approximation and Projection (UMAP) diagram of cells in the sham and MCAO groups, where red dots represent cells derived from the sham group, and green dots show cells derived from the MCAO group (n = 1). (D) Uniform Manifold Approximation and Projection (UMAP) diagram of 4871 cells in the sham group and Uniform Manifold Approximation and Projection (UMAP) diagram of 6836 cells in the MCAO group. (E) t-SNE plot of 21 cell clusters in the sham and MCAO groups. (F) t-SNE plot of 11 major cell lineages. t-SNE, t-distributed stochastic neighbor embedding.

We performed marker gene analysis of the 21-cell clusters to gain further insight into gene expression heterogeneity in each cell subpopulation in the MCAO rat cortex. The Wilcoxon rank-sum test was used to analyze genes in all cell clusters and then scored using One-vs-Res [27]. Genes with high specific expression, logFC $>$0.25, and expression in at least 20% of cells in each cell cluster were selected as significant marker genes of cell subpopulations. Subsequently, the top two marker genes in each cell cluster were drawn into a heatmap to display the marker gene expression in different cell clusters more intuitively (Fig. 4A). Fig. 4A depicts that each cell cluster was separated according to marker gene expression, indicating that the selected marker genes could accurately represent the cell clusters. The marker genes for each cell cluster are listed in Table 1. Additionally, t-SNE plots and violin plots of the marker genes of clusters 4, 6, 11, 12, and 17, with significant differences between sham and MCAO groups in the 21 cell clusters, were drawn to demonstrate the expression differences of the marker genes in the different cell clusters (Fig. 4B–K). In conclusion, 42 marker genes were identified in the rat cortex, elucidating gene expression heterogeneity in each cell subpopulation.

Fig. 4.

Identification of marker genes. (A) Heatmap plots of the top two marker genes in each cell cluster. (B–F) The t-SNE plot depicted selected marker gene expression in clusters 4, 6, 11, 12, and 17. (G–K) Violin plots of selected marker gene expression of clusters 4, 6, 11, 12, and 17 in various cell clusters.

The number of cells in each cell cluster was analyzed in the cortex of rats, and it was found that the number of cells in clusters 0–3 increased significantly after MCAO (Fig. 5A). Moreover, the proportion of cells in clusters 0–3 increased dramatically after MCAO (Fig. 5B). The heatmap of the proportion of 21 cell clusters revealed that clusters 0–3 were significantly different in the sham and MCAO groups. Cluster 20, an ependymal cell subset, was not observed in MCAO group (Fig. 5C). These results indicated that microglia, astrocytes, and oligodendrocytes, especially microglia, proliferated in the rat cerebral cortex after MCAO. As the most abundant resident immune effector cells in the central nervous system (CNS), significant microglial proliferation demonstrates the importance of microglia-based immunomodulatory interventions in preventing and treating ischemic stroke. A comparison of the number of genes in each cell cluster between the sham and MCAO groups revealed a significant difference in nine cell subsets between the two groups (Fig. 5D), including clusters 1 (microglia), 2 (astrocytes), 3 (oligodendrocytes), 4 (endothelial cells), 5 (macrophages), 6 (endothelial cells), 7 (ependymal cells), 10 (oligodendrocytes), and 19 (fibroblasts).

Fig. 5.

Analysis of differences between sham and MCAO groups. (A) The number of cells per cluster in the sham and MCAO groups. (B) Pie charts depicting the cell proportion of each cell cluster in the sham and MCAO groups. (C) Heatmap presenting the cell proportion of each cell cluster in the sham and MCAO groups. (D) Gene numbers of cell clusters in rat cerebral cortex.

We selected the top 40 DEGs with the smallest p-value in the cell subpopulations for GO analysis to elucidate the biological processes involved in the six cell subpopulations that changed significantly after MCAO. Our analysis revealed that the biological enrichment process of the top 40 DEGs in cluster 1 was primarily related to the stress response after hypoxia (Fig. 6A), while cluster 2 may maintain protein homeostasis after MCAO, mainly through autophagy (Fig. 6B). Cell clusters 3 and 6 displayed similar biological processes, such as small GTPase-mediated signal transduction (Fig. 6C,D). Cluster 7 may be related to synapse formation (Fig. 6E), whereas cluster 10 mainly involved cell proliferation (Fig. 6F). In conclusion, GO analysis demonstrated that the biological processes involved in cell subpopulations were significantly altered in the rat cerebral cortex following MCAO.

Fig. 6.

Go analysis of DEGs in cell subpopulations. We performed GO analysis of DEG in cell clusters 1, 2, 3, 6, 7, and 10 to further clarify the biological function of significantly changed cell subsets after MCAO. (A) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 1. (B) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 2. (C) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 3. (D) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 6. (E) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 7. (F) Bubble plot of GO analysis of the top 40 DEGs in cell cluster 10. The bubble size represents the number of DEGs enriched in the biological process, and the color represents the p-value. DEGs, differentially expressed genes; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.

In addition, to GO analysis, we performed KEGG analysis on the top 40 DEGs in the cell subpopulations that were significantly altered after MCAO to further indicate the involved signaling pathways. Fig. 7 illustrates that KEGG analysis revealed the signaling pathways involved in Clusters 1, 2, 3, 6, 7, and 10 after MCAO. Notably, GO analysis revealed that the stress response after hypoxia in cluster 1 may be related to the HIF-1 signaling pathway, while autophagy in cluster 2 may be related to the mTOR signaling pathway. In conclusion, KEGG analysis further clarified the signaling pathways involved in DEGs in the significantly altered cell subpopulations after MCAO.

Fig. 7.

KEGG analysis of DEGs in cell subpopulations. KEGG analysis was performed on DEGs in cell clusters 1, 2, 3, 6, 7, and 10. (A) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 1. (B) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 2. (C) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 3. (D) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 6. (E) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 7. (F) Bubble plot of KEGG analysis of the top 40 DEGs in cell cluster 10. The bubble size represents the number of DEGs enriched in the biological process, and the color represents the p-value.