Selectively breeding rats for high or low locomotor activity in a novel environment (LocoScore) produced the bHR line (Wakil:bHR, RRID:RGD_405847397) and bLR line (Wakil:bLR, RRID:RGD_405847400), respectively (Stead et al., 2006). After 37 generations, 12 bHRs and 12 bLRs (F₀) were chosen from 24 distinct families to crossbreed. The F₁ offspring with similarly high or low LocoScores were then bred with each other to produce a re-emergence of diverse behavioral phenotypes in the F₂ generation (Figure 1). These 48 F₂ litters generated 540 rats [n = 216 behaviorally phenotyped as juveniles (1 month old), n = 323 behaviorally phenotyped as young adults (2–4 months old)]. Our current study sampled a subset of the adults (F₀: n = 24, n = 6/phenotype per sex; F₂: n = 250, n = 125/sex) that overlapped with previous genetic studies (Chitre et al., 2023; Zhou et al., 2019) but was distinct from our previous male-only meta-analysis (Birt et al., 2021).

Behavioral phenotyping for the F₀ and F₂ rats used in our current study was performed in adulthood, in the morning, with each test occurring on separate days. Testing order was counterbalanced for phenotype (bHR/bLR), with males tested before females on separate days. All F₀ and F₂ rats were tested for locomotor activity in a novel environment (protocol: Stead et al., 2006). During testing, rats were placed in a box akin to their home cage but located in a different room with novel cues. Over 60 min, lateral and rearing movements were counted via beam breaks, and the cumulative total defined as locomotor score (LocoScore). All F₀ and F₂ were tested for exploratory and anxiety-related behaviors on an EPM [dimly lit (40 lux), protocol: Chitre et al., 2023]. Rats began the five-minute test at the intersection of the arms. A video tracking system (Ethovision, Noldus Information Technology) recorded the percent time spent in the open and closed arms, distance traveled (cm), and time immobile (sec). A subset of F₂s (n = 209) subsequently underwent seven sessions of PavCA training to measure their bias in favor of reward cues over the reward itself (protocol: Meyer et al., 2012). To create a summary PavCA Index, three behavioral variables were averaged: (1) Probability difference: the probability of lever contact minus the probability of food magazine entries, (2) Response bias: the total conditioned stimuli (CS: lever) contacts minus the total food magazine entries, divided by the sum of the two behaviors. (3) Latency score: the average latency to enter the food magazine minus the average latency to contact the lever, divided by the length of the CS duration (8 s). The PavCA Index for the last 2 days of testing (days six and seven) was used to classify rats as “sign trackers” (ST: values > 0.5), “intermediate” (IN: values −0.5 to 0.5), or “goal trackers” (GT: values < −0.5).

For each of the continuous behavioral variables, the interacting effects of sex and phenotype (F₀ bHR vs. F₀ bLR vs. F₂) were examined using analysis of variance (ANOVA type 3, contrasts = “contr.sum”). Exploratory analyses were also performed to examine the potential effect of batch-related variables, testing order, maternal lineage, and paternal lineage. Correlations between behaviors were characterized in the F₀ and F₂ datasets separately using parametric methods (Pearson’s R, simple linear model). For the F₀s, Lineage was included as a dummy variable (bHR as reference 0, bLR coded as 1). For PavCA, a Fisher’s Exact Test was performed on the ratios of male to female animals classified as ST, IN, or GT.

Adults (postnatal days 113–132) were decapitated without anesthesia and brains rapidly extracted (<2 min). For the F₀s, whole hippocampus was immediately dissected and flash frozen. For the F₂s, whole brains were flash frozen, and later hole punches from the dorsal and ventral hippocampus were pooled from four hemisected coronal slabs per rat (−2.12 to −6.04 mm Bregma; Paxinos and Watson, 2013). RNA was extracted using the Qiagen RNeasy Plus Mini Kit. A NEB PolyA RNA-seq library was produced and sequenced using a NovaSeq S4 101PE flowcell (targeting 25 million reads/sample).

RNA-Seq data preprocessing was performed using a standard pipeline including alignment (STAR 2.7.3a: genome assembly Rnor6), quantification of gene level counts (Ensembl v103, Subread version 2.0.0), and basic quality control. All downstream analyses were performed in Rstudio (v.1.4.1717, R v. 4.1.1). Transcripts with low-level expression (<1 read in 75% of subjects) were removed. Normalization included the trimmed mean of M-values (TMM) method (Robinson and Oshlack, 2010, edgeR v.3.34.1; Robinson et al., 2010), and transformation to Log2 counts per million [Log2 cpm (Law et al., 2016), org.Rn.eg.db annotation v.3.13.0; Carlson, 2019]. Following quality control, the F₀ dataset contained n = 23 subjects (subgroups: n = 5 bHR females, n = 6 for each of the other subgroups: bHR males, bLR females, bLR males) with Log2 cpm data for 13,786 transcripts, and the F₂ dataset contained n = 245 subjects (subgroups: n = 122 males, n = 123 females) with Log2 cpm data for 14,056 transcripts.

Differential expression was calculated using the limma/voom method (Law et al., 2014, package: limma v.3.48.3) with empirical Bayes moderation of standard error and FDR correction. For the F₂s, the same differential expression model was used for each variable of interest (LocoScore, EPM % time in open arms, EPM distance traveled, EPM time immobile, PavCA Index). Technical co-variates were included if they were strongly related to the top principal components of variation identified using Principal Components Analysis (PCA) or had confounding collinearity with variables of interest [covariates: percent of reads that were intergenic (%intergenic) or ribosomal RNA (%rRNA), dissector, sequencing batch, and PavCA exposure (“STGT_Experience”)].

Equation 1: F₀ differential expression model:

Equation 2: F₂ differential expression model:

