In the present study, we systematically mapped brain region activation following 24 h social isolation by labeling c-Fos protein in the SDMN as well as in other regions involved in social stress processing, such as Hb-IPN axis (Sugama et al., 2002; Benekareddy et al., 2018; Xu et al., 2018; Ogawa and Parhar, 2022). Previous studies have shown that even a single day of social isolation significantly elevates c-Fos expression in the Hb system of rodents (van Kerkhof et al., 2013; Stanisavljević et al., 2020). This increased activity is likely accompanied by distinct patterns of brain activation in several regions of the SDMN, which may serve as a basis for differentiating VPA-treated mice from control mice. Beyond its role in acute social stress processing, growing evidence links the habenula to the pathophysiology of autism spectrum disorder. In ASD models, altered habenular excitability and disrupted glutamatergic signaling within the lateral habenula have been shown to contribute to social deficits and motivational impairments (Murru et al., 2021; Xu et al., 2023). Moreover, structural alterations of the habenula have been reported in individuals with ASD, further underscoring its relevance in social behavior regulation and autism-related neurocircuitry (Germann et al., 2021). Since idiopathic autism and embryonic VPA exposure affect multiple neural pathways rather than an isolated brain region, we applied a network-based approach to examine the broader functional organization of affected circuits (van Kerkhof et al., 2013; Nagode et al., 2017; Ádám et al., 2020; Dorsey et al., 2024). To achieve this, we constructed a functional connectome based on correlational activation patterns across brain regions (Wheeler et al., 2013; Stefaniuk et al., 2023; Kotlarz et al., 2024). This method enabled us to explore alterations in the broader neural network structure, rather than focusing solely on localized changes.

To support region-specific interpretation of neuronal activation, four functionally distinct brain regions were included as controls (with presumably no direct influence on social behavior): the dorsomedial hypothalamus (DMH), the arcuate nucleus (ARC), the Edinger–Westphal nucleus (EW), and the dentate gyrus (GD). The DMH was selected based on its broad integrative role in autonomic and behavioral regulation within the hypothalamus (DiMicco et al., 2002; Kataoka et al., 2013). The ARC serves as a key neuroendocrine and metabolic center, providing a reference for activation patterns in nuclei involved in homeostatic regulation (Song and Choi, 2023). The EW, particularly its centrally projecting subpopulation, represents a functionally divergent midbrain structure, allowing comparison to regions with different neurochemical and anatomical profiles (Kozicz et al., 2011). Finally, GD, the primary gateway of cortical input to the hippocampus, was included to represent a general input structure with a well-defined role in hippocampal information processing (Hainmueller and Bartos, 2020).

We hypothesize that VPA-treated mice, modeling the autistic phenotype, exhibit diminished stress-related neural activation and altered SDMN network-wide connectivity in response to social separation compared to typically developing mice. By analyzing these circuits, we aim to deepen our understanding of the mechanisms underlying ASD and identify potential targets for therapeutic intervention.

C57BL/6 mice were maintained on a 12:12 light–dark cycle, with the dark phase beginning at 8:00 AM. Food and water were provided ad libitum. All experiments were conducted at least 2 h after lights-off to ensure alignment with the animals’ active phase. Female mice were housed in pairs, with timed mating achieved by introducing a male for 24 h. Pregnancy was confirmed by monitoring weight gain on day 13.5 of gestation (E13.5) when the VPA and control treatment took place (Kim et al., 2011; Mehta et al., 2011; Kang and Kim, 2015; Zarate-Lopez et al., 2024). Pregnant mice were injected subcutaneously with either Na-VPA (Depakine® 400 mg/4 mL 500 mg/bwkg) (VPA-treated) or volume-matched saline for control groups (CTR). Male pups from 14 different litters were weaned on postnatal day 24 and used for the experiment on postnatal day 31 (n_CTR = 10, n_VPA = 11). The animal study was approved by the Food Chain Safety and Animal Health Directorate of the Government Office for Pest County, Hungary (PE/EA/926–7/2020). The study was conducted in accordance with the local legislation and institutional requirements.

Group level behavioral differences caused by VPA treatment were validated using the three-chamber sociability test on the day of weaning (day 24). The apparatus consisted of three connected transparent Plexiglass boxes (19 × 45 cm) (Kaidanovich-Beilin et al., 2011). In the initial session, mice were habituated to the environment for 5 min. During the social preference session, an unfamiliar control mouse was placed in one terminal chamber, while the other remained empty. Each stimulus mice were used in two tests (one for a VPA and one for a control animal) in randomized order. Such stimulus mice were later incorporated into the social housing groups as companion animals, however, they were never used as experimental individuals. The test lasted 10 min, with data from the first 5 min analyzed, as exploratory behavior and social interactions primarily occurred during this initial period. After this time, animals generally became less active, and their chamber exploration plateaued. The entire session was video recorded and analyzed using Solomon Coder (2024) software. Sociability was evaluated by measuring both the total time spent in each chamber and the average visit duration. The average visit duration was defined as the mean duration of individual visits to a chamber, calculated as the time elapsed between an entry and an exit. The total time spent in a chamber represented the cumulative duration of all visits within the session. The data were statistically compared using the Wilcoxon signed-rank test.

