Genomic responses to selection for tame/aggressive behaviors in the silver fox (Vulpes vulpes)

The profound behavior differences appeared rapidly after selection, and changes in brain gene expression levels might have played an important role in the response. To investigate this, Illumina RNA-seq experiments were performed on brain tissues from 12 aggressive and 12 tame individuals (SI Appendix, Figs. S1–S3), including the right prefrontal cortex and right basal forebrain (SI Appendix, Fig. S4). These experiments yielded a total of 1.57 billion RNA-seq reads, with an average of 30 million reads per sample (SI Appendix, Tables S1 and S2). These reads were aligned to both the fox draft genome scaffolds and de novo brain transcriptome assembly (Materials and Methods), producing high-quality read count data on the 48 samples for 12,808 annotated genes in the transcriptome. Among these genes, 146 are differentially expressed in prefrontal cortex of tame and aggressive individuals at a 5% false discovery rate (FDR; q <0.05) (Fig. 1B and SI Appendix, Fig. S5 and Table S3). In addition, there were 33 differentially expressed genes in the basal forebrain (SI Appendix, Fig. S6 and Table S4).

Among these hits, the two most significant genes are DKKL1 and PCDHGA1 (P < 10⁻⁸ in prefrontal cortex and P < 10⁻¹¹ in basal forebrain; Fig. 1B), and their up-regulation in tame fox was confirmed using qRT-PCR in the same RNA-seq samples (SI Appendix, Fig. S7 and Table S5; Materials and Methods). DKKL1, Dickkopf-like protein 1, has signal transducer activity and interacts with the noncanonical Wnt pathway. In the mouse brain, DKKL1 displays region-specific expression, with the highest expression level in the cortical neurons of the adult cortex (7). Little is known about the function of DKKL1 in the brain except that it bears sequence similarity to DKK1, an antagonist of canonical Wnt signaling implicated in a wide spectrum of physiological processes, including neurogenesis, neuronal connectivity, and synapse formation. Overexpression of DKKL1 in the ventral hippocampus, but not in the prefrontal cortex, has been associated with increased susceptibility to social defeat stress in mice (23). PCDHGA1, the Protocadherin Gamma Subfamily A1 gene, encodes a neural cadherin-like cell adhesion protein. Protocadherins are known to play critical roles in the establishment and function of specific cell–cell connections in the brain, such as synapse development (24), as well as dendrite arborization and self-avoidance in the central nervous system (25, 26). Pcdhga1 expression was down-regulated in a learned helpless rat model, suggesting that its expression might affect behavioral phenotypes (27). Twenty-eight of the 146 significant genes in the prefrontal cortex (SI Appendix, Table S6) and 15 of 33 genes in the forebrain (SI Appendix, Table S7) are under known fox behavior QTL peaks (17, 20), and more than one-half of these genes overlap with the two significant QTL peaks on fox chromosome 12. The RNA-seq experiments identified 163 differentially expressed genes that might be responsible for the behavioral phenotype changes after selection.

Previous studies of pathological aggression and anxiety in humans and other animals strongly suggest that genes involved in several neurologic receptor pathways may have altered expression levels in tame foxes. Serotonin (5-HT) is a neurotransmitter known to play a role in feelings of contentment/happiness in humans, and deficiency has been linked with many mood disorders, including anxiety and depression (28). Altered expression levels of serotonin receptors have been documented in patients with schizophrenia and bipolar disorder (29). Serotonin and serotonin metabolite (5-HIAA) levels have been found to be significantly elevated in tame foxes compared with aggressive foxes (7), similar to other mammals and invertebrates (19, 20). In this study, we examined genes in the serotonin receptor pathways based on the KEGG database (30, 31) and found significantly differentially expressed genes, including serotonin receptors 5A, 3A, and 7 and a pair of downstream signaling genes: DUSP1 in the cAMP/PKA pathway and AKT1 in the PI3K/AKT pathway (Fig. 2A and SI Appendix, Figs. S8 and S9). Nearly all the changes are in the direction of increased serotonin signaling in the tame animals.

