Humanized substitutions of Vmat1 in mice alter amygdala-dependent behaviors associated with the evolution of anxiety.

The human vesicular monoamine transporter 1 (VMAT1) harbors unique substitutions (Asn136Thr/Ile) that affect monoamine uptake into synaptic vesicles. These substitutions are absent in all known mammals, suggesting their contributions to distinct aspects of human behavior modulated by monoaminergic transmissions, such as emotion and cognition. To directly test the impact of these human-specific mutations, we introduced the humanized residues into mouse Vmat1 via CRISPR/Cas9-mediated genome editing and examined changes at the behavioral, neurophysiological, and molecular levels. Behavioral tests revealed reduced anxiety-related traits of Vmat1 Ile mice, consistent with human studies, and electrophysiological recordings showed altered oscillatory activity in the amygdala under anxiogenic conditions. Transcriptome analyses further identified changes in gene expressions in the amygdala involved in neurodevelopment and emotional regulation, which may corroborate the observed phenotypes. This knock-in mouse model hence provides compelling evidence that the mutations affecting monoaminergic signaling and amygdala circuits have contributed to the evolution of human socio-emotional behaviors.

Introduction

Recent years have seen a surge in the study revealing distinct human brain characteristics at the molecular, cellular, and circuit levels in both cortex and subcortical structures (Pryluk et al., 2019; Raghanti et al., 2018; Sousa et al., 2017a; Zhu et al., 2018). Differences in brain monoaminergic signaling are likely to be key factors conferring human attributes, as various monoaminergic pathways are involved in emotional regulation, memory consolidation, cognitive flexibility, and complex social behaviors (Cools et al., 2008; Ebstein et al., 2010; Shohamy and Adcock, 2010). Comparative studies have indeed demonstrated the distinct monoaminergic mechanisms in the human brain (Raghanti et al., 2018; Sousa et al., 2017b). Moreover, the diversification of monoaminergic mechanisms has been linked to the evolution of social behavior in non-human primates (Bergey et al., 2016; Staes et al., 2019). Thus, the genetic differentiation of monoaminergic genes may underlie our unique emotional and social characteristics, such as empathy and altruism (Raghanti et al., 2018).

In line with these studies, we recently reported that the vesicular monoamine transporter 1 (VMAT1, also known as SLC18A1) gene has been under positive selection in the human lineage as evidenced by two amino acid substitutions (Glu130Gly and Asn136Thr) not shared by other primates (Sato and Kawata, 2018). The VMATs transport monoamine neurotransmitters into synaptic vesicles, thereby regulating the kinetics of synaptic transmission. There are two isoforms of VMAT, VMAT1 and 2, with the VMAT2 being highly expressed in the brain and widely studied (Erickson et al., 1996; Lohr et al., 2014). Although VMAT1 had previously been thought to be expressed solely in the peripheral nervous system and chromaffin cells, its functional importance in the central nervous system (CNS) has recently become evident (Lohoff et al., 2014). Our previous study showed that a new mutation at the 136th amino acid (136Ile) of VMAT1 has emerged relatively recently in the course of human evolution and that the human-specific polymorphism Thr136Ile has likely been maintained by balancing selection within non-African populations (Sato and Kawata, 2018). Fluorometric assays additionally demonstrated that these human-specific residues enhance monoamine neurotransmitter uptake, suggesting the tendency of lower monoamine uptake in early hominids (Sato et al., 2019). The 136Thr, the dominant variant across populations, is known to be marginally associated with psychiatric diseases, such as bipolar disorder, depression, and anxiety, as well as with neuroticism, a personality trait conferring a greater risk of psychiatric disease (Lohoff et al., 2006, 2008a; Vaht et al., 2016). Given the functional changes and behavioral outcomes associated with the VMAT1 sequence substitutions, this gene also might have affected a certain portion of our emotional evolution.

Although the evolution of VMAT1 and its basic functions in the mammalian CNS have been clarified, their contributions to unique human attributes are still unclear. Here we examine the molecular, neurological, and behavioral changes in mice caused by human-specific VMAT1 mutations. We show that humanized Vmat1 (Vmat1Ile) induces changes in gene expression and neuronal activity in the amygdala, reduces anxiety-like behaviors under anxiogenic conditions, and even enhances preference for social novelty. These findings collectively suggest the functional importance of human-specific mutations in VMAT1, and the present work provides a new experimental strategy for evolutionary studies of human emotions.

Results

Conservation of the 136th residue of VMAT1 across non-human vertebrates

To investigate the evolutionary significance of the 136th residue in human VMAT1 (hVMAT1), we first examined the conservation of this residue in a wide range of taxa in vertebrates. Multiple sequence alignments of VMAT1 among 236 vertebrate species revealed the unique evolution of this residue in the human lineage (Figure 1A). Notably, we confirmed that almost all vertebrates except bicolor damselfish (Stegastes partitus) retain asparagine (Asn) at the 136th amino acid, whereas humans are the only species that carries the unique Thr136Ile (rs1390938) polymorphism.

Figure 1

An ultra-conserved residue in VMAT1 exhibits functional variants unique to humans

(A) A phylogenetic tree constructed by multiple sequence alignment of the VMAT1 gene. Almost all genes across 236 vertebrate species retain asparagine (Asn) on the 136th residue, whereas humans are the only vertebrate species except for bicolor damselfish (Stegastes partitus; shown in yellow) with a unique polymorphism (Thr136Ile; rs1390938). A gap in the aligned sequence is shown in gray. Note that the phylogenetic relationship presented here is not necessarily consistent with the known species tree.

(B) hVMAT1 and (C) mVMAT1 protein structure predicted by homology modeling and the effects of human-type mutagenesis (corresponding to 133Thr and 133Ile in mVMAT1). (Top) When introduced in silico, 133Ile, and not 133Thr, exhibits hydrophobic interactions (shown by green dotted lines) with surrounding sites, which likely influences the folding and/or stability of mVMAT1 protein. Blue, red, and orange dotted lines represent amide bonds, hydrogen bonds, and weak van der Waals interactions, respectively. (Bottom) Δ Vibrational entropy energy between WT (133Asn) and mutants. Amino acids are colored according to the vibrational entropy change conferred by the given mutation. Red represents a gain of flexibility and blue represents a rigidification of the structure. The 133Ile mutation leads to the increased flexibility of the first luminal loop, a receptor-like domain affecting the affinity of ligands.

Predicted effects of humanized residues on mouse VMAT1 structure and function

In silico predictions of wild-type mouse VMAT1 (mVMAT1) protein structure revealed that like the hVMAT1 residue 136, the corresponding residue in mVMAT1 (133Asn) was located between the first luminal loop and the second transmembrane domain (Figures 1B and 1C), suggesting a similar function. Given the evolutionary conservation across vertebrates and the difference in chemical properties between amino acids, the humanized mutations to mVMAT1, especially 133Asn to 133Ile, were predicted to have a significant impact on protein function (Table 1). Furthermore, the Asn133Ile mutation was also predicted to stabilize the protein structure as evidenced by greater folding free energy (ΔΔG: 1.106 kcal/mol; Table 1) and to increase the flexibility of the loop domain (Figure 1C), which possibly increases the efficiency of monoamine uptake.

Table 1

In silico prediction of mVMAT1 tolerance to the humanized mutations Asn133Thr and Asn133Ile

Ref.	Alt.	Provean		SIFT		DynaMut
		Score	Prediction	Score	Prediction	ΔΔG (kcal/mol)	Prediction
Asn	Thr	–3.27	Deleterious	0.243	Tolerated	–0.062	Destabilizing
Asn	Ile	–6.55	Deleterious	0.014	Damaging	1.106	Stabilizing

The bold text indicates significant effects of the substitutions as predicted by Provean score < –2.5 or SIFT score < 0.05. The Provean term “Deleterious” and SIFT term “Damaging” refer to significant effects on protein function but do not necessarily indicate that they deteriorate protein function.

Generation of Vmat1-humanized mouse models by CRISPR/Cas9 genome editing

To replace the 133Asn of mouse Vmat1 with Thr or Ile, we employed the CRISPR/Cas9-mediated genome editing strategy (Figure 2A). We selected a guide RNA targeting the 133rd Asn codon in Vmat1 exon 4 by using the web-based CRISPR design tool to minimize off-target cleavage risk (Tables S1 and S2, see STAR Methods for details). By electroporating CRISPR components (Vmat1 exon 4 crRNA, tracrRNA, and Cas9 nuclease) along with the ssDNA donor containing the desired substitutions into C57BL/6J fertilized eggs (Figures 2A and S1A), we successfully generated 133Thr founder No.1 and 133Ile founder No.10 (Figures S1B–S1D). After confirming the absence of off-target cleavages (Figure S2), we crossed 133Thr founder No.1 and 133Ile founder No.10 with wild-type (WT) C57BL/6J mice to obtain heterozygous F1 generations that were again verified to possess the designed substitutions by Sanger sequencing. The homozygous mice were further crossed after the F5 generations to obtain three humanized genotypes, Thr/Thr, Thr/Ile, and Ile/Ile, in addition to WT, for behavioral, neurophysiological, and transcriptomic analyses (Figure S3, see STAR Methods).

Figure 2

Generation of the Vmat1-humanized mouse models by CRISPR/Cas9 genome editing

(A) Targeting strategy for mVMAT1 133Asn humanization. The genetic configuration of the mouse Vmat1 gene is shown above. Exon 4 encoding 133Asn is enlarged, and the primers used for genotyping (Exon4_F and Exon4_R) are depicted. To replace the mouse 133Asn with 133Thr or 133Ile by CRISPR/Cas9-mediated genome engineering, a guide RNA with minimum off-target effects was designed. GGG (gray) represents the PAM sequence. In addition to 133Asn humanization, restriction enzyme recognition sites (EcoRI and FspI) were synonymously incorporated to avoid unwanted re-editing and to simplify genotyping. Sanger sequencing profiles of 133Asn/Asn (WT), 133Thr/Thr, 133Ile/Ile, and 133Thr/Ile are shown on the right.

(B) PCR-RFLP assay, in which PCR products amplified using Exon4_F and Exon4_R were digested by EcoRI and FspI, respectively, could be used to distinguish four genotypes without sequencing.

The human Ile substitution in mVMAT1 reduced behavioral anxiety