Following the three-chamber test, animals were cross-fostered into mixed groups of five, including both VPA-treated and control (CTR) mice. The animals in each group were from different litters, and each group consisted of a mix of VPA-treated (2 or 3) and CTR individuals. On postnatal day 31, one mouse from each group was separated for 24 h. While 24 h of social deprivation might appear too long to induce specific c-Fos activation, van Kerkhof et al. (2013) showed that 24 h of isolation caused a higher increase in c-Fos expression in some brain areas than 3.5 h of isolation. The separated mice were placed in a new home cage with fresh bedding, supplemented with a small portion of bedding from their previous cage. The cages of separated mice were kept in the same room but positioned far from their previous cage mates. On postnatal day 32, after 24 h of isolation and at least 2 h into the dark phase, the animals’ active period, mice were anesthetized using isoflurane (Forane®) and transcardially perfused with physiological saline, followed by 4% paraformaldehyde (PFA).

Brains were stored in 4% PFA and then transferred to 25% sucrose in PBS 1 day before sectioning. After sucrose infiltration, brains were frozen and sectioned at 50 μm using a freezing microtome (Leica SM 2000R). For c-Fos immunohistochemistry. Free-floating sections were rinsed in PBS, and endogenous peroxidase activity was blocked for 20 min at room temperature (RT) using PBS containing 2% H₂O₂. Blocking was performed for 1 h at RT in PBS containing 1% normal goat serum and 0.2% Tween-20, followed by incubation with a polyclonal rabbit anti-c-Fos antibody (1:5000, Abcam, ab190289) in PBS with 1% normal goat serum and 0.02% Tween-20 for 1 h at RT and then overnight at 4°C. On the following day, the slices were incubated for 2 h at RT in a biotinylated anti-rabbit IgG secondary antibody (1:200, Biotinylated IgG, Vector Laboratories, BA-1000), followed by overnight incubation at 4°C in avidin-biotin solution (1:500, ELITE ABC kit, Vector Laboratories) containing 0.01% Tween-20. After rinsing in PBS and Tris buffer (pH 8.0), sections were incubated in NiDAB solution (0.15 mg/mL DAB + 4.25 mg/mL Ammonium nickel sulfate Hexahydrate) for 10 min at 4°C, then developed for 5 min at RT with the addition of 1 μL dense H₂O₂. Following the staining, the sections were mounted onto chrome-gelatin-coated slides, dehydrated through an ascending ethanol series [50, 70, 90, 96, 100% (v/v)], cleared in xylene, and coverslipped using DePeX mounting medium.

Images were captured digitally using a Nikon Eclipse E800 microscope equipped with an Olympus DP74 camera. Brain regions were identified based on the SDMN framework (O’Connell and Hofmann, 2011), with the following modifications: The Decision-Making Network was expanded to include ACC, and some of the SDMN regions were further subdivided according to functional differences. Regions associated with social stress and autism were also analyzed. Regions were imaged at the same rostrocaudal coordinates, determined by anatomical landmarks, using 4 × magnification to ensure that the entire region of interest (ROI) was fully visible within a single image. Each ROI was imaged bilaterally, with at least one image taken from each hemisphere per subject, allowing for the analysis of laterality effects. For larger nuclei spanning multiple sections, such as the BNST, multiple images were taken from different slides within the same rostrocaudal coordinate range, ensuring comprehensive anatomical representation (Atlas Thumbnails, 2024). For cell quantification, the entire ROI was outlined based on anatomical references from the mouse brain atlas, and c-Fos + cells were counted within the fully delineated area to maintain consistency across subjects.

c-Fos positive cells were quantified using the particle analysis function of ImageJ (2024), following background subtraction and Gaussian blur to enhance detection accuracy. A size filter (70–200 pixels) and a circularity threshold (0.85–1.00) were applied, and signals meeting these criteria were considered c-Fos-positive cell nuclei. These validated counts were then used for subsequent statistical analyses. Data were analyzed using the statistical software R. For each brain region, c-Fos-positive cell counts were compared between the experimental groups using a negative binomial mixed-effects model (model parameters: response variable: c-Fos positive cell number; main effect: Embryonic treatment + nested random effects litter, individual) + offset [log(Area of ROI)] (Yu et al., 2022). False discovery rate (FDR) correction was applied to p-values using the Holm-Bonferroni method, setting the significance threshold at 0.05. Within subject pairwise comparisons failed to detect any significant inter-hemispheric differences, therefore laterality as a factor was not included in the final model. In cases where the p-value ranged between 0.05 and 0.1, Cohen’s D values were calculated using the “effsize” R-library (Torchiano, 2013). Comparisons with an effect size (Cohen’s D) <0.8 were visualized separately to highlight near significant differences (Lakens, 2013) in order to screen for potential trends and identify regions with suggestive group differences that may warrant further consideration.

