All experimental procedures were performed in accordance with the National Institutes of Health Guide for Care and Use of Laboratory Animals and were approved by the University of New Mexico Health Sciences Center Institutional Animal Care and Use Committee (IACUC). All studies were conducted in in healthy 7-8-week-old female mice, anesthetized with 1.5 mg/g (i.p.) urethane (Sigma-Aldrich) full dose injection followed by 0.75 mg/g (i.p.) half dose after 20 min. Animals were then placed in a stereotaxic frame and remained under anesthesia for the entirety of the experiment; body temperature was maintained 37 ± 0.5°C. A skin flap was created to expose the skull, and mineral oil placed to facilitate 2-D intrinsic optical imaging (IOS). The exposed cortex was illuminated with white light collected through a CCD camera; the changes in absorption and reflection of white light due to SD propagation allow SDs detection with non-fluorescence imaging (Hillman, 2007; Pinkowski et al., 2023). We used a Nikon SLR camera NIKKOR 50 mm f/1.8–16 lens (f-stop for imaging = 5.6) and a 12 mm extension tube was attached to a Mightex USB CCD (CXE-B013-U). IOS signals were then monitored in real time and recorded for 2.5 min. Light intensities were converted to a 16-bit grayscale image by the CCD camera, and SDs propagation images were generated after normalization and ΔF/F₀ analysis (Image J 1.53e)¹ to ensure SD generation (Supplementary Figure 1 and Supplementary Videos 1, 2). In some animals, a burr hole was prepared to allow for subsequent KCl application to induce SD (see below).

Spreading depolarizations were then induced repetitively (4 SDs in 2 h at 30 min intervals) SDs were induced in 6 mice with either focal application of KCl (1M) through a burr hole in C57Bl6J mice (n = 3) or with optogenetic stimulation (3 mW for 20 s) in Thy1-ChR2-YFP positive mice (n = 3). Thy1-ChR2-YFP breeders were initially purchased from The Jackson Laboratory (B6.Cg-Tg(Thy1-COP4/EYFP)18Gfng/J, Strain #:007612, RRID:IMSR_JAX:007612) and crossed to maintain the line, as previously described (Arenkiel et al., 2007). With either stimulus, evoked SDs propagated throughout the neocortex ipsilateral to the stimulus, without crossing to the contralateral hemisphere (Supplementary Figure 1 and Supplementary Videos 1, 2). A matched number of sham controls was used with focal NaCl administration through a burr hole, in C57Bl/6 mice (n = 3), or 0.4 mW for 20 s in Thy1-ChR2-YFP positive mice (n = 3). Thirty min after the last SD, the anesthetized animal was decapitated and the brain extracted, washed with ice cold PBS and kept on ice.

The contralateral and ipsilateral cortices were rapidly dissected from underlying striatum and the hippocampus and total RNA was extracted using TRIzol (Invitrogen, Thermo Fisher, Waltham, MA, USA). Extractions were either from intact ipsilateral or contralateral cortices (Figures 2–6) or from regions ∼18 mm² rapidly micro-dissected on ice (Figure 6). Where indicated, microdissections were made from 3 different regions of the ipsilateral hemisphere: (1) SD initiation site; (2) intermediate site >3 mm from initiation, and (3) remote site >5 mm from initiation (Supplementary Figure 3). All collected tissues were then flash frozen and kept at –80°C until RNA extraction. RNA was quantified using Qubit (Bio-Rad) and its quality was determined using a NanoDrop 1000 (Thermo Fisher), using absorbance at 260, 280 and 230 nm. Aliquots of 2 μg RNA were sent to Arraystar, Inc. for Illumina paired-end RNAseq as described below.

RNAseq analysis was performed by Arraystar Inc. (Rockville, MD, USA). Total RNA (2 μg) was used to prepare the sequencing library in the following steps: total RNA was enriched by oligo (dT) magnetic beads (rRNA removed); RNAseq library was prepared using KAPA Stranded RNASeq Library Prep Kit (Illumina), which incorporates dUTP into the second cDNA strand and renders the RNAseq library strand-specific. The completed libraries were qualified with an Agilent 2100 Bioanalyzer and quantified by using the absolute quantification qPCR method. To sequence the libraries on an Illumina NovaSeq 6000 instrument, the barcoded libraries were mixed, denatured to single stranded DNA in NaOH, captured on Illumina flow cell, amplified in situ, and subsequently sequenced for 150 cycles for both ends on an Illumina NovaSeq 6000 instrument. Sequence quality was examined using the FastQC software. The trimmed reads (trimmed 5′, 3′-adaptor bases using cutadapt) were aligned to a reference genome (Mouse genome GRCm38) using Hisat2 software.

The expression levels of known genes and transcripts were calculated using the StringTie and Ballgown software packages (Pertea et al., 2016) and expressed as fragments per kilobase of transcript per million mapped reads (FPKM). The number of identified genes and transcripts per group was calculated based on the mean of FPKM in group >0.5. Prior to performing DESeq2 analysis, genes were filtered to eliminate genes with less than 0.5. counts per million. Differentially expressed gene (DEG) and transcript analyses were performed using unpaired t-tests followed by Welch’s correction; a fold change (FC) threshold of 1.25 and p-value < 0.05 were used together to classify genes and transcripts as differentially expressed. We subsequently performed DESeq2 analysis (Ge et al., 2018) with a false discovery rate (FDR) adjusted p-value cutoff of < 0.01 and FC > 1.5. We then performed a post hoc power analysis: we estimated that with a n of 6 and alpha error of 0.05, the effect size is equivalent to a Cohen’s d of 5.47 (or a rho of 0.939) (Rothstein et al., 1992).

To validate gene expression from our RNAseq data set, separate aliquots containing 1 μg of total RNA were reverse transcribed using SuperScript II reverse transcriptase (Life Technologies) following the manufacturer’s protocol and qPCR carried out with a CFX96 Touch Real-Time PCR Detection System using SYBR Green mix (Life Technologies). Primer sequences for qRT-PCR for mouse mRNAs were obtained from PrimerBank.² The GAPDH gene (FW TGTGATGGGTGTGAACCACGAGAA, RV GAGCCCTTCCACAATGCCAAAGTT) was used to normalize values. As shown in Supplementary Figure 3, GAPDH expression itself was unchanged by SD with either KCl or optogenetic induction of SD. The relative expression of target genes was determined using the comparative 2-ΔCt method (Livak and Schmittgen, 2001).

Values were expressed as means ± SEM. Statistical analysis was performed using GraphPad Prism version 9 (La Jolla, CA, USA). The data were analyzed by unpaired t-test followed by Welch’s correction or analysis of variance (ANOVA) followed by Dunnett’s Multiple Comparison tests. Differences were considered statistically significant when p-values were < 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001). We performed a post hoc power analysis: we estimated that with a n of 12, and alpha error of 0.05, the effect size is equivalent to Cohen’s d of 2.37, or a rho value of 0.99 (Rothstein et al., 1992).