The full F₀ and F₂ differential expression results were compared to our previous meta-analysis of hippocampal RNA-Seq studies from late generation bHR/bLR males (Supplementary Table S2 in Birt et al., 2021) using parametric and non-parametric methods (Pearson’s and Spearman’s correlation of Log2FC values) and visualized using gene rank-rank hypergeometric overlap [RRHO v. 1.38.0 (Plaisier et al., 2010; Rosenblatt and Stein, 2014), ranking by t-statistics] and VennDiagram [v.1.7.3 (Chen, 2022)]. For downstream analyses, we defined bHR/bLR differentially expressed genes as the 1,063 genes with FDR < 0.10 in either the F₀s or late generation meta-analysis, or nominal replication (p < 0.05) in both with consistent direction of effect.

As an exploratory analysis, we compared bHR/bLR hippocampal differential expression to the pattern of differential expression in other brain regions in previous small transcriptional profiling studies of bHR/bLR adults, including publicly available data from the amygdala [GSE88874: n = 5/group, generation F31 (Cohen et al., 2015; McCoy et al., 2017), GSE86893: n = 6/group, generation F34-F36 (Cohen et al., 2017)] and dorsal raphe [GSE86893, n = 6/group (Cohen et al., 2017)], and unpublished data from the cortex (GSE286181, n = 6/group) and hypothalamus (GSE286181, n = 6/group) from an early generation of selective breeding (F4). To run this comparison, differential expression was calculated for each dataset using the limma pipeline (Supplementary methods). To reduce noise and increase statistical power, we used a standardized pipeline (Hagenauer M. et al., 2024) to perform a simple random effects meta-analysis (Viechtbauer, 2010, package: metafor) to summarize the amygdala differential expression results (collective sample size of n = 11/group for 7,133 genes), and—to potentially identify bHR/bLR differences that might exist brain-wide—all non-hippocampal differential expression results (collective n = 29/group for 11,509 genes). These results were compared to our hippocampal findings using non-parametric methods [Spearman’s correlation of Log2FC values and RRHO v. 1.38.0 (Plaisier et al., 2010; Rosenblatt and Stein, 2014), ranking by t-statistics].

To elucidate functional patterns, we ran Gene Set Enrichment Analysis (fgsea v.1.2.1, nperm = 10,000, minSize = 10, maxSize = 1,000, FDR < 0.05) using a custom gene set database (Brain.GMT v.1, Hagenauer M. H. et al., 2024) that included standard gene ontology, brain cell types, regional signatures, and differential expression results from public databases. We created a continuous variable representing bLR-like vs. bHR-like differential expression for each gene by averaging the t-statistics for bLR vs. bHR comparisons in the F₀ dataset and former late generation RNA-Seq meta-analysis (Birt et al., 2021) and for each of the F₂ behaviors (with bHR-like phenotype set as reference). A second non-directional analysis used the absolute value of the average t-statistic.

Hippocampal cis-eQTL mapping was performed using published methods Munro et al. (2022, unpublished)¹. Quality-controlled F₂ RNA-Seq data (Log2 CPM, n = 245 following quality control) was corrected for technical covariates (Equation 2, residualized), followed by rank-based inverse normal transformation ((see footnote 1) Ongen et al., 2016). F₂ genotypes were generated by low coverage whole genome sequencing followed by imputation (from data release for (Chitre et al., 2023): doi: 10.6075/J0K074G9, n = 4,425,349 SNPs). Principal Components Analysis was run on the gene expression matrix and genotype matrix [following pruning for linkage disequilibrium, Plink2 v.2.00a2.3 (Chang et al., 2015)], and principal components 1–5 from both analyses included as covariates within the single-SNP linear regression for cis-eQTL mapping [tensorQTL v.1.0.6 (Taylor-Weiner et al., 2019)]. We tested SNPs with minor allele frequency (MAF) >0.01 within ±1 Mb of each gene’s transcription start site (tss). A significant eVariant-eGene relationship (cis-eQTL) was defined using empirical beta-approximated p-values calculated using permutations for each gene, and false discovery corrected (FDR < 0.05) using results from the top SNP for all genes. When SNPs were in perfect linkage disequilibrium, a single SNP was selected randomly. Additional, conditionally independent cis-eQTLs for each eGene were identified using stepwise regression (tensorQTL: default settings). We estimated cis-eQTL effect size (allelic fold change, aFC) using an additive cis-regulatory model (package aFC.py) with raw expression read counts and the same covariates as cis-eQTL mapping (Mohammadi et al., 2017).

We extracted F₀ genotype information for each eVariant (n = 10 bHR/n = 10 bLR in Chitre et al., 2023, release: doi: 10.6075/J0K074G9) using VcfR (v1.14.0, Knaus and Grünwald, 2017).² We defined partial bHR/bLR segregation using myDiff() G_st’ > 0.27 [(Hedrick, 2005), akin to all 0/0 vs. all 0/1 in our dataset]. We assigned the direction of effect for the Log2 aFC to reflect the bLR-enriched allele vs. bHR-enriched allele and compared these predictions to both the F₀ differential expression results (Log2FC) and bHR/bLR late generation RNA-Seq meta-analysis results (estimated d) using parametric (linear regression) and non-parametric (Spearman’s rho) methods.

We used Summary Data-based Mendelian Randomization (SMR; Zhu et al., 2016) to test for colocalization of cis-eQTLs with QTLs from the full F₂ cohort [GWAS results: (Chitre et al., 2023)] for adult behaviors included in our differential expression analysis (LocoScore, EPM time immobile, EPM distance traveled, EPM % time in open arms, PavCA Index) and juvenile behaviors targeting analogous traits (open field (OF) time immobile, OF distance traveled, OF % time in center). Z-scores for the cis-eQTL and GWAS associations with each top eVariant were used to calculate the SMR approximate chi-squared test statistic, with p-values determined using the chi-squared distribution’s upper tail [df = 1, FDR correction: mt.rawp2adjp() (proc = “BH”) in multtest v.2.26.0 (Pollard et al., 2005)]. Results were visualized using the manhattan() plot function in the qqman package [v.0.1.8 (Turner, 2018)].

