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

Abstract

Animal domestication efforts have led to a shared spectrum of striking behavioral and morphological changes. To recapitulate this process, silver foxes have been selectively bred for tame and aggressive behaviors for more than 50 generations at the Institute for Cytology and Genetics in Novosibirsk, Russia. To understand the genetic basis and molecular mechanisms underlying the phenotypic changes, we profiled gene expression levels and coding SNP allele frequencies in two brain tissue specimens from 12 aggressive foxes and 12 tame foxes. Expression analysis revealed 146 genes in the prefrontal cortex and 33 genes in the basal forebrain that were differentially expressed, with a 5% false discovery rate (FDR). These candidates include genes in key pathways known to be critical to neurologic processing, including the serotonin and glutamate receptor pathways. In addition, 295 of the 31,000 exonic SNPs show significant allele frequency differences between the tame and aggressive populations (1% FDR), including genes with a role in neural crest cell fate determination.

Keywords: domestication; fox; selection; transcriptome.

Fig. 1.

RNA-seq analysis identified differentially expressed genes in brain tissues between the tame and aggressive fox populations. (A) Artificial selection scheme for tameness and aggression in foxes. The conventional population of farm-bred foxes (blue arrow) was a founding population for both tame and aggressive fox populations. The population of conventional farm-bred foxes is still maintained in Novosibirsk. Starting in 1959, the selection experiment for tame foxes has been carried out to recreate the evolution of canine domestication. In 1970, an aggressive population was also selected to compare with the tame population. (B) A volcano plot showing differentially expressed genes detected in 12 tame and 12 aggressive fox prefrontal cortex samples. Plotted on the x-axis is the log2 fold difference between tame and aggressive samples. Plotted on the y-axis is the –log10 (P value) calculated with the R package edgeR. Significant differentially expressed genes (FDR <0.05) are indicated in red, and nonsignificant genes are shown in gray. (C) Bar plot of qRT-PCR validation results in prefrontal cortex and forebrain samples for the top two significant candidate genes, PCDHGA1 and DKKL1.

Fig. 2.

Genes that are differentially expressed between tame and aggressive fox populations in serotonin and glutamate receptor pathways. Diagrams of a serotonergic synapse (A) and a glutamatergic synapse (B) show the presynaptic and postsynaptic terminals (adapted from KEGG pathway database). The RNA-seq expression levels in both tissues are plotted in individual bar plots for significantly differentially expressed genes (q < 0.10 in at least one tissue) between tame and aggressive foxes. Differentially expressed receptors and genes involved in downstream signaling pathways (assigned by KEGG; SI Appendix, Fig. S8) are in blue boxes. (A) In tame individuals, serotonin receptors HTR5A-like is up-regulated in both tissues. HTR3A is up-regulated only in the prefrontal cortex, and HTR7 is down-regulated in the cortex. DUSP1 is in the cAMP/PKA pathway (middle right part of the figure), and AKT1 is a major component of the PI3K/AKT pathway (bottom right of the figure). They are both up-regulated in tame foxes. (B) A subclass of glutamate receptors, NMDA receptor 2D (GRIN2D; glutamate receptor, ionotropic, N-methyl-d-aspartate 2D), and downstream signaling genes ITPR3 and ADCY7 (pathways in red boxes in the middle right and bottom right parts of the figure, respectively) are differentially expressed between tame and aggressive foxes, with up-regulation in the tame animals.

Fig. 3.

GRM3, a metabotropic glutamate receptor gene with significant allele frequency changes in the tame population. (A) Gene dropping simulation scheme to determine the adjusted P value under genetic drift, inbreeding, and founder effect. A null distribution assuming no association between SNP genotypes and behavioral phenotypes was generated by simulating all founder genotypes under a grid of starting founder allele frequencies (0.01∼0.99 in increments of 0.01). Then alleles were dropped down the observed tame and aggressive pedigree structures (SI Appendix, Figs. S2 and S3) based on Mendelian inheritance. This was repeated to produce a null distribution of the magnitude of allele frequency changes. From this, we obtained P values for the observed allele frequency difference between tame and aggressive RNA-seq samples. A total of 295 SNPs were significant across all starting allele frequencies at a 1% level based on 10,000 simulations. (B) A volcano plot showing allele frequency differences between tame and aggressive RNA-seq samples on the x-axis and the –log10 P values on the y-axis. The 295 significant SNPs are labeled in red. (C) GRM3 (metabotropic glutamate receptor 3) has a C → G nonsynonymous SNP change causing a Thr to Ser missense mutation (T52S). In the RNA-seq data, aggressive foxes have 100% C alleles, and tame foxes only have 30% C alleles (P = 4 × 10⁻⁷; adjusted P < 0.01). PBP1_mGluR_groupII, ligand-binding domain of the group II metabotropic glutamate receptor; NCD3G, nine cysteines domain of family 3 GPCR; 7tm_3, 7 transmembrane sweet-taste receptor of 3 GCPR. Annotation from RCSB Protein Data Bank (UniProt ID code Q14832). (D) Crystal structure of the GRM3 extracellular region (PDB ID code 3SM9) viewed by jmol software. T52S (labeled in blue) is near the ligand-binding site, suggesting that it might alter the protein function. (E) Integrative genomics viewer screen shot at the GRM3 SNP position in pooled gDNA-seq samples (SI Appendix, Figs. S10 and S11). In independently selected gDNA resequencing samples, the tame G allele frequency (67%) is confirmed in the tame population and is missing in the aggressive population. (F) The C allele is conserved in dogs, other mammals, and chickens. The tame G allele is the derived allele.