Besides the critical role of serotonin, dopamine and glutamate also have been linked to aggression (32). In our dataset, no genes in the dopamine receptor pathway were identified as significantly differentially expressed. For the glutamate receptor pathway, the NMDA receptor 2D subunit and downstream signaling genes ITPR3 and ADCY7 were significantly up-regulated in the tame animals (Fig. 2B and SI Appendix, Figs. S8 and S9). NMDA receptors are a subclass of glutamate receptors important for synaptic plasticity, learning, and memory. This pathway also plays a key role in fear conditioning (33). Up-regulation of NMDA signaling might be consistent with increased responsiveness to keepers by the tame foxes. These results suggest that the gene expression response to selection for tameness in silver foxes impacts neurotransmitter receptor pathways and shed light on the biological basis of affiliative and aggressive behaviors by relating to neurologic and pharmacologic correlates with those behaviors.

In addition to the expression response, other genes may manifest changes in coding sequences that could affect protein function. Such genes often show allele frequency changes in their coding SNPs. In the RNA-seq data, we identified 31,025 high-quality exonic SNPs (Materials and Methods) and tested allele frequency differences at these positions. Founder effect, inbreeding, and random genetic drift can result in allele frequency changes, and these factors must be controlled to allow accurate assessment of the role of selection. The tame and aggressive fox populations were selected solely for specific behavioral traits, and full pedigree data for the tame (6,670 individuals) and aggressive (1,863 individuals) populations were maintained during the entire breeding program (SI Appendix, Figs. S2 and S3) (11).

Using this information, we directly simulated the precise effect of genetic drift and inbreeding on allele frequency changes by “gene dropping”, a method that uses the known pedigree structures for an ascertained sample of genotypes drawn from the population (in this case, the 24 RNA-seq individuals) (Fig. 3A and SI Appendix, Fig. S10). At an adjusted P value of 0.01, 295 SNPs in 176 genes had significantly different allele frequencies between the tame and aggressive populations (Fig. 3B and SI Appendix, Table S8), with a mean allele frequency difference of 0.79. Nonsynonymous SNPs were slightly enriched in the significance of allele frequency changes compared with all exonic SNPs (25.9% vs. 23.9%), but the difference did not reach statistical significance (SI Appendix, Fig. S11 C and D). One-third of the 176 significant genes overlapped with significant fox behavior QTL peaks (SI Appendix, Table S9).

Ten whole-genome sequences were obtained for each of the tame, aggressive (SI Appendix, Figs. S12 and S13), and conventional farm-bred fox populations (34) at 25× coverage per population. These genome sequencing samples were from different individuals, allowing independent cross-validation of allele frequency changes. Overall, the SNP allele frequency changes were significantly correlated (Spearman ρ = 0.73; q value < 0.01) between our RNA-seq results and these whole-genome sequences (SI Appendix, Fig. S9B). Among 176 genes containing SNPs with significant allele frequency differences, 73 are located in the regions highlighted in the genome paper (34). Almost all of these genes (n = 71) showed extreme F_ST values in comparisons of tame and aggressive foxes from the genome paper (SI Appendix, Table S12). SorCS1, a transporter important for trafficking AMPA glutamate receptors to the cell surface, is one of the QTL positional candidate genes with decreased heterozygosity and increased divergence between populations that was identified in the analysis of resequenced genomes (34). Six SorCS1 coding SNPs are among the 295 SNPs with significant tame vs. aggressive allele frequency difference, including the third most significant SNP on the list (SI Appendix, Table S8), highlighting the consistency of allele frequency divergence. SorCS1 is also among the 52 genes under the behavior QTL peaks (SI Appendix, Table S9). Despite the consistent allele frequency changes occurring in SorCS1 due to selection, no changes in expression level were detected, perhaps because the parts of the brain that we used in this study are not among the brain regions showing intense SorCS1 expression signals in mouse brain (35).