To determine whether the strength of cis-eQTL/QTL co-localization (SMR t-statistic) correlated with F₂ behavioral differential expression, we assigned a predicted direction of effect based on the relationship between genotype and behavior within the larger F₂ sample (adults: n = 323 adults, juveniles: n = 216) and genotype and expression within the cis-eQTL analysis (n = 245). We then examined the correlation between the F₂ Log2FCs for each adult behavior and the “directional” SMR T-statistics for the same adult behavior or analogous juvenile behavior (OF distance traveled, OF time immobile, OF % time in center), both in the full dataset (all 5,937 cis-eQTLs) and within the subset of cis-eQTLs that we had already confirmed were segregated in bHR/bLRs in a manner predictive of bHR/bLR differential expression (492 cis-eQTLs representing 456 eGenes).

To narrow down our pool of top candidate genes for mediating the effect of genetic variation on behavior, we required that our final top candidate genes have expression strongly related to genetic variation (cis-eQTLs) that is segregated in bHR/bLR (G_st’ > 0.27) in a manner that correctly predicts bHR/bLR differential expression and is co-localized with a QTL for behavior (SMR FDR < 0.10) in a manner that correctly predicts F₂ differential expression. Using a conservative estimate (Supplementary methods), this convergence of results should only be observable once, at most, in our dataset due to random chance.

The bHR/bLR (F₀) crossbreeding scheme produced F₂ animals with behaviors ranging between the more extreme bHR and bLR phenotypes [all behaviors: p < 1.5e-06 for effect of group (F₀ bHR vs. F₀ bLR vs. F₂); examples: Figures 3A,B; Supplementary Figure S4, full statistics: Supplementary Table S1]. F₂ behavior sometimes appeared more similar to bLRs than bHRs (e.g., Figure 3B) suggesting a floor effect or that genetic contributions to internalizing-like behavior may be more dominant.

Exploratory locomotion, anxiety-like behavior, and PavCA behavior remained strongly correlated within the F₂s, as previously observed in the bHR/bLR lines (examples: Figures 3C,D; Supplementary Figure S4, full statistics: Supplementary Table S2). Although these behaviors often differed by sex (example: Figure 3E, Supplementary Figure S4, full statistics: Supplementary Table S1), sex differences were not responsible for driving the correlation between different behaviors (with sex in the model: all behavior–behavior relationships still p < 0.0284). These findings are consistent with results using the full F₂ cohort (Chitre et al., 2023) and imply that locomotor activity in a novel environment echoes a broader behavioral temperament, reflecting genetic and environmental influences shared across anxiety, mood, and reward-related behaviors.

The hippocampus plays an important role in many processes relevant to bHR/bLR behavioral phenotype, including novelty processing, exploratory behavior, behavioral inhibition, emotional regulation, environmental reactivity, and stress-related responses. We observed robust bHR/bLR differential expression in the hippocampus. Within the F₀ RNA-Seq dataset, there were 131 differentially expressed genes with elevated expression in bLRs versus bHRs, and 86 differentially expressed genes with higher expression in bHRs (False Detection Rate (FDR) < 0.10, Figures 4A,C, Supplementary Table S3). In contrast, despite the observed sex differences in behavior, there were only 21 genes upregulated in females (versus males) and 22 genes upregulated in males (versus females) (Figures 4B,C, Supplementary Tables S3, S4). The effect sizes (Log(2) Fold Changes, or Log2FC) for bHR/bLR differentially expressed genes were also larger than those for sex, with the exception of a few X and Y chromosome genes (Figures 4A,B, Supplementary Tables S3, S4). There were no significant interactions between the effects of Lineage and Sex on gene expression (FDR > 0.10), but our sample size was underpowered to detect these effects (n = 5−6/subgroup). The presence of robust bHR/bLR hippocampal differential expression in both male and female F₀s replicated previous male-only studies (Birt et al., 2021).

The bHR/bLR hippocampal differential expression in our current study replicated many effects observed in our previous meta-analysis of hippocampal transcriptional profiling studies in bHR/bLR males (Birt et al., 2021), with the F₀ Log2FC correlating positively with the bLR versus bHR estimated effect size (d) observed in RNA-Seq data from later generations (Figures 5A,B). Sixty-two of the 984 bHR/bLR differentially expressed genes in either dataset (FDR < 0.10) were significant (FDR < 0.10) in both datasets (Figure 5C, Supplementary Table S3). More genes showed replication of nominal bHR/bLR effects (p < 0.05) with consistent direction of effect in both datasets, so that, in total, 1,063 genes had evidence of bHR/bLR differential expression in the hippocampus (Figures 5D,E, Supplementary Table S3). As many generations of selective breeding for a behavioral phenotype are likely to produce an enrichment of eQTL alleles influencing the phenotype, these 1,063 genes were prioritized in downstream analyses.

As an exploratory analysis, we also compared the pattern of bHR/bLR differential expression identified in the hippocampus in our current study to bHR/bLR differential expression in other brain regions using data from previous transcriptional profiling studies in bHR/bLR adults, including the amygdala [n = 5/group (Cohen et al., 2015; McCoy et al., 2017), n = 6/group (Cohen et al., 2017)], dorsal raphe [n = 6/group (Cohen et al., 2017)], and unpublished data from the cortex (n = 6/group, Supplementary Table S5) and hypothalamus (n = 6/group, Supplementary Table S5). To increase power, we performed a meta-analysis of the two amygdala datasets (collective n = 11/group, 7,133 genes, Supplementary Table S5) and a meta-analysis encompassing all of the non-hippocampal data to identify bHR/bLR differences that might exist brain-wide (n = 29/group, 11,503 genes, Supplementary Table S5). These comparisons suggested that at least some of the bHR/bLR differential expression identified in the hippocampus may also be present in other brain regions, whereas other differential expression may be hippocampal specific (Supplementary results).

