Authors: Patrick B. Chen, Rachel Chen, Nathan LaPierre, Zeyuan Chen, Joel Mefford, Emilie Marcus, Matthew G. Heffel, Daniela C. Soto, Jason Ernst, Chongyuan Luo, Jonathan Flint
Categories: Article, QTL mapping, fear conditioning, inbred mouse strains, quantitative complementation, single-nucleus ATAC-seq, single-nucleus RNA-seq
Source: Cell Genomics
Knowing the genes involved in quantitative traits provides an entry point to understanding the biological bases of behavior, but there are very few examples where the pathway from genetic locus to behavioral change is known. To explore the role of specific genes in fear behavior, we mapped three fear-related traits, tested fourteen genes at six quantitative trait loci (QTLs) by quantitative complementation, and identified six genes. Four genes, Lamp, Ptprd, Nptx2, and Sh3gl, have known roles in synapse function; the fifth, Psip1, was not previously implicated in behavior; and the sixth is a long non-coding RNA, 4933413L06Rik, of unknown function. Variation in transcriptome and epigenetic modalities occurred preferentially in excitatory neurons, suggesting that genetic variation is more permissible in excitatory than inhibitory neuronal circuits. Our results relieve a bottleneck in using genetic mapping of QTLs to uncover biology underlying behavior and prompt a reconsideration of expected relationships between genetic and functional variation.
Keywords: inbred mouse strains, fear conditioning, QTL mapping, quantitative complementation, single-nucleus RNA-seq, single-nucleus ATAC-seq
A major challenge in behavior genetics is to turn genetic information into mechanistic understanding of the sort that would, for example, be useful in designing new treatments for psychiatric disorders and, more generally, understanding how genetic variation leads to behavioral variation. A key step is progressing from quantitative trait locus (QTL) to gene, which has only been achieved for any complex trait, in any species, in a small number of cases.^1^^,^^2^ Of the 5,000 QTLs identified in rodents, less than 100 genes have been identified, almost all on the basis of correlative evidence, such as proximity of a gene to the QTL or alterations in transcript abundance, rather than by a causal test of a gene’s candidacy.^2^ Currently, there is no consensus on how to proceed from QTL to gene.
Here, we demonstrate the power of a quantitative complementation (QC) test, first applied in Drosophila^3^ and later shown to work in rodents,^4^ to directly query the causal gene impacted by the QTL. Construction and phenotyping of F1 hybrids with and without a knockout (KO) of a candidate gene, and of inbred strains with and without the KO of a candidate gene, test whether the QTL operates through the gene under investigation. By separately assaying the joint effect of QTLs (from the phenotypic difference between strains) and the effect of the mutation (from the phenotypic difference between KOs and wild type [WT]), the QC test reveals a QTL whose effect depends on the presence of the candidate gene as a significant interaction between the effect of mutation and the effect of strain. Applying the test requires access to inbred animals carrying a KO on the same genetic background as the QTL, which has been difficult to achieve when KOs were generated using homologous recombination in 129 strains or C57BL/6N. The development of the CRISPR-Cas9 technology now lifts that restriction by enabling genetic manipulation in any strain.^5^^,^^6^
We set out to apply the QC test to fear-related behaviors. We chose these traits because we could use extensive information available about the brain regions involved (ventral hippocampus and amygdala^7^^,^^8^), and their underlying circuitry,^9^^,^^10^^,^^11^^,^^12^ to explore where and how in the brain genetic variation results in behavioral change. We reasoned that by knowing the causal genes, we could identify changes in their expression and regulation at a single-cell level, revealing the molecular mechanisms that mediate the impact of genetic variants on traits.
Our first task was to identify QTLs involved in fear-related behaviors that ideally indicated one or a small number of genes for QC testing. We chose a hybrid mouse diversity panel (HMDP) for mapping because of its potential to deliver high resolution^2^ (thus increasing the chance of finding QTLs containing a small number of candidate genes) and because it includes recombinant inbreds derived from C57BL/6J, hence providing many QTL alleles on a strain readily amenable to CRISPR-Cas9 modification (the HMDP strains and numbers of animals used are listed in Table S1). We mapped variation in a conditioned fear assay in which animals were exposed to an auditory cue associated with an aversive shock. Re-exposure to the auditory cue elicits a fear response, measured as the amount of time spent freezing^13^^,^^14^ and referred to in this paper as FC-cue. The context in which conditioning occurs also elicits freezing (FC-context),^15^ and we assayed this too. We also mapped variation in an unconditioned fear assay using the elevated plus maze (EPM), whose open arms act as anxiogenic environment. We counted the number of entries into the EPM’s open arms as a measure of fear, referred to henceforth as EPM-open.
To enable us to perform a well-powered meta-analysis and robustly identify loci at high resolution, we combined our data with results from 29 publications.^16^^,^^17^^,^^18^^,^^19^^,^^20^^,^^21^^,^^22^^,^^23^^,^^24^^,^^25^^,^^26^^,^^27^^,^^28^^,^^29^^,^^30^^,^^31^^,^^32^^,^^33^^,^^34^^,^^35^^,^^36^^,^^37^^,^^38^^,^^39^^,^^40^^,^^41^^,^^42^^,^^43^^,^^44^^,^^45^ We implemented a rigorous pipeline to harmonize phenotypes across the multiple datasets (STAR Methods: harmonizing phenotyping data for meta-analysis), generating phenotypes on a total of 6,544 mice from seven cohorts (Tables S1 and S2).
Our next requirement was a set of genotypes at markers sufficiently dense to capture most causative variants. Short-read sequencing provides almost complete catalogs of genetic variants for sixteen inbred strains,^46^ which are also a resource for imputing variants in the other mapping strains.^47^ We carried out imputation using a hidden Markov model-based technique with proven effectiveness in inbred model organisms.^48^ After imputation, we obtained genotypes on 16,767,664 markers (the imputation strategy and validation results are described in STAR Methods: imputation of genotypes). For mapping, marker sets were chosen that are polymorphic within the strains of each dataset, a number varying from 4 to 16 million SNPs segregating between the inbred strains, sufficient to identify causative variants at QTLs.^49^^,^^50^
The phenotypic dataset we obtained is highly structured, consisting of closely related individuals (recombinant inbreds share half their genome, like full siblings) as well as distantly related individuals. Genome-wide association was carried out with GEMMA, which implements a linear mixed-model marker association test,^51^^,^^52^ including a genomic relationship matrix to control for the population structure arising from the inclusion of individuals with different degrees of relatedness. We mapped each phenotype in each cohort and carried out a meta-analysis of the results using a random effects model implemented in METASOFT.^53^^,^^54^ We used Genome Reference Consortium Mouse Build 38 (GCA_000001635.2, mm10), and all coordinates in this paper refer to that build.
To assess significance, we used GCTA^55^ to simulate data for each cohort and meta-analyzed the results in the same way as we did for the real datasets. From 100,000 simulations, we obtained a 5% threshold of negative logarithm (base 10) of the association p value (logp) of 4.50, consistent with thresholds estimated by others using the same experimental design.^35^^,^^36^^,^^37^^,^^38^^,^^39^^,^^40^^,^^41^^,^^42^^,^^43^^,^^44^^,^^45^^,^^46^^,^^47^^,^^48^^,^^49^^,^^50^^,^^51^^,^^52^^,^^53^^,^^54^^,^^55^^,^^56^ Applying this significance threshold, we identified 14 loci for FC-context, 15 for FC-cue, and 72 for EPM-open (Table S3). At a threshold corrected for testing three phenotypes (logp = 4.97) the number of loci for EPM-open fell to 60 while remaining unchanged for the other two phenotypes. Four loci were common to FC-context and FC-cue and four were common to EPM-open and FC-cue, giving 93 unique loci. The Manhattan plots in Figure 1 show all loci that exceeded the logp threshold of 4.5.
Figure 1 Manhattan plots for three fear-related behaviors in mice(A–C) The data show the results of a meta-analysis of (A) 1,931 mice for cue conditioning (FC-cue), (B) 2,671 for contextual conditioning (FC-context), and (C) 1,942 for entries into the open arms of the elevated plus maze (EPM-open arms). Chromosome numbers are listed on the horizontal axis. The vertical scale is the negative logarithm (base 10) of the association p value. A horizontal dotted red line gives the location of the 5% significance threshold (from simulation).
While the genetic architecture of fear-related behaviors was highly polygenic, we found evidence for a few large-effect loci, notably on chromosomes 7 and 13 for FC-cue (Figure 1A) and chromosomes 4 and 13 for FC-context (Figure 1B). More than five times as many QTLs were identified underlying EPM-open than either of the fear-conditioning phenotypes, likely reflecting higher polygenicity of the trait (Figure 1C). The size of the QTL intervals varied considerably from 0.02 to 19.6 Mb (confidence intervals estimated by simulation^57^) with a median of 2.5 Mb, and each locus contained a median of 23 genes (range from 0 to 329; Table S3).
Which of the genes at a locus is causal for the phenotype? To answer this for specific loci, we first chose a 2 Mb QTL at the end of chromosome 13 associated with FC-cue (logp = 7.6), previously identified in a panel of BXD recombinant inbreds,^27^ that includes a hyperpolarization-activated cyclic nucleotide-gated channel 1 (Hcn1). Since pharmacological blockade of HCN1 reduces freezing, Hcn1 has been proposed as the causal gene at this locus^27^ (although deletion of Hcn1 does not alter conditioned fear^58^). We used an exon-excision CRISPR-Cas9 strategy to knock out out all four annotated and one unannotated gene lying within the 95% confidence intervals of the QTL (Figure 2A) (we omitted the single-exon unannotated transcript Gm6416). We used a quantitative polymerase chain reaction to confirm that the engineered mutations altered RNA abundance. The CRISPR-Cas9 strategy and characterization of the mutants is described in STAR Methods and Tables S4 and S5.
Figure 2 Quantitative complementation of five genes at a locus on chromosome 13 for cue fear conditioning(A) QTL regional information from LocusZoom.^59^ The top part shows the association results from a meta-analysis, with the position of the highest-scoring variant annotated in purple. The vertical scale is the negative logarithm (base 10) of the association p value. The bottom section gives the location and orientation of genes at the locus. A red box identifies each gene used in the quantitative complementation tests.(B) Design of the quantitative complementation test. Black boxes indicate the knockout (KO) allele, and strain is indicated by color, where blue is C57BL/6J and red is DBA/2J. An “X” indicates a cross between the named groups above and below the X, and an arrow points to the progeny of each cross. The four groups used in the quantitative complementation test are wild-type C57BL/6J (B6 WT), heterozygote KOs (B6 KOs), an F1 from crossing DBA/2J to C57BL/6J (D2 WT), and the F1 KO by DBA/2J (D2 KO).(C) Results of quantitative complementation testing of four annotated genes in the region. The vertical axis is the quantile normalized duration of freezing to a cue from the fear-conditioning test. The horizontal axis lists the groups of mices used. Groups are B6.WT: WT C57BL/6J; B6.KO: heterozygote KOs on C57BL/6J; D2.WT: F1 from crossing DBA/2J to C57BL/6J; D2.KO: the F1 from crossing the KO onto DBA/2J. These genes were not significant by quantitative complementation (Table 1).(D) Results of the quantitative complementation test for the long non-coding RNA 4933413L06Rik. The axes labels are the same as in (C). The difference between these groups yielded a significant interaction result in the quantitative complementation test (Table 1).
Each KO was subject to a QC test (Figure 2B) to find out which gene was mediating the effect of the locus on FC-cue. Four genes tested negative (Figure 2C), while one was positive (Figure 2D), which to our surprise was an unannotated long non-coding RNA (lncRNA), 4933413L06Rik. To interpret Figures 2C and 2D, consider that the first two groups in Figure 2B (C57BL/6J wild type [B6.WT] and DBA/2J wild type [DBA.WT]) measure a strain effect arising from a single copy of the D2 genome in the F1 animals (assuming strict additivity, the strain effect should be about half that found in inbred strain comparisons, consistent with our results). Any difference between the B6.WT and B6.KO groups is attributable to the presence of the KO.
The pattern of results in Figure 2C for the four genes tested looks a strain effect can be seen but no discernable difference between groups carrying the KO and the WTs. By contrast, Figure 2D shows a different pattern. Animals with the KO do differ from their respective WTs, indicating that the KO influences the phenotype but with differences between the two strains. The B6.KO animals spend less time freezing than their WT siblings (B6.WT), while the DBA.KO group freezes more than their WT siblings (DBA.WT). In other words, the effect of the mutation depends on the strain background, which we assume to be due to a nearby QTL.
We tested this relationship between strain and KO by analysis of variance and found a significant interaction (p = 0.0015; Table 1), indicating that the effect of the QTL is mediated by the unannotated gene 4933413L06Rik. We found no evidence of the involvement of Hcn1, demonstrating that QC testing can unambiguously identify genes mediating the effects of QTLs.
We next tackled the most significant region of association for FC-context, a 15-Mb region on chromosome 4 (between 70 and 85). There were two QTLs here (70–78 and 79–85 Mb; Figure 3), together containing 46 genes. To narrow the choice of genes, we tested genes close to the most highly associated markers and with confirmed gene expression in the hippocampus and amygdala (from published sources^60^^,^^61^^,^^62^^,^^63^^,^^64^^,^^65^^,^^66^^,^^67^). At the first locus, the interaction p value of Ptprd exceeded a multiple correction testing threshold (p < 0.0036 for testing 14 genes listed in Table 1), while at the second two genes, Psip1 and Sh3gl3, exceeded the threshold (results shown in Figure 3 and Table 1).
Figure 3 Quantitative complementation of four genes at two QTLs on chromosome 4 for contextual fear conditioning(A and C) QTL regional information from LocusZoom.^59^ The top part shows the association results from a meta-analysis, with the position of the highest-scoring variant annotated in purple. The vertical scale is the negative logarithm (base 10) of the association p value. The bottom section gives the location and orientation of genes at the locus. A red box identifies each gene used in the quantitative complementation tests.(B and D) Results of quantitative complementation testing of six annotated genes. The vertical axis is the quantile normalized duration of freezing to the context from the fear-conditioning test. The horizontal axis lists the groups of mice used. Groups are B6.WT: WT C57BL/6J; B6.KO: heterozygote KOs on C57BL/6J; D2.WT: F1 from crossing DBA/2J to C57BL/6J; D2.KO: the F1 from crossing the KO onto DBA/2J. Genes significant by quantitative complementation are shown in gray (results in Table 1).
Finally, we chose three loci contributing to variation in EPM-open behavior (Figure 4). The first, on chromosome 16, had a highly significant association (logp = 9.48) and lies in an intron of a cell adhesion molecule, Lsamp (limbic system-associated membrane protein), whose deletion is known to alter behavior in the EPM.^68^^,^^69^ We found a significant interaction (p = 0.001; Table 1), indicating that the effect of the QTL is mediated by the Lsamp gene. We tested Stim2 and Nptx2 at two other loci using the same CRISPR-Cas9 exon deletion approach. While the QC test confirmed the candidacy of Nptx2, it was negative for Stim2. Table 1 summarizes the results of the QC tests for all genes at the six loci we examined.
Figure 4 Quantitative complementation of three genes at three QTLs for the number of entries into the open arms of the EPM(A, C, and E) QTL regional information from LocusZoom.^59^ The top part shows the association results from a meta-analysis, with the position of the highest-scoring variant annotated in purple. The vertical scale is the negative logarithm (base 10) of the association p value. The bottom section gives the location and orientation of genes at the locus. A red box identifies each gene used in the quantitative complementation tests.(B, D, and F) Results of quantitative complementation testing. The vertical axis is the quantile normalized duration of the number of entries into the open arms of the elevated plus maze. The horizontal axis lists the groups of mice used. Groups are B6.WT: WT C57BL/6J; B6.KO: heterozygote KOs on C57BL/6J; D2.WT: F1 from crossing DBA/2J to C57BL/6J; D2.KO: the F1 from crossing the KO onto DBA/2J.
One hypothesis about how QTLs act is that they alter gene expression, which in our case would manifest as a difference in the abundance of transcripts from C57BL/6J and DBA/2J genomes in relevant cell types. With six genes in hand, we set out to test this assumption. We generated single-nucleus RNA-sequencing (snRNA-seq) datasets in C57BL/6J and DBA/2J animals from the ventral hippocampus and amygdala. We identified 17 cell types in the ventral hippocampus and 22 in the amygdala with two replicates per strain per region (Figure 5, joint uniform manifold approximation and projection). The clusters included excitatory (glutamatergic) and inhibitory (GABAergic) neurons, as well as several types of non-neuronal cells that were annotated based on known marker genes. We define individual neuronal cell types based on neurotransmitter expression and region or top marker gene (for example, Exc-DG in the ventral hippocampus refers to excitatory cells in the dentate gyrus, and Exc-Tafa1 in the amygdala refers to excitatory cells marked by the gene Tafa1), and we define neuronal cell class as either excitatory or inhibitory through combining all excitatory or all inhibitory cell types, respectively (one amygdala neuronal type was excluded from both cell classes due to the expression of both excitatory and inhibitory markers). Non-neuronal cell types were not further defined due to an enrichment step for neuronal nuclei during sample processing limiting the number of non-neuronal nuclei captured.
Figure 5 Cell-type uniform manifold approximation and projection for ventral hippocampus and amygdala in B6 and DBA from single-nucleus RNA, ATAC, and methylation sequencing. In all panels, vertical and horizontal axes represent UMAP space.(A–C) Ventral hippocampal clusters plotted by (A) modality, (B) strain, and (C) cell-type identity.(D–F) Amygdala neuronal clusters plotted by (D) modality, (E) strain, and (F) cell-type identity.
Given that QTL action is believed to have its effect on transcription through the alteration of the sequence at regulatory regions, we also generated single-nucleus data for CpG methylation^70^ and from a single-nucleus assay of transposase-accessible chromatin (snATAC-seq).^71^ Both epigenetic modalities can be used to identify regulatory elements, such as promoters and enhancers.^72^ Cell-type clusters generated from snRNA profiles were propagated to snATAC and methylation datasets through joint embedding (Figure 5). We confirmed that cell types were integrated by modality (Figures 5A and 5D), that there was no obvious strain bias for cell types (Figures 5B and 5E), and that cell types were well defined (Figures 5C and 5F). Importantly for interpreting our analyses, the three epigenetic datasets were obtained independently from different animals, with two independent replicates per modality.
Five out of six genes showed significant differences in expression between C57BL/6J and DBA/2J (significance defined as exceeding a 10% p value [logp > 1.0], adjusted for the total number of transcripts tested, from DEseq2^73^ output; negative log2 fold change [log2FC] indicates higher expression in B6): 4933413L06Rik (in eleven cell types; median log2FC = −0.74, range = −1.63 to −0.53), Lsamp (in eleven cell types; median log2FCc = −0.21, range = −0.78 to 1.1), Psip1 (in five cell types; median log2FC = 0.41, range = 0.30 to 0.49), Nptx2 (in one cell type; log2FC = 1.82), and Ptprd (in eight cell types; median log2FC = −0.31, range = −1.24. to 1.85). Of the 36 cell types exhibiting differential expression for these genes, 31 were within the hippocampus. Different hippocampal cell types showed different degrees of overlap for differential gene expression of genes; for example, four of the five genes were differentially expressed within Exc-CA1-Galntl6 cells, while only 4933413L06Rik was differentially expressed within Exc-CA1-Camk2d cells. We noted that 22 of the 31 significant differences in the hippocampus and three out of five significant differences in the amygdala occurred in cell types categorized as excitatory neurons. Results of the DESeq2 analyses are given in Table S6.
We asked if the pattern of gene expression also held for snATAC-seq and methylation profiles, namely showing more strain differences in the hippocampus with an enrichment in excitatory neurons. To do this, we made an inventory of variable snATAC-seq and methylation sites. We counted differential ATAC sites (using a 5% threshold from adjusted p values derived from DESeq2^73^) that contained sequence variants under the assumption that differences at these sites were more likely genetic in origin than those without such differences (more than 90% of the significantly different ATAC sites contain a sequence variant compared to 35% of non-variable sites). For the methylation data, we counted the number of fully methylated CpG sites in each cell type and divided them into those with and without mutations that disrupted methylation.
Figure 6 plots differential methylation and ATAC sites in the hippocampus and amygdala at the five QTLs where we identified causal genes. Three general observations arose from the data presented in Figure 6. First, there were more strain differences in the hippocampus than amygdala, as reflected in the greater number of dots compared to crosses. Second, at each locus, the strain differences were found primarily in excitatory neurons, in both the hippocampal and amygdala, as shown by the preponderance of red. Third, the pattern of strain differences for each modality varied between loci and genes, providing no unequivocal location or molecular signature to identify causal variants. The Ptprd locus (Figure 6A) contained a concentration of methylation, ATAC sites, and sequence differences at the 3′ end of the gene, suggesting a location for one or more causal variants here. At the Psip1 and Sh3gl2 locus (Figure 6B), strain differences were spread across the region with no indication of which might be relevant. The QTLs containing 4933413L06Rik (Figure 6C) and Nptx2 (Figure 6D) included a concentration of strain differences in methylation sites around the genes but with no obvious candidates for causative variants. One possible exception was the single ATAC site that exceeded a corrected significance threshold (logp = 6) in an intron of the Lsamp gene (Figure 6E) coinciding with a 2-kb deletion in DBA/2J. The function of this site and the consequences of its deletion are not known.
Figure 6 ATAC-seq and methylation data at five QTLs in the ventral hippocampus and amygdala(A–E) Data from the hippocampus are marked by dots (⋅) and data from the amygdala by crosses (+). In each image, the lowest row shows the distribution of genetic variants (black). The next two tracks, immediately below the genes and indicated by the blue background, show the location of methylation sites that coincide with a sequence variant and are thus present in only one of the two strains. These sites are colored according to the class of cell type in which they red for excitatory neurons, blue for inhibitory neurons, and gray for non-neuronal cell types. Above the methylation tracks are shown ATAC sites, which differ significantly between the strains. The vertical axis shows the absolute value of log2 fold change for these ATAC sites. Sites are again colored according to the class of cell type in which they occur using the same color scheme as for methylation. The vertical positioning of the genes is only for ease of plotting and carries no meaning.(F) Comparison of the percentage of significantly different RNA species in QTL regions and in the genome for the hippocampus. The vertical axis is the percentage of transcripts that are significantly different between the two strains, and the horizontal axis lists the cell types. Data from QTL intervals are indicated by dots colored according to excitatory (red), inhibitory (blue), and non-neuronal (gray). Data from the entire genome are indicated by black crosses.
The apparent enrichment of genetically mediated variation in excitatory neurons at the five QTLs containing the causal genes led us to examine whether the same was true for all the QTLs we had identified. We compared gene expression and epigenetic variation at 93 QTLs with the rest of the genome, expecting that any enrichment would be confined to, or more prominent at, the QTLs. Surprisingly, although the QTLs occupy about 0.01% of the genome, the percentage of significantly different transcripts (again defined as exceeding a 5% threshold from adjusted p values derived from DESeq2^73^) in QTL regions was indistinguishable from the percentage of significantly different transcripts in the genome (Figure 6F), a finding replicated for methylation and snATAC-seq data (Tables S7 and S8). In other words, there was no enrichment in QTLs for genetically mediated variation in excitatory neurons; all regions of the genome contribute to this phenomenon.
Is the enrichment of genetic variants in excitatory neurons significant? It certainly seems so, but at least two confounds could produce this result. For the snATAC-seq and snRNA-seq data, identification of significant strain differences is more likely with higher read coverage, and for methylation, the more sites we identify, the higher the likelihood that some will coincide with a sequence variant. To assess the contribution of sequence coverage, we compared a null model, in which sequence coverage predicted the number of differentially expressed genes, ATAC peaks, or methylation sites with mutations, to one that additionally included cell type and found highly significant improvements in fit (all p values < 2.2E−16; Table S9), demonstrating that cell type was an important contributor, over and above the contribution of sequence coverage. A generalized linear model to predict differences also showed that sequence coverage does not account for the differences between excitatory and inhibitory neurons (Table S10).
Second, the inclusion of more excitatory cell types than inhibitory types could bias findings (eight excitatory cell types compared to three inhibitory in the hippocampus) because the addition of extra cell types increases the chances of finding variant sites. Against this, we note that in the amygdala, there were more inhibitory than excitatory cell types (eight versus six). We ran a test of this hypothesis, by comparing the contribution of class (excitatory and inhibitory) to cell type, considering that cell type was completely nested within class. We resorted to a Bayesian analysis, which revealed a significant effect of class (the 95% confidence intervals do not overlap zero) (Table S11).
Is the enrichment in excitatory cells a peculiarity of the difference between C57BL/6J and DBA/2J? We answered this question for the hippocampus using methylation data from six other strains. Using CAST/EiJ as an outgroup, we compared the number of mutations in CpG sites at excitatory and inhibitory neurons in the ventral hippocampus. Table S12 shows that for each strain, compared to CAST/EiJ, the proportion of mutated CpG sites was higher in excitatory than inhibitory neurons, supporting our finding that at least for methylation, alterations in functional elements directly attributable to genetic variation were enriched in excitatory neurons compared to other cell types. Together, these results suggest a model whereby genetic effects act on specific hubs of circuits/cell types.
We have established the feasibility of combining mapping with QC testing for gene identification, breaking open a bottleneck in the genetic dissection of complex traits in mice, and provided a set of six genes for further mechanistic understanding of fear-related behavior. Unexpectedly, genetically mediated variation in three independent molecular analyses, methylation, snRNA-seq, and snATAC-seq, was found to occur preferentially in excitatory neurons and to be enriched in the hippocampus compared to the amygdala. These observations demonstrate the possibility of moving quickly from locus to gene, identify unexplored biology underlying fear-related behavior, and raise questions about expected relationships between genetic and functional variation.
Our findings highlight the power of unbiased genetic examination of the biological basis of behavioral variation. At a locus on chromosome 13, we made the unexpected discovery that an lncRNA (4933413L06Rik) is involved in cue fear conditioning. Little is known about 4933413L06Rik (one study found it to be differentially expressed during development of the auditory forebrain^74^). Our finding opens the possibility for clarifying the function of this lncRNA in fear-related behavior. Two other causal genes (Psip1, Sh3gl2) have also never been studied in the context of fear-related behaviors. We note that of the five genes with known functions, four (Lsamp, Ptprd, Nptx2, and Sh3gl2) are implicated in synapse development or function.^75^^,^^76^^,^^77^^,^^78^^,^^79^ Ptprd belongs to the type IIA receptor-type protein tyrosine phosphatase family of phosphatases, involved in synaptic formation and structure during development and neurogenesis.^76^ Lsamp is a cell surface adhesion molecule that is involved in neuritogenesis and axon guidance (and is also implicated in the maturation of serotonergic^80^ and thalamocortical circuits^81^). Nptx2 is a secreted glycoprotein that localizes to both pre- and post-synaptic compartments of excitatory synapses containing AMPA receptors and is involved with synaptogenesis^82^ and adult neurogenesis.^83^ Sh3gl2 (endophilin1), a cytoplasmic Src homology 3 domain-containing protein,^84^ localizes to presynaptic nerve terminals, where it functions in synaptic vesicle endocytosis^85^ and the regulation of exocytosis.^86^ The significance of these functional observations is difficult to assess, but it at least suggests that variation in fear-related behavior in part originates in synaptic development.
The QC test could serve as a gold standard for gene identification following QTL mapping in inbred mice. Following mapping, the construction of KOs is relatively straightforward, and the test requires a simple breeding protocol. We recommend powering the QC tests based on results from the mapping experiment (the effect size of the locus in the mapping population will be lower than that in QC test^87^), but note that our failure to detect a candidate gene at a locus on chromosome 5 could be due to low power or because we tested the wrong gene.
Our instantiation of the QC test is not the only one; alternative designs are possible. For example, in an experiment to test the candidacy of a gene at a QTL for resistance to obesity in mice, a QC test was constructed from a KO, a subcongenic containing the QTL, and their corresponding background strains (a cross giving rise to four genotypes on a uniform genetic background).^88^ However, it is essential to use KOs on an isogenic background so as to avoid the complication of the effect of additional QTLs that may be present on the non-isogenic DNA.
Interpretation of the QC results is subject to some qualifications. First, while a positive QC test result reveals a gene through which a QTL has its effect, it does not necessarily mean that the QTL is the one containing the candidate gene. The QC test will detect the effect of any QTL that requires the gene to be present. Typically, that QTL will be the one local to the gene tested, but that assumption is not essential for the functioning of the QC test. Second, the QC test makes no assumptions about how the QTL and KO interact. Does this matter? It may, because in the QC test, non-additive interactions are typically easier to detect than additive ones. Under an additive model, the KO removes the effect of one of the QTLs so that the phenotypic effect is reduced by the allelic effect of that QTL. In non-additive models, the presence of the KO could reduce one QTL effect to zero (recessive effect) or increase it, potentially in a multiplicative fashion. In these cases, the effect of the interaction is to augment the allelic effect over that of the additive interaction. We have some evidence that non-additive interactions predominate in our results.
QC tests results from four genes (Lsamp, Ptprd, 4933413L06Rik, and Psip1) are consistent with a non-additive interaction between QTL and KO. In each case, the phenotype of the KO on the hybrid background is larger than either the WT DBA/2J or the KO on the C57BL/6J (Figures 2, 3, and 4). One interpretation is that non-additivity is common. Another interpretation is that we failed to detect additive interactions because non-additivity is easier to the larger effect attributable to the interaction confers greater power for the QC test. This may contribute to false negative findings, but it does not undermine the conclusion that the effect of the QTL is mediated by the gene. Under both additive and non-additive models, a significant result means that the effect of the QTL requires the gene.
One tantalizing corollary emerging from the analysis of gene causality in fear-related behaviors is that genetic variation may act preferentially in excitatory, rather than inhibitory, neurons. We note that any patterns of strain differences could not be attributed solely to coverage, as even downsampling our RNA-seq and ATAC-seq data to as low as 25% of the total unique molecular identifiers still resulted in substantial detection of significant strain differences (see https://dx.doi.org/10.6084/m9.figshare.25521304). We had expected that the pattern of strain differences might point to the involvement of a particular cell type in one region of the brain. Instead of finding cell-type enrichment, we observed that strain differences in three different modalities, snRNA-seq, methylation, and snATAC-seq, were more prevalent in excitatory than inhibitory neurons, particularly in the hippocampus compared to the amygdala. Moreover, we found this enrichment pattern to be true not just for the QTL regions but for the entire genome. There was little evidence of strain differences in non-neuronal tissues, although we can be less confident in this assertion because our sample was enriched for neuronal cells and may not be representative of non-neuronal variation.
What might explain this finding? One possibility is that it is specific to a comparison between C57BL/6J and DBA/2J. We think that this is unlikely and provide some evidence against this idea by running a similar analysis for methylation data in other strain comparisons, though we have not extended this observation to other modalities. Alternatively, the finding represents a preference for genetic effects in excitatory neurons. Why might this be so?
The idea that genetic effects operate preferentially in excitatory neurons fits with a neuronal circuit model where inhibitory neurons sculpt behaviors driven by excitatory neurons. Excitatory neurons are extremely diverse, resulting in high-dimensional and sparse coding,^89^ consistent with the expectation that such an architecture permits highly efficient information transfer.^90^ Consequently, encoding within excitatory networks cannot be predicted from their transcriptome alone. In contrast, activity within inhibitory neuronal networks is lower dimensional, with correlations between inhibitory cell types determined primarily by cell type, so that cells belonging to a single transcriptomic class are highly correlated.^91^^,^^92^^,^^93^ In this model, genetic variation is more permissible in excitatory neuronal networks than inhibitory, suggesting that disturbances of inhibitory networks are more likely to be pathological. If this interpretation is correct, then identifying genetic effects that disturb function at a circuit level may be more fruitful than looking for cell-type-specific genetic effects.
Our results are subject to some limitations. First, the results from the meta-analysis have not been replicated, so it is possible that some of the loci may be false positive. This may explain why the QC test was negative for Stim2 on chromosome 5. Second, some of the negative QC test results may be due to relatively low power to detect additive effects compared to non-additive effects. Third, we did not test every expressed transcript at every QTL. For instance, it is possible that single-exon unannotated transcripts (e.g., Gm6416, which we did not test at the chromosome 13 locus) might also contribute. However, the failure to detect the involvement of additional genes does not mitigate the main result of our experiment, namely, that the QC test can be used in a systematic fashion to identify genes mediating the QTL effect.
Further information and requests for reagents and resources should be directed to and will be fulfilled by the lead contact, Jonathan Flint (JFlint@mednet.ucla.edu).
All mutant KO animals are deposited as frozen sperm at the Transgenic and Genome Editing Core Facility at Augusta University. Ordering information can be found at https://doi.org/10.6084/m9.figshare.25521283.
All experimental procedures using live animals were approved by UCLA’s Animal Care and Use Committee (protocol number ARC-2018-026). Male mice from eight inbred strains A/J, C57BL/6J, BALB/cJ, FVB/J, DBA/2J, WSB/EiJ, PWK/PhJ, and CAST/EiJ were purchased from Jackson Laboratories at 8 weeks of age and transferred to UCLA where they were kept for at least 7 days before tissue extraction. Animals were housed with ad libitum food and water in a 12 h light-dark cycle. Mutant KO animals were generated on a C57BL/6J background.
Published data were obtained from three the Jackson Laboratory website (https://phenome.jax.org), the Gene Network website (https://genenetwork.org) and PubMed. PubMed articles were identified using search terms including ‘fear’, ‘anxiety’, ‘elevated plus maze’, and ‘fear conditioning’. This identified 29 potentially useful datasets.^16^^,^^17^^,^^18^^,^^19^^,^^20^^,^^21^^,^^22^^,^^23^^,^^24^^,^^25^^,^^26^^,^^27^^,^^28^^,^^29^^,^^30^^,^^31^^,^^32^^,^^33^^,^^34^^,^^35^^,^^36^^,^^37^^,^^38^^,^^39^^,^^40^^,^^41^^,^^42^^,^^43^^,^^44^^,^^45^
To harmonize the phenotypic data for the meta-analysis, we had to deal with the heterogeneity of the datasets collected. Each experiment used different protocols, different testing equipment and often reported different, though related, assays. Different studies used different names for what was supposedly the same phenotype. For example, for fear conditioning some studies reported the time freezing,^26^ other reported its opposite, the amount of activity.^19^ Some recorded 3 min of freezing, some 5 min, or more. EPM behavior was recorded in apparatus of different dimensions, under different lighting conditions, at different times of day, all of which are known to affect behavior.^104^
To deal with this, after downloading data and examining what was most collected, we decided to use the time freezing during a 5-min period after exposure to an altered context (contextual freezing), and the mean time spent freezing during exposure to a conditioned stimulus (most usually a tone). For the elevated plus maze we used the number of entries into the open arms. We created a data dictionary to convert the different phenotype names to a common set, and this is provided in Table S2.
Some studies also provided information about sex and age, and we tested, using a linear model, whether there was an effect on phenotypes of either. Since we found no significant effects (at a 5% threshold after correcting for multiple testing) and since some studies did not provide this information, we did not include covariates in downstream analyses.
We examined the relationship between phenotypes from different studies in two ways. First, for each strain for each study we obtained the trait mean and then calculated the correlation between the same strains in different studies. Not every study tested the same strains, so although there were 42 possible pairwise comparisons for FC-context, only 30 comparisons could be carried out. One study (Palmer5) was dropped because of insufficient overlap in the strains used with other studies. There was a considerable range in the correlations between studies (from −0.04 to 0.93), though only 13 comparisons were significant (at a p < 0.05, uncorrected for multiple testing), in part reflecting differences in the variation in the number of strains in each comparison (from 3 to 31). In some cases, a pair of studies were correlated, but not with other studies (for example for the FC-cue phenotype, study Bolivar2 is correlated with Bothe-2005 (r = 0.85) but neither study was highly correlated with others).
As a second measure of the relationship between studies, we ordered the means of each strain, and ranked them. We then applied the order for one study to another to examine the extent to which studies agreed on the order of strain means. We decided to use the pattern of correlations to select studies for the genetic analysis, requiring the presence of more than 5 strains in common, and a positive correlation greater than 0.2 (though we did not require this to be significant since). The studies we chose, together with the number of animals and the phenotype names, is given in Table S1.
Mice from the Mouse Diversity Panel (HMDP) were used for the behavioral analyses. Mice (n = 700) were obtained through Jackson Laboratory at approximately 60 days old and housed for a 14-day acclimation period prior to testing. Mice were housed in groups (3–4 per cage) under a 12h/12h day/night cycle with ad libitum access to food and water. Testing was carried out between 10 a.m. and 4 p.m. Auditory background stimulus in the form of white noise (80db) was delivered through overhead speakers. All protocols conformed to NIH Care and Use Guidelines and were approved by UCLA’s Animal Care and Use Committee (protocol number ARC-2018-026).
All tests were completed at the Behavioral Testing Core at UCLA. Animals were handled for 5 days prior to experiments. The order of experiments for genetic mapping was EPM -> FC cue/context. EPM -> 7 days-> FC training -> 24 h -> FC context -> 24 h -> FC tone.
Animals were placed into an elevated plus maze apparatus (arms 5 cm × 30 cm) for 5 min with 15 lux at the center. We measured time spent in open/closed arms, and total entries into open/closed.
Mice were placed in the conditioning chamber (context A) for 3 min before the onset of the discrete conditioned stimulus. The parameters were as
Tone: pure tone, 2700Hz, 80dB, rise/fall 50, duration 30s.
Shock: 0.75 mAmps, duration 1s, inter-trial interval 180s.
Animals were exposed to three pairings of the CS and US. After the CS-US pairings, the mice were left in the conditioning chamber for another 60 s and then placed back in their home cages. FC boxes and software to control the boxes were purchased from MedAssociates.
*Fear-conditioning context * 24 h after conditioning, mice were returned to the same context (context A) with no shock or tone and freezing was measured for 8 min.
Fear-conditioning cued test: 24 h after the context test, mice were placed in a novel context (context B) where, they are exposed to the same procedure as the training but without any shocks (3 min exposure followed by 3 CS-US pairings, with parameters equal to fear acquisition). Freezing was measured during the initial 3 min exposure, during tone, and during the inter-trial interval.
Behavior was recorded digitally from a camera mounted above each test chamber at 30 FPS. Following the completion of experiments, each recorded EPM video was analyzed using ANY-maze for positional tracking of the animal along the apparatus. For EPM mapping we measured open arm entries, time spent in open arms, and closed arm. VideoFreeze (MedAssociates) was used to measure animal freezing in QC testing from recorded videos. Freezing was determined as an absence of all visible movement except that required for respiration; in VideoFreeze, we set the motion threshold to 18 au with a minimum freeze duration of 1 s (30 frames). Contextual fear was assessed by total freezing time over the first 3 min of the test. Cued fear was assessed by total freezing time in response to the presentation of the first, second, and third CS during ITI. One exception was with freezing measurements during genetic mapping with BXD lines. Due to the variable coat color across lines, we used ANY-maze for animal tracking and freezing measurements as it performed better than VideoFreeze.
Our starting point for genotyping each mouse strain was variant calls from sequencing reads for those strains that have been sequenced,^46^^,^^105^ and genotypes from high density arrays for those that have not^106^ (https://www.informatics.jax.org/). We confirmed genotypes in the BxD recombinant inbred lines by randomly choosing 24 lines and subjecting to genotyping them on the GigaMUGA mouse array (140 thousand SNPs).^107^ This genotyping was done commercially (Neogen GeneSeek Operations, 4131 N 48th St. Lincoln, NE 68504). All lines were of the expected genotype.
We used imputation to generate a near complete set of genotyped markers for all strains. The highest density set of genotyped markers is that derived from 16 strains for which we have complete sequence data^46^ but this does not include any recombinant inbreds and lacks some of the strains used for mapping. Since the missing strains are descendants of a small number of founder haplotypes,^106^ we can impute the genotypes in the non-sequenced strains. Using genotypes from arrays as a scaffold^56^^,^^57^^,^^58^^,^^59^^,^^60^^,^^61^^,^^62^^,^^63^^,^^64^^,^^65^^,^^66^^,^^67^^,^^68^^,^^69^^,^^70^^,^^71^^,^^72^^,^^73^^,^^74^^,^^75^^,^^76^^,^^77^^,^^78^^,^^79^^,^^80^^,^^81^^,^^82^^,^^83^^,^^84^^,^^85^^,^^86^^,^^87^^,^^88^^,^^89^^,^^90^^,^^91^^,^^92^^,^^93^^,^^104^^,^^105^^,^^106^^,^^107^^,^^108^ we imputed missing genotypes using data from sequenced animals. Most imputation algorithms are designed for fully outbred populations, where a key problem is correctly phasing haplotypes (e.g.,^109^). Alternative methods are needed for efficient imputation in mouse strains and so we used the software tool EMINIM,^48^ which requires a reference file, containing strains with known haplotypes, and a target file for imputation, containing strains genotyped at some SNPs. We ran this with the recommended parameter settings for inbred mice given according to the software documentation.
We prepared these two files by placing strains that were genotyped at all SNPs in the reference file, while strains with one or more missing SNPs were placed in the target file. For the BXD mice, SNPs between genotype markers from the same founder haplotype were pre-filled due to the negligible likelihood of multiple recombination events between consecutive genotyped SNPs. This resulted in 27 reference strains and 227 target strains. We then ran EMINIM imputation on each chromosome separately and obtained calls on 16,767,664 SNPs (triallelic SNPs, constituting 1% of the total, were discarded). Since EMINIM outputs high confidence values for one of the alleles, we translated a 90% or higher confidence for an allele into a hard call for that allele, with lower confidence calls being called as missing. We combined these imputation results with the genotyped data. The rate of missing calls was less than 0.8% for 18 of 19 autosomal chromosomes (chromosome 19 had a missing call rate of 1.03%).
We validated the imputation results by comparing imputed calls to genotypes from sequenced BXD strains obtained from.^110^ The sequence-derived set contained 4,325,552 SNPs, of which 3,812,095 variant sites overlapped with our imputation. Across all chromosomes, each BXD strain was imputed with between 93% and 99% accuracy. Imputation accuracy for all strains combined varied between 89% and 99% per chromosome, with 96.29% overall. Due to the pre-filling procedure described previously, this accuracy reflects only the most challenging SNPs to those between genotyped SNPs from different founder haplotypes. When including the genotyped and pre-filled SNPs, the accuracy for each chromosome was greater than 99.64%.
We performed a second assessment of imputation accuracy by removing the C57BL/6J and DBA/2J strains from our reference panel, randomly masking 50% of the bases of those strains, and using the remaining inbred strains in the reference panel to impute the masked bases. Once again, we used EMINIM with its default inbred imputation parameter settings and translated a 90% or higher confidence call for an allele into a hard genotype call. Calls with less than 90% confidence were translated to a missing call, leading to a missing genotype rate of 0.09%–0.59% per chromosome. The accuracy of imputation on the masked bases was 95.6–99.7% per chromosome, and 97.9% overall.
All phenotypes were mapped in a two-stage approach. First, we used a mixed model, implemented in GEMMA^51^^,^^52^ (obtained from https://github.com/genetics-statistics/GEMMA), to map each study separately. Input files for GEMMA were generated using PLINK,^103^ obtained from https://zzz.bwh.harvard.edu/plink. We used the PLINK binary format of the required files, including all alleles with a frequency greater than 10%. A genetic relationship matrix was generated using the
gemma -bfile input.file -gk -o kinship.file.
Mapping was carried out using the
gemma -bfile pl input.file -k kinship.file -lmm -o output.file.
After mapping with GEMMA, the results files were combined with purpose written perl scripts into a format suitable for meta-analysis with METASOFT.^53^^,^^54^ METASOFT (downloaded from http://genetics.cs.ucla.edu/meta_jemdoc/) was invoked with the following
java -jar Metasoft.jar -input input.file -mvalue -output output.file.
GEMMA assigns effects to the minor allele, as determined by the.bim files from PLINK which means there is considerable variation between studies as to which allele gets an effect (due to allele frequency variation). After running each study though GEMMA we swapped alleles to the most common arrangement, before proceeding to the meta-analysis.
METASOFT calculates a genomic-control inflation factor for the mean effect (lambda_Mean Effect) and a genomic-control inflation factor for heterogeneity (lambda_heterogeneity). These values were calculated for each run, and the software was run a second time, to include the effect of these variables. For
java -jar Metasoft.jar -input input.file -mvalue -output output.file -lambda_hetero 2.738124 -lambda_mean 2.430179.
Simulations to estimate the appropriate significance threshold were carried out using GCTA,^55^ obtained from https://yanglab.westlake.edu.cn/software/gcta. Phenotypes were simulated using real genotype data under a simple additive genetic model. A set SNPs was chosen as the causal variants that were seen to have effects in the mapping studies. We obtained estimated effect sizes for these SNPs from GEMMA output files and converted the beta estimate from the linear model to the values required by GCTA using
Beta_gcta=beta_orig∗sqrt(2∗p∗(1−p))
(where p is the allele frequency from the GEMMA output file)
Input SNP files were those used for the mapping with GEMMA. 1000 simulations were carried out using the
gcta64 --bfile input.SNP.file --simu-qt --simu-causal-loci causal.snplist --simu-hsq 0.2 --simu-rep 3 --out.
After simulating data for each component study, files were processed through GEMMA and METASOFT exactly as for the real data.
To investigate the contribution of genotype structure to the inflation seen in quantile-quantile plots we simulated data in the same way, but with no causal SNPs.
We identified genes lying within QTLs, and determine whether they are expressed in brain tissue from in situ hybridization data (Allen Brain Atlas,^60^ the mouse ENCODE transcription dataset^61^), and RNA-seq data obtained from the following ^62^^,^^63^^,^^64^^,^^65^^,^^66^^,^^67^ For ISH data from the Allen dataset, we looked for presence of expression. For ENCODE transcriptome data we set a cutoff for an RPKM >5 in at least one brain region category, and for RNA-seq studies we ranked expression of each gene and averaged across studies; any gene in the top 25% of all genes (or >rank 6000) was considered expressed. For the Allen institute scRNA-seq dataset, we only included sequencing data from cell type clusters from CA1-CA3 of the hippocampus. Genes that passed expression thresholds for all three datasets were included for further consideration.
Where an existing knockout had been successfully created for the gene, usually on a C57BL/6N background, we adopted the same design, targeting the same part of the gene but in the C57BL/6J strain. To generate each KO, we used a CRISPR-Cas9 genome editing approach. For each gene, two single guide RNA (sgRNA) were generated by Synthego Corp. (Redwood City, California, USA) for targeting, except for 4933413L06Rik where we carried out a knock-in with a single-stranded DNA synthesized by IDT, Inc (Coralville, Iowa, USA). Guide RNAs (gRNAs) were designed to maximize both efficiency and specificity scores.^111^ The sgRNA and/or ssDNA, and Cas9 protein (Alt-R SpCas9 Nuclease V3 from IDT) were co-injected into zygotes of C57BL/6J mice (Jackson Laboratory, Stock#009086). After microinjection, zygotes were transferred into the oviduct of pseudo-pregnant Swiss Outbred mice (Jackson Laboratory, Stock#034608) to generate founder mice. Founder mice were obtained and confirmed by PCR genotyping and Sanger sequencing. The correctly targeted founders were bred with C57BL/6J mice, and PCR genotyping and Sanger sequencing were again performed to confirm germline transmission. We used quantitative PCR to confirm that the deletion affected abundance of RNA transcript in hippocampal tissue from F1 animals. Full descriptions of each knockout design, and sgRNA sequences are given at dx.doi.org/10.6084/m9.figshare.25521283. Genotyping primers are listed in Table S4.
To confirm the effect of the deletion on RNA production at each locus we ran quantitative polymerase chain reactions using a standard RT–qPCR assay (Taqpath COVID-19 Combo Kit). We analyzed RNA extracted from the hippocampus of heterozygous knockouts on C57BL6/J animals. Primer sequences can be found in Table S5.
We generated offspring from four crosses.^4^ C57BL/6J animals were mated to DBA/2J, and to heterozygote C57BL/6J animals where one allele was the knockout. DBA/2J animals were also mated to the heterozygote knockout mice. This generated four wildtype C57BL/6J, F1 DBA/2J/C57BL/6J (wildtype), heterozygote C57BL/6J knockouts, and F1 DBA/2J/C57BL/6J (knockout). All four groups were assayed for the relevant phenotypes associated with the QTL, using the phenotyping protocols described above.
Evidence that the gene is the QTL gene was then detected as a statistical ‘cross’ (mutant or wildtype) by ‘strain’ (DBA/2J or C57BL/6J) interaction in a linear model, using the R statistical programming language.^112^ Batch and sex were included as covariates, using the R command
summary(lm(phenotype∼sex+batch+cross∗strain,data=data))
Adult male animals (Jackson Laboratories) were euthanized at 10–16 weeks old in an isoflurane chamber and decapitated. The brain was removed and the ventral region of the hippocampus was microdissected, snap frozen in dry ice, and stored at −80°C until processing. Tissue from ∼2 animals were combined into a single tube and considered a replicate, with 2 replicates per strain for snmC-seq2, snRNA-seq, and snATAC-seq experiments.
Adult male animals (Jackson Laboratories) were euthanized at 10–16 weeks old in an isoflurane chamber and decapitated. The brain was removed and coronal brain slices containing amygdala tissue were generated on a 1mm brain matrix (World Precision Instruments). Amygdala tissue was micro-dissected from these slices under a dissecting scope in cold PBS, snap frozen in dry ice, and stored at −80 until processing. Tissue from ∼2 to 3 animals were combined into a single tube and considered a replicate, with 2 replicates per strain for snmC-seq2, snRNA-seq, and snATAC-seq experiments.
We carried out snmC-seq2 on microdissected tissue as previously described.^70^ Briefly, frozen tissue was homogenized into single nuclei suspensions with Dounce homogenization, then immediately sorted on into a 384-well plate with a FACSAria sorter (BD Biosciences) at the UCLA Flow Cytometry Core. We selected for a 75-25 enrichment of neuronal vs. non-neuronal nuclei during FACS sorting using NeuN-488/DAPI counterstains (Millipore Sigma MAB377X). Bisulfite conversion and single-cell methylome libraries were generated following this step.
Single nuclei suspension and library generation were completed at the Cedars Sinai Applied Genomics, Computation and Translational Core and followed the 10X protocol for the Chromium Next GEM Automated Single Cell 3′ Library and Gel Bead Kit v3.1 (cat# PN-100014) as described except for the following
Suspensions from cell nuclei were generated using the recommended method from the 10X scMultiome protocol (CG000375 Rev C) to lyse cells and obtain nuclei. Following single nuclei suspension generation, nuclei were counterstained for 7-AAD and NeuN-405 antibody (Novus Biologicals, 200) and sorted on a MACSQuant Tyto (Miltenyi Biotech) prior to GEM generation. We selected for a 75-25 split of NeuN+/7-AAD+ nuclei for neurons and NeuN-/7-AAD+ for non-neuronal nuclei respectively.
We captured ∼10,000 nuclei per genotype per region per replicate on a single 10X GEM chip. All downstream library preparation was done according to the 10X Genomics protocol (CG000286) and sequenced on a Novaseq 6000 with a target of ∼40-50k reads per nucleus.
Single nuclei suspension and library generation were completed at the Cedars Sinai AGCT core and followed the 10X protocol for Next GEM scATAC-Seq v1.1 (PN-1000175) as described except for the following
Nuclei suspensions were generated using the recommended method from the 10X scMultiome protocol (CG000375 Rev C) to lyse cells and obtain nuclei.
Following single nuclei suspension generation, nuclei were counterstained for 7-AAD and sorted on a MACSQuant Tyto prior to GEM generation. NeuN was not used for neuronal enrichment due to dye incompatibility between our NeuN antibody and a nuclear counterstain. After the sort, we carried out permeabilization of nuclei as per the protocol. We aimed to capture 10,000 nuclei per well x 8 wells, for a total of 80,000 nuclei over 8 total samples (∼10,000 nuclei per genotype per region per replicate). All downstream library preparation was done according to the 10X Genomics protocol (CG000209) and sequenced on a Novaseq 6000 with a target of >35k reads per nucleus.
All quality control and pre-processing were done under the Seurat package framework.^96^ Per biological sample, we filtered out cells that (1) fall below the 5th percentile of the total UMI counts (nCount_RNA) or the 5th percentile total number of unique genes expressed (nFeature_RNA) or 700 unique genes expressed, whichever was more stringent; (2) are over the 95th percentile quantile in terms of either the total UMI counts or the total number of unique genes expressed; (3) have larger than 5% mitochondria fraction (percent.mt).
Global coverage counts per million (CPM) was applied to each cell followed by log transformation (“LogNormalize”). We then projected cells from each biological sample to low dimensional space using principal components analysis (PCA) on highly variable features selected by Seurat. Potential doublets were identified and subsequently removed from the downstream analysis by DoubletFinder,^99^ ran in the top 15 principal components space with the expected doublet rate set to the recommended amount from 10X genomics based on high quality yield volume.
All quality control and pre-processing were done under the Seurat and Signac package framework.^97^ Per sample, we first used Model-based Analysis for ChIP-Seq (MACS) to call sample specific de-novo peaks from its fragments file.^98^ We then merged sample-specific sets of peaks to a unified peaks set while removing peaks with length larger than 10000bp or smaller than 20bp. A unified peaks by cells count matrix was constructed from the fragments file while removing cells with lower than 200 peaks detected and peaks only present in less than 10 cells. Cells were filtered based on the following (1) appropriate number of non-duplicate, usable read-pairs (passed_filters from Cell Ranger’s output singlecell.csv). Specifically, we set it to larger than 3000, 4000, 2500, and 5000 for the 2 BL6 and 2 DBA samples collected in the hippocampus region. Similarly, we set it to be larger than 7500, 5000, 7500, and 7500 for their amygdala counterpart. (2) number of fragments overlapping peaks (peak_region_fragments from Cell Ranger’s output singlecell.csv) falls within the 5th percentile and the 95th percentile. (3) ratio of fragments overlapping peaks over the total number of non-duplicate, usable read-pairs falls within the 5th percentile and the 95th percentile. (4) nucleosome_signal: the ratio of fragments between 147 bp and 294 bp (mononucleosome) to fragments <147 bp (nucleosome-free) is smaller than 4 (5) TSS enrichment score is larger than 2. We did not include a filter for ratio of peaks in black list regions over the total number of non-duplicate, usable read-pairs as this was removed during the construction of the DBA SNP-swapped reference genome. We normalized the count data with Text Frequency Inverse Document Frequency (RunTFIDF) and performed Singular Value decomposition (RunSVD) on top 90% informative features selected by Signac (FindTopFeatures). The first low dimensional embedding was excluded from downstream doublet detection and clustering analysis due to high correlation with sequencing depth. Potential doublets were identified and subsequently removed from downstream analysis by DoubletFinder,^99^ ran on the 2^nd^ - 11^th^ low dimensional embedding with the expected doublet rate set to the recommended amount from 10X genomics based on high quality yield volume. Finally, we built a gene-by-cell transcriptional activity matrix that counts per cell, at the gene body and 2000bp upstream to capture the promoter region, the total number of ATAC-seq counts.
Cells were filtered on the basis of several metadata (1) mCCC level <0.03; (2) global mCG level >0.5; (3) global mCH level <0.2; (4) total mapped reads >100,000; (5) Bismarck mapping rate >0.5; and 6 (percent genome covered >2). Methylation features were calculated as fractions of methylcytosine over total cytosine across gene bodies ± 2kb flanking regions and 100kb bins spanning the entire genome. Methylation features were further split into CG and CH methylation types. Features overlapping our methylation mm10 blacklist were removed. 100kb bin features were then filtered on minimum mean coverage >500 and maximum mean coverage <3000. Gene body features were filtered on minimum coverage >5 and all remaining features were normalized per cell using the beta binomial normalization technique in allcools.^72^ Individual CpG sites were only counted when they had ≥5 reads covering the site.
Gene counts were normalized using SCTransform,^102^ and regressed out percentage of reads from mitochondrial genes. We then integrated cells from all samples using reciprocal principal components analysis (rPCA) implemented in Seurat 4.0.5^96^ on the top 5000 integration genes and in the top 30 reciprocal principal components space. For clustering, we standardized the integrated data, performed PCA on all integrated genes and ran de-novo Louvain clustering algorithm in the top 15 principal components space with resolution set to 0.1. Cluster markers that are conserved between the strains were called using non-parametric Wilcoxon rank-sum test and subsequently used for annotation. We annotated hippocampals clusters by manually checking conserved markers against the ALLEN BRAIN MAP’s Mouse Whole Cortex and Hippocampus dataset. For cells collected in the amygdala, we performed a second round of analysis restricted to those identified to be neuronal in the initial annotation to increase our annotation resolution. We used the same integration pipeline and parameters described above except when clustering we ran de-novo Louvain clustering algorithm in the top 25 principal components space with resolution set to 0.75, due to the increased expected complexity of cell types in this brain region. All clusters were first grouped into broad types of inhibitory or excitatory neurons except for one cluster that simultaneously expressed Gad2 and Slc17a7. We thus simply annotated it as “Neuronal” followed by its by cluster marker gene names. We also manually compared our clusters to those that were previously reported in amygdala,^113^ and adapted their annotation for clusters with high concordance in top marker genes.
We first jointly projected all cells’ ATAC peak profile to uncorrected Latent Semantic Indexing (LSI) embeddings with TFIDF transformation followed by calculating SVD on the top 90% most informative peaks. Peak profile embeddings were then integrated in shared low dimensional space via integration anchors identified in the 2^nd^ to 30^th^ reciprocal LSI space as implemented in Signac 1.5.0. We then integrated cell transcriptional activity profiles by performing SCTransform after regressing out percentage of activity from mitochondrial genes and carried out rPCA integration on integration genes identified from the single-nucleus RNA experiment and in the top 10 reciprocal principal components. We transferred the single-nucleus RNA annotation onto the single-nucleus ATAC cells by linking the RNA’s expression profile with ATAC’s transcriptional activity profile through canonical correlation analysis (CCA) described in Seurat. Pairs of cells from each modality that are mutual nearest neighbors in the top 15 canonical components space were identified as “transfer anchors”. “Anchors” were further filtered and weighted by distances in the integrated peak embeddings prior to imputing ATAC cells’ annotation of cell type. For amygdala cells, after the initial assignment to cell type, we repeated the above-mentioned pipeline only on those cells predicted to be neuronal to obtain a set of higher resolution annotation. The second round of analysis was performed on integration genes identified from the subset of single-nucleus RNA data also categorized as neuronal.
We used the negative of the average mCH fraction of the gene body ± 2kb as the proxy of methylation cells’ transcriptional activity as described previously.^72^ We first integrated, across strains, single-nucleus methylation gene-body mCH profiles on RNA integration genes via “integration anchors” identified in the top 30 reciprocal PCA space. Gene body mCH profiles were then linked to single-nucleus RNA expression profiles via “transfer anchors” identified in the top 15 canonical components space. “Transfer anchors” were further filtered and weighted by distances in the integrated mCH embeddings prior to imputing methylation cells’ annotation using their RNA counterpart as reference. For amygdala, we directly linked to the higher resolution annotation used on the neuronal cells in RNA. These served as our final annotation for all the neuronal cells in both hippocampus and amygdala. Since mCH methylation, used to construct the “transfer anchors”, is largely not present in the non-neuronal population,^114^ we further performed de-novo annotation on methylation cells using both the genome wide mCG and mCH fraction at 100kb as described previously.^115^ Non-neuronal clusters were annotated by their canonical marker genes (gene-body hypo-methylation). These served as our final annotation for all the non-neuronal cells in both hippocampus and amygdala.
Both single-nucleus methylation and single-nucleus ATAC cells’ expression profiles were imputed with previously computed “transfer anchors”. All three modalities were merged on their integrated or imputed expression profiles, projected to low dimensional space via PCA, and visualized by UMAP (performed on the top 15 principal components).
We used DESeq2^73^ for cell-type specific pseudobulk level differential expression analysis and differential accessibility analysis. Per cell type, raw counts were aggregated to replicate level and DESeq2 was run under default parameters to detect statistical evidence of strain differences. Any genes with no coverage in either or both strains were excluded.
We used the brms package in R to test for the effect of class analysis.^116^ We defined a Bayesian model using brms with a negative binomial distribution for the number of significant differences (for RNA and ATAC data) or the number of mutations for methylation. We used up to eight chains and visually inspected the output to ensure chains were mixed. We used the following code to define the
This work was funded in part through NIH grants R01MH115979, R01MH125252, U01MH130995, and DP1DA044371, the UCLA Jonsson Comprehensive Cancer Center, the Eli and Edythe Broad Center of Regenerative Medicine, and the Stem Cell Research Ablon Scholars Program. We would like to gratefully acknowledge the UCLA Behavioral Testing Core for assisting with the behavioral experiments, the Cedars-Sinai Applied Genomics, Computation & Translational Core for assisting with the sequencing experiments, and the Genome Editing Core at Augusta University for assisting with the generation of the mouse models necessary to support this research.
R.C. performed the mouse breeding and behavioral experiments. C.L. performed the single-cell methylation analyses, data from which were processed by M.G.H. and Z.C. P.B.C. performed the snRNA and snATAC analyses, with data processing by Z.C. J.F., P.B.C., and D.C.S. analyzed data from the single-cell experiments. J.F. and N.L. analyzed the mouse mapping data. J.M. provided statistical support. J.F. and P.B.C. designed the study. P.B.C., J.F., and E.M. co-wrote the manuscript, which was subsequently reviewed and edited by the rest of the authors.
The authors declare no competing interests.
Published: May 1, 2024