One of the 176 genes exhibiting a significant SNP frequency change is GRM3, the metabotropic glutamate receptor 3. This glutamate receptor has been associated with schizophrenia, bipolar, mood disorders, and delayed sexual maturity in human studies (36, 37). In our exonic SNP data, GRM3 exhibited a C-to-G change, causing a threonine-to-serine missense mutation (T52S) in the coding region, with a C frequency of 100% in the aggressive foxes but only 30% in the tame foxes (P = 4 × 10⁻⁷; adjusted P < 0.01) (Fig. 3C). The altered amino acid is in the extracellular region near the glutamate-binding site, which might affect binding affinity (Fig. 3D). The allele frequencies were validated in independently selected tame, aggressive, and unselected individuals (Fig. 3E and SI Appendix, Figs. S10 and S11). The tame allele (G) is missing in both aggressive and unselected foxes. Evolutionarily, the ligand binding region is highly conserved, with all genome-sequenced mammals and chicken having the C allele (Fig. 3F). The increased G allele frequency might be a direct response to the artificial selection for tameness in the farm fox experiment.

Our results show that both gene expression and allele frequency responses in the tame foxes occurred in the glutamate receptor signaling pathway; the expression of GRIN2D was elevated in tame individuals in the forebrain, and GRM3 showed significant allele frequency changes. This same pathway also experienced significant changes in both ancient domestication events and recent selection experiments in other mammals. The parallels with the domestic dog are particularly noteworthy, with genes in glutamate receptor signaling (i.e., GRIA1 and GRIN2A) also showing significant changes during the course of domestication (38). Similarly, in domestication of the cat, two glutamate receptor genes, GRIA1 and GRIA2 were also found to be under positive selection (39). A recent selective sweep was also found in GRIK2 in domestic rabbits (6). This convergence of selection signals on glutamate receptor signaling strongly motivates additional experimental confirmation of a functional role for glutamate signaling in behavioral differences of domesticated mammals.

Similarly, genes in the protocadherin family also display both expression and allele frequency changes during selection for tameness in foxes. Three protocadherins, PCDH9, PCDH17, and PCDH20 all have multiple SNPs with significant allele frequency changes (adjusted P < 0.05) (SI Appendix, Table S10). PCDHGA1, a protocadherin gamma gene, is the second most significant differentially expressed gene between tame and aggressive fox brains (Fig. 1B). Remarkably, another member of the same protocadherin gamma subfamily A, Pcdhga11, was among the 20 genes with the highest correlations between gene expression value and tameness score in a F2 rat population (40). Comparative genomic analysis between domestic and wild cats also identified protocadherin A1 and B4 (PCDHA1 and PCDHB4) under the selection peaks (39), suggesting a shared role of protocadherins in tame phenotypes across multiple mammalian species.

Among the candidate domestication genes identified in studies of dogs and other mammals (41, 42), ITGA6 was also found to have significant allele frequency changes in three informative SNPs in our fox study, which also implicates a role for cell surface adhesion and signaling molecules in behavior. A QTL and transcriptome study using an F2 population of two outbred rat lines selected for tameness and aggression identified four candidate genes for the behavioral difference (40). Five criteria were applied to determine these genes’ involvement in tameness. The top two candidate genes, Gltscr2 and Lgi4, were significant for all five criteria, and both have informative SNPs in our fox data. Two synonymous coding SNPs in Lgi4 both showed significant allele frequency differences at an adjusted P < 0.05 (SI Appendix, Table S11). Two nonsynonymous and three synonymous SNPs were found in Gltscr2, and they were marginally significant, with an allele frequency difference of 0.375 (adjusted P = 0.10) (SI Appendix, Table S11). In summary, selection for tame/aggressive phenotypes in different mammals can lead to expression and genetic changes in genes in the same pathways (43).

Charles Darwin, along with many others, observed that selection for domestication in mammals often leads to a collection of phenotypes including shortened snout, curly tail, white spotting of fur on the chest, and floppy ears, often referred to as the “domestication syndrome.” These features all seem to occur in tissues that are derived from neural crest cells, suggesting that the process of selection for domestication impacts neural crest cell function (44). Intriguingly, several of the genes that manifested significant allele frequency changes in our tame foxes may play a role in neural crest cell fate (45). Wnt signaling plays a key role in initial neutral crest cell differentiation, and both Wnt3 and Wnt4 in the fox have more than one SNP with significant allele frequency changes. Protocadherins are also important in neural crest cell function. Direct assessment of whether these genes play roles in neural crest cell function in the fox presents an interesting experimental challenge.