Mouse brain-wide mitochondrial connectivity anchored in gene, brain and behavior

1 Division of Behavioral Medicine, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 2 Department of Psychiatry and Behavioral Sciences, Stanford University 3 The Sackler Institute, Department of Psychiatry, Columbia University Irving Medical Center, New York NY Division of Systems Neuroscience, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 5 Brain Mind Institute, EPFL, Switzerland 6 Division of Molecular Therapeutics, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 7 Department of Pediatrics, Columbia University Irving Medical Center, New York NY 8 Division of Developmental Neuroscience, Department of Psychiatry, Columbia University Irving Medical Center, New York NY 9 New York State Psychiatric Institute, New York NY 10 Department of Neurology, Columbia University Irving Medical Center, New York NY


Introduction
The shaping of behaviors by life experiences is driven by energetically demanding circuitry across the brain 1 . The brain's enormous energetic demand is mainly met by ATP produced through oxidative phosphorylation (OXPHOS) within mitochondria 2 . As a result, mitochondria influence multiple aspects of brain development and function ranging from dendritic and axonal branching 3,4 , remodeling of gene expression 5 , regulation of neurotransmitter release and excitability in mature synapses [6][7][8] , neurogenesis 9 , and inflammation 10 . Thus, the mounting molecular and functional evidence that the brain is under energetic constraints 11,12 suggests that a fundamental understanding of brain mitochondrial biology is essential for uncovering the basis of brain function and behavior.
Mitochondria are small, dynamic, multifunctional organelles with their own genome that generate ATP via the respiratory chain (RC) 13 . But not all mitochondria are created equal.
Mitochondria subserving different cellular demands (i.e., in different cell types) have different relative molecular compositions and functional phenotypes [14][15][16][17] . Therefore, developing a comprehensive understanding of the association between mitochondria and animal behavior calls for assessments of multiple mitochondrial features, across multiple brain regions. By analyzing mitochondrial features across multiple brain regions simultaneously, we can also potentially uncover unknown brain "mitochondrial circuits". This idea aligns with our evolving understanding of large-scale brain circuitry and metabolism 11,18 , and of the network distribution of neural activity achieved using brain-wide, high spatial resolution methods (e.g., MRI-based functional and structural connectivity maps) 19,20 . Although a similar degree of resolution for mitochondrial phenotyping is not feasible with current technologies, miniaturization of biochemical and molecular mitochondrial assays 21 open new possibilities to map mitochondriato-behavior across multiple cortical and sub-cortical regions in the same animal.
Both animal behaviors and mitochondrial health exhibit naturally occurring and acquired (i.e., experience-dependent) variation, providing an opportunity to map their associations. For example, behaviorally, exposure to social stress such as chronic social defeat predictably increases potential threat vigilance and social avoidance 22 , but animals differ in their vulnerability/resilience to these stress-induced changes in behavior 23,24 . Experimental challenges aimed at modeling neuroendocrine disturbances resulting from chronic stress, such as chronic exposure to corticosterone (CORT), also induce avoidance behaviors associated with recalibrations in specific brain circuitry 25 , gene expression 26, 27 , and anatomical plasticity (i.e., atrophy) in stress-sensitive brain regions like the hippocampus 28 . In relation to mitochondria, a separate body of literature similarly documents naturally-occurring mitochondrial variation 21,29 , as well as stress-induced functional recalibrations in mitochondria that occur within days to weeks of stress exposure [30][31][32] (meta-analysis in 33 ), potentially linking mitochondrial health to behavior. More direct causal experiments show that mitochondrial RC function directly influence the brain and specific behavioral domains including working memory 34,35 , social dominance 29,36,37 and anxiety-related behavior 4 . Targeted mitochondrial defects even cause mood disorder-like phenotypes in animals 38 and are likely implicated in the etiology of psychiatric, neurological and degenerative disorders in humans 39 . Moreover, in vivo human brain metabolic imaging studies show that energy metabolism in specific brain areas (e.g., nucleus accumbens, NAc) predict performance and anxiety [40][41][42] . Thus, current evidence positions mitochondria as an upstream modulator of brain function and behavior.
Although the importance of mitochondria for brain structure and function is unequivocal, we lack an understanding of potential differences in mitochondria across different brain regions, therefore calling for a more complete mapping of region-specific brain mitochondrial function in relation to behavior. To address the hypothesis that mitochondrial health in specific brain regions is associated with behaviors, we have further miniaturized existing biochemical and 5 molecular assays of mitochondrial RC function for sub-milligram tissue samples and deployed them across 17 cortical and sub-cortical brain regions in mice with a wide range of stressenhanced behavioral phenotypes, quantified through four behavioral tests. Using network-based connectivity analysis, we also sought to explore mitochondrial properties across brain regions, finding initial evidence that mouse brain mitochondria specialize as distinct functional networks linked to behavior. Further, we use gene co-expression and structural/anatomical connectivity data, which anchor the newly observed mitochondrial circuits into other modalities.