Since some bHR/bLR differential expression may be due to either linkage disequilibrium with causal variants or genetic drift specific to our bred lines, we performed RNA-Seq on hippocampal tissue from a large F₂ intercross sample (n = 250) to identify differential expression that continued to independently correlate with exploratory locomotion, anxiety-like behavior, and reward-related behavior. Hippocampal gene expression associated with bLR lineage resembled the expression associated with lower F₂ LocoScore, as indicated by the negative correlation between the F₀ bLR vs. bHR Log2FCs for all genes and the F₂ LocoScore Log2FCs for all genes (R = −0.20, p < 2.2e-16, Figure 6A, Supplementary Table S6). Similarly, gene expression associated with bLR lineage partially resembled expression in F₂s exhibiting lower exploration (distance traveled) and greater anxiety (greater time immobile, less time in the open arms) on the elevated plus maze (EPM) task (Supplementary Table S6), and greater goal-tracking behavior on the PavCA task (lower PavCA Index; Supplementary Table S6). This pattern of correlations confirmed that a portion of the hippocampal differential expression that emerged following selective breeding was related to behavioral temperament.

These correlations strengthened when focusing specifically on bHR/bLR differentially expressed genes (1,063 genes in Figures 5C–E), of which 1,045 were present in the F₂ dataset (Supplementary Table S6). Of these genes, 111 showed both upregulation in bLR rats and at least one nominal (p < 0.05) association with bLR-like behavior in the F₂s (i.e., expression correlated with decreased locomotor activity, decreased EPM distance traveled, increased EPM time immobile, decreased EPM % time in open arms, or decreased PavCA Index), a 3.27X enrichment beyond random chance (Fisher’s exact test: 95%CI: 2.61–4.08, p < 2.2e-16, Figure 6B), and 81 genes showed both upregulation in bHR rats and at least one nominal (p < 0.05) association with bHR-like expression in the F₂s, a 2.72X enrichment beyond random chance (Fisher’s exact test: 95%CI: 2.10–3.51, p = 7.14e-13, Figure 6C).

However, we were unable to identify differentially expressed genes for F₂ behavior with strong enough effects to survive false discovery rate correction (FDR < 0.10). This was also true when using a model that included sex-specific differential expression for F₂ behaviors (sex*behavior interaction: all FDR < 0.10). This inability to detect significant differential expression related to F₂ behavior was particularly striking because the F₂ sample size was much larger than the sample sizes used to detect differential expression in our bred model (F₂: n = 250, F₀: n = 24), and this greater statistical power lead to the expected increase in the detection of more subtle differential expression related to sex (1,679 genes with FDR < 0.10, Supplementary Figure S5, full results: Supplementary Table S4).

These findings drive home the role of cumulative small, polygenic effects in generating complex behavior, and suggest a need for larger sample sizes to reliably detect these polygenic effects on gene expression in a heterogeneous population. These findings also demonstrate the utility of selective breeding in behavioral genetics: the highly divergent phenotype and minimized within-group variability made it possible to detect relevant differential expression in a much smaller size. For downstream analyses, we chose to focus on the differential expression with the strongest converging evidence supporting its potential to mediate behavioral temperament from both the selectively bred lines and F₂ rats (the 111 genes that were upregulated in bLRs and nominally with bLR-like behavior in the F₂s and 81 genes upregulated in bHRs and nominally with bHR-like behavior in the F₂s).

To determine generalizability, we compared our list of differentially expressed genes implicated in behavioral temperament by converging evidence from the bred lines and the F₂s (identified in Figures 6B,C) to a database of 2,581 genes previously identified as differentially expressed in the hippocampus of other bLR-like and bHR-like rat models targeting hereditary behavioral traits resembling extremes on the internalizing/externalizing spectrum [database from Birt et al. (2021), summarized in Figure 7A; Andrus et al., 2012; Blaveri et al., 2010; Díaz-Morán et al., 2013; Garafola and Henn, 2014; Meckes et al., 2018; Raghavan et al., 2017; Sabariego et al., 2013; Wilhelm et al., 2013; Zhang et al., 2005]. Sixteen of 111 genes that were upregulated in bLRs and with bLR-like behavior in our study were also upregulated in other bLR-like models (Figure 7B, enrichment OR: 2.50 (95%CI: 1.37–4.29), Fisher’s exact test: p = 0.00242) and 14/81 genes that were upregulated in bHRs and with bHR-like behavior in our study were downregulated in other bLR-like models [Figure 7C, enrichment OR: 2.07 (95%CI: 1.07–3.74), p = 0.0189]. Notably, Tmem144 had elevated hippocampal expression in three other bLR-like rat models (Figure 7D, Blaveri et al., 2010; Meckes et al., 2018; Wilhelm et al., 2013). Five other genes were differentially expressed in two other rat models (Figure 7D, Bphl, Ist1, RGD1359508, Nqo2, Fcrl2). We expect less than one gene (0.39) in our dataset to show this degree of convergence due to random chance (Supplementary methods).

To ascribe functional trends to the differential expression associated with behavioral temperament, we performed Gene Set Enrichment Analysis using a combined score for each gene summarizing bLR-like vs. bHR-like expression across the bHR/bLR and F₂ analyses. Sixty-three gene sets were upregulated (FDR < 0.05) with a bHR-like phenotype (i.e., in bHRs and with bHR-like F₂ behavior; Supplementary Table S7). Nineteen of these implicated hippocampal subregions or cell types, mostly neuronal (n = 10), emphasizing GABA-ergic cells (n = 3) and dendrites (n = 3). The dentate gyrus was implicated (n = 1), epithelial cells (n = 4, including gene C2cd3) and vasculature (n = 3, including Mfge8). Fourteen gene sets were derived from previous differential expression experiments (Baker et al., 2012; Lim et al., 2021), with most related to stress or fear conditioning (n = 12, upregulated: n = 9, including gene C2cd3). Other implicated functions included nervous system development, proliferation, and cell fate (n = 13, including genes Mfge8, Nqo2, Ucp2, and C2cd3) and transcription regulation (n = 9, including Ucp2).