Predicted molecular functions and biological pathways enriched in genes differentially expressed after SD were identified using Ingenuity Pathway software [IPA, Content version: 62089861 (Release Date: 2021-02-17), Qiagen, (Krämer et al., 2014)].³ The analysis was set to filter nervous system expressed genes only; reference set: Ingenuity Knowledge Base (Genes Only); and including direct and indirect relationships.

Healthy female mice were subjected to multiple SDs (4 over a period of 2 h) which were generated to match clusters of SDs in clinical settings (Dreier et al., 2017). This procedure has previously been shown to induce robust transcriptional changes (Urbach et al., 2006; Sintas et al., 2017; Takizawa et al., 2020). Thirty minutes after the last SD or sham stimulus, the neocortex was removed, and RNA extracted from tissue subjected to SD or shams for subsequent RNA sequencing (see section “2 Materials and methods and Supplementary Figure 1).

Both focal KCl and optogenetic stimulation are utilized in studies of SD in brain injury, and thus to increase generalizability of this initial study of SD regulation of target genes, data from these two SD-stimulation methods were combined for differential expression (DE) analysis. For the combined DE analysis, we performed unpaired t-tests followed by Welch’s correction on the mapped 12,944 genes. We identified 102 significantly differentially expressed genes (DEGs) (Figure 1A, black, FC > 1.25, p-value < 0.05, n = 6, Supplementary Table 1). We subsequently performed DESeq2 analysis (Ge et al., 2018) and identified 16 genes that were also significant after false discovery rate correction (FDR adjusted p-value < 0.01, FC > 1.5) when comparing the SD group with the sham tissues (Figure 1A, red, FDR < 0.01, FC > 1.5, n = 6, Supplementary Table 1).