Miniaturization of mitochondrial assays for sub-milligram resolution
Although enzymatic activity assays directly reflect the capacity of RC complexes and therefore energy production capacity, currently available assays suffer from low throughput that preclude a multi-region, brain-wide analysis in dozens of animals. Building from our previous efforts to scale throughput of mitochondrial biochemical assays 21 , here we miniaturized and optimized assays for 96-well plate format to establish the lowest detection limit for mouse brain tissue. With this optimized platform, we can quantify the activities of RC complex I (CI, NADH dehydrogenase), complex II (CII, succinate dehydrogenase), and complex IV (CIV, cytochrome c oxidase) enzymatic activities, as well as the Krebs cycle enzyme citrate synthase (CS, a marker of mitochondrial content) in two 1mm-diameter, 200μm deep tissue punches. This represents <1mg of tissue (estimated 0.33mg), over an order of magnitude more sensitive than currently available methods.
In our pipeline, the same biological sample used for enzymatic activities is also used to quantify mtDNA abundance, both by i) the classical metric relative to the nuclear genome -mtDNA/nDNA ratio, termed mtDNA copy number (mtDNAcn) 43 , and ii) per tissue volume (mtDNA copies per mg), termed mtDNA density. From our qPCR data of nuclear genome abundance across the brain, independent of the mtDNA, we find that cellular density varies by 6 up to 8.5-fold between brain regions (highest: cerebellum, lowest: visual cortex). Unlike in peripheral tissues, this remarkably large variation in neuronal and glial cells somata density significantly skews mtDNAcn estimates. Figure S1 shows the inter-correlations of mtDNAcn and mtDNA density in relation to enzyme activities, establishing that mtDNA density is a more appropriate and generalizable estimate of mitochondrial genome abundance across mouse brain regions.
We also combine these five primary mitochondrial features into a mitochondrial health index (MHI) by dividing RC function (CI+CII+CIV) by mitochondrial content (CS+mtDNA density), thus creating an index of energy production capacity on a per-mitochondrion basis 21 (See Methods for details). Although each of the six resulting features are partially correlated with other features, the proportion of shared variance between individual features across brain regions is 31-64%, indicating that they each contribute some non-redundant information about mitochondrial phenotypes, which can be deployed at scale.
Our study design aimed to profile animal-to-animal differences in mitochondrial health across a broad set of brain regions that are associated with anxiety-related behavior, social behaviors, cognition, and mitochondrial disorders, and to relate these measures to each animal's behavioral phenotypes. In total, we enzymatically and molecularly phenotyped 571 samples covering 17 cortical, sub-cortical and brainstem brain regions isolated by bilateral punches at defined stereotaxic coordinates (Table 1, Figure S2B). To eventually compare the specificity of our findings related to brain and behavior, 5 peripheral tissues were also collected and analyzed from each animal, which we expected to show somewhat related but more modest associations with behavioral outcomes.

Protein levels do not reflect enzymatic activity
We initially explored if RC protein abundance is a viable surrogate for mitochondrial RC activity 29 , which could theoretically allow high spatial resolution imaging of the entire brain. We 7 focused on the cerebellum due to its well-defined cellular composition and laminar organization, where the Purkinje cell layer is flanked by molecular and granular layers 14 (Figure S3A).
Compared to protein levels, enzymatic activity ultimately determines mitochondrial RC function and energy production capacity and consequently should be regarded as the most representative measure of mitochondrial function. In consecutive cerebellar slices, we compared SDH (CII) enzymatic activity to the protein abundance of a RC complex II subunit, SDHA (succinate dehydrogenase, subunit A), for which a validated high-affinity antibody is available. Enzyme activity did not correlate with protein abundance for any of the three cerebellar layers (proportion of shared variance, r 2 =0.02-0.07), indicating that protein abundance is not an appropriate surrogate for enzymatic activity ( Figure S3B). Reasons for this finding could include the action of post-translational modifications, stoichiometry of the four SDH subunits, or the biochemical context that drive activity independent of protein content (e.g., 44 ). Therefore, we focus downstream analyses on direct measures of mitochondrial RC enzymatic activity and mtDNA density.

Diversity of mitochondrial RC activities between animals
We sought to examine mitochondria-behavior associations in a cohort of adult male mice exhibiting naturally-occurring behavioral and mitochondrial variation. Our goal was to identify robust and generalizable associations between mitochondrial function and behavior, rather than potential correlations that would exist among only subgroups of animals (naïve or stressed). Therefore, to further extend the spectrum of mito-behavior variation using well-characterized rodent stress models, subgroups of the cohort were either chronically administered corticosterone (CORT) 45 , or exposed to chronic social defeat stress (CSDS) 46 ( Figure S2A).
To create further diversity among groups, half of the CSDS mice were allowed to recover for two months. Based on previous work 23,47 , we also expected naturally-occurring behavioral and brain molecular phenotype differences between animals that are resilient or susceptible to CSDS 8 (based on the social interaction test; see Methods) 48 . The naturally occurring variation in mitochondria and behavior plus the effects of various exposures provides a strong test of robustness for our hypothesis, which is that across a diverse population of mice, mitochondrial health in specific brain regions are consistently associated with behaviors.
Across our cohort of inbred mice with a range of exposures, a wide range of mitochondrial health was observed. The average variation for mitochondrial activities (excluding MHI), across all animals and all 17 brain regions was 36% (coefficient of variation, C.V.) ( Figure   S4B). For any given brain region, the variation in mitochondrial health between mice reached up to 2.9-fold between the animals with the lowest and highest activities. This means that for a given brain region, there are large mouse-to-mouse differences in mitochondrial content and RC activity, even among inbred control mice not exposed to stressors ( Figure S4). Relative to brain regions, peripheral tissues showed about a third less variation between animals, with an average 25% C.V. (Figure S4D).
To verify that interventions used to enhance naturally-occurring variation in mitochondria were effective, we compared mitochondrial features in mice exposed to CORT and CSDS relative to naïve mice. The effects of CORT and CSDS on mitochondrial health were quantified as standardized effect sizes (Hedges g) and shown in (Figure 1A, detailed in Figure S5). Both interventions significantly altered mitochondrial features in a region-specific manner. CORTtreated mice had higher mitochondrial activities than non-stressed animals in ~60% of brain regions, whereas CSDS animals had lower activities in ~82% of brain regions (Figure 1D), suggesting opposing effects of these two different stress models on brain mitochondria. The amygdala showed the greatest CORT-induced increase in CII enzymatic activity (+49%, p<0.05), whereas the periaqueductal gray (PAG) showed the largest CSDS-induced decrease in CI activity (-42%, p<0.05). 9 In order to determine if brain regions were co-regulated in their stress-induced mitochondrial recalibrations, we employed a topological data analysis (TDA)-based Mapper approach. Mapper is a variant of nonlinear dimensionality reduction (manifold learning) methods that produces a graph or network embedding of the high-dimensional data (a.k.a. shape graph) while recovering projection loss using an additional partial clustering step 49,50 . When applied to mitochondrial features, separately for each group, the shape graphs revealed differences in regional mitochondrial recalibrations across the two groups ( Figure 1E-F). This was quantified using a graph-theoretical measure of participation coefficient (PC) 51 , where higher values of PC indicate uniformly distributed connectivity across regions (integrated network) and lower values indicate a segregated connectivity pattern ( Figure 1G). CORT-induced mitochondrial recalibrations were relatively more region-specific (or segregated) whereas CSDS caused a more coherent response across the brain (or integrated).
Although we were not powered to make statistical comparisons between susceptible and resilient mice, we did observe relatively large differences in effect sizes (g > 0.8) in mitochondrial recalibrations between the groups ( Figure S6). Additionally, compared to nonrecovered CSDS mice, animals who were allowed to recover from CSDS exhibited only marginal mitochondrial health changes in a limited number of brain regions, as visualized by the Mapper-generated shape graphs ( Figure S7). Overall, these results confirm the existence of naturally-occurring and acquired variation in brain mitochondrial health between animals, providing a strong basis to test the existence of conserved associations with behaviors.