The c-Fos density for each counted nucleus was calculated by summing the cell counts from both hemispheres in bilateral regions, along with their respective area measurements, then dividing the cell number by the area (Stark et al., 2006; Stefaniuk et al., 2023). For the nuclei spanning multiple slides, the mean density was used.

Within each of the two experimental groups (VPA, CTRL), all possible pairwise correlations between the Fos signal in the 36 brain regions were determined by computing Pearson’s correlation coefficients. The resulting p-values and R² values were compiled into matrices. The R² matrix was visualized as heatmaps (Figures 1, 2) using the “heatmaply” R-library (Galili et al., 2018).

For functional connectome network analysis, we used a methodology based on the work of Wheeler et al. (2013). To reconstruct the functional connectome, we included only R² values that overlapped with connectivity data from the Allen Brain Atlas, ensuring that the heatmaps reflected physiologically plausible connections (Oh et al., 2014). This filtering step was crucial to eliminate indirect correlations, where two nuclei might exhibit similar activity not because of a direct anatomical connection, but due to a common upstream hub that synchronizes their activity. Our approach specifically aimed to identify these hub regions, as they play a key role in mediating large-scale network coordination. These values were thresholded on the basis of two criteria: an upper limit of p < 0.05 and a lower limit of R² > 0.7. These correlations were considered as functional connections between brain regions and were visualized as edges in subsequent network visualizations using the igraph R-package (Csárdi et al., 2025). Black edges indicated positive correlations, while red edges indicated negative correlations (Figure 3).

For functional network analysis, degree and betweenness centrality values were calculated using the igraph package. Degree centrality quantifies how many direct connections (edges) a node has. A region with a high degree centrality is strongly connected to many other regions but does not necessarily act as a key relay in communication. Betweenness centrality, on the other hand, measures how often a node lies on the shortest path between other nodes. Betweenness centrality was calculated as B(v) = ∑[σ_uw(v)/σ_uw], where σ_uw is the number of shortest paths between nodes u and w, and σ_uw(v) is the number of shortest paths between u and w that pass through node v. The sum in the expression ranges over all pairs of distinct nodes u and w (Nguyen, 2014). A high betweenness score means that a node serves as a bridge, influencing how efficiently information flows through the network (Zhang and Luo, 2017; Kotlarz et al., 2024). For centrality analysis, a limit was defined at the 80% quantile of the sum of both centrality values. This criterion ensures that a region is not only highly connected (high degree) but also strategically positioned within the network, allowing it to act as a relay for indirect interactions. Nuclei ranking above the 80% quantile in both betweenness and degree centrality were classified as hubs.

In the control (CTR) group, animals spent significantly more time in the chamber containing an unfamiliar conspecific compared to the empty chamber [Z₍₁₀₎ = −1.988, p = 0.047]. In contrast, there was no significant difference in chamber preference (Figure 4A) for the VPA-treated group [Z₍₁₁₎ = −0.978, p = 0.328]. Regarding the average duration of visits (durations between the entry and exit into a chamber) in each chamber, CTR animals showed a more pronounced difference. The visits (Figure 4B) to the conspecific lasted significantly longer than those to the empty chamber [Z₍₁₀₎ = −2.666, p = 0.008]. However, no such difference was observed in VPA-treated group [Z₍₁₁₎ = −0.889, p = 0.374].

Overall, the CTR group exhibited greater neural activity in response to social isolation compared to the VPA-treated group in the observed brain regions (Figure 5 and Supplementary Table S1). In all regions where group differences were detected, activity was consistently higher in the CTR group (Figure 5), and this trend was also apparent across regions qualitatively. Using a multivariate (multiple brain regions) mixed model on the separate networks c-Fos activity appeared significantly higher in the control animals in the MRS [χ²_{(1, 279)} = 4.463, p = 0.035]. A subsignificant trend was visible also in the SBN [χ²_{(1, 191)} = 3.719, p = 0.054]. A similar trend cannot be excluded even in the chosen control regions [χ²_{(1, 90)} = 2.83, p = 0.09], but there was no significant difference in the stress-regulating network regions [χ²_{(1, 171)} = 0.837, p = 0.360] (Table 1).