We accordingly conducted comprehensive behavioral tests among the 20 sets of mice including one of each genotype (Vmat1WT, Vmat1Thr/Thr, Vmat1Thr/Ile, and Vmat1Ile/Ile) from the first batch at 9–55 weeks (see Table S3 for detailed information on age and the number of individuals used in each test). Detailed results were shown in Figures S4–S19. Tests of social interaction (SI) between stranger mice of the same genotype in a novel environment (one-chamber SI test) and preference for a novel stranger over a previously exposed stranger (now familiar) in Crawley’s three-chamber social interaction (CSI) test revealed marginal differences in behavioral phenotypes between mice with and without Vmat1Ile (P = 0.062 for interactive effects of genotype (Ile) and place, generalized additive model (GAM) with quasi-Poisson distribution; Figures 3A, S9, and S10). In the CSI, Vmat1Thr/Thr mice significantly preferred familiar to stranger mice (P = 3.26×10-4), whereas no such tendency was observed in Vmat1Ile mice (P = 0.759 for Vmat1Thr/Ile mice and (P = 0.195 for Vmat1Ile/Ile mice; Figure 3A). Moreover, we found consistently lower levels of anxiety among Vmat1Ile mice across tests, including a relatively greater preference for the (anxiogenic) light box during the two-chamber light/dark transition (LD) test and the center of an open field (OF) (Figures S5 and S6).

Figure 3

Comprehensive behavioral tests reveal distinct behavioral changes in Vmat1Thr/Thr and Vmat1Ile mice, including reduced anxiety

(A) In Crawley’s social interaction test, Vmat1Thr/Thr mice preferred a familiar (previously exposed) mouse over a novel stranger.

(B) Structural Equation Modeling (SEM) revealed the best fit model to explain the effects of the Vmat1 genotype on locomotor activity (Act) and anxiety (Anx). According to the model, Vmat1Ile (Vmat1Thr/Ile and Vmat1Ile/Ile) significantly reduces anxiety-like behavior. The boxes and circles represent measured and latent variables, respectively. The paths represent causal relationships on which numerical values indicate standardized coefficient and the gray-scale intensity of the paths indicates statistical significance tested by t-tests. See STAR Methods for the detailed modeling procedure.

(C) Vmat1Ile mice exhibited lower anxiety than WT and Vmat1Thr/Thr genotypes. Male, 10 to 19-week-old mice were used in the behavioral experiments shown here (see the detailed information in Table S3). For the composite anxiety score, a higher value indicates lower anxiety. Statistical significance was evaluated by paired t-test for (A), and by pair-wise t-test with FDR correction by the Benjamini–Hochberg method for (C). Interactive effects of genotype (Ile allele labeled by 1, and 0 otherwise) and cage place were also assessed by the generalized additive model with quasi-Poisson distribution in (A). †: 0.05 < P < 0.1, ∗0.01 < P < 0.05, ∗∗0.001 < P < 0.01, ∗∗∗P < 0.001.

To confirm a pervasive effect on specific behavioral domains by humanized Vmat1, we constructed composite behavioral phenotype scores by summating standardized scores on related tests (LD, OF, elevated plus maze (EP), and SI tests) and compared the results among genotypes. We further conducted Structural Equation Modeling (SEM) to investigate the effects of genotype on these composite scores. This analysis revealed a significant effect of 133Ile on anxiety scores (R2 = 0.15, P = 8.7×10-3) but only marginal effects on locomotor activity (R2 = 0.062, P = 0.053; Figure 3B). Anxiety levels were generally lower in Vmat1Ile mice (P = 0.016, one-way ANOVA; Figure 3C).

In addition to standard tests of anxiety, exploratory tendency, and locomotor activity, we also performed a novel delayed reward task to evaluate the differences in impulsivity between Vmat1Thr/Thr and Vmat1Ile/Ile mice (n = 10 for both genotypes, see STAR Methods for the detailed experimental procedure; Figures S19). Throughout the test, there was no statistically significant difference in food preference between genotypes (P = 0.50 for the first term, GAM with quasi-Poisson distribution; Figure S19D). However, Vmat1Ile/Ile showed a lower preference for the higher calorie but delayed food reward than Vmat1Thr/Thr mice every test day from 5 to 9 except for day 6 compared with the basal level after habituation (day 4) (Figure S19D).

Differential regulation of dopaminergic and neurodevelopmental genes in the amygdala by Vmat1 genotype

We then conducted RNA-Seq analysis to investigate the effects of the four Vmat1 genotypes on gene expression patterns in three brain regions likely associated with the behavioral phenotypes (prefrontal cortex, amygdala, and striatum) (n = 4 mice for each genotype; see Table S4 for detailed information on samples and reads). Principal component analysis revealed distinct expression patterns among brain areas, but no obvious differences between genotypes and ages (4 vs. 10 months) (Figure S20).

We thus screened for differentially expressed genes (DEGs) by pair-wise comparisons between genotypes, which yielded a total of 80 DEGs exclusively in the amygdala (Figure 4A). A large proportion of these genes (56 out of 80) were differentially expressed between Vmat1WT and Vmat1Ile/Ile mice (Figure 4A). The expression levels of DEGs for each sample were shown in Figure S21. We then utilized the correlations between individual composite anxiety scores and expression levels of the detected amygdalar DEGs among a subset of eight mice (n = 2 for each genotype). Particularly strong correlations were found for nine genes, Baz1a, Camk2n1, Ccdc114, Cldn1, Fsbp, Gdnf, Gm42608, Gsg1l, and KI (P < 0.1; Figures 4B and S22). The DEGs were further examined for gene ontology (GO) enrichment, which revealed a significant overrepresentation of genes involved in “behavior,” “fear response,” “response to amphetamine,” “cardiac muscle tissue development,” or “cAMP catabolic process” (Table 2) between WT and Ile/Ile mice, whereas no GO terms were enriched in DEGs for other genotype comparisons (as there were few such genes). Moreover, gene expression patterns in the Ile mice amygdala were significantly positively correlated with those of striata of transgenic mouse models of Huntington’s disease (P < 1.0 × 10-28; Table S5 and Figure S23).

Figure 4

Differentially expressed genes (DEGs) in the brain among Vmat1 genotypes and predicted co-expressing modules

(A) The number of DEGs detected by pair-wise comparisons among the four genotypes. All DEGs were found in the amygdala (with none in the prefrontal cortex or striatum).

(B) Correlations between individual DEG expression levels for the WT vs. and Vmat1Ile/Ile comparison and composite anxiety scores from LD, EP, OF, and SI tests (see STAR Methods). Only genes with strong Spearman’s correlations (P < 0.1) are shown. The bands are 95% confidence intervals.

(C) Network dendrogram from co-expression modules based on the expression data of all 47 regional brain samples. Each branch represents an individual gene, and the colors below represent the module, correlation (ρ) with the behavioral phenotype (locomotor activity and anxiety), and the relative expression level in the amygdala across genotypes. The samples with asterisks are from 10-month-old male mice with behavioral data and were used to calculate the correlations between expression levels and behavioral phenotypes, and the others are from four-month-old male mice. The M1 module (shown in turquoise), showing negative correlations with the anxiety score, exhibited significant overrepresentation of the DEGs detected between WT and Vmat1Ile/Ile mice.

(D) Protein-protein interaction networks among the genes in M1. The DEGs detected between WT and Vmat1Ile/Ile mice are shown in red. The thickness of the line indicates the strength of data supports analyzed by STRING.

Table 2

List of gene ontology (GO) terms significantly overrepresented in the set of differentially expressed genes (DEGs) between Vmat1WT and Vmat1Ile/Ile mice

Direction	Pathways	Genes	Adj. Pval
Upregulated	Regulation of cell junction assembly	Ace, Cldn1, Enpp2	5.5 × 10-3
Downregulated	Locomotory behavior	Gng7, Gpr88, Ido1, Pde1b, Penk, Ppp1r1b, Rasd2	7.2 × 10-5
Downregulated	Behavior	Asic4, Foxp1, Gng7, Gpr88, Ido1, Pde1b, Penk, Ppp1r1b, Rasd2	4.9 × 10-4
Downregulated	Mononuclear cell differentiation	Cd4, Foxp1, Pde1b	2.8 × 10-3
Downregulated	Response to amphetamine	Pde1b, Ppp1r1b, Rgs9	2.8 × 10-3
Downregulated	Multicellular organismal response to stress	Asic4, Gng7, Ido1, Penk	2.8 × 10-3
Downregulated	Cardiac muscle tissue development	Foxp1, Myom3, Rarb, Rbp4, Wnt2	3.9 × 10-3
Downregulated	Regulation of cardiac muscle cell proliferation	Foxp1, Rbp4, Wnt2	4.0 × 10-3
Downregulated	Behavioral fear response	Asic4, Gng7, Penk	4.6 × 10-3
Downregulated	cAMP catabolic process	Pde7b, Pde10a	4.8 × 10-3

P-values were corrected by the Benjamini–Hochberg method. Only terms with adjusted P < 0.05 are presented. Downregulation indicates Vmat1WT > Vmat1Ile/Ile, whereas upregulation indicates Vmat1Ile/Ile > Vmat1WT.

To further characterize the biological pathways affected by the Vmat1 genotype, we performed a weighted gene correlation network analysis (WGCNA; Figure 4C). A co-expressing gene module (M1; Figure 4C) stood out as it was significantly overrepresented with DEGs detected in the WT vs. Ile/Ile comparison (odds ratio = 3.79, P = 2.88×10-8 by Fisher’s exact test) and was negatively correlated with composite anxiety score. These genes include many involved in adrenergic (Adra2b), dopaminergic (Adcy5, Adora2a, Drd1, Drd2, Gnal, Gng7, Gpr6, Gpr88, Pde1b, Pde1c, Pde7b, Pde10a, Pdyn, Penk, Ppp1r1b, Rasgrp2, Rgs9 and Tac1), glutamatergic (Grm4), and serotonergic (Htr1b, Htr1d, Htr1f, and Htr4) signaling pathways (Figure S24). Sub-network analysis further detected another functional gene cluster involved in neural development within the co-expressing module (cluster 2 in Figure 4D). Strikingly, three genes, Foxp1, Foxp2, or Six3, all of which are important regulators of neural growth, were in hubs connecting genes involved in postsynaptic signaling (cluster 1) as mentioned above and neural development (cluster 2; e.g., Dlx6, Ebf1, Isl1, and Wnt2).

Quantitative reverse transcription-PCR (qRT-PCR) was additionally conducted to examine the expression levels of amygdala marker genes and eliminate the possibility that a certain proportion of amygdala samples was derived from different nuclei or microstructures, which could cause the observed differential expression of genes. We confirmed that the expression levels of four genes (i.e., Cplx1, Negr1, and Nts and Prkcd, each representing basal (BA), lateral (LA), and central (CeA) nucleus of the amygdala, respectively) were not significantly different among samples (Figure S25A). Moreover, their expression levels were overall correlated with those revealed by RNA-seq (Figure S25B). These complement our results in transcriptomic analysis and indicate that not anatomical or structural but regulatory differences derived from Vmat1 genotypes were reflected in differential expression of genes detected above.

The humanized mVMAT1 disrupted anxiety-associated theta band activity in the basolateral amygdala

Behavioral and transcriptomic results suggest that the Vmat1 genotype influences anxiety-related neuronal activity patterns. Inspired by a previous report that theta (4–12 Hz) oscillations in the medial prefrontal cortex (mPFC) and basolateral amygdala (BLA) are specifically enhanced under anxiogenic environments such as the open arms of an elevated plus maze (EP) and the center area of an OF (Adhikari et al., 2010; Kuga et al., 2022; Likhtik et al., 2014), we next performed simultaneous local field potential (LFP) recordings from the dorsomedial prefrontal cortex (dmPFC) and BLA in freely moving mice during elevated plus maze (EP) exploration (Figure 5A; n = 5 Vmat1WT and n = 5 Vmat1Ile/Ile male mice, 16–20 weeks). The histological confirmation of electrode positions and raw LFP signals were shown in Figures 5B and 5C, respectively. Power spectral analysis of WT mice revealed lower 4–7 Hz power in both dmPFC and BLA during exploration of the EP open arms compared with close arms, suggesting a relationship with anxiety-like behavior (Figure 5D). Further, the Granger causality test revealed the directionality of 4–7 Hz oscillations from the mPFC to BLA (Figure 5E; n = 5 Vmat1WT mice), consistent with anatomical evidence that the dmPFC preferentially projects to the BLA (Bukalo et al., 2015; Gabbott et al., 2005) and that mPFC theta oscillations drive BLA neuronal activity (Likhtik et al., 2014). Based on these observations, we compared 4–7 Hz power in dmPFC and BLA between the WT and Ile/Ile genotypes. First, we applied a two-way ANOVA analysis on each brain region and found a significant difference in 4–7 Hz signal power between the two genotypes in both regions (main effect of genotypes, dmPFC: F1,19 = 4.95, P = 0.041; BLA: F1,19 = 8.34, P = 0.011). In addition, the power changes associated with open/closed arms were further tested in each genotype. WT mice exhibited significantly lower LFP 4–7 Hz signal power in the dmPFC and BLA under the anxiogenic condition of open arm exploration compared with closed arm exploration (dmPFC: P = 0.031; BLA: P = 0.031, Mann–Whitney U test followed by Bonferroni correction), whereas Vmat1Ile/Ile mice showed no such consistent activity pattern (dmPFC: P > 0.99; BLA: P > 0.99, Mann–Whitney U test followed by Bonferroni correction) (Figure 5F), suggesting that anxiety-related neuronal mechanisms in the dmPFC–BLA circuit are disrupted by the Vmat1Ile/Ile genotype.

Figure 5

Wild-type and Vmat1Thr/Thr mice, but not Vmat1Ile/Ile mice, exhibit a reduction in amygdalar 4–7 Hz local field potential (LFP) power under anxiogenic conditions

(A) An elevated plus maze (EP) test.

(B) (Left) LFP recordings were simultaneously performed from the dorsomedial prefrontal cortex (dmPFC) and basolateral amygdala (BLA). (Right) Histological confirmation of electrode locations in the dmPFC and BLA.

(C) Typical LFP signals from the dmPFC and BLA.

(D) Spectrograms of dmPFC (left) and BLA (right) LFP power in (anxiogenic) open arms relative to closed arms. Data were averaged from all Vmat1WT mice. The bar above indicates the 4–7 Hz band, showing pronounced decreases in LFP power in both regions.

(E) Spectral Granger causality averaged over 20 dmPFC–BLA electrode pairs.

(F) Comparison of dmPFC and BLA LFP 4–7 Hz power (z-scored) between open and closed arms. ∗P < 0.05, Mann–Whitney U test followed by Bonferroni correction. Each line represents one mouse.

Localization of Vmat1-expressing neurons and Vmat1 proteins in the CNS

Whereas the presence and functional involvement of Vmat1 in the CNS have been recently reported (Lohoff et al., 2014), its detailed expression profiles have not been fully elucidated. Given the transcriptomic and neurophysiological changes in the amygdala caused by the Vmat1 mutation, we hypothesized that the synaptic modulation via Vmat1-expressing neurons projecting to the basolateral amygdala led to the observed regulatory changes. To this aim, we sought to examine the localized expression of Vmat1 in those neurons. Four commercially available antibodies were first used to detect Vmat1 proteins in mouse brain tissues, but none of them were good enough to obtain feasible signals for synaptically localized Vmat1, probably owing to the general difficulty in recognizing highly lipophilic twelve-transmembrane transporters like Vmat1. To address this lack of specificity of antibodies, we generated two new knock-in mouse lines using genome editing techniques (Figure S26A).

To locate cell bodies of Vmat1-expressing neurons, we combined the Vmat1-T2A-iCre knock-in line with a prevailing Cre-reporter, Ai9 (Figures S26B and S26C). In addition, we employed the epitope-tagging strategy, in which proteins with no antibodies available can be detected by a highly specific monoclonal antibody against an epitope-tag, generating the Vmat1-3×HA knock-in line (Figures S26C and S26D). After confirming the reliability of these two knock-in mouse lines by using the adrenal medulla and hippocampal dentate gyrus as the positive controls (Figures S26E and S26F), we analyzed brain tissues from 6-week-old males. As a result, we detected some dopaminergic neurons in the ventral tegmental area (VTA) expressing Vmat1, visualized by Cre recombinase-dependent tdTomato reporters (Figure S26G). We also found a certain proportion of serotonergic neurons in the dorsal raphe nucleus (DR) expressing Vmat1 (Figure S26H). By utilizing the spinning-disk confocal microscopy technology at super-resolution, we further detected the co-localization of HA-tag signals with a synaptic vesicle marker, synaptophysin, in the BLA (Figure S26I). This indicates the possibility that the terminals of Vmat1-expressing monoaminergic neurons reach the BLA and could affect synaptic modulations.

Discussion

Whereas the potential importance of VMAT1 in neuroendocrine signaling has been suggested, few studies have focused on this gene until recently, largely owing to its relatively low expression levels in the CNS (Erickson et al., 1996; Peter et al., 1995; Schütz et al., 1998). As we have shown in this study, however, a certain proportion of monoaminergic neurons markedly express Vmat1 in the CNS. In addition, a series of studies have demonstrated the psychopathological effects of this gene, suggesting its importance in emotional regulation. The Thr136Ile polymorphism (rs1390938) is allegedly associated with bipolar disorder (Lohoff et al., 2006), autism spectrum disorder (Noroozi et al., 2017), cognitive function related to schizophrenia (Need et al., 2009), anxiety, depressiveness, neuroticism, or maladaptive impulsivity (Vaht et al., 2016). Our previous study uncovered from an evolutionary perspective that this variant is a human-specific polymorphism (as shown in Figure 1A) and under positive selection in the human lineage (Sato and Kawata, 2018). Two substitutions (i.e., Glu130Gly and Asn136Thr) occurred in the human lineage after the divergence from the common ancestor of chimpanzees and humans, and these substitutions combined have been shown to decrease the monoamine uptake of VMAT1 (Sato et al., 2019). On the other hand, the new “hyper-function” allele (Lohoff et al., 2014), which increases the monoamine uptake of VMAT1 and is associated with fewer psychopathological symptoms, by contrast, emerged just before the Out-of-Africa event of modern humans and is currently 20–61% in frequency in non-African populations (Sato and Kawata, 2018). We also found that the Thr136Ile polymorphism has been under balancing selection in non-African populations (Sato and Kawata, 2018). Given such evolutionary significance and psychopathological effects of the Thr136Ile polymorphism in VMAT1, this variant could play a very important role in exploring the evolutionary origins of human psychological traits and their diversity.

Phenotypic changes in mice caused by the humanized Vmat1

Our comprehensive behavioral tests show the protective effects of the Vmat1Ile allele against anxiety-like behavior (Figures 3B and 3C), which is well consistent with the human phenotype (Vaht et al., 2016). Compared with Vmat1Ile (Vmat1Thr/Ile and Vmat1Ile/Ile) mice, both wild-type and Vmat1Thr/Thr mice showed similar behavioral phenotypes and have higher levels of anxiety in both solitary (LD, EP, and OF) and social (SI and CSI) environment. This is expected from the molecular evidence that 136Asn and 136Thr are comparable in the level of monoamine uptake, which is lower than that of 136Ile (Sato et al., 2019). It is noteworthy that the genotypic effect was not observed for a fear response of mice (Figure S16). Given the potential difference in the neurological mechanism behind fear and anxiety (Davis et al., 2010; Tovote et al., 2015), Vmat1 may specifically regulate anxiety-related circuits of the brain as discussed later. The results of the delayed reward task may also suggest that the Vmat1Ile/Ile mice are more impulsive, at least initially, under this condition but contradicts an epidemiological study reporting an association of 136Thr with maladaptive impulsivity (Vaht et al., 2016). Kim et al. (2018) reported that the deletion of the D2 dopamine receptor (D2R) increased impulsive behavior in mice, whereas restoration of D2R in the amygdala alone normalized impulsive behavior (Kim et al., 2018). This is compatible with our results, as dopaminergic genes including Drd2 were downregulated in Vmat1Ile/Ile mice (Figure S24). The possible difference in the association of VMAT1 genotype with impulsivity between humans and mice may stem from functional connectivity changes in related neural circuits and warrants further investigation. Overall, these observations suggest that the 136Ile substitution may contribute to uniquely human behavioral traits, namely high exploratory tendency and relatively lower anxiety under threatening or novel conditions. Indeed, 136Ile allele frequency increases along the migration route of modern humans (Sato and Kawata, 2018), possibly suggesting that the exploratory tendency and/or boldness conferred by this mutation was advantageous for survival. In this regard, Vmat1Ile is analogous to the D4R polymorphism associated with novelty-seeking behavior (Matthews and Butler, 2011).

A large proportion of differentially expressed genes was detected between Vmat1WT and Vmat1Ile/Ile, and most of them were downregulated in Vmat1Ile/Ile mice. Moreover, these genes were significantly overlapped and highly correlated in expression patterns with those differentially expressed in Huntington’s disease (HD) model mice (Table S5, Figure S23). It is possible that the Ile mice recapitulate some dimensions of HD pathology, given that dopamine (DA) signaling likely plays a significant role in HD, and the only current FDA-approved drug for HD is a VMAT2 inhibitor, tetrabenazine (Chen et al., 2013). Notably, the expression levels of DEGs in Vmat1Ile/Ile mice were negatively correlated with those in Snca-KO mice known to be a model of Parkinson’s disease (PD) (Table S5). As up- and downregulation of dopaminergic neurotransmission have been considered as one of the key mechanisms of HD and PD, respectively (Höglinger et al., 2004; Jahanshahi et al., 2010), our disease correlative implications might provide the basis for pharmaceutical and medical application of Vmat1-humanized mice.

WGCNA detected a co-expressing module (M1; Figure 4C) in which the DEGs were significantly overrepresented and revealed that many genes involved in monoaminergic signaling pathways were co-expressed with the detected DEGs (Figure 4D). Lohoff et al. reported that Vmat1 KO affected DA signaling in particular, with upregulation of D2R and downregulation of tyrosine hydroxylase (TH) observed in both the frontal cortex and striatum (Lohoff et al., 2019). Given the association between Vmat1 KO and decreased DA in the frontal cortex, we speculated that Vmat1Ile/Ile mice would also exhibit downregulation of D2R and upregulation of TH compared with WT mice. Indeed, D2R (along with many other dopaminergic pathway genes) was downregulated in Vmat1Ile/Ile mice compared to WT mice. However, there were no significant differences in TH expression among genotypes (Figure S24), suggesting that the Vmat1 Ile allele alters dopaminergic transmission but not DA synthesis. The sub-network analysis further detected another functional gene cluster involved in neural development within the co-expressing module (cluster 2 in Figure 4D), such as Foxp1, Foxp2, Six3, Dlx6, Ebf1, Isl1, and Wnt2. Foxp1 (forkhead box protein P1) is a transcription factor that regulates the development of various tissues including the brain. Foxp1 heterodimerizes with Foxp2, and a human-specific Foxp2 substitution is associated with language impairment (Enard et al., 2009, 2002). Foxp2 also regulates D2R expression (Spiteri et al., 2007) and extracellular DA levels (Enard et al., 2009), potentially involved in neural development (Vernes et al., 2011). Dlx6, Ebf1, Isl1, and Six3 are all implicated in the development of striatal neurons that express DA receptors (Ehrman et al., 2013; Wang et al., 2011; Xu et al., 2018). The DEG Wnt2 (also known as Irp) has been found to increase dopaminergic neurons by inducing the proliferation of progenitor cells in the developing midbrain (Sousa et al., 2010). The differential expression of these genes may also be related to the regulatory role of VMAT1 in hippocampal neurogenesis (Multani et al., 2013), which emerging evidence suggests is associated with anxiety-like behavior (Anacker and Hen, 2017; Revest et al., 2009). Taken together, the differential regulation of DA signaling by Vmat1 genotypes may affect widespread-cascading postsynaptic machinery on neural development and plasticity.

The electrophysiological analysis revealed abnormal neuronal activity in the amygdala of Vmat1Ile/Ile mice under fearful conditions (i.e., open arms in EP test; Figure 5F). This was not seen in the medial prefrontal cortex, and thus suggests that amygdalar neuronal disturbance underlies the reduced anxiety-like behavior associated with this genotype. This is consistent with a functional imaging study showing that 136Ile was associated with increased reactivity and decreased habituation of the amygdala toward threat-related stimuli in humans (Lohoff et al., 2014). Additional studies are required to reveal the precise association between these LFP changes and amygdalar output in Vmat1Ile mice. In light of our behavioral results, we suggest that Vmat1Ile may promote the activity of neurons that dampen anxiety, such as those projecting from the BLA to the central amygdala (Janak and Tye, 2015; Tye et al., 2011).

Transcriptomic and neuronal changes in the amygdala associated with anxiety imply intra- and inter-species differences in emotion and social behavior

One of the major findings of this study is that the Vmat1 genotype affected transcriptomic regulation and neuronal activity primarily in the amygdala. Given that amygdalar DA signaling is a powerful modulator of fear and anxiety (de la Mora et al., 2010), the alterations in gene expression and oscillatory neuronal activity possibly contribute to the observed reduction in anxiety-like behaviors among Vmat1Ile mice compared with WT mice. The abnormal neuronal activity pattern significantly detected in the amygdala of Vmat1Ile mice may also be involved in a difference in VMAT1 expression levels across the brain. Previous studies verified that VMAT1 mRNA and protein expression levels are relatively high in the amygdala compared with other brain regions in both rats and humans (Hansson et al., 1998; Lohoff et al., 2006). Although it is necessary to consider measurement variability owing to its low expression, Vmat1 was expressed highly in the amygdala compared with the striatum (P = 0.0097, Dunnett’s test; Figure S27), whereas Vmat2 expression was comparable across regions (Figure S27). This may indicate distinctive regional regulation between subtypes as shown in a previous study (Multani et al., 2013). Our additional knock-in mouse lines validated the localization of Vmat1-expressing serotonergic and dopaminergic neurons in the DR and VTA, respectively (Figure S26). As these neurons are considered to project to the BLA, altered monoaminergic transmission via these neurons could lead to the differential transcriptomic or neuronal regulation observed in the amygdala. Incidentally, we detected no transcriptional changes in dopamine synthesis as represented by the expression of TH, which is different from the phenotype of Vmat1-KO mice (Lohoff et al., 2019). On the other hand, both WT and Vmat1-KO mice in the present and previous study showed relative upregulation in downstream DA pathway genes (e.g., Drd2; Figure S24). This can happen because Vmat1WT proteins are still functional unlike Vmat1-KO and do not necessarily lead to a complete lack of transported dopamine, minimizing the feedback effects activated by the null allele.

Neuronal circuits within the amygdala store associations between environmental cues and adverse stimuli as changes in synaptic strength, and these modified circuits drive the expression of emotion-related behaviors and physiological responses. Thus, the amygdala is a critical regulator of emotional and social behaviors (Murray, 2007), and amygdala dysfunction is indeed implicated in many neuropsychiatric disorders (Avino et al., 2018; Chang et al., 2017; Posner et al., 2016). Although it is an evolutionarily primitive brain region with relatively well conserved gross structure and general function across species (Janak and Tye, 2015), recent studies have provided evidence for more subtle structural and/or functional differences within and among species. The size of the amygdala has been shown to correlate with creativity (Asari et al., 2010), mental state inference (Rice et al., 2014), and the size and complexity of social networks in humans (Bickart et al., 2011) and non-human primates (Sallet et al., 2011). The DA response in the medial amygdala network is associated with mother-infant bonding (Atzil et al., 2017). Collectively, these findings strongly suggest a substantial contribution of the amygdala to the evolution of the human social brain by acting as a hub among brain networks associated with emotion, cognition, and communication (Bickart et al., 2014).

A single substitution in humanized mouse models highlights molecular and neurological evolution underlying human emotional traits

Genetic deletion (KO) and overexpression have become predominant strategies for examining the neurological, behavioral, and pathogenic functions of specific genes. However, gene deletion and overexpression may alter the expression levels of many additional genes, and thus the observed phenotype may be distinct from that induced by target-specific pharmacological manipulations. The present study employed an alternative “knock-in” strategy for non-synonymous Vmat1 mutations and succeeded in capturing moderate differences in transcriptomic, neuronal, and behavioral phenotypes. As discussed earlier, our knock-in alleles confer quantitative differences in transported monoamine and exhibited different regulatory effects on the transcriptomic landscape than the null allele. These results exemplify the strength of this strategy and provide keen insight into the functional significance of natural mutations.

Noticeably, the variant evaluated is a human-specific substitution under a selection not possessed by other mammals (Figure 1A). This knock-in model thus provides a unique opportunity to examine the molecular and neurological mechanisms that distinguish human behavior from that of other primates. Such extensive functional analysis of the “humanized” brain has been limited to a few genes such as FOXP2 (Enard et al., 2009; Schreiweis et al., 2014) or ARHGAP11B (Xing et al., 2021) and never for effects on emotional traits (Zhu et al., 2019). Therefore, the present study is the first to suggest a new experimental strategy for studying the evolution of human brain mechanisms underlying emotion and related behavioral traits.

Limitations of the study

Lastly, it needs to be noted that recent genome-wide association studies (GWAS) on psychiatric disorders such as depression (Wray et al., 2018) and schizophrenia (Ripke et al., 2014), and specific personality traits (Lo et al., 2017) have not detected VMAT1 as a top hit despite its alleged associations with these phenotypes under study (Lohoff et al., 2008b; Vaht et al., 2016). This may indicate that VMAT1 has a weaker influence on phenotype than other significant loci, or that this gene must interact with other loci (epistasis) and/or environmental factors (G × E interaction) for measurable effects on phenotype. Such dependence on other genes or external factors could obscure phenotypic associations of VMAT1 evaluated by conventional approaches of GWAS. We found that a single amino acid substitution (Asn133Ile) altered the expression of 56 genes within the amygdala, any or all of which may contribute to the observed changes in behavioral phenotype. The effects of various environmental factors (which were largely controlled in this experiment through group housing and shared dams) and DEGs on these behavioral and neurological phenotypes would warrant further study.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources used in this study should be directed to and will be fulfilled by the lead contact, Masakado Kawata (kawata@tohoku.ac.jp).

Materials availability

Mouse lines generated in this study will be made available from TI (tinoue@ncnp.go.jp) with a material transfer agreement.

Experimental model and subject details

C57BL/6J, B6C3F1, and B6.Cg-Gt(ROSA)26Sortm9(CAG-tdTomato)Hze/J (also referred to as Ai9) mice were purchased from Charles River, Japan, SLC Japan, and The Jackson Laboratory, respectively. One-cell-stage embryos were used for genome editing procedures. Animals were kept in clean conditions, and males were used for experiments. For Vmat1-humanized mice, behavioral tests were conducted between 9 and 55 weeks of age (Table S3). A subset of these mice (1st batch) was also used for transcriptome analyses at 42–45 weeks. Additional sets of mice (2nd batch) were generated and used for transcriptome analysis at 15 weeks and electrophysiological recordings at 16–24 weeks (Figure S3). For Vmat1-T2A-iCre and Vmat1-3×HA knock-in mice, immunostaining of adrenal glands and brains was carried out at 6 weeks. All animal care protocols and experiments were performed according to guidelines approved by the Animal Care and Use Committee of the National Institute of Neuroscience, National Center of Neurology and Psychiatry (approval numbers: 2017005, 2020007), the Animal Experiment Committee of Fujita Health University (approval number: APU19063), and the Experimental Animal Ethics Committees of Graduate School of Medicine, and the University of Tokyo (approval number: P29-14).

Method details

Multiple sequence alignments and phylogenetic tree of VMAT1

Two hundred sixty-three orthologous sequences of vertebrate SLC18A1 (VMAT1) were obtained from Ensembl101. After aligning the sequences using webPRANK (Löytynoja and Goldman, 2010) with default parameters, we calculated the proportions of amino acids identical to the human VMAT1 within aligned protein regions. Sequences with < 50% identity were then excluded. In the case of multiple sequences for a given species, that with the highest proportion of identical amino acids to human VMAT1 was selected. A phylogenetic tree was then constructed from the 236 species sequences with > 50% identity to human VMAT1 using IQTREE v2.1.1 (Minh et al., 2020) with the “-nt AUTO” setting, including the Neanderthal sequence, which we constructed by replacing 136Ile in the human reference sequence with 136Thr. The phylogenetic relationship was visualized using iTOL v5.6.3 (Letunic and Bork, 2019) through a web interface (https://itol.embl.de). Animal silhouettes were obtained from phylopic (http://phylopic.org/) under a Public Domain license.

In silico prediction of mVMAT1 structure and tolerance of mutated residues

Homology modeling of the mouse VMAT1 (mVMAT1) protein structure was performed using the SWISS-MODEL (Waterhouse et al., 2018) web server (http://swissmodel.expasy.org). The human organic anion transporter MFSD10 (PDB: 6s4m.1.A) was selected as the best template, with a sequence identity of 17.04 (verified June 2020) and QMEAN score of –8.40 according to SWISS-MODEL. The predicted 3D structure of mVMAT1 was visualized using PyMOL 2.4.1 (DeLanoScientific, San Carlos, CA). Provean v1.1.3 (Choi and Chan, 2015) and SIFT (Kumar et al., 2009) were used to estimate the intolerance for individual amino acid mutations introduced in mVMAT1 based on the evolutionary conservation and the chemical properties of the exchanged residues. The computations were made using the Provean web server (http://provean.jcvi.org/protein_batch_submit.php?species=mouse; last accessed June 2020). A mutation is predicted as deleterious (i.e., likely to affect its protein function) when the Provean score is smaller than −2.5 or the SIFT score is smaller than 0.05. Furthermore, we used DynaMut (Rodrigues et al., 2018) to evaluate the effects of humanized mutations on the stability and flexibility of VMAT1 protein structure. Briefly, DynaMut enables an accurate assessment of mutation impact on protein stability by implementing and integrating well established normal mode approaches with graph-based signatures in a consensus predictor (Rodrigues et al., 2018). Here we focused on two measurements: 1) folding free energy (ΔΔG), an index of the difference in stability between WT and mutant proteins in which positive values represent increased stability of the mutant, and 2) Δ vibrational entropy energy (ΔΔSVib ENCoM), a per-site index of the difference in flexibility between WT and mutant proteins. Calculations were conducted using a web server (http://biosig.unimelb.edu.au/dynamut/).

Generation of Vmat1-humanized mouse models by CRISPR/Cas9 genome editing

Preparation of the CRISPR components and donor DNAs

The target region within the mouse Vmat1 exon 4 was first analyzed by using the web-based CRISPR design tool CRISPOR (Concordet and Haeussler, 2018) (http://crispor.org), to select a guide RNA sequence for replacing mouse 133Asn with 133Thr or 133Ile via the CRISPR/Cas9 system with minimal risk of off-target cleavages. Then, two parts of the CRISPR guide RNA, crRNA and tracrRNA (shown in Figures 2A, S1A, and Table S1), were chemically synthesized and purified by reversed-phase column chromatography (FASMAC, Japan). Recombinant Cas9 protein (EnGen Cas9 NLS from Streptococcus pyogenes) was purchased from New England Biolabs (Ipswich, MA). The single-strand DNA (ssDNA) donors containing the humanized substitution and the restriction enzyme recognition sites flanked by 42–50 bases homology arms (Table S1) were subsequently designed and chemically synthesized by Eurofins Genomics (Tokyo, Japan). In addition to the intended humanization, the synonymous substitutions depicted in Figure 2A were introduced near the 133rd site to prevent unwanted re-editing. For 133Thr substitution, an EcoRI site was synonymously integrated into the protospacer adjacent motif (PAM) sequence to be destroyed. An FspI site was similarly incorporated within the guide RNA target sequence for the 133Ile allele. These synonymous substitutions also made the genotyping procedure much easier by eliminating the need for laborious sequencing (Figure 2B).

Electroporation of mouse one-cell embryos

One-cell-stage embryos were collected by in vitro fertilization using C57BL/6J super-ovulated females (Charles River, Japan) and stud males. The guide RNAs, Cas9 proteins, and the donor ssDNAs were electroporated into mouse zygotes following the standard protocol (Hashimoto et al., 2016) to obtain founder knock-in mice (Figure S1A).

PCR Genotyping and sequencing analyses

Genomic DNAs were prepared from the newborn mouse tail by treatment with proteinase K in Lysis Buffer. The 133Thr and 133Ile knock-in founders were screened by PCR-RFLP assay (Figure S1C). Briefly, PCR products were amplified using the primer pairs depicted in Figure 2A (also listed in Table S1) and then digested by EcoRI or FspI, respectively. For founders and F1 generations, PCR products were also analyzed by Sanger sequencing to confirm the correct substitution allele. After that, PCR-RFLP could be used reliably for genotyping of subsequent generations.

Off-target analyses of founder mice

To exclude the possible side effects from off-target cleavages, the genomic DNAs extracted from the founder mice (133Thr No.1 and 133Ile No.10, shown in Figures S1C and S1D) were analyzed by Sanger sequencing. We employed CRISPOR (Concordet and Haeussler, 2018) (http://crispor.org) to predict potential off-target candidate loci, and the 12 loci identified (Table S2) were amplified by PCR using the primers listed in Table S1 and sequenced (results are shown in Figure S2).

Preparation of humanized mouse models for behavioral tests, transcriptome analyses, and neurophysiological analyses

After confirming the absence of off-target cleavages, 133Thr founder No.1 and 133Ile founder No.10 were crossed with WT C57BL/6J mice to obtain heterozygous F1 generations, which were again verified to possess the designed substitutions by Sanger sequencing. Homozygous 133Thr and 133Ile mice after F2 generations were then selected and maintained. For behavioral tests, F5 or F6 homozygous mice were crossed to obtain Thr/Thr, Thr/Ile, and Ile/Ile genotypes. In addition, C57BL/6J WT males and females carrying the Asn/Asn genotype were crossed to supply control mice. To eliminate differences in rearing environment, newborn males of the 4 genotypes (Asn/Asn, Thr/Thr, Thr/Ile, and Ile/Ile) were grouped in sets and nursed by the same mothers. At 4 weeks old, 20 sets of 4 males, one of each genotype, were weaned and housed in 20 separate cages. Behavioral tests were conducted between 9 and 55 weeks of age as indicated (Table S3). A subset of these mice (1st batch) was also used for transcriptome analyses at 42–45 weeks. Likewise, additional sets of mice (2nd batch) including all 4 genotypes were generated and used for transcriptome analysis at 15 weeks and electrophysiological recordings at 16–24 weeks (Figure S3).

Generation of Vmat1-T2A-iCre and Vmat1-3×HA knock-in mice by CRISPR/Cas9 genome editing

To insert functional gene cassette sequences in-frame and 5’ to the stop codon of the Vmat1 gene’s 17th exon via the CRISPR/Cas9 system (depicted in Figure S26A), the guide RNA sequences closest to the stop codon TAG were selected. Two parts of the CRISPR guide RNA, crRNA and tracrRNA (listed in Table S1), were chemically synthesized (Integrated DNA Technologies; IDT) and recombinant Cas9 protein (Alt-R S.p. HiFi Cas9 V3) was purchased from IDT.

For Vmat1-T2A-iCre knock-in mice, T2A self-cleaving peptide sequences were employed to bicistronically express codon improved Cre recombinase (iCre) in Vmat1-expressing cells. Targeting vector containing T2A-iCre sequences flanked by 1 kb homology arms on both sides was constructed via the In-Fusion cloning system (Takara Bio). The mixture of annealed crRNA/tracrRNA, Cas9 protein, and targeting vector was injected into mouse zygotes (B6C3F1, SLC Japan) at concentrations of 50, 100, and 15 ng/μl, respectively, and the zygotes were transferred into the oviduct of surrogate mothers to raise pups. Knock-in mice were screened by PCR using the four primers listed in Table S1. For knock-in candidates, PCR products were analyzed by Sanger sequencing to confirm the correct insertions. Obtained founders were crossed with homozygous B6.Cg-Gt(ROSA)26Sortm9(CAG-tdTomato)Hze/J (also referred to as Ai9) to select double heterozygous F1 males (Vmat1 ::Gt(ROSA)26Sor ), and their brains and adrenal glands were sampled at 6 weeks for immunostaining.

To insert three tandem copies of HA-tag sequences just upstream from the stop codon, the ssDNA donors containing 3×HA-tag sequences flanked by 36 bases homology arms on both sides were artificially synthesized (eurofins Genomics, listed in Table S1). The guide RNAs, Cas9 proteins, and ssDNA donors were electroporated into mouse zygotes (B6C3F1, SLC Japan) following the standard protocol (Hashimoto et al., 2016) to obtain founders. Knock-in mice were screened by PCR using the two primers listed in Table S1. For knock-in candidates, PCR products were analyzed by Sanger sequencing to confirm carrying designed insertions. The founders were intercrossed to obtain homozygous Vmat1-3×HA males (Vmat1 ), and their brains and adrenal glands were sampled at 6 weeks for immunostaining.

Immunohistochemistry

Primary antibodies were used at the following dilutions: rabbit anti-HA-tag (1:1000), chicken anti-tyrosine hydroxylase (1:2000), sheep anti-tryptophan hydroxylase (1:2000), and mouse anti-synaptophysin (1:1000). Alexa Fluor-conjugated secondary antibodies were used at dilutions of 1:600.

Brains and adrenal glands were intracardially perfused with 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS; pH 7.4), post-fixed with the same fixative for 16-24 h, and cryo-protected with 30% sucrose in PBS, then embedded in an OCT compound (Sakura Finetek Japan) to be frozen. Brains were sliced into 40 μm coronal sections and adrenal glands were sliced into 12 μm sections using a Leica cryostat (CM3050). Sections were treated in blocking buffer (10% normal donkey serum, 0.2% Triton X-100 in PBS) for 1 h at room temperature with gentle agitation. Sections were then incubated with primary antibodies in a blocking solution (3% normal donkey serum, 0.2% Triton X-100 in PBS) overnight at 4°C. After being washed with PBS, sections were treated with appropriate Alexa Fluor-conjugated secondary antibodies and DAPI for 1 h at room temperature. After washing, sections were mounted onto glass slides with VECTASHIELD Vibrance antifade mounting medium (Vector Laboratories). Fluorescent images were obtained with a KEYENCE BZ-X710. For confocal and super-resolution imaging, an OLYMPUS SpinSR10 was used to take z-stacked images.

Behavioral battery tests

Animal care and experimental conditions for behavioral tests

Animals were housed in groups of four individuals of each genotype under a 12 h light/dark cycle (7:00 AM to 7:00 PM) with ad libitum access to food and water. Adult male mice were used in all tests to eliminate the behavioral effects of the estrus cycle. The composition of cohorts and ages of the mice for every experiment are summarized in Table S3. All behavioral tests were conducted in a soundproof room, and as much effort as possible was spent on controlling the effects of confounding factors such as light intensity, temperature, and humidity for the tested group of cage mates. Furthermore, to minimize the impact of one test on subsequent tests in individual mice, behavioral assessments were conducted in the following order after general health and neurological screens (GHNS): light/dark transition (LD), open field (OF), elevated plus maze (EP), hot plate (HP), social interaction in a novel environment (SI), rotarod (RR), Crawley’s 3-chamber social interaction test (CSI), startle response/prepulse inhibition (PPI), Porsolt forced swim (PS), T-maze (TM), Barnes maze (BM), tail suspension (TS), fear conditioning (FZ), and home cage social interaction tests (HCSI).

General health and neurological screens (GHNS)

The presence of whiskers and bald patches was checked daily to ensure the health status of mice (Miyakawa et al., 2001). The righting, whisker touch, and ear twitch reflexes were also evaluated to ensure age-appropriate neurological function. Body weight and rectal temperature were measured. In addition, the neuromuscular function was assessed by grip strength and wire hanging tests. Briefly, grip strength was measured using a grip strength meter (O’Hara & Co., Tokyo, Japan). The mouse was positioned to spontaneously grasp a wire grid by the forelimbs and then pulled backward by the tail until wire release. The peak force was recorded in Newtons (N). Each mouse was tested three times, and the greatest value obtained was used for further analyses. In the wire hanging test, a mouse was placed on a wire mesh at the top of the apparatus (O’Hara & Co.), and the wire mesh was then gently turned upside down. The mouse gripped the wire in order not to fall off, and the latency to fall was recorded with a 60 s cut-off time.

Light/dark transition (LD)

The light/dark transition test was conducted to measure anxiety levels as previously described (Takao and Miyakawa, 2006). The apparatus consisted of two equal-sized plastic boxes (20 × 20 cm), one illuminated (390 ± 20 lx) and the other dark (< 2 lx), separated by a central partition plate with a small 3 × 5 cm opening allowing the mouse to transit from one box to the other. A mouse was first placed in the dark box and allowed to freely explore the apparatus for 10 min. The time to first entry into the light box, time spent in the light box, number of transitions between boxes, and distances traveled in light and dark boxes were measured using Image LD software (see “Image analysis of behavioral tests”).

Open field (OF)

Locomotor activity and explorative tendency were measured using an open field apparatus (40 × 40 × 30 cm; Accuscan Instruments, Columbus, OH) as described previously (Takao et al., 2010). The test chamber was illuminated at 100 ± 5 lx. A mouse was placed in the corner of the apparatus and total distance traveled (cm), time spent in the center area (20 × 20 cm), vertical activity, and stereotypic counts over 120 min were recorded by the VersaMax system and analyzed using Image OF software (see “Image analysis of behavioral tests”).

Elevated plus maze (EP)

An elevated plus-maze test was conducted as previously described (Komada et al., 2008). The apparatus consisted of two open arms (25 cm long × 5 cm wide) crossing two enclosed arms of identical dimensions but with 15-cm high transparent walls, all connected by a central platform. The entire apparatus was elevated 50 cm above the floor and illuminated uniformly at 100 lx. A mouse was placed at the center of the maze, facing one of the closed arms, and allowed to freely explore for 10 min. The number of open and closed arm entries, distances traveled, proportions of entries into open arms, and proportions of time spent in open arms were measured using Image EP software (see “Image analysis of behavioral tests”).

Hot plate (HP)

Pain sensitivity was measured by the hot plate test. A mouse was placed on a 55°C hot plate (Columbus Instruments, Columbus, OH), and latency to escape from the plate was recorded with a 15-s cut-off time.

Social interaction in a novel environment (SI)

The social interaction test in a novel environment test (one-chamber social interaction test) was performed as described previously (Miyakawa et al., 2001). Two mice of the same genotype but reared in different cages were placed together in a box (40 × 40 × 30 cm) and allowed to explore freely for 10 min. Mouse behaviors were recorded with a CCD camera. The total duration of contact, the total number of contacts, total duration of active contact, mean duration per contact, and total distance traveled were measured automatically using Image SI software (see “Image analysis of behavioral tests”). Active contact was defined as maintained contact while either mouse traveled more than 10 cm between two successive image frames acquired once per second.

Rotarod (RR)

Motor coordination and motor learning were examined using an accelerating rotarod (UGO Basile, Comerio, VA, Italy). The latency to fall from the rod was recorded during three daily trials conducted on 2 consecutive days (3 × 2 trials). During each trial, rod speed was increased from 4 to 40 rpm over the 5-min test period.

Crawley’s 3-chamber social interaction (CSI)

Crawley’s 3-chamber social interaction test was performed to measure sociability and social novelty preference (Takao et al., 2010). The experimental apparatus was a 41 × 62 cm rectangular non-transparent gray Plexiglas box separated into three equal-sized chambers (20 × 40 cm) by transparent Plexiglas plates with small openings allowing mice to freely transit from one chamber to another. One wire cylinder-shaped cage was placed in the corner of each side chamber. Tested mice (aged 16–19 weeks, see Table S2) were first allowed to freely explore the chambers for 10 min as a habituation period. A stranger mouse was then placed randomly in one of the side-chamber wire cages and the tested mouse was again allowed to freely explore the chambers for 10 min. In a final step, a second stranger mouse was placed in the wire cage in the opposite side chamber, and the tested mouse was again allowed to freely explore the chambers for 10 min. The number of contacts and interaction times with strangers 1 and 2 were recorded using Image CSI software (see “Image analysis of behavioral tests”).

Startle response / prepulse inhibition (PPI)

Sensory-motor integration was tested by a prepulse inhibition test (Miyakawa et al., 2001), which measures the effect of an audible pre-tone (prepulse) of variable intensity on the startle response to a louder pulse. The system used for the detection of startle reflexes (O’Hara & Co.) consisted of four standard cages placed on a movement sensor within sound-attenuated chambers with fan ventilation. Before each testing session, acoustic stimuli and mechanical responses were calibrated using devices supplied by the manufacturer. A test session was composed of 6 different trial types: two startle stimulus-only trials (110 or 120 dB) and four prepulse inhibition trials (74 or 78 dB prepulses delivered prior to a 110 or 120 dB startle). Six blocks of the six trial types (i.e., 36 trials in total) were conducted, with each trial type presented once in pseudorandom order within each block. The PPI (%) was calculated for each trial type according to the following equation: {[(startle amplitude of trial without prepulse) – (startle amplitude of trial with prepulse)]/(startle amplitude of trial without prepulse)} × 100.

Porsolt forced swim (PS)

The Porsolt swim test was performed as previously described (Miyakawa et al., 2001). The apparatus consisted of four Plexiglas cylinders (12 cm diameter × 22 cm height) filled with water (room temperature) up to a height of 7.5 cm. Four mice (a set) were placed individually in each cylinder, and images were captured at two frames per second for 10 min. Swimming distance and immobility time (% of total) were recorded automatically using Image TS software (see “Image analysis of behavioral tests”).

Tail suspension (TS)

Mice were suspended 30 cm above the floor in a visually isolated area by adhesive tape placed approximately 1 cm from the base of the tail. Behavior during suspension was recorded for 10 min using Image TS software (see “Image analysis of behavioral tests”).

T-maze (TM) forced alteration

A T-maze forced alternation task was performed as previously described (Shoji et al., 2012). The T maze was constructed of three white plastic runways with walls 25-cm high (O’Hara & Co.) partitioned into 6 areas by sliding doors. The stem of the “T” was divided into a start compartment area S1 at the intersection with the arms and the main area S2 (13 × 24 cm). Similarly, the T arms were divided into passageway areas P1 and P2 with S1 and main areas A1 and A2 (11 × 20.5 cm). Each trial consisted of one forced run followed by one free run. First, the mouse was placed in S2 and forced to proceed in one direction (S2 → A2 → P2 → S1, or S2 → A1 → P1 → S1) by opening the doors successively. After returning to S1 and a set inter-trial delay, mice were then allowed to choose to enter A1 or A2. If the mouse chose the opposite direction from the previous forced run, it was counted as a correct response, while if the mouse chose the same direction as the previous run, it was counted as an incorrect response. After 30 practice trials (10 per day for 3 consecutive days), trials with delay (4 trials with 3-s delays and 2 trials each with 30-, 60-, and 120-s delays) were conducted for 3 consecutive days. Data were acquired and sliding doors were controlled automatically using Image TM software (see “Image analysis of behavioral tests”).

Barnes maze (BM)

Spatial learning, spatial memory, and behavioral flexibility were tested using the Barnes maze as described previously (Takao et al., 2010). The maze was a circular whiteboard (1m in diameter) with 12 holes equally spaced around the margin and elevated 81.5 cm from the floor (O’Hara & Co.). A black plastic box (17 × 13 × 7 cm) lined with paper cage bedding was positioned under one of the holes (the target). The board position of the target (1–12) was constant for a given mouse but randomly assigned across individuals. The apparatus was illuminated at 800 lx or more. The board was rotated daily so that the spatial location of the target changed relative to distal visual room cues while proximal cues were held constant. A training trial was conducted per day over 20 days until 24 trials were completed in total before the 1st probe test. Each trial lasted a maximum of 5 min and was completed when the mouse entered the black plastic box. Probe tests were conducted for 3 min without the black plastic box one day or one month after the final training trial. The number of errors, latency and distance traveled to reach the target hole for the first time, the number of omission errors, and time spent around each hole were recorded by Image BM software (see “Image analysis of behavioral tests”).

Fear conditioning (FZ)

Contextual and cued fear conditioning tests were conducted and analyzed as described in a previous study (Miyakawa et al., 2001). Mice were placed in a transparent acrylic chamber (26 × 34 × 33 cm, O’Hara & Co.) with 55 dB ambient white noise and 100 lx illumination. In the conditioning session, the conditioned acoustic stimulus (CS, 55 dB) was presented for 30 s at three times (2, 4, and 6 min), and a mild foot shock (2 s, 0.3 mA) was presented as the unconditional stimulus (US) at the end of each CS. One day after the conditioning session, contextual fear was measured in the same chamber, while cued fear was measured thereafter in a triangular box (33 × 33 × 33 cm) constructed of white opaque Plexiglas. Mice were placed in the chamber for 3 min with neither CS nor US presented. Thereafter, the CS (55 dB) was presented for the last 3 min and images were captured at 1 frame per second. Each pair of successive frames in which the mouse moved was measured. When this measure was below a certain threshold (20 or 30 pixels), the mouse was considered “freezing”; alternatively, when this measure equaled or exceeded the threshold, the behavior was considered “non-freezing.” “Freezing” that lasted less than 2 s was not included in the analysis. Data acquisition, control of stimuli (i.e., tones and shocks), and data analysis were conducted automatically using Image FZ software (see “Image analysis of behavioral tests”).

Home cage social interaction test (HCSI)

The 24-hour home cage locomotor activity and social interaction test was conducted as previously described (Miyakawa et al., 2003) for one continuous week. The monitoring system included a standard home cage (29 × 19 × 13 cm) with a filtered cage top and an infrared video camera (O’Hara & Co.). Two mice of the same genotype that had been housed separately were placed together in the home cage, and locomotor activity and social behavior were monitored by video. Social interaction was measured by counting the number of “particles” (tracers for each mouse) detected in each frame, with two separate particles indicating no contact and one particle indicating physical contact (i.e., because the tracking software could not distinguish two separate bodies). Data acquisition and analysis were conducted using Image HA software (see “Image analysis of behavioral tests”).

Delayed reward task

Following these comprehensive behavioral tests, we conducted an original delayed reward task to evaluate the impulsivity of Vmat1Thr/Thr and Vmat1Ile/Ile mice (n = 10 for both genotypes). In brief, a mouse was housed in a cage (22.5 × 32.5 × 21 cm; FDB-300W/FDL-8D, Melquest, Toyama, Japan) with controlled access to two separate feeding trays for 16 days. Consumption of normal food (3.57 kcal/g; CRF-1, Oriental Yeast Co., Ltd., Tokyo, Japan) and high-calorie food (Banana chip containing 5.42 kcal/g; 4515996091582, KALDI COFFEE FARM, Tokyo, Japan), considered small and large rewards, respectively, were measured daily by weighing the feeding trays. After habituating mice to this system for 4 days with ad libitum access to both foods, we programmed feeding tray access so that the normal food was provided for the first 45 min and high-calorie food was provided only for the last 15 min of the three 1-h daily feeding time (7:00–8:00 PM, 10:00–11:00 PM, and 1:00–2:00 AM) for 6 consecutive days. Under such delayed reward conditions, consumption of the lower-calorie food is considered an impulsive behavior (i.e., does not maximize reward). Furthermore, to maximize the impulsivity, we set a fasting day after the 6 days of delayed feeding while monitoring animal weights (so as not to jeopardize the animal’s health) and conducted the same test again for 5 days.

Image analysis of behavioral tests

We developed in-house programs to analyze the imaging data acquired for many of these behavioral tests (Image LD, Image OF, Image EP, Image SI, Image CSI, Image TS, Image TM, Image BM, Image FZ, and Image HA) based on ImageJ (U.S. National Institutes of Health; available at https://imagej.nih.gov/ij/) and modified for each test by Tsuyoshi Miyakawa (available through O’Hara & Co.).

Structural Equation Modeling (SEM)

We obtained comprehensive metrics for two major domains of mouse behavior: activity and anxiety, by standardizing, normalizing, and combining sets of related behavioral scores. Total distances traveled in LD, OF, EP, SI, and CSI were used as indices for activity, while time in the light compartment of the LD shuttle box, duration in the central area of the OF, time spent in open arms of the EP, and total duration of active contacts in the SI were used as indices of anxiety. To further confirm the effects of the Vmat1 genotype on the behavioral composites and to enhance statistical power and improve reliability, we used Structural Equation Modeling (SEM) (Kline, 2016). The SEM model included a measurement model and a regression model. The measurement model consisted of two latent factors grouping tests of locomotor activity and anxiety-like behavior, which were assessed by the same metrics above. The regression model evaluated the effect of genotype on the two latent factors. In the best fitted model, genotype was labeled by the presence of Ile allele as follows: Vmat1WT = 0, Vmat1Thr/Thr = 0, Vmat1 Thr/Ile = 1, and Vmat1 Ile/Ile = 1. The goodness of model fit was evaluated by normed Comparative Fit Index (CFI), Root Mean Square Error of Approximation (RMSEA), and the Standardized Root Mean Square Residual (SRMR), and the significance of paths was tested by t-test. Normalization of behavioral scores and SEM were performed using the bestNormalize and lavaan packages in R, respectively.

In vivo electrophysiological recordings

Subjects

Five male WT C57BL/6J mice (16–20 weeks old, SLC Shizuoka, Japan) with preoperative weights of 20–30 g as well as four male Vmat1Ile/Ile mice (16–20 weeks old) and four male Vmat1Thr/Thr mice (20–24 weeks old) obtained from the 2nd batch were implanted with intracranial electrodes for in vivo electrophysiological recordings. The animals were housed under a 12 h/12 h light/dark schedule with lights on at 7:00 AM prior to surgery.

Surgery

Animals were anesthetized with isoflurane gas (1%–3%) and 1-mm diameter circular craniotomies were made using a high-speed drill at the target coordinates. For local field potential (LFP) recordings, an array of 3 immobile tetrodes was stereotaxically implanted above the dmPFC (2.00 mm anterior and 0.50 mm lateral to bregma) at a depth of 1.40 mm and an array of 4 tetrodes was implanted in the amygdala (0.80 mm posterior and 3.00 mm lateral to bregma) at a depth of 4.40 mm using guide cannulae. The tetrodes were constructed from 17-μm diameter polyimide-coated platinum–iridium (90/10%) wire (California Fine Wire, Grover Beach, CA), and the individual electrode tips were plated with platinum to lower the impedance to 200–250 kΩ. Stainless steel screws were implanted on the skull and attached to the brain surface to serve as ground/reference electrodes. To prevent brain tissue drying at the implant sites, a dummy cannula was inserted through the guide cannula, and both were covered by a cap. The entire recording device was secured to the skull using stainless steel screws and dental cement. After all surgical procedures were completed, anesthesia was discontinued, and the animals were allowed to awaken spontaneously. Following surgery, each animal was housed in a separate transparent Plexiglas cage with free access to water and food for at least 7 days before recordings.

Electrophysiological recording

The mouse was connected to the recording equipment via a digitally programmable amplifier (Cereplex M, Blackrock Microsystems, Salt Lake City, UT) placed close to the animal’s head. For recording electrophysiological signals, the EIB of the microdrive array was connected to a Cereplex M digital headstage, and the digitized signals were transferred to a Cereplex Direct data acquisition system. Electrical signals were sampled at 2 kHz and low-pass filtered at 500 Hz. The animal’s moment-to-moment position was tracked at 15 Hz using a video camera attached to the ceiling.

The elevated plus maze used during electrophysiological recording was made of ABS resin and consisted of a central square (7.6 × 7.6 cm) and four arms (each 28 cm long × 7.6 cm wide), two open arms with no railing and two closed arms enclosed by transverse walls (15 cm in height). The maze was elevated 30 cm from the floor and illuminated with two 32-W fluorescent overhead lights, which produced light intensities of 250 ± 20 lx and 170 ± 20 lx in the open and closed arms, respectively. In a recording session, a mouse was placed in the center of the central square facing an open arm and allowed to explore the maze apparatus for 10 min.

Confirmation of electrode locations by histological analysis

The mice were overdosed with isoflurane, perfused intracardially with 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS, pH 7.4), and decapitated. After dissection, the brain was fixed overnight in 4% PFA/PBS and then cryoprotected by successive overnight incubations in 20% sucrose and 30% sucrose in PBS. Frozen coronal sections (100 μm) were cut using a microtome, mounted, and processed for cresyl violet staining. Briefly, slices were rinsed in water, stained with cresyl violet, and coverslipped with Permount. The positions of all electrodes were confirmed by identifying the corresponding electrode tracks in histological tissue sections.

LFP analysis

LFP data were analyzed using MATLAB. To compute the time-frequency representation of LFP power, LFP signals were convolved using complex Morlet wavelet transformation at frequencies ranging from 1 to 250 Hz and z-scores were computed for each frequency band across an entire period. When data were obtained from multiple electrodes in a mouse, they were averaged to single values within each mouse. To compare LFP power changes between open and closed arms, we first evaluated data normality by the F test and confirmed that some datasets used in this study were not normally distributed. We thus compared two-sample data by Mann-Whitney U test, a non-parametric test. Multiple comparisons were performed by post hoc Bonferroni corrections. The null hypothesis was rejected at the P < 0.05 level.

Extraction of RNA and RNA sequencing

Total RNA was extracted from the prefrontal cortex, basolateral amygdala, and striatum of four mice per genotype (two 4-month-old and two 10-months-old individuals from batches 1 and 2, respectively), yielding 48 samples in total. Cage mates were selected for sequencing when possible (see Figure S3 and Table S4 for detailed sampling information). Briefly, mice were habituated to a novel cage for 30 minutes and then decapitated under anesthesia. The indicated brain regions were then isolated using a mouse brain atlas (Palkovits and Brownstein, 1988) by the same experimenter (TS). The extracted tissues were homogenized in BioMasher Standard (TaKaRa, Shiga, Japan) and total RNA was isolated using the RNeasy® Plus Mini Kit (Qiagen, Hilden, Germany). RNA concentration and purity were assessed using a Nano-Drop® ND-1000 spectrophotometer (Thermo Scientific, Waltham, MA), and total RNA integrity was quantified by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). From the 48 RNA samples, 48 cDNA libraries were prepared and 100 bp paired-end reads were sequenced on a DNBSEQ platform at the Beijing Genomics Institute (BGI, Hong Kong, China). A total of 2,344,530,700 raw sequencing reads were generated (24–25 million read pairs per sample; Table S4) and used for gene expression profiling.

Data processing

The raw reads were quality checked, trimmed, and filtered using fastp 0.20.0 (Chen et al., 2018), yielding 20–21 million filtered reads (Table S4). STAR v2.7.5c (Dobin et al., 2013) was then used for mapping the reads to the reference mouse genome (GRCm38.p6). We assigned reads to the Ensembl100 annotation and generated fragment counts by the featureCounts utility of the Subread package (v1.6.3) in R (Liao et al., 2014). We obtained reads with mapping quality above 20 (essentially uniquely mapped reads) and only pairs that were properly aligned on the same chromosome. The fragment counts were then used to perform differential expression analysis with iDEP v0.91 workflow (Ge et al., 2018), a web application that includes several conventional software applications such as DESeq2 (Love et al., 2014) and edgeR (Robinson et al., 2010), and can perform a comprehensive analysis of RNA-Seq data. After normalizing all counts according to effective library size, we retained genes with > 1 count per million (CPM) in half of the samples (24 samples). A total of 13,608 genes were then assessed in downstream analyses. Principal component analysis (PCA) was first performed to reveal the overall expression pattern of each sample. While most samples demonstrated distinct expression patterns, one supposed amygdala sample (derived from a 4-month-old Vmat1Thr/Thr mouse) exhibited an expression pattern like the other striatal samples (Figure S20). Given the inconsistency in the expression pattern, we excluded this sample and conducted the normalization step again for the remaining 47 samples.

Differential expression and weighted gene correlation network analysis

We conducted pair-wise between-genotype comparisons of the same brain regions to identify differentially expressed genes (DEGs), which were defined by |log2(fold-change)| > 1 and false discovery rate (FDR; corrected by the Benjamini and Hochberg method) < 0.05. The DEGs were further characterized by enrichment analysis for GO terms in the iDEP pipeline. An Illumina BaseSpace application (http://basespace.illumina.com/apps) was used to investigate the correlation in expression levels of DEGs among currently and previously collected datasets. In this analysis, we set FDR < 0.05 as the cutoff to select genes for input, yielding 131 genes recognized in total. To identify specific biological pathways affected by Vmat1 genotypes, we conducted a weighted gene correlation network analysis (WGCNA) (Langfelder and Horvath, 2008) of 15 samples from the amygdala. First, we calculated the transcripts per million (TPM) for each gene, and the expression matrix was then adjusted for the mouse batch as a potential confounding factor. The 1000 genes with the most variable expression levels (largest coefficients of variation) among samples were retained. We chose a soft thresholding power = 21, as it yielded a scale-free topology fitting index (R squared) greater than 0.8. The correlations among gene expression levels across samples were then calculated to detect co-expressing modules. We also calculated the correlations between the expression levels of a given gene in the amygdala and the behavioral composite (locomotor activity and anxiety), based on the data of 8 individuals (two for each genotype). Genes belonging to modules of interest were further included in protein-protein interaction networks constructed using STRING (Szklarczyk et al., 2021), and functional gene clusters were detected by the MCODE module of Cytoscape 3.8.2 (Shannon et al., 2003).

Quantitative RT-PCR

To further validate the results obtained in RNA-Seq and to eliminate the possibility that a large proportion of amygdala samples of each Vmat1 genotype was derived from different nuclei (i.e., basal, lateral, or central amygdala), quantitative reverse transcription-PCR (qRT-PCR) was performed with the same samples used in RNA-Seq (n = 4, 3, 4, and 4 for Vmat1WT, Vmat1Thr/Thr, Vmat1Thr/Ile, Vmat1Ile/Ile genotypes, respectively). cDNA was prepared from the total RNA of each sample using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems), following the manufacturer’s instructions. PCR was performed using the forward and reverse primers specific to the flanking sequences of Cplx1, Negr1, and Nts and Prkcd, which are marker genes of basal (BA), lateral (LA), and central (CeA) nucleus of the amygdala, respectively (Table S6). The expression levels of each transcript of the four genes were quantified using the StepOne Real-Time PCR System (Thermo Fisher) with Power SYBR Green PCR Master Mix (Applied Biosystems). Normalized gene expression was estimated using Gapdh as an internal control, and the expression levels were visualized relative to the mean expression across samples. Differences in expression levels among genotypes were evaluated with the generalized linear model (GLM) assuming a gamma distribution.

Quantification and statistical analysis

Information on sample sizes, statistical tests, and P values are provided in the respective figure and its legend. The difference in mean values among genotypes in behavioral tests were assessed with one-way ANOVA and pair-wise t-tests (alpha = 0.05). Two-way repeated measures ANOVA was used to assess the difference among genotypes in time-series behavioral data. Paired t-tests were used to evaluate the difference in values between conditions within each genotype. Generalized additive model with a quasi-Poisson distribution was used to test the interactive effect of genotype and cage place in Figure 3A. The significance of paths in the SEM model was assessed by t test in Figure 3B. In addition to one-way ANOVA, pair-wise t-tests with FDR correction by the Benjamini-Hochberg was used to assess the variance among genotypes in Figure 3C. Two-way ANOVA on each brain region was first applied to evaluate the variance between arm conditions and genotypes in Figure 5F, followed by Mann-Whitney U test with Bonferroni correction to assess the difference in values between arm conditions within each genotype. Generalized additive model with a gaussian and quasi-binomial distribution was used to assess the difference in values between genotypes in Figure S19. Generalized linear model with a Gamma distribution was used to evaluate difference in expression levels among genotypes in Figure S25A. Dunnett’s test was used to assess the difference in fold changes of expression levels among genotypes in Figure S27. Statistical significance was evaluated with alpha = 0.05 in all the tests. All statistical analyses were performed with R 4.0.2.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Rabbit monoclonal anti-HA-tag	Cell Signaling TECHNOLOGY	C29F4, #3724; RRID: AB_1549585
Chicken polyclonal anti-Tyrosine hydroxylase	abcam	ab76442; RRID: AB_1524535
Sheep polyclonal anti-Tryptophan hydroxylase	Millipore	AB1541; RRID: AB_90754
Mouse monoclonal anti-Synaptophysin	Boehringer Mannheim	902314; RRID: AB_2313789
Donkey anti-Rabbit IgG H&L Alexa Fluor 568	abcam	ab175470; RRID: AB_2783823
Donkey anti-Chicken IgY H&L Alexa Fluor 488	Jackson ImmunoResearch LABOLATORIES	703-545-155; RRID: AB_2340375
Donkey anti-Sheep IgG H&L Alexa Fluor 488	abcam	ab150177; RRID: AB_2801320
Donkey anti-Mouse IgG H&L Alexa Fluor 488 preadsorbed	abcam	ab150109; RRID: AB_2571721
Chemicals, peptides, and recombinant proteins
EnGen Cas9 NLS, S.pyogenes	New England Biolabs	M0646
Alt-R S.p. HiFi Cas9 Nuclease V3	Integrated DNA Technologies	1081058
Critical commercial assays
In-Fusion Snap Assembly Master Mix	Takara Bio	638947
RNeasy® Plus Mini Kit	Qiagen	74904
High Capacity cDNA Reverse Transcription Kit	Applied Biosystems	4368814
Power SYBR Green PCR Master Mix	Applied Biosystems	4367659
Deposited data
RNAseq fastq files deposited to NCBI Sequence Read Archive (SRA)	NCBI	SRA: PRJNA660500
Experimental models: Organisms/strains
C57BL/6J	Charles River Japan	000664
B6C3F1	SLC Japan	B6C3F1/Slc
B6.Cg-Gt(ROSA)26Sortm9(CAG-tdTomato)Hze/J (Ai9)	The Jackson Laboratory	007909
Vmat1 humanized (Vmat1Thr/Thr, Vmat1Thr/Ile, Vmat1Ile/Ile) mice	This paper	N/A
Vmat1-T2A-iCre mice	This paper	N/A
Vmat1-3×HA mice	This paper	N/A
Oligonucleotides
CRISPR RNAs and donor DNAs, see Table S1	This paper	N/A
Primers for PCR-RFLP assay, see Table S1	This paper	N/A
Primers for off-target analysis, see Table S1	This paper	N/A
Primers for knock-in mice genotyping, see Table S1	This paper	N/A
Primers for qRT-PCR, see Table S6	This paper	N/A
Recombinant DNA
T2A self-cleaving peptide sequences	Trichas et al. 2008	https://bmcbiol.biomedcentral.com/articles/10.1186/1741-7007-6-40
Codon improved Cre recombinase sequences	Shimshek et al. 2002	https://onlinelibrary.wiley.com/doi/abs/10.1002/gene.10023
HA-tag sequences, see Table S1	This paper	N/A
Software and algorithms
MATLAB R2018b	MathWorks	https://www.mathworks.com/products/matlab.html
CerePlex Direct Software Suite	Blackrock Microsystems	https://www.blackrockmicro.com/technical-support/software-downloads/
WebPRANK	Löytynoja and Goldman, 2010	https://www.ebi.ac.uk/goldman-srv/webprank/
IQTREE v2.1.1	Minh et al. 2020	http://www.iqtree.org
iTOL v5.6.3	Letunic and Bork, 2019	https://itol.embl.de
SWISS-MODEL	Waterhouse et al. 2018	http://swissmodel.expasy.org
PyMOL 2.4.1	DeLanoScientific	https://pymol.org/2/
Provean v1.1.3	Choi and Chan, 2015	http://provean.jcvi.org/index.php
SIFT	Kumar et al. 2009	http://provean.jcvi.org/index.php
DynaMut	Rodrigues et al. 2018	http://biosig.unimelb.edu.au/dynamut/
CRISPOR	Concordet and Haeussler, 2018	http://crispor.tefor.net
R 4.0.2	R Core Team	https://www.r-project.org
fastp 0.20.0	Chen et al. 2018	https://github.com/OpenGene/fastp
STAR v2.7.5c	Dobin et al. 2013	https://github.com/alexdobin/STAR
featureCounts	Liao et al. 2014	http://subread.sourceforge.net
iDEP v0.91	Ge et al. 2018	http://bioinformatics.sdstate.edu/idep/
BaseSpace	Illumina	http://basespace.illumina.com/apps
WGCNA	Langfelder and Horvath, 2008	https://cran.r-project.org/web/packages/WGCNA/index.html
STRING	Szklarczyk et al. 2021	https://string-db.org