Diversity of anxiety and depressive like behaviors
Just as there was a spectrum of mitochondrial health across mice, as expected from previous work, animals also naturally exhibited large variation in their behavioral phenotypes.
Behavioral tests included potential threat ("anxiety") monitored by the open field test (OFT) and elevated plus maze (EPM), hyponeophagia monitored by novelty suppressed feeding test 1 0 (NSF), and approach-avoidance conflict using the social interaction test (SI). Because specific behavioral tests are generally administered in conjunction with specific interventions (CORT, CSDS), some animals were only tested on some and not all behaviors. Both CORT and CSDS produced the expected elevation in anxiety-related behavior compared to naïve mice ( Figure   S8). Results from the behavioral tests either correlated moderately with each other (EPM and OFT, r=0.60), were not correlated (OFT and NSF, r=-0.02), or were negatively correlated (OFT and SI, r=-0.46). This indicated that each test captures different aspects of behavior, thereby providing a robust basis to examine how different aspects of behavior might relate to brain mitochondrial health.

Brain mitochondrial health correlates with specific behaviors
We next evaluated the extent to which mitochondrial health in different brain regions was associated with behavior ( Figure 2A). Behavioral scores were transformed so that higher scores on each test indicate higher avoidance/anxiety-like behaviors (as in 23 ) (See Methods, Figure S8). Behavioral scores were then correlated with mitochondrial function across the 6 mitochondrial features, where higher values indicate higher mitochondrial content or functioning.
A frequency distribution of the effect sizes for all mitochondrial-behavior pairs revealed a significant non-zero correlation between MHI in the 17 brain regions and behavior on OFT (p<0.01), EPM (p<0.0001), and SI (p<0.0001), but not for NSF ( Figure 2B, frequency distributions for the other mitochondrial features are shown in Figure S9). Time to feed on the NSF was capped at 600 seconds, so the correlations are less precise than for other tests. To better understand which regions were driving the direction and magnitude in the distributions, we then examined the patterns of correlations for all 17 brain regions, for all six mitochondrial features, across the four behavioral tests ( Figure 2C).
For both OFT and EPM, higher mitochondrial health in the majority of brain regions was 1 1 correlations were of 0.51 for M1 CII and OFT (p=0.025, linear regression), and 0.92 for NAc MHI and EPM (p=0.0005). For NSF there was more heterogeneity between brain regions, with a similar number of positive and negative correlations. Finally, for SI, we observed the opposite relationship; animals with higher mitochondrial functioning generally displayed greater sociability (lower social avoidance) (average r=-0.21, highest for MHI in SN and SI, r=-0.78, p=0.0035).
Interestingly, of the six mitochondrial features, the strongest mito-behavior correlations (for 3 out of 4 behaviors) were for MHI, which may reflect the superiority of MHI as an integrative measure of mitochondrial health over individual enzymatic and molecular features.
As expected, the average correlation between mitochondrial function and behaviors was significantly more consistent for brain mitochondria than for mitochondria in peripheral tissues ( Figure 2D). This finding has two implications. First, it reinforces the specificity of these mitobehavior findings for the brain. Second, it implies that mitochondria across the brain and other tissues within an individual mouse are not equivalent, and differentially regulated. This naturally raised the question whether specific brain regions within an animal could also exhibit independently regulated mitochondrial properties. In other words, whether there exist a functional organization of the brain into separate networks based on mitochondrial properties.