Consistent with previous studies (Kokaia et al., 1993; Dietrich et al., 2000; Rangel et al., 2001; Kaido et al., 2012), we found significant changes in immediate early genes (IEG) ARC, FOS, and JUN proto-oncogenes, and early growth response (EGR) genes. Other DEGs with significantly increased levels included adrenoceptor beta 1 (ADRB1), dual specificity phosphatase (DUSP) gene family, the nuclear receptor subfamily 4 group A member 1 (NR4A1), and prostaglandin-endoperoxide synthase 2 (PTGS2).

We next used Ingenuity Pathway Analysis (IPA) to identify molecular networks and biological pathways regulated by DEGs. This information was then used to obtain a prediction of physiological and pathological brain processes altered by SD. Enriched pathway analysis on FDR significant genes showed that genes upregulated by SD are involved in nervous system development and function, tissue morphology and development. We also found significant enrichment in genes involved in cell death and survival and cell morphology (Supplementary Table 2).

To achieve a more extensive enrichment analysis we then expanded IPA pathway analysis to genes with p-value < 0.05 and FC > 1.25. Figures 1B, C show the top two scored networks (score 63 and 35) generated by IPA with the combined DEGs upregulated in the SD group as compared to sham animals. The intensity of the color red for each molecule is proportional to the fold change (FC) of the labeled gene. Genes with FDR < 0.01 are also highlighted with thicker red outlines. Full arrows indicate a direct relationship between molecules, while dotted arrows indicate a predicted indirect relationship. Orange arrows and molecules indicate the predicted activation of a gene, while yellow indicates predicted inhibition. These networks implicated SD target genes associated with cell proliferation and synaptic activity, including cAMP responsive element binding protein 1 (CREB1), mitogen-activated protein kinase 1 and 3 (MAPK1 and MAPK3), and HOMER1 (Figure 1C).

Predictions of SD-induced changes associated with known disease states and functions are shown in Table 1 and they include: neurological diseases, psychological disorders, nervous system development and function, tissue morphology, and tissue development. The complete list of diseases and functions and the genes associated with each pathway are reported in Supplementary Table 2. Figure 1D also shows the number of genes associated with each pathway and associated p-value range. They are shown together so that different strengths of the network elements can be compared. A complete list of networks and associated DEGs is reported in Supplementary Table 3.

Figure 2 illustrates pathways related to cellular functions induced by repetitive spreading depolarization. The top common pathways identified include: cellular assembly and organization, cellular function and maintenance, and cellular development and proliferation (Supplementary Table 4). The heatmaps in Figures 2A, B show the p-value and the activation z score for each pathway, and a graphical representation of the cellular functions and genes involved, respectively. As shown in Figure 2, each gene expression FC is proportional to their color intensity; orange arrows and pathways indicate a predicted activation, whereas blue indicates increased levels of genes which lead to inhibition of the pathway; genes with FDR < 0.01 are also highlighted with thicker red outlines. IPA analysis predicted an activation of axonogenesis, axon branching, dendritic growth and branching and growth of neurites, based on increased levels of expression of genes encoding VGF, FOS, PTGS2, JUN, BDNF, gamma-aminobutyric acid type A receptor subunit beta3 (GABRB3), activating transcription factor 3 (ATF3), and cyclin dependent kinase like 3 (CDKL3). We also found predicted activation of pathways regulating neuronal development, long-term potentiation, and neurotransmission due to increased expression of FOS, PTGS2, EGR2, BDNF, JUN, KCNJ2, ARC, and HOMER1.

Interestingly, we also found predicted inhibition of neuronal cell death (z score −0.903), apoptosis (z score −1.197) (Figure 2A), together with increased development of neurons (z score + 2.224) survival of cerebral cortex cells (z score + 1.774), and cell viability (z score + 1.413) (Figure 2B), due to increased expression of genes such as FOS, BDNF, ATF3, JUNB, DUSP1, neuronal PAS domain protein 4 (NPAS4), nuclear factor interleukin 3 regulated (NFIL3), and SERPINE1.