Thirty-seven gene sets were upregulated (FDR < 0.05) with a bLR-like phenotype (i.e., in bLRs and with bLR-like F₂ behavior). Eleven of these implicated hippocampal subregions or cell types, especially microglia (n = 8, including gene Tmem144). Other emphasized pathways included mitochondrial function, oxidative phosphorylation, and cellular respiration (n = 6, including genes Wdr93 and Idh1), metabolism (n = 5, including Pex11a, Lsr, Ist1, and Idh1), and immune response (n = 4). A non-directional analysis produced weaker results (12 gene sets with FDR < 0.10) highlighting similar functions (metabolism: n = 3, including Pex11a, Lsr, Ist1, Mcee, and Idh1; and microglia: n = 4, including Fcrl2 and Tmem144). Gene sets related to a bLR-like model, Flinders Sensitive Line, were also highlighted (n = 3).

To identify hippocampal gene expression that might be influenced by proximal genetic variation, we integrated our current F₂ RNA-Seq data (n = 245) with previous genotyping results [n = 4,425,349 single nucleotide polymorphisms (SNPs) (Chitre et al., 2023)] to identify 5,351 genes (eGenes) with hippocampal expression tightly correlated (FDR < 0.05) with nearby genetic variation [cis-eQTLs: within ±1 MB of the transcription start site (TSS)]. Using stepwise regression, we identified additional conditionally-independent cis-eQTLs beyond the strongest cis-eQTL for each eGene (Supplementary Figure S6A), distinguishing a final total of 5,937 cis-eQTLs representing 5,836 unique eVariants. Like previous cis-eQTL analyses, these eVariants were predominantly located within ±400 kB of the TSS of their respective eGene (Supplementary Figure S6B). A comparison with existing rat cis-eQTL databases [RatGTEx: Hong-Le et al. (2023), and other tissues in RatGTEx: https://dx.doi.org/10.17504/protocols.io.rm7vzyk92lx1/v1], Supplementary Figure S6C, indicated that most hippocampal eGenes were also significant eGenes within at least four other tissues (out of 11 tissues characterized, Supplementary Figure S6D), and confirmed that previously-identified brain ciseQTLs showed a similar direction of effect on gene expression within the hippocampus (R = 0.67–0.75, rho = 0.65–0.77, Supplementary Figure S8, S9) when there was at least a nominal (p < 0.05) relationship in our dataset, although many cis-eQTLs remained region specific. As our hippocampal cis-eQTL database represents a valuable resource for the interpretation of rat genomic results, we have shared it on RatGTEx.³

We used our cis-eQTL database to predict the effect of genetic variation that segregates the bHR/bLR lines on gene expression. Many eVariants (2,452) showed at least partial bHR/bLR segregation in the F₀ rats (n = 10 bHR/n = 10 bLR sequenced in Chitre et al., 2023), such that if all subjects from one phenotype (e.g., bHRs) had 2 reference alleles (0/0), all subjects from the other phenotype had at least 1 alternate allele (0/1); population segregation statistic G_st’ > 0.27 (Hedrick, 2005). To predict the effect of these bHR/bLR segregated eVariants on gene expression, we calculated the allelic Log2FC (aFC) for each eVariant and assigned the direction of effect based on the allele frequency within the bLR vs. bHR F₀ rats (Figure 8A). These predictions correlated strongly with the F₀ differential expression results (Figures 8B,C, 2,500 eGene/eVariant combinations: R = 0.77, rho = 0.63, p < 2e-16) and our previous bHR/bLR late generation meta-analysis effect sizes (2,114 eGene/eVariant combinations: R = 0.52, rho = 0.61, p < 2e-16, Supplementary Figure S10). These results validated our hippocampal cis-eQTL database and confirmed that bHR/bLR differential expression of eGenes is likely driven by bHR/bLR genetic segregation.

To identify cis-eQTLs that might mediate the influence of genetic variation on behavior, we determined which hippocampal cis-eQTLs co-localized with regions of the genome associated with bHR/bLR-like behavior (QTLs) within the larger F₂ sample [adults: n = 323 adults, juveniles: n = 216 (Chitre et al., 2023)] using Summary Data-based Mendelian Randomization (SMR; Zhu et al., 2016). We focused on QTLs for behaviors measured in F₂ adults that were included in our differential expression analysis (LocoScore, EPM time immobile, EPM distance traveled, EPM % time in open arms, PavCA Index), and for analogous behaviors measured in an independent sample of F₂ juveniles (open field (OF) time immobile, OF distance traveled, OF % time in center). This analysis identified 79 cis-eQTLs that were co-localized with QTLs for LocoScore (FDR < 0.10), including 1 cis-eQTL that was also co-localized with a QTL for EPM distance traveled (FDR < 0.10), and 13 of the 14 cis-eQTLs that were co-localized with QTLs for OF distance traveled (FDR < 0.10). Most cis-eQTLs that strongly co-localized with behavioral QTLs were on chromosome 1, as expected due to the strength of the QTLs on this chromosome (Figures 9A,B).

To narrow down our pool of top candidate genes for mediating the effect of genetic variation on behavioral temperament, we used converging information from our different samples and analyses. First, we narrowed our scope to cis-eQTLs that we had confirmed are segregated in bHR/bLRs with differential expression matching predictions based on the distribution of alleles in the two lines (Figure 8C; 492 cis-eQTLs representing 456 eGenes). Within this subset of cis-eQTLs, the strongest co-localization with QTLs tended to predict F₂ differential expression with behavior, especially when considering the predicted direction of effect based on the relationship between genotype and behavior within the larger F₂ sample (adults: n = 323 adults) and genotype and expression within the cis-eQTL analysis (n = 245) (Figure 10A). This was particularly true for LocoScore (Figure 10B, R = 0.56, p < 2.2e-16), but also other F₂ adult behaviors (Supplementary Figure S11, R = 0.33–0.57, all p < 3.07e-14). It was also true when comparing F₂ differential expression to SMR co-localization results with QTLs for two analogous juvenile behaviors (Figure 10C, Supplementary Figure S12, OF distance traveled: R = 0.35, p = 1.48e-15, OF time immobile: R = 0.26, p = 4.44e-09).

The most compelling candidates for mediating the effect of genetic variation on behavioral temperament should have expression strongly related to genetic variation (cis-eQTLs) that is segregated in bHR/bLR, correctly predicts bHR/bLR differential expression, and co-localizes with a QTL for behavior that correctly predicts F₂ differential expression associated with that behavior (Figure 2). Among the SMR results, 16 of the 80 genes surviving FDR correction (FDR < 0.10) met all these criteria (Figure 11, examples: Supplementary Figures S13–S16). By conservative estimate (Supplementary methods), one gene or less in our dataset is expected to show this degree of convergence due to random chance. These 16 genes were clustered within seven genomic regions on chromosomes 1, 7, and 19, suggesting that there still remained some false discovery due to linkage disequilibrium (Figures 9A,B). That said, when cross-referencing with functional annotation, eight of these candidate genes—representing five of the identified regions—were clearly related to mitochondrial function and bioenergetics (Figures 11, 12), hinting at the relevant genes in each region.

Among these 16 top candidate genes, eight are directly involved in bioenergetics (Figure 12). The differential expression results overall similarly showed upregulation in gene sets related to mitochondria, oxidative phosphorylation, and metabolism in bLR-like vs. bHR-like animals. These findings complement previous evidence that adult bLRs have elevated oxidative phosphorylation in the hippocampus, as indicated by increased mitochondrial oxygen consumption and elevated electron transport chain activity (McCoy et al., 2019). Since the expression of our top candidate genes was strongly correlated with genetic variation tied to behavioral phenotype in both bHR/bLR and F₂ samples, our results imply that variation in energy production may mediate the effect of heredity on temperament and provide insight into the responsible mechanisms.

In particular, our results suggest that bLR-like animals have enhanced fatty acid oxidation, which is a pathway that is particularly important during times of high energy usage (Nsiah-Sefaa and McKenzie, 2016) because it can release twice as much energy as glucose metabolism (Schönfeld and Reiser, 2013). bLR-like animals had upregulation of multiple fatty acid oxidation gatekeepers [Lsr, Ist1, and Pex11a (Chang et al., 2019; Mattiazzi Ušaj et al., 2015; Mindthoff et al., 2016; Renne and Ernst, 2023)]. Downstream, there was also upregulation that could facilitate the tricarboxylic acid (TCA) cycle (Idh1) (Gherardi et al., 2020) and electron transport chain (Wdr93) (GeneCards, n.d.; InterPro, n.d.; Meyer et al., 2009) to increase energy production. These findings have widespread functional implications, as the brain consumes disproportionate energy to maintain neurotransmission and synaptic repolarization (van Rensburg et al., 2022), especially during times of heightened activity, such as stress (Zalachoras et al., 2020).

In contrast, energy production in bHR-like animals may be kept under tight regulation by upregulation of Spg7, Afg3l1, Ucp2, and Lipt2. Ucp2 encodes a mitochondrial transporter and anion carrier that promotes homeostasis by serving as a metabolic switch, decreasing TCA cycle function (Vozza et al., 2014) and mitochondrial proton gradient (Ardalan et al., 2022; Berardi and Chou, 2014; Hass and Barnstable, 2021; Keita et al., 2007). Lipt2 plays a similar feedback role, coupling TCA cycle enzyme activity to its input via the mitochondrial fatty acid synthesis pathway (Habarou et al., 2017; Nowinski et al., 2020; Solmonson and DeBerardinis, 2018). Spg7 and Afg3l1 encode subunits of the m-AAA complex, which tailor mitochondrial protein levels to cellular need (Opalińska and Jańska, 2018). Moreover, Ucp2, Spg7 and Afg3l all regulate mitochondrial calcium intake (König et al., 2016; Koshenov et al., 2020; Patron et al., 2018), which couples energy production to synaptic activity by stimulating TCA cycle enzymes (Gherardi et al., 2020; Stoler et al., 2022). As discussed below, this tight regulation may limit oxidative phosphorylation under some conditions, but also reduce reactive oxygen species production and allow for greater biosynthesis.

Our results bolster growing evidence that bioenergetic genes and pathways regulate behaviors like exploratory activity, anxiety, and reward learning. The mitochondrial m-AAA complex, fatty acid oxidation pathway, and fatty acid synthesis feedback pathway are all critical for movement and motor activity in animals and humans (Bernardinelli et al., 2017; Martinelli et al., 2009; Murru et al., 2019; Nsiah-Sefaa and McKenzie, 2016; Patron et al., 2018), with severe, pathogenic mutations in Spg7, Afg3l1, and Ist1 producing hereditary paraplegia and ataxia (Lallemant-Dudek and Durr, 2021; Nsiah-Sefaa and McKenzie, 2016; Schüle and Schöls, 2011), sometimes with altered cognition, executive function, and social/emotional function (Hedera et al., 2002; Lupo et al., 2020; Ringman et al., 2020; Schüle and Schöls, 2011; Zhang et al., 2017). Logically, more subtle changes within these pathways could alter exploratory activity.

Energy production is also theorized to critically modulate the energy-demanding circuitry necessary for behavioral inhibition (Killeen et al., 2013; Russell et al., 2006), and some of our candidate bioenergetic genes are broadly implicated in behavioral temperament. Ucp2 knock-out animals consistently demonstrate bLR-like behaviors, including decreased exploration, and anxiety- and depressive-like behaviors, especially following stress (Andrews et al., 2006; Du et al., 2016; Gimsa et al., 2011; Hermes et al., 2016; Sun et al., 2011; Wang et al., 2014; Yasumoto et al., 2021). Human GWAS also link UCP2, SPG7, and WDR93 to the stress response, psychiatric disorders, externalizing behavior, and substance use disorders (Karlsson Linnér et al., 2021; Li et al., 2022; Orhan et al., 2012; Pehlivan et al., 2020; Rimpelä et al., 2019; Russell et al., 2020; Saunders et al., 2022; Yasuno et al., 2007).

Metabolic differences have been observed in humans and animal models with anxiety and internalizing-like behavior (Ait Tayeb et al., 2023; Filiou et al., 2011, 2014; Filiou and Sandi, 2019; Liu et al., 2022) and hyperactivity and externalizing-like behavior (Chang et al., 2018; Dimatelis et al., 2015; Dupuy et al., 2021; Zametkin et al., 1990). Our findings suggest that genetic vulnerability may contribute to these metabolic differences, bolstering support for metabolic interventions in psychiatry (e.g., Chang et al., 2018; Danan et al., 2022; Filiou and Sandi, 2019; Liu et al., 2022). That said, the evidence linking energy production to internalizing-like vs. externalizing-like behavior is inconsistent across measurements and models, suggesting that the critical vulnerability may lie downstream in bioenergetically regulated functions like apoptosis, oxidative stress, and biogenesis (Filiou and Sandi, 2019). We have evidence supporting each of these possibilities.

During fatty acid oxidation and oxidative phosphorylation, reactive oxygen species are produced as a byproduct (Schönfeld and Reiser, 2021; van Rensburg et al., 2022). Both energy production and reactive oxygen species increase with elevated synaptic activity and environmental stress (Salim, 2017; van Rensburg et al., 2022; Zalachoras et al., 2020). Thus, many metabolic genes are regulators of oxidative stress, with upregulation in bLR-like animals linked to greater oxidative stress and upregulation in bHR-like animals sometimes appearing protective [e.g., Ucp2, Spg7/Afg3l1, Pex11a (Arsenijevic et al., 2000; Atorino et al., 2003; Du et al., 2016; Gimsa et al., 2011; Hass and Barnstable, 2021; Rodríguez-Serrano et al., 2016)]. bHR-like animals also had upregulation of protective Mfge8 (Liu et al., 2014) and upregulation of Nqo2, which can enhance reactive oxygen species production or reduce oxidative stress (Janda et al., 2015; Rashid et al., 2021; Vella et al., 2005) in a manner important for encoding novelty in hippocampal interneurons (Gould et al., 2020) and potentially stress-related disorders (Bainomugisa et al., 2021).

These results bolster evidence that natural and genetically-selected variation in anxiety is consistently associated with markers of oxidative damage in animals and humans (Filiou and Sandi, 2019). Reactive oxygen species are also implicated in the development of anxiety and depressive-like behavior following chronic stress (Schiavone and Trabace, 2016; van Rensburg et al., 2022). The hippocampus is particularly vulnerable to oxidative stress (Salim, 2017) and accumulating evidence implicates oxidative stress in psychiatric disorders, including internalizing disorders and comorbid substance abuse (Bouayed et al., 2009; Cecerska-Heryć et al., 2022; Hovatta et al., 2010; Schiavone and Trabace, 2016; Tobore, 2019; van Rensburg et al., 2022; Zalachoras et al., 2020).

Gene sets related to immune activation and microglia were upregulated in bLR-like animals. This upregulation may be driven by bLR/bHR bioenergetic differences: both the ATP and reactive oxygen species produced by fatty acid oxidation and oxidative phosphorylation can cause microglial activation (Illes et al., 2020; Rojo et al., 2014) and promote microglial release of pro-inflammatory factors (Guevara et al., 2020; Guo et al., 2022). Notably, two of the top candidates upregulated in bHR-like animals, Ucp2 and Mfge8, are also master regulators of microglial activation, promoting an anti-inflammatory and pro-repair state (De Simone et al., 2015; Fang and Zhang, 2021; Gao et al., 2021). Both Mfge8 and Ucp2 encourage microglial engulfment of damaged cells and unwanted synapses. Disrupting this process causes hippocampal dysfunction, inflammation, anxiety-like behavior, insomnia, and depressive-like behavior (Choudhury et al., 2022; Fuller and Van Eldik, 2008; Yasumoto et al., 2021; Yi, 2016), mirroring a bLR-like behavioral phenotype. Fcrl2 was also upregulated in bHR-like animals and is likely abundant in microglia, dampening immune responses (Hammond et al., 2019; Matos et al., 2020). In contrast, Tmem144 was upregulated in bLR-like animals in our study and three others (Blaveri et al., 2010; Meckes et al., 2018; Wilhelm et al., 2013) and is highly expressed in microglia during development (Cao et al., 2020; La Manno et al., 2016), but with unknown function.

These findings complement previous findings that bLR microglia exhibit an “intermediate activation” hyper-ramified morphology (Maras et al., 2022) resembling that observed following chronic stress (McEwen and Akil, 2020), when reactive oxygen species and microglial activation are critical for the development of anxiety-like behavior (Guevara et al., 2020; Lehmann et al., 2019). Moreover, inhibiting microglial activity reduced bLR-like behavior (Maras et al., 2022). Microglial activation has also been implicated in affective and substance use-related behaviors (Choudhury et al., 2022; Northcutt et al., 2015).

Neuroimmune activation could also be caused by mitochondrial regulation of apoptosis and cell death. Excessive mitochondrial calcium intake, decreased m-AAA complex function, decreased Spg7, and decreased Lipt2 can all trigger the mitochondrial membrane potential collapse that drives apoptosis (Bernardinelli et al., 2017; Patron et al., 2018; Shanmughapriya et al., 2015). m-AAA complex deficiencies can also cause dysfunctional mitochondrial protein synthesis, respiration, transport, and fragmentation (Patron et al., 2018) and are linked to neurodegeneration (König et al., 2016; Patron et al., 2018) whereas Ucp2 is considered neuroprotective (Hass and Barnstable, 2016; Kumar et al., 2022). Therefore, down-regulation of Spg7, Afg3l1, Lipt2, and Ucp2 in bLR-like animals might increase risk for cell loss and neuroimmune activation, especially after periods of intense neuronal activity, such as occurs during stress (Zalachoras et al., 2020).

bHR/bLR bioenergetic differences may also contribute to the upregulation of gene sets related to nervous system development and proliferation in bHR-like animals. Both energy availability and use exert control over proliferation, cell differentiation, and growth-related processes, and biosynthesis using glucose-derived products directly competes with oxidative phosphorylation for essential substrates (Beckervordersandforth, 2017). Therefore, many of the candidate metabolic genes also influence proliferation and growth (e.g., Ucp2, Lsr, Lipt2: Esteves et al., 2014; Pecqueur et al., 2008; Takahashi et al., 2021; Wang et al., 2023; Zhang and Ma, 2021). Other top candidates regulate growth-related processes, including Mfge8 and Fzd6 (Wang et al., 2012; Yli-Karjanmaa et al., 2019; Zhou et al., 2018). Fzd6 has also been linked to anxiety and depressive-like behavior (Sani et al., 2012; Voleti et al., 2012). These results are noteworthy due to known bHR/bLR differences in neurogenesis, proliferation, and growth factor response (Birt et al., 2021; Perez et al., 2009; Turner et al., 2019), and extensive literature implicating both hippocampal atrophy in internalizing disorders and growth-related processes in antidepressant function (Duman and Monteggia, 2006).

A vast literature exists documenting sex differences in anxiety and reward-related behaviors and circuitry (Bangasser and Cuarenta, 2021; Becker and Koob, 2016). Our study was not specifically designed to study these differences, but sex differences were observed for several behaviors, some of which have appeared consistently over many generations of bHR/bLR rats (anxiety-like behavior: Hebda-Bauer et al., 2017) or in previous studies. For example, more females than males demonstrated sign-tracking vs. goal-tracking PavCA behavior (Hughson et al., 2019; Pitchers et al., 2015), reflecting sex differences in reward processing (Davis et al., 2008; Flagel et al., 2011). Other sex differences may represent batch effects (e.g., LocoScore), as males and females were by necessity tested separately.

Our hippocampal differential expression analyses could shed light on these sex differences. Beyond the x- and y-chromosomes, the effect sizes for sex differences in hippocampal expression tended to be small, making them difficult to measure in the F₀ sample, but in the larger F₂ sample 1,679 genes had sex differences that survived false discovery rate correction (FDR < 0.10). We did not find evidence that sex modulated the relationship between gene expression and bHR/bLR phenotype or F₂ behavior, but there was a notable overlap (10%) of the genes with sex differences in expression in either the F₀ or F₂ datasets with genes with bHR/bLR differential expression (Supplementary Table S4). In future studies, these genes may be excellent candidates for mediating sex differences in behavior.

We focused on a single brain region because of the need to generate a large sample size (n = 250) from a heterogeneous F2 population, however, each of the measured behaviors depends on the activity of broader brain circuitry. Can our hippocampal results provide insight into the functioning of other brain regions? To address this question, we ran an exploratory analysis comparing bHR/bLR hippocampal differential expression to the pattern of differential expression in other brain regions in previous transcriptional profiling datasets from bHR/bLR adults, including the amygdala (Cohen et al., 2015, 2017; McCoy et al., 2017), dorsal raphe (Cohen et al., 2017), and unpublished data from the cortex and hypothalamus. These comparisons suggest that at least some of the bHR/bLR differential expression that we identified in the hippocampus may also be present in other brain regions, whereas other differential expression may be hippocampal specific. Our ongoing studies using spatial transcriptomics and fluorescent in situ hybridization (FISH) should provide further insight into bHR/bLR differential expression in other brain regions, as well as illuminate the specificity of our findings to particular hippocampal cell types and subregions. We also have ongoing work characterizing bHR/bLR gene expression and chromatin accessibility in the nucleus accumbens, another region noted for its role in sensation-seeking, reward processing, and addiction.

The behaviors quantified in our F₂ rats did not encompass all hippocampal-dependent behaviors that differ in the two lines. Notably, bHR and bLR rats show differences in contextual fear conditioning (Prater et al., 2017; Widman et al., 2019) that could be influenced by many of the pathways implicated in our results, including mitochondrial function, oxidative stress, microglial function, and neurogenesis (Gao et al., 2018; Olsen et al., 2013; Villasana et al., 2016; Yu et al., 2022). Future work following up on these findings could provide insight into the heritable contributions underlying internalizing disorders like post-traumatic stress disorder (Banerjee et al., 2017).

One limitation of our approach is that loci which influence behavior may still be part of haplotypes that include multiple eQTLs, despite the added resolution provided by a F₀-F₁-F₂ cross. We have addressed this limitation by integrating individual genes into higher order biological concepts. This approach improves the translatability of our results and should be robust to the presence of some false positives. That said, we may be missing genetic variation contributing to our phenotype by only focusing on single nucleotide variants that could mediate effects on behavior via basal hippocampal gene expression levels. For example, one of the most compelling candidates that we identified was AABR07071904.1, with a cis-eQTL near the strongest LocoScore QTL peak (Chitre et al., 2023). According to genome assembly Rnor6 (Ensembl v103), AABR07071904.1 generates long non-coding RNA, but in mRatBN7.2 (Ensembl v106) the gene was retired, potentially mapping to Zfp939-201. In either form, it could play some important regulatory role, but it is noteworthy that the implicated cis-eQTL is in linkage disequilibrium with a missense coding variant for Plekhf1 (Supplementary Figure S13, Chitre et al., 2023). Plekhf1 was not differentially expressed in our study, but has been linked to stress and mood (Chitre et al., 2023). Future work will explore other mechanisms that may contribute to our phenotype, including coding variants, structural variants, epigenetic modifications, epistatic interactions, and context-dependent activity.