Mitochondrial function-based organization of the brain
To address this question, we first asked whether mitochondrial features in each tissue were independent or correlated with other tissues. Using the same rationale that underlie functional connectivity analysis in fMRI data 52-54 , we first generated a correlation-based similarity matrix of all mitochondrial features across examined brain regions and non-brain tissues ( Figure 3A). Within the 17 brain regions, mitochondrial features were generally positively correlated (average r=0.22, p<1 -100 , two sample t-test), with a few exceptions. Higher mitochondrial activities in one brain region generally also implies higher activities in other regions (p<1 -82 , 1-sample t-test). Further, a modular structure was also apparent, within brain 1 2 and peripheral regions, such that mitochondrial features within each region was more similar than with others regions (p<0.0001, permutation test); Figure 3B also shows the average correlations of each mitochondrial feature with other features among all brain regions -RC enzymes were most strongly co-regulated. To compute how similar each region is to other regions, we computed nodal degree (i.e., average inter-regional correlation) for all regions.
Nodal degree, or global connectivity, was highest in the cerebellum (average r=0.34) and lowest in the brainstem vestibular nucleus (r= 0.05) ( Figure 3C).
As suggested above from the divergence between brain regions and peripheral tissues in the mito-behavior correlation patterns, the mitochondrial features in peripheral tissues were not correlated with brain mitochondria (average r=0.02) or with other peripheral tissues mitochondria (average r=-0.03) ( Figure 3D). The lack of association between brain and nonbrain tissues of the same animals suggest that MHI (as well as content and RC activities) is not a ubiquitous animal-level feature, but rather relatively independently defined in the brain, and further specified within different brain regions.
Given the positive connectivity across the brain and recent evidence of circuit-level metabolic coupling in the Drosophila brain 12 , we examined whether certain regions exhibited particularly strong mitochondrial feature co-regulation. Functional brain networks in the living brain are defined by synchronous activity patterns (e.g., default-mode network, fronto-parietal network). Similarly, we reasoned that connectivity patterns of mitochondrial features (based on single measurements) could reflect networks of brain regions with similar bioenergetic properties. To examine this hypothesis with a framework agnostic to anatomical categorization and inclusive of all mitochondrial features, we performed multi-slice community detection analysis 55 , with mitochondrial features represented in six separate layers ( Figure 3E).
Categorical multi-slice community detection allowed for detecting cohesive groups of brain regions, or communities, that: i) are more similar to each other than they are to the rest of the 1 3 regions, and ii) have cohesion converging across the six layers of mitochondrial features.
To examine whether this functional organization of the brain revealed using mitochondrial features is also evident cross-modally, we used Allen Brain Atlas datasets of brain-wide gene co-expression 56 as well as EYFP-labeled axonal projections that define the structural connectome 57 in the same mouse strain. Specifically, we examined whether gene expression correlations and structural connectivity within brain regions that are functionally grouped together as a network, based on mitochondrial features, is higher than expected by chance. We used two independent graph theoretical metrics: strength fraction (S.F.) 58 and quality of modularity (Q_mod) 59 , to examine whether similar communities exist in gene and structural connectivity data. Using permutation testing, we randomly shuffled community structure of brain regions 10,000 times to determine whether the community structure derived from mitochondrial features is also evident in gene co-expression and structural connectivity data. Both strength fraction and quality modularity statistics suggest that mitochondria-derived networks also have higher similarity in gene co-expression and structural connectome data than expected by chance (S.F.-gene-coexpression, p=0.020; S.F.-structural connectivity, p=0.029; Q_mod-gene-coexpression, p=0.008); and Q_mod-structural connectivity, p=0.015) ( Figure   S10). Hence, this provides convergent multimodal evidence of mitochondrial functional organization overlapping with gene expression and the structural connectome.
Brain region clusters' contributions to the mitochondrial-behavioral correlation 1 4 Finally, to understand the potential significance of this effective brain-wide mitochondrial connectivity in relation to animal behavior, we used our mitochondria-derived network results to partition the brain into three networks, and then analyzed the ability of each network to linearly predict behaviors. The mitochondrial-behavior correlations varied both in strength and direction between the three networks. The average mitochondrial activity in the cortico-striatal network consistently showed the strongest average correlations with behaviors measured in the OFT, EPM, and SI tests (r=0.20, 0.58, -0.41, respectively), whereas the spatial navigation network exhibited the weakest average correlations for the same three behaviors (r=-0.03, 0.18, -0.14) ( Figure 3G). Together, these results suggest the presence and functional significance of brainwide mitochondrial networks, embedded within existing neural circuitry in the mouse brain.