Many studies utilize the contralateral neocortex as a convenient comparison tissue for DEG studies (Shen and Gundlach, 1998; Volobueva et al., 2023), as SD does not directly propagate across the corpus callosum. However, even if SD does not directly propagate through the contralateral cortex, the profound and widespread disruption of neuronal network activity produced by SD might be expected to have some impact on neuronal activity and gene expression in the contralateral hemisphere. The assumption that contralateral tissue is equivalent to tissue from sham animals for DE analysis has not previously been tested.

We first performed DE analysis between these 2 tissues (i.e., contralateral hemispheres in animals exposed to SD in the ipsilateral hemisphere vs. sham tissues). As shown in Supplementary Figure 2 only 9 genes were significantly differently expressed after unpaired t-test followed by Welch’s correction, but none of were still significant after FDR correction, confirming a high degree of similarity between gene expression in these tissues.

We next performed the same analyses described above (Figure 1), but comparing tissue exposed to repetitive SDs with the contralateral cortex in the same animals, rather than sham animals. Of the 12,893 mapped genes 149 were significant after unpaired t-test followed by Welch’s correction (Figure 3A, black, FC > 1.25, p-value < 0.05, n = 6, Supplementary Table 1), of these 21 were also significant after FDR correction (FDR adjusted p-value < 0.01, FC > 1.5, n = 6) (Figure 3A, red). IPA analysis confirmed a high similarity between the sham and contralateral cortex when used as control compared to the SD group. Figures 3B, C show the top scored networks generated from the SD vs. contralateral DEGs list, genes with FDR < 0.01 are also highlighted with thicker red outlines. The top diseases and function are reported in Table 2 and Supplementary Table 2 (compare with Figure 1). Likewise, Figure 4 shows cellular and molecular pathways that are comparable with DEG analysis with sham tissues (compared with Figure 2). These results imply that, at least at the level of gene expression changes after 2 h, the contralateral hemisphere that is not directly invaded by SD can likely serve as an unaffected control.

To validate the predicted activation of neuroprotective pathways, we analyzed the expression levels of a subset of SD target genes by qRT-PCR. The genes selected for validation where based on the IPA analysis; indeed, enrichment pathways analysis was used to identify genes functionally related to SD target genes that were not found significant in DE analysis. Based on the IPA results above, we focused on genes regulating 3 cellular functions: (1) cell differentiation and proliferation, (2) synaptic plasticity, and (3) stress responses. qPCR was performed on the same RNA that was subjected to the RNAseq analyses, as well as 3 additional animals for each group. As shown in Figure 5, multiple SDs resulted in significant increases in the expression of genes activating cell proliferation and differentiation such as BNDF, DUPS6, FOS, and JUNB (Figure 5). Considering the predicted activation of neurons and long-term potentiation, as well the increased axonogenesis, growth of neurites, and axon and neurites branching, we then analyzed the expression of genes known for their regulation of synaptic plasticity. As shown in Figure 5 we found significantly increased expression of Homer1a and ARC. Finally, we examined levels of PTGS2, NR4A1, and EGR2 as genes regulating stress and inflammatory responses, showing a significant upregulation.

Lastly, we tested whether the levels of gene expression varied with distance from the SD initiation site. In this set of studies, KCl was used as the stimulus, as the site of initiation is clearly restricted to small zone of the application of KCl through the burr hole. Thus, 2 h after onset of the initial SD, RNA was extracted from either total contralateral cortex, or from three different regions of the ipsilateral hemisphere: (1) SD initiation site; (2) intermediate site, and (3) remote site (see section “2 Materials and methods” and Supplementary Figure 3). By qRT-PCR we analyzed the levels of the previously identified target genes (see section “3 Results” above) and compared these to the levels in the contralateral hemisphere of SHAM animals. For the majority of the genes examined (BDNF, DUSP6, PTGS2 ARC, and Homer1a), highest expression was detected in the intermediate area (>3 mm from the initiation site; Figure 6) while FOS and JUNB genes showed significantly increased expression in the remote site (>6 mm from the initiation site). Interestingly, when analyzing Homer1a levels, we found a decrease in their expression at the SD initiation site (Figure 6) when compared to the contralateral and intermediate sites (Figure 6). The difference between the contralateral and remote sites for NR4A1 did not reach statistical significance (p-value = 0.09) and EGR2 did not show differential expression across regions (Supplementary Figure 3).