To ensure that the correlative activation between brain regions (Figure 1), truly represents meaningful connections, we tested whether the region pairs with confirmed direct anatomical connections (based on the Allen Brain Atlas connectome [Oh et al., 2014)] exhibited higher correlations than those pairs without proven axonal connections (t = 2.601, df = 1952.1, p = 0.009) in the control mice. The correlations were higher between connected pairs than between unconnected pairs in the VPA animals, too (t = 19.157, df = 41.881, p < 0.001). The correlation matrices, based on known anatomical connections, and high (R² > 0.7) and significant (p < 0.05) correlations (Figure 2), revealed a generally stronger co-activation of brain areas in VPA-treated animals as compared with CTR animals. When the network was reconstructed based on the matrices shown in Figure 2, distinct organizational differences were observed between VPA-treated and CTR groups when the correlation matrices were visualized as graphs (Figure 3). The VPA network exhibited a greater number of strong correlations, indicating broader functional connectivity across multiple brain regions. This suggests that VPA treatment alters network organization, leading to increased co-activation among several regions.

In CTR animals, the network was characterized by a structured and modular organization, with distinct clusters of nuclei and fewer interconnections (Figure 3A). A central group of stress-related regions was evident, while peripheral nodes maintained more selective connections. In contrast, the network in VPA-treated animals displayed a denser and highly interconnected structure, with more widespread links between regions (Figure 3B). This shift suggests a reorganization of functional connectivity, potentially reflecting altered information processing and integration across brain areas following VPA exposure.

Quantitative network analysis revealed that, during social isolation, VPA-treated animals exhibited greater regional connectivity (degree centrality) compared to CTR animals (W = 376.5, p = 0.002) (Figure 6). Despite this increase in connectivity, the betweenness centrality, a measure of how crucial a node is in facilitating communication between different parts of the network (Wheeler et al., 2013; Zhang and Luo, 2017), did not differ significantly between the groups (W = 694.5, p = 0.572) (Figure 7). This suggests that, while VPA-treated animals had more connections overall, these connections did not increase the prominence of specific nodes as central hubs in the network.

Regions that displayed high ranks (above 80% of all ROIs) in both degree (Figure 6) and betweenness (Figure 7) centrality were considered as hubs within the observed brain network, as these regions were both highly connected and played a key role in mediating network-wide communication. In the CTR animals only IPN matched these criteria, ranking above the 80% threshold in both degree (Figure 6A) and betweenness centrality (Figure 7A). In contrast in VPA-treated mice, only the vGD exceeded the 80% threshold for both degree (Figure 6B) and betweenness centrality (Figure 7B). Overall, regions within the functional brain network of VPA individuals exhibited higher degree centrality (Figure 6) but similar—or even lower, considering the highest quantiles—betweenness centrality (Figure 7) compared to CTR mice.

While the method of reconstructing functional networks based on cFos activation does not replace direct circuit-tracing approaches, such as fMRI, optogenetics or electrophysiological recordings, recent research supports the use of c-Fos expression patterns to infer large-scale functional networks. Stefaniuk et al. (2023) demonstrated that c-Fos-based network analyses can accurately capture behaviorally relevant functional connectivity, while Kotlarz et al. (2024) further confirmed its usefulness in identifying experience-dependent neural network changes. Additionally, Stark et al. (2006) compared c-Fos activity patterns with fMRI-based functional connectivity, concluding that c-Fos mapping provides highly reliable spatial resolution of functionally active circuits. The possibility cannot be excluded that brain regions without direct anatomical connections could exhibit similar Fos activation patterns due to shared functional roles or parallel processing of social stress responses. By including only neuroanatomically verified connections, we aimed to reduce spurious correlations and highlight the differences in the network structure. This approach helps mitigate the risk of falsely interpreting activity synchronization between two nuclei as a direct functional link when, in reality, their activity is driven by a third common hub. Furthermore, our methodology was specifically intended to identify these hub regions, as they are the critical mediators of network-wide coordination.

Research suggests that females may be more resilient to prenatal VPA exposure due to compensatory mechanisms, including hormonal modulation and alternative social coping strategies, which make ASD-like traits less evident in early developmental stages (Lai et al., 2013; Hull et al., 2020). In juvenile rodents, females often display more adaptive social behaviors, potentially masking the full impact of prenatal VPA exposure (Hiller et al., 2014; Ratto et al., 2018). This aligns with human studies showing that autistic females are often diagnosed later than males, as their behavioral phenotypes diverge from the classical, male-based diagnostic framework (Lai et al., 2011; Belcher et al., 2023). Given these sex differences, our study specifically aimed to investigate the classical ASD-like phenotype, which is more reliably expressed in male subjects during juvenile stages. While the question of sex differences in ASD is a worthy topic of research, it was beyond the scope of the present study. Future work incorporating both sexes, particularly at later developmental stages, may provide further insights into sex-specific ASD manifestations.