Discussion
Using a high-throughput approach to functionally phenotype hundreds of brain samples from a heterogenous mouse cohort, we have defined brain-wide associations between mitochondrial health and behaviors. Combined with previous findings 29,32 , the diverging mitobehavior associations between brain regions, and between brain and non-brain tissues, brought to light the possibility that different brain regions might exhibit different mitochondrial phenotypes. Examining the network characteristics of mitochondrial features across the brain provides evidence for the modular distribution of mitochondria across cortical and sub-cortical regions, and their relevance to depression-like and anxiety-like behaviors. Based on these data, we conclude that mouse brain mitochondria may exist as behaviorally-relevant networks overlapping with, but distinct from, other modalities such as gene expression and structural connectivity.
Utilizing CORT and CSDS interventions to induce behavioral variation in the population also allowed us to compare the effects of these two interventions on mitochondrial functioning. 1 5 Different brain regions were found to respond, in some cases, in opposite directions, particularly after exposure to CORT. This could be related, among other factors, to regional differences in glucocorticoid and mineralocorticoid receptor density. In contrast to CORT, the recalibrations of brain mitochondria to CSDS was more uniform, with the majority of brain regions exhibiting a coordinated reduction in most mitochondrial features. These results highlight a significant difference between the effects of behavioral stress and endocrine interventions on mouse brain bioenergetics, which should be taken into account when selecting and interpreting rodent models of chronic stress.
Previous research in rodents showed that mitochondrial function in specific brain regions, such as the NAc, is linked to complex behaviors such as social dominance and anxiety 4, 29 . Here we confirm the strong correlation between NAc mitochondrial health and anxious behavior, as the strongest correlation between EPM-based anxiety-like behavior and MHI is in the NAc. We also extend this finding to show that mitochondrial health is also linked to behaviors across distributed brain networks. Similar to the conceptual shift from regional and cellular perspectives towards distributed brain networks and neuronal ensembles 60,61 , our discovery advances the notion that mitochondria may regulate brain function and behavior through distributed mito-networks. This notion is consistent with recent evidence of metabolic coupling among distributed patterns of behavior-related neural activity in Drosophila 12 .
While we cannot directly explain the apparent harmonization of mitochondrial function across brain regions, it may arise from three main factors. First, mitochondria could respond to differences in neuronal circuit functioning, in agreement with observed coupling of neuronal and metabolic activities 12 , such that the cellular infrastructure of neural circuits that fire together not only wire together, but also generate similar mitochondria. A second possibility is that brain regions which are regularly co-activated, i.e., within the same functional networks, harbor similar levels of receptors for neuroendocrine factors (stress and sex steroids) known to influence 1 6 mitochondrial biogenesis and/or functional specialization 62,63 . A third possibility would involve differences in metabolites and substrates that contribute to mitochondrial functioning and regulation, arising similarly among co-activated brain regions. These and other potential factors underlying the modularity of mitochondria within the mouse brain remain to be investigated.
From a mitochondrial biology perspective, our data highlights a potentially important distinction between mitochondrial content (abundance of mitochondria) and RC activity normalized to mitochondrial content, such as the MHI (respiratory chain activity per mitochondrion). These frozen tissue measurements naturally reflect the maximal functional capacity of the mitochondria RC, rather than their actual in vivo rates, which are likely driven by neural activity. Interestingly, the behavioral correlations with content features (mtDNA density and CS activity) for social avoidance behavior were similar to one another, and differed somewhat from RC enzymes. Experimentally, the specificity of these findings highlights the value of parallel assessments of multiple mitochondrial features reflecting unique mitochondrial properties (i.e., mitochondria are not all created equal). In neuroimaging terms, the composite MHI can be understood as the mitochondrial analog to fMRI-based multi-variate pattern analysis (MVPA) 64,65 , where multiple features (mitochondrial enzymes for MHI, voxels for MVPA) are combined to create a more stable and statistically accurate metric of the desired outcome (mitochondrial health for MHI, or brain activation for MVPA). Moreover, mtDNAcn has previously been assessed across multiple brain regions 66 , but our results here go beyond in showing that quantifying mtDNA content on a per-cell basis (mtDNA:nDNA ratio) is heavily skewed by cellularity. As such, it is driven by how many cell bodies are present in the tissue, and correlates poorly with both mitochondrial content and RC activity. Therefore, our data reinforce the notion that mtDNAcn on its own is not a valid measure of mitochondrial health 67,68 . In the mouse brain, our findings shows that mtDNA density per unit of volume (rather than per cell) is a more biologically meaningful mitochondrial feature. 1 8 regions. Interestingly, we observe some parallels in our mitochondria-derived networks with established large-scale networks. For example, the cortico-striatal network identified includes the CPu, NAc, mOFC, mPFC, and motor and visual cortex, which are implicated in decision making and executing actions. The second identified network (Cereb, VN, VTA, Thal, CA3, DGv and DGd) is the most heterogenous but comprised of regions that are known to be connected and involved in salience and spatial navigation 70 . Lastly, the final circuit is made up of limbic and limbic-associated regions (Amyg, Hypoth, PAG and SN) involved in threat responses 71,72 .
In comparing these mitochondria-derived networks to gene expression and structural connectome data we found that the communities significantly overlap across modalities, but are not identical. Finally, each network's mitochondrial health exhibited different associations with behavioral responses. Unsurprisingly, the cortico-striatal network explained the greatest proportion of variance in behavior, increasing the likelihood that the identified large-sclae mitochondrial networks are functionally relevant. In our view, this novel perspective should provide a basis to further delineate the functional brain-wide organization of mitochondria.
Until now, the methods used to reveal neurobiological and metabolic networks in mammals have typically been through indirect functional and structural connectivity analysis 18 .
Here we have developed a scalable approach to examine bioenergetic, and more specifically mitochondrial features across a large number of brain regions in mice with pre-defined behavioral phenotypes. Further, we showed that bioenergetic connectivity was tightly linked with gene co-expression and structural connectivity, thereby providing converging evidence of mitochondria-based networks across modalities. This study synergizes with recent work in Drosophila 12 , providing the technical and empirical foundation to bring mitochondrial biology into brain-wide, network-based models of neural systems in mammals 19,20 . Applied to mitochondria, network-based analytics will contribute to develop more accurate maps, and eventually an understanding of what determines mitochondrial and metabolic properties across 1 9 interacting neural circuits. Developing a spatially-resolved understanding of brain mitochondrial biology in diverse populations of animals and humans should help to resolve the energetic basis of brain function and behavior.

Animals
This study was carried out in accordance with NIH Guidelines, and was approved by the Institutional Animal Care and Use Committee (IACUC) at New York State Psychiatric Institute.
Three month-old CD1 retired breeder male mice were obtained from Charles River, and were used as the aggressors in the social defeat model.

Chronic Social Defeat Stress (CSDS)
Aggressors Prescreening: All CD1 mice used in the experiment were pre-screened for aggressive behaviors as previously described 47 . During a three-day screening procedure, a novel C57BL/6J mouse was placed in the cage of the CD1 mouse for 3 min. C57BL/6J screener mice were not further used in the study. The latency of the CD1 mouse to attack the C57BL/6J screener mouse was recorded. CD1 mice that attacked in less than 1 minute on at least the last two consecutive screening days were considered to be aggressive.
Experimental Groups: Thirteen C57BL/6J mice underwent social defeat stress, while 6 remained in their cages, to serve as the control group. Animals were randomly assigned to these groups. One mouse subjected to social defeat died of unknown causes during the duration of the experiment. Of the remaining 12 social defeat mice, half (n=6) were sacrificed two days after the completion of the stressor, and they are referred to as the 'stressed' CSDS group. Three control mice were sacrificed at the same time. The other half of the CSDS mice (n=6) along with 3 control mice, were allowed an 8.5 weeks (59 days) stress recovery period prior to sacrifice. This group is referred to as the 'recovered' group. Animals were selected for the recovery group from the CSDS group so that there was an even ratio of susceptible to resilient mice in the 'stressed" and 'recovered' groups. 2 1 Social Defeat Paradigm: Adult experimental male C57BL/6J mice (n = 13) were exposed to a CSDS paradigm daily, for 10 days. The experimental C57BL/6J mice were placed in the cage of a new CD1 aggressor mouse for 5 minutes every day, for 10 consecutive days. After the 5 minutes of physical defeat, the C57BL/6J mice were housed in the same cage as the CD1 aggressor, with a perforated plexiglass divider to separate them for 24 hours. After the 10 days of defeats, experimental C57BL/6J mice and CD1 mice were singly housed. Adult control male C57BL/6J mice (n = 6) were housed 2 mice per cage separated by a perforated plexiglass divider. To control for the effects of experimental handling, each control mouse was paired with a new control mouse every day for 10 days.

Corticosterone administration (CORT)
Corticosterone (Sigma Aldrich, St Louis, MO) was dissolved in vehicle (0.45% ßcyclodextrin) at a concentration of 35 µg/ml, equivalent to administration of ~5 mg/kg/day per mouse 45 . C57BL/6J mice (n=5) were group-housed and administered corticosterone in their drinking water. Control mice (n=5) were group-housed and received only vehicle. Animals were randomly allocated to the two groups. Water bottles were prepared twice a week, and the mice never had access to other water. Behavioral testing began on day 56 of CORT administration. Experimental and control mice were sacrificed on day 63 following the completion of behavioral testing. The duration of administration was 9 consecutive weeks.

Behavioral tests
Social Interaction (SI) Test: Social avoidance was measured 24 hrs after the completion of the last day of defeat (day 11). During the first trial, experimental mice were allowed to explore an open field arena (40 cm × 40 cm × 40 cm) containing an empty wire enclosure for 2.5 minutes. During the second trial, a CD1 mouse was placed into the wire enclosure, and the C57BL/J6 mouse was reintroduced for 2.5 minutes. Time spent in the SI zone and time spent in the corner zones during the first and second trial were recorded. SI ratios were calculated as 2 2 time spent in the SI zone during the second trial divided by time spent in the interaction zone during the first trial 23 . Mice with SI ratios of <1 were considered "susceptible" (n = 6), while mice with SI ratios >1 were considered "resilient" (n = 7) 73  the amount of food consumed in 5 min were measured (home cage consumption). NSF latency was capped at 10 min, with animals that did not consume any portion of the pellet receiving a score of 600 sec.
Behavioral Z-Scores: Social avoidance scores were determined by averaging the z scores for 4 measures as previously described 23 ; SI ratio, corner ratio, time spent in SI zone, and time spent in corner zones. The z scores for SI ratio and time spent in SI zone were multiplied by -1, so that higher scores across all 4 measures indicated higher avoidance.
Similarly, the OFT score was determined by averaging the z scores for the 2 measures. The z scores were multiplied by -1, so that higher OFT score indicate higher avoidance of the brightly lit center of the OF. EPM scores were also multiplied by -1, so that higher EPM scores indicate higher avoidance of the open arms of the EPM. High NSF scores already indicate higher avoidance, so they were not inverted. Therefore, across all 4 tests, higher scores indicate higher anxiety-like/avoidant behavior.

Tissue Collection
Animals were sacrificed by rapid decapitation to maintain mitochondrial integrity in brain tissue. Brains were rapidly flash frozen in ice-cold isopentane, stored at -80°C, and later transferred to -170°C (liquid nitrogen vapor) until mitochondrial measures were performed.
Brains were transferred back to -80°C during the week prior to being sectioned and were then transferred to -30°C the night before sectioning. Brains were sectioned coronally on a Leica Model CM3050 S cryostat. The cryostat internal temperature and blade temperature were set to -22°C during sectioning. Brains were mounted onto a specimen disk using optimal cutting temperature (OCT) compound. The brain was sectioned coronally, alternating between two 200μm thick slices and then two 20μm, and were deposited onto microscope slides. Brain sections were kept at -80°C until brain region-specific tissue sample collection. One of the 2 4 recovered CDSD mouse brain contained blood and could not be reliably sectioned, so this animal was excluded from analysis. A second CSDS mouse brain and a control brain both cracked during slicing, so some of the brain regions (n= 7 for CSDS mouse, n=8 for control) could not be obtained from those brains, but the regions that were obtained were included in analysis.
Tissue biopsy punches on frozen brain sections: The scalable Allen Mouse Brain volumetric atlas 2012 75 was used to determine the location of each brain region of interest, and their distance from bregma. The atlas' Nissl-stained images were used as a reference for estimating each section's distance from bregma, which determined which brain slices would be used for punch biopsy collection of each region of interest ( Figure S2B). Brain slices were punched using 1.00mm diameter Robbins True-Cut Disposable Biopsy. Two bilateral punches were collected for each brain region over dry ice. For brain regions in the midline, punches were taken from two consecutive slices. All tissue punch locations were approximated by using landmarks on the slice and by comparing to the atlas. It is important to note that the 1mm puncher was larger than the actual brain region in some instances, and so parts of neighboring brain regions may have been included in the punches. Therefore, the punching technique may include some error. Punches were stored at -80°C until they were ready to be used for enzymatic activity assays.
The tissue punches were too light to be accurately weighed, so we had to estimate the weight based on the reported mouse brain density of 1.04g/cm 376 . The punches were 0.5mm  2 5 Tissue Preparation: The punches from each brain region were homogenized in 0.2ml of homogenization buffer (1mM EDTA and 50mM Triethanolamine) (2x 1mm tissue punch/ 0.2mL of homogenization buffer), with 2 tungsten beads to disrupt the tissues' cells and release the mitochondria. Tissues were homogenized using a Tissue Lyser (Qiagen cat# 85300), which was run at 30 cycles/sec for 1 minute. The tissues were then incubated in ice for 5 minutes, and were then re-homogenized for 1 minute. Tissues were vortexed to ensure homogeneity.

Mitochondrial Measurements
Peripheral tissues were cut over dry ice, weighed, and were then homogenized 1:180 (weight:volume, mg: μ l), except for heart samples that were further diluted to 1:720 to be in the dynamic range of the assays.

Enzymatic activities
Enzymatic activities were quantified spectrophotometrically for Citrate Synthase (CS), Complex I, Succinate Dehydrogenase (SDH, complex II), Cytochrome C Oxidase (COX, complex IV) and and were expressed per mg of tissue, as described previously with some modifications 21 . All miniaturized assay measurements were performed in 96-well plates and enzymatic activity assays recorded on a Spectramax M2 (Spectramax Pro 6, Molecular Devices). Linear slopes reflecting changes in absorbance of the reporter dye were exported to Microsoft excel and converted into enzymatic activities using the molar extinction coefficient and dilution factor for each assay. The assays were optimized to determine the minimal amount of tissue required to obtain reliable results assessed by the C.V. between duplicates. The assays were then further optimized to determine the minimal amount of brain tissue homogenate required for each individual assay. Assay validation involved regressing increasing tissue amounts (number of punches) with observed activities, which confirmed that an increase or decrease in tissue used produced proportional changes in total activity.
Samples were run in duplicates for each enzyme, along with a nonspecific activity control, and every plate had a positive control (heart homogenate). The 96-well plates were 2 6 designed so that each brain region/tissue from all animals were run on a single plate, which prevents potential batch variation for comparisons between the animals. However, due to the size of the plates, no more than 2 types of tissues could be run together on a single plate. This Mitochondrial DNA quantification 2 8 mtDNA density and mtDNA copy number (mtDNAcn) were measured as previously described 77  nDNA Ct. mtDNAcn was calculated by 2 (ΔCt) x 2. For measures of mtDNA density, the Ct value was linearized as 2 Ct / (1/10 -12 ) to derive relative mtDNA abundance per unit of tissue.
In tissues of similar cellular density (number of cell nuclei per area or mass of tissue), mtDNAcn (mtDNA:nDNA) provides an accurate reflection of mtDNA genome density per cell.
However, different brain regions vary widely in their cellularity (up tp 8.5-fold), mostly because some defined areas such as the granular layer of the cerebellum are populated with numerous small cell bodies, whereas other regions such as the molecular layer of the DG are completely acellular and devoid of cell bodies/nuclei. Nevertheless, acellular tissue compartments filled with dentrites can be rich in mitochondria and mtDNA, and therefore the number of mtDNA copies per unit of tissue (μm 3 or mg of tissue) is a more generalizable and accurate estimate of mtDNA density between brain regions.

Mitochondrial Health Index
The mitochondrial health index (MHI) integrates the 5 primary mitochondrial features, yielding an overall score of mitochondrial respiratory chain activity on a per mitochondrion basis 21

Confocal imaging
The immunofluorescent staining was imaged on a Leica SP8 confocal microscope 29.58nm and a z-stack with a step size 0.15μm was taken through the stained area. A 405nm laser was used for exciting the DAPI channel with a power of 8.0% and gain of 60%. A 561nm laser was used for the AlexaFluor 546 channel with a power of 4.0% and a gain of 20%. 3 1 Final images were deconvolved in the Lightning software using the default settings.

SDHA staining analysis
Images were analyzed in ImageJ version 1.52p. Fifty consecutive slices from each stack were selected for a section thickness of 7.5μm, which roughly covered the depth of penetration of the antibodies. The 16bit images were first thresholded using the Triangle algorithm within the integrated Auto-thresholding plugin to create a binary image. The particle analyzer was then used to quantify the percent of the area stained for each slice. The percent area stained from each slice was summated to calculate the percent of volume stained. For compartmentalized data, the molecular and granular layers were separated using a rectangular selection of approximately 50px width through each stack. For the purkinje layer, a rectangle was drawn that narrowly captured the stained area of each slice. The volumetric staining for each region of interest was calculated from the percent of area stained quantified using the particle analyzer.

TDA-based mapper analysis
After creating a delta of mitochondrial features between stressed (CORT or CSDS) and control mice, the data matrix for each group was process through the TDA-based Mapper pipeline 49 . The input data matrix contained 102 concatenated rows for brain regions (17 brain a large number of studies using animal models and computational research suggest that interregional interactions in the brain are multivariate and nonlinear [78][79][80] . Thus, to better capture the 3 2 intrinsic geometry of the data, a nonlinear filter function based on neighborhood embedding was used 49 . The second step of Mapper performs overlapping n-dimensional binning to allow for compression and reducing the effect of noisy data points. Here, we divided the data into lower dimensional space into 64 bins with 70% overlap. Similar results were observed for different number of bins (e.g., 49, and 81). Third, partial clustering within each bin is performed, where the original high dimensional information is used for coalescing (or separating) data points into nodes in the low-dimensional space. Partial clustering allows to recover the loss of information incurred due to dimensional reduction in step one 81 . Lastly, to generate a graphical representation of the shape of the data, nodes from different bins are connected if any data points are shared between them.
The Mapper-generated graphs can be annotated (or colored) using meta-information that was not used to construct the graphs. Here, we annotated these graphs using region-labels to examine whether mitochondrial features were similarity expressed across all regions or whether regional specificity was observed across the two groups. To quantify the extent of segregation (high regional specificity) or integration (low regional specificity) across the six mitochondrial features, we used a graph theoretical measurement of participation coefficient 51 .
Participation coefficient ܲ of a node ݅ is defined as: is the number of links of node to 0 if its links are mostly within its own community (and hence segregated). 3 3 One of the most commonly studied mesoscale aspect of a graph is modularity, where highly modular graphs consist of cohesive groups of nodes (or communities) that are more strongly connected to each other than they are to the rest of the network 82 . Examination of modularity has been recently extended to multi-slice networks that are defined by coupling multiple adjacency matrices across time or modality 55 . Here, we used six slices, each derived from one of the mitochondrial features, where each slice represented weighted adjacency matrix of pair-wise Pearson's correlation between brain regions. We use categorical multi-slice community detection algorithm with the presence of all-to-all identity arcs between slices 55 . The generalized modularity function for multi-slice community detection is given by defined by mitochondrial features. The iterative community detection 85 was run multiple times (1,000), followed by consensus clustering to get stable results for identifying large-scale networks of brain regions 86 .

Comparison of mito-based networks with transcriptomic-and structural connectivitybased networks
To compare our mitochondrial communities against other modalities, we examined the organization of gene co-expression in the mouse brain as well as mouse structural connectome data. We utilized the Allen Brain Atlas' ISH (in situ hybridization) feature, which maps each gene onto a reference standard coordinate atlas image, providing a spatial estimate of transcript levels representing gene expression 56 . Specifically, we used the Anatomic Gene Expression Atlas (AGEA) (https://mouse.brain-map.org/agea), which is a a three-dimensional male adult C57BL/6J mouse brain atlas based on the ISH gene expression images. The AGEA feature allows for selecting both a 'seed voxel' and a target 'selected voxel' from exact brain coordinates coordinates, yielding the transcriptome-wide correlation between the seed and target regions. The correlation is a measure of the average co-expression between the two voxels. The co-expression values were obtained for all possible pair of brain regions (17x17 matrix), yielding a gene co-expression correlation matrix to which the structure of our mitobased networks could be compared using a permutation text.
The structural connectome data was obtained from The Allen Mouse Brain Connectivity Atlas 57 , which is a mesoscale connectome of the adult mouse brain. We utilized the normalized projection strength values for the all brain regions of interest. Because this atlas does not distinguish between dorsal and ventral dentate gyrus, data was obtained for 16 brain regions.
For three brain regions where the Atlas provides connectivity values for multiple subregions (OFC, VN, Cerebellum), the connectivity values were averaged across the subregions to yield a global measure for each region, thus matching the dimensionality of our mitochondrial dataset.

Statistical analysis
Standardized effect sizes (Hedge's g) were computed to quantify the effect of stress conditions on mitochondrial features. Significant effect sizes were determined by the 95% confidence intervals. Two-way ANOVAs were used to compare the effects of the two stressors on the mitochondrial measures across all tissues, and to compare the average mitochondrialbehavior correlations between brain regions and tissues for all behaviors. A binomial test was used to determine if the proportion of brain regions exhibiting an increase or decrease in mitochondrial features differed from chance level (50%). Correlations between behavioral Zscores and mitochondrial activities were estimated using Spearman's r (r) to guard against inflation. Frequency distributions of the effect sizes for the two stressors and of the mitochondrial-behavior correlations were fitted with Gaussian curves and analyzed by one sample t tests compared to the null distribution. Correlations between tissues' mitochondrial activities were measured as Pearson's r. A survival curve using Mantel-Cox log-rank test was generated for novelty suppressed feeding test latencies because the test is capped at 600 sec.
Hierarchical clustering was performed on the data using Euclidian distance with Ward's clustering algorithm. Statistical analyses were performed with Prism 9 (GraphPad) and Metaboanalyst version 4 87 .

Funding
Work of the authors is supported by GM119793, MH122706, MH108719, MH090964, MH111918, MH119735.    Fig. 8). (C) Individual correlations for the 17 brain region, across the 6 mitochondrial features, for each behavioral score (see Extended Data Figure 9 for additional details) quantified as Spearman's r. The strongest correlations for each behavioral test are plotted below. An adjusted p-value of <0.002 was applied (false-discovery rate 1%). All tests have been adjusted so that a higher score indicates higher anxiety-like behavior (see Methods for details). (D) Average correlation of each behavior for the brain (B) and tissue (T) mitochondrial features; 2-way ANOVA. (F) three distinct communities or brain networks. Modular structure confirmed by permutation test, p<0.0001). (G) Average mito-behavior correlation by module for each behavioral test (top). The bottom panel shows networks with each region color-coded by its average correlation with behaviors; for OFT and EPM (left) and SI (right). Modularity metrics anchored in whole-brain transcriptome and the structural connectome data are shown in Extended Data Fig. 10.

brain regions
Behavioral Tests

Connectivity matrix based on mitochondrial features
Brain regions (17) Peripheral tissues (5)  Average Inter-regional Correlation

Global connectivity by brain region
Average inter-regional correlation