Fox samples and history of the fox experimental populations

Samples were collected from adult foxes maintained at the experimental farm of the Institute of Cytology and Genetics (ICG) (Novosibirsk, Russia).

The samples from three populations maintained at the ICG farm were used in this study:

1. The conventional farm-bred population is a standard farm bred population that is outbred and has not been deliberately selected for behaviour. The conventional farm-bred population originated from foxes from eastern Canada16 where fox farm breeding began in the second part of the nineteenth century. 2. The tame population was developed through selection of conventional farm-bred foxes for a tame response to humans beginning in 1959 at the ICG. The population began with 198 individuals that were selected from several fox farms across the former Soviet Union due to their less aggressive and fearful behaviour towards humans. A description of the selective breeding programme was published previously20,22,23,88. Pedigree records were carefully maintained and a significant effort was made to avoid inbreeding throughout the breeding programme. A representative video of tame fox behaviour is available online: https://www.youtube.com/watch?v = vrqOSgEh0fQ 3. The aggressive population was developed by selecting conventional farm-bred foxes for an aggressive response towards humans, beginning in the late 1960s at the ICG. The population started with approximately 150 initial founders, but an additional 70 conventional farm-bred foxes were introduced into the aggressive population in 1990s. This introduction aimed to increase the population size, which had been reduced shortly after the dissolution of the Soviet Union (1993). A description of the selective breeding programme was published previously20,22,23,88. Pedigree records were carefully maintained and a strong effort was made to avoid inbreeding during the entire breeding programme. A representative video of behavior of aggressive foxes is available online: https://www.youtube.com/watch?v=GeAWbLLNesY

Sample used for whole-genome sequencing

A blood sample from an F 1 male produced by cross-breeding a female from the aggressive strain and a male from the tame strain was used for whole-genome sequencing. DNA from blood was extracted using the phenol–chloroform method89.

Samples used for re-sequencing

Blood samples from 30 individuals, corresponding to 10 from each of the tame, aggressive and conventional farm-bred populations, were collected for re-sequencing. Samples were chosen so as not to share any parents or grandparents, and each population sample included an equal number of males and females (Supplementary Table 4). DNA was extracted using Qiagen Maxi Blood Kits, as per the manufacturer’s instructions.

Samples used for RNA-seq

Brain samples were collected from 24 male foxes (12 from the tame and 12 from the aggressive populations) into RNAlater and then stored at −80 °C. RNA was extracted from three brain regions: the right basal forebrain, the right prefrontal cortex, and the right part of the hypothalamus. Sequencing was performed on an Illumina HiSeq2000. The basal forebrain and prefrontal cortex samples were sequenced using single-end 50 bp reads, and the hypothalamus samples were sequenced using single-end 100 bp reads. In total, 37.2, 41.3 and 72.6 Gb of data were produced for samples from the basal forebrain, right prefrontal cortex and hypothalamus, respectively. The RNA-seq reads were quality filtered and used for annotation of the fox assembly.

RNA-seq quality filtering included several steps. Data quality, GC content and distribution of sequence length were initially assessed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and then reads were processed with flexbar90 in two passes: the first to trim adapters, remove low-quality reads and remove reads less than 35 bp in length, and the second to remove polyA tails. Third, reads that mapped to fox mitochondrial DNA sequences from NCBI (accession numbers JN711443.1, GQ374180.1, NC_008434.1 and AM181037.1) using Bowtie291,92 were discarded, and finally, any remaining reads that mapped to ribosomal DNA sequences were discarded.

Samples used for genotyping

Samples from 64 tame, 70 aggressive, 109 F 1 and 537 F 2 foxes were used for genotyping. Fox F 2 pedigrees were produced by cross-breeding tame and aggressive foxes to produce F 1 and then breeding F 1 foxes to each other to produce F 2 pedigrees. The same set of F 2 pedigrees was previously used for QTL mapping46.

Sequencing and assembly of the fox genome

Fox paired-end and mate-pair DNA libraries with nine different insert size lengths (from 170 bp to 20 kb) were constructed (Supplementary Table 1). The libraries were sequenced on an Illumina HiSeq2000, with the short insert size libraries yielding read-lengths of 100 and 150 bp and the long insert size, mate-pair libraries yielding 49 bp ends (Supplementary Table 1). In total, 366 Gb of raw reads were produced. A series of strict filtering steps was performed to remove artificial duplications, adapter contamination and low-quality reads93. The program SOAPdenovo v.2.04.434 was used for de novo assembly (Supplementary Table 1). Briefly, reads from the short-insert libraries (<2,000 bp) were first assembled into contigs on the basis of k-mer overlap information. Then, reads from the long-insert libraries (≥2,000 bp) were aligned onto the contigs to construct scaffolds. Finally, we used the paired-end information to retrieve read-pairs and then performed a local assembly of the collected reads to fill gaps between the scaffolds. The program SSPACE v.2.094 was used to extend the pre-assembled scaffolds with reads from all long-insert (2–20 kb) libraries (9 libraries, in total). SSPACE v.2.0 was run with the following parameters: -x 0 -k 5 -n 20. Genome assembly quality was evaluated using GC content and the sequencing depth distribution by mapping all the reads back to reference genome using SOAP295.

The fox genome was assembled into 676,878 scaffolds with a total length of 2,495,544,672 bp, contig N50 of 20,012 bp and scaffold N50 of 11,799,617 bp (Supplementary Table 1). The raw reads and the longest 82,429 scaffolds, which are all scaffolds at least 200 bp in size, were deposited in NCBI (BioProject PRJNA378561).

Annotation of the fox genome

Fox RNA-seq data, de novo gene prediction and homology with canine and human proteins were used to annotate the protein-coding genes in the fox assembly (Supplementary Table 2).

Homologue-based prediction

Protein sequences available for the dog and human from Ensembl release-70 were mapped to the fox genome assembly using TBLASTN (BLASTall 2.2.23) with an e-value cutoff 1 × 10–5. The aligned sequences were then analysed with GeneWise (v.2.2.0)96 to search for accurate spliced alignments.

De novo prediction

Repetitive sequences were masked in the fox genome assembly using RepeatMasker (v.3.3.0) (http://www.repeatmasker.org/). De novo gene prediction was then performed with AUGUSTUS (v.2.5.5)97. The parameters were optimized using the gene models with high GeneWise scores from the homologue-based prediction.

RNA-Seq prediction

Filtered RNA-Seq reads from three tissues were aligned against the fox genome assembly using TopHat98. The candidate exon regions identified by TopHat were then used by Cufflinks99 to construct transcripts. Finally, the Cufflinks assemblies for the three tissues were merged using the Cuffmerge option in Cufflinks.

The three gene sets obtained by each of the three approaches (homologue-based prediction, de novo prediction and RNA-Seq prediction) were integrated based on gene structures. Finally, all gene evidence was merged to form a comprehensive and non-redundant gene set. In total, 21,418 protein-coding genes were identified in the fox genome.

Gene annotation

In order to assign gene symbols to the fox genes with high confidence, a reciprocal blast method was applied. Fox protein sequences and the dog protein sequences that are located on the dog chromosomes, not chromosome fragments (as downloaded from Ensembl release-73), were analysed with BLASTP in both directions. The BLASTP-aligned results were filtered using an e-value cutoff of 1 × 10–5, and reciprocal best hit (RBH) pairs were determined using the following condition: for two genes (for example A and B) from the fox gene set and the dog gene set, respectively, they would be accepted as an RBH pair if and only if they were reciprocally each other’s top-BLASTP-score hits, meaning there was no gene in the fox gene set with a higher score than A to B, and there was no gene in the dog gene set with a higher score than B to A. This analysis of the fox predicted genes against the dog Ensembl database identified 16,620 dog Ensembl IDs, 14,419 with gene symbols available. Fox protein sequences and human protein sequences on human chromosomes, not chromosome fragments (downloaded from Ensembl release-73), were analysed with BLASTP using the same protocol. This analysis idendified 15,826 human Ensembl ID’s, all having an associated gene symbol. These 15,826 high confidence gene symbols were assigned to the associated fox genes and were used in downstream analysis.

The 21,418 predicted protein coding genes were compared against several databases to produce a preliminary annotation. Genes were aligned using BLASTP to the SwissProt and TrEMBL databases100and assgined to the best match of their alignments. Motifs and domains of genes were determined by InterProScan101 against protein databases including ProDom, PRINTS, Pfam, SMART, PANTHER and PROSITE. Furthermore, all genes were aligned against the KEGG102 proteins, and the pathway in which the gene might be involved was derived from the matched genes in KEGG. The fox genome annotaion statistics are presented in Supplementary Table 2.

Alignment of the fox scaffolds against the dog genome

The top 500 longest scaffolds (size range: 47,686 bp to 55,683,013 bp), which contain 94% of the fox genome by length, were aligned against the CanFam3.1 assembly (autosomes, mitochondrial DNA and X-chromosome) and the dog Y-chromosome assembly103 using LAST104. Because each scaffold mapped to multiple locations in the dog genome, we sought to identify the dog chromosome(s) to which it was most likely syntenic. For each scaffold, the maximum LAST score corresponding to each dog chromosome was identified. These scores were Z-transformed using the formula \(\left( {x_{i} - \bar{x}} \right)/{\sigma}\), and the dog chromosome(s) with Z-scores significant at P < 0.05 to a particular scaffold were considered syntenic to that scaffold (Supplementary Table 3).

To confirm the accuracy of the assignment of each fox scaffold to one or more syntenic positions in the dog, the LAST mapping results were then scanned with a Python script to determine the best hit at each nucleotide along each scaffold. The LAST mapping data were imported into a MySQL database to identify which dog chromosome corresponded to the highest-scoring mapped segment overlapping each nucleotide along the fox scaffold. Regions mapping to an individual chromosome were plotted as lines using MatPlotLib, with the position on the scaffold as the x-axis and the position on the dog chromosome as the y-axis. Dog chromosomes to which the scaffold mapped robustly are identified in the legend on the plot. Robust mapping was defined as cases where the best mapping score for the scaffold against that chromosome was at least one standard deviation above the average highest score across all chromosomes. This strategy allowed for visualization of the relationship between each scaffold and the dog genome based on this high score alone, and the fact that it showed an overwhelming consensus with the Z-score data supported the assignment of dog syntenic fragments using the first approach (Supplementary Fig. 2).

Re-sequencing of fox samples from three populations

DNA samples from 30 foxes (10 foxes from tame, 10 foxes from aggressive, and 10 foxes from conventional farm-bred population) were sequenced using individual libraries. The libraries were constructed using the Nextera DNA Sample Preparation kit V2 (Illumina) and included individual barcodes. The libraries were quantified by qPCR and pooled by combining five individuals from a single population (six pools in total). Each pool was sequenced on one lane of Illumina HiSeq2000 using a TruSeq SBS sequencing kit version 3 (Illumina) for 100 cycles from each end of the fragments. Reads were analysed with Casava1.8.2 (Illumina). The genome of each individual was sequenced at approximately 2.5×. Nine samples, four from the tame and five from the conventional populations, that received lower sequence coverage were re-sequenced on a part of a lane to balance the total amount of sequencing data obtained for all individuals. In total, 75.9 Gb, 81.8 Gb and 67.5 Gb of sequencing were obtained for the tame, aggressive and conventional samples, respectively (Supplementary Table 4). The sequencing data was deposited to NCBI (BioProject PRJNA376561).

Read alignment and SNP calling

The reads obtained for each sample were mapped, for each individual, with Bowtie291,92 to the 676,878 scaffolds of the fox assembly. Reads that mapped to more than one location or that mapped with a quality lower than a Phred score of 20 were removed using SAMtools105. The MarkDuplicates tool of Picard (http://picard.sourceforge.net) was utilized to remove duplicated reads. The ten samples from each population were then combined into population pools, and GATK106,107,108 was used to re-align indels. Fox SNPs were identified using two SNP-calling programs, UnifiedGenotyper and ANGSD (Supplementary Table 5):

1. SNPs were called by GATK UnifiedGenotyper with the pooled data from each of the three populations (three pools of 10 individuals each). SNPs with more than 2 alleles and with extremely high or low read coverage (more than 3 × the average depth across all samples, or less than 1/3 the average depth across all samples) were removed using vcftools (–min-alleles 2 –max-alleles 2 –min-meanDP 8.60543 –max-meanDP 77.44887). 2. SNPs were also called using ANGSD109 for individual samples from each of the three populations (30 individual samples). SNPs were called using the parameters: -doMajorMinor 1 -GL 2 -doMaf 2 -doGeno 7 -realSFS 1 -doSNP 1 -doPost 1 -doCounts 1 -dumpCounts 4 -doHWE 1 and then filtered with parameter –lrt 50.

SNPs called by both programs were identified using scaffold locations, and a total of 8,458,133 SNPs identified by both programs were retained for further analysis (Supplementary Table 5).

Principal component analysis

PCA was performed using the genotypes of all individuals across the set of 8,458,133 SNPs without providing any information about the populations of origin of the re-sequenced samples (tame, aggressive or conventional). A covariance matrix for the SNP data was calculated using the EIGENSOFT software110. The eigenvectors from the covariance matrix were generated with the R function ‘eigen,’ and significance was determined with a Tracy–Widom test to evaluate the statistical significance of each principal component (P < 0.01 for both the first and the second principal components). The results of PC analysis were visualized using R.

Construction of the individual tree

A tree of relationships among the sequenced individuals from the tame, aggressive and conventional farm-bred populations was constructed using the neighbour-joining method36. Individual genotypes for the 8,458,133 SNPs were used. The distances (D ij ) between each pair of individuals (i and j), were calculated using the formula:

$$D_{ij} = \mathop {\sum }\limits_{m = 1}^{M} d_{ij}/L$$

Where M is the number of segregating sites in i and j; L is the length of regions and d ij is the distance between individuals i and j at site m. We set d ij equal to 0 when individuals i and j were both homozygous for the same allele (AA/AA); 0.5, when at least one of the genotypes of an individual i or j was heterozygous (Aa/AA, AA/Aa or Aa/Aa); and 1, when individuals i and j were both homozygous but for different alleles (AA/aa or aa/AA). We used the distance matrix of d ij to construct a phylogenetic tree using the neighbour-joining method and the program fneighbor111.

STRUCTURE analysis

Clustering analysis was performed using the Bayesian inference program STRUCTURE 2.3.437,38,39,40. Individual genotypes for 680,000 SNPs randomly chosen from the 8,458,133-SNP set were used. Four independent runs were performed at each level of k from 1 to 5 with a burn-in of 100,000 and 100,000 Markov-chain Monte Carlo replicates using the admixture model without prior information about populations. The values for estimated log probability of data, L(K), were used to calculate delta k for the levels of k from 2 to 4 in order to find the optimal number of subpopulations following a typical procedure112 (Supplementary Table 21). The value for both delta k and the mean of the estimated log probability of the data were highest at k = 3 (Supplementary Table 21).

Analysis of allele frequency differences

Pooled heterozygosity

Pooled heterozygosity (H p ) is a measure of heterozygosity in a set of samples across a region containing multiple SNPs113. Re-sequenced samples from each population (10 samples per population) were combined, and H p was estimated for each of the three populations separately. Because each individual was sequenced with low coverage (~2.5×) we used allelic read depth in pooled data (~25 coverage) for H p estimation in each population. The depth of each individual allele was counted using the SNP data from the GATK/UnifiedGenotyper run and used to determine the major and minor allele frequencies for each SNP in each population. H p was calculated using a sliding window approach. The selection of window size has considered several factors, including the estimated linkage disequilibrium (LD) length in tame and aggressive populations55, simulations of the allelic fixation rate (Supplementary Note 1; Supplementary Fig. 5), and the results of a pilot analysis with smaller window sizes. The 500 kb windows were moved along the fox scaffolds in 250 kb steps. Only scaffolds of 500,000 bases and longer were included in this analysis, corresponding to the largest 309 scaffolds. Within the scaffolds, only windows containing 20 or more SNPs were considered. The average number of SNPs per window was 1,784 (median: 1,739; standard deviation: 1,084; max: 6,730). In total there were 9,151 windows in the analysis. The average read depth per window is presented in Supplementary Table 7. H p was calculated separately for each population using the formula: H p = 2Σn MAJ Σn MIN /(Σn MAJ + Σn MIN )2, where n MAJ and n MIN are the number of reads for major and minor alleles for each SNP, respectively, Σn MAJ is the sum of the reads of the major alleles for all SNPs in that window, and Σn MIN is the same for the minor alleles113. Calculations were performed using in-house scripts written in R. Because the window H p values were not normally distributed (Supplementary Fig. 11), the significance threshold was established in each population by 10,000 permutations following a previous study114. The allele depth data were permutated using the complete set of 8,458,133 SNPs. SNP positions were held constant, and H p was calculated for all windows with over 20 SNPs in every permutation run. 10,000 permutations were conducted in R, and the minimum H p values and values at multiple percentile levels were recorded from each permutation.

For a threshold P-value of <0.0001, the 0.0001 percentile of the minimum values from the 10,000 permutations was calculated in R for each population. All windows in a population with H p values at that calculated value or lower were considered to be significant at P < 0.0001 (Supplementary Table 6). The P-value threshold of 0.0001 (1/10,000) was chosen because there were 9,151 (just under 10,000) windows analysed. This criterion represents a stringent threshold with an expected false positive rate of less than one window per population.

For the window corresponding to the SorCS1 gene, region 94, we estimated the probability of observing the H p values in the tame and aggressive populations compared to a null distribution estimated using 10,000 permutations. We compared the tame and aggressive H p values in the region to the minimum H p value for various percentiles recorded while running the permutations, that is, if the lowest H p value at percentile 0.01 for all of the 10,000 permutations for that population was higher than the observed H p value, the P-value for the observed value is <0.01. The lowest possible percentile for which this is true was reported.

Combined H p windows

The significant H p windows that were identified on the same scaffold and in the same population when the gap between them was not larger than 1 Mb were merged into combined H p windows (Supplementary Table 7). Our reasons for combining these windows were twofold: (1) uneven distribution of reads among windows could impact our analysis; (2) evaluation of the H p values in gap windows (windows located in the 1 Mb interval between H p windows significant within a single population) showed low heterozygosity although these windows did not meet the population’s significance cut-off.

Fixation index

The fixation index (F ST ) was calculated in R using the estimator formula reported previously115, following an earlier publication116, which allows for the use of pooled data in windows. F ST was calculated for the same 9,151 windows that were used in the pooled heterozygosity analysis. For each SNP the following estimators were calculated:

\(\hat N^{\left[ k \right]} = \left( {\frac{{a_{1}}}{{n_{1}}} - \frac{{a_{2}}}{{n_{2}}}} \right)^{2} - \frac{{h_{1}}}{{n_{1}}} - \frac{{h_{2}}}{{n_{2}}}\)

$$\hat D^{\left[ k \right]} = \hat N^{\left[ k \right]} + h_1 + h_2$$

$$h_i = \frac{{a_i\left( {n_i - a_i} \right)}}{{n_i\left( {n_i - 1} \right)}}$$

Where k is the individual SNP; a 1 is the number of reads for allele 1 in population 1; n 1 is the depth of reads for that SNP; a 2 is the number of reads for allele 1 in population 2; and n 2 is the depth of reads for that SNP.

For each window F ST was estimated using the formula:

$$\hat F = \frac{{\mathop {\sum }

olimits_{k = 1}^K \hat N^{\left[ k \right]}}}{{\mathop {\sum }

olimits_{k = 1}^K \hat D^{\left[ k \right]}}}$$

Combined F ST windows

The significant F ST windows that were identified on the same scaffold and in the same type of analysis when the gap between them was not larger than 1 Mb were merged into combined F ST windows (Supplementary Table 7).

Identification of 103 regions of interest

The positions of all significant windows identified in the fox genome were analysed and used to establish regions where either a single significant window was identified, or any combination of classes of significant windows (H p or F ST in any population(s)) were located on a single scaffold within 1 Mb of each other (Supplementary Table 7).

Simulations

Simulations were conducted in forqs117 (see Supplementary Note 1 for more details). Population parameters were selected for the simulation based on pedigree information and breeding records from 1959 (when the population was founded) through 2010, as the DNA samples used in the current study were collected no later than 2010. A base simulation with fifty generations of breeding and 240 animals was conducted and three parameters were varied:

1. To evaluate the effect of population size, the population was simulated with population sizes of 120, 480 and 960 individuals. Each of these scenarios assumed that every founder had two unique haplotypes and that the population was bred for 50 generations. 2. To evaluate the effect of the relatedness of the founding animals, two alternate levels of relatedness were simulated. The populations were set to have either 50 or 100 founding haplotypes distributed evenly in the first generation, in contrast to the 480 in the base simulation. In these scenarios, populations of 240 individuals were bred for 50 generations. 3. To evaluate the effect of the number of generations, breeding of the base population (240 unrelated individuals) was simulated over 100, 250 and 500 generations.

The simulations were run using fox chromosome 1 (VVU1) as a proxy for the fox genome. The chromosomal length (220 Mb) and recombination map (120 cM) were approximated using a meiotic linkage map of VVU1 aligned against the dog genome35 (Supplementary Note 1, Supplementary Fig. 3).

Haplotype frequencies were calculated at 100,000 bp intervals in each simulation scenario. The distribution of the haplotype frequencies (Supplementary Fig. 4) included all non-zero haplotype frequencies across all 100 replications of each scenario. The length of haplotypes that were identical-by-descent with founder haplotypes was calculated for every haplotype in every individual in the final generation. The haplotype lengths were recorded for all 100 replicates of each simulation scenario. The proportion of the genome represented by haplotypes of a given size or shorter was calculated and is shown in Supplementary Fig. 5. The distribution of the average haplotype lengths along chromosome 1 was calculated by dividing the chromosome into one hundred 2.2 Mb windows and averaging the lengths of all haplotypes that have a midpoint falling in the window (Supplementary Fig. 6).

Mapping the fox windows against the dog genome

The 9,151 windows used in the H p and F ST analyses were mapped against the dog genome (CanFam3.1) using LASTZ (v.1.03.66)118 to identify the window order on the fox chromosomes. The ‘multiple’ option of LASTZ was used to map to the entire dog genome in one run, and the alignments were then chained using the ‘–chain’ option. All other parameters were set to default. LASTZ computed alignments separately for the forward and reverse sequence of each window and produced a separate list of alignments for each strand. To identify the best match and the secondary best match for each window, the LASTZ alignments were then filtered using the following protocol:

1. The mapped window segments were sorted by their starting nucleotide positions in the window. The alignments of the first two mapped segments in each window were compared, and if they overlapped by more than 50% of the length of either (after chaining by LASTZ, so this only happened when the same region mapped in different directions), the segment with the lower mapping score was removed, and the one with the higher mapping score was compared again to the next mapped segment in the window for overlap. All mapped window segments that did not overlap with other mapped segments were also retained. 2. Segments that mapped sequentially and in the same direction to the same dog chromosome were combined into a single segment if the ratio of the length of the combined dog segment to the length of the combined fox segment was between 0.8 and 1.2. This step allowed the identification of extended regions where fox segments were mapped to the same dog chromosome in the expected order and without large gaps. When segments were combined, the mapping score of the new, longer segment was calculated as the sum of the mapping scores of the two combined segments. 3. Short mapping segments (<1,000 bp) remaining after the joining of sequential segments were removed. 4. The second filtering step (combining segments mapped to the same dog chromosome, in the same orientation, and of similar length between dog and fox) was run again to combine any segments that were previously separated by a short segment. 5. Medium size mapping segments (<10,000 bp) were removed. 6. The second filtering step was run again to combine any segments that were previously separated by a medium segment.

When there was one filtered result for a window, this result was considered to be the main hit. When there were two hits for a window, the hit with the higher mapping score was reported as the main hit and the lower score was reported as the secondary hit. When there were three or more remaining hits, the window was examined manually and if two or more non-adjacent mapping segments were on the same dog chromosome, in the same direction, and were located close to each other, they were combined to a single extended segment. The top score is used as the primary mapping location and the second highest is reported as the secondary hit. All subsequent matches are not reported.

Out of 9,151 windows analysed, 8,715 (95.3%) mapped to one location in the dog genome, 402 to two locations, 18 to more than two locations and 6 did not receive a location after filtering. The order of windows in the fox genome (Fig. 2) was established using the alignment of the fox scaffolds against the dog genome and the known synteny between dog and fox chromosomes29,30,31,35.

Gene enrichment analysis

The human gene symbols assigned by reciprocal blast in the course of the gene annotation of the fox genome were used in this analysis. Fox orthologues of human genes located inside of, or overlapping with, windows used in the pooled heterozygosity (H p ) and F ST analyses are listed in Supplementary Table 7. To determine the genes overlapping with each window, the intersect tool of bedtools was used with the options –wa and –wb with the windows as the ‘a’ file and the genes as the ‘b’ file.

GO term over-representation analysis

GO term over-representation analysis was performed for the significant windows identified in the H p and F ST analyses using the PANTHER (protein analysis through evolutionary relationships) classification system (PANTHER, v.13.0)41. The six data sets (genes identified in significant H p T, H p A, H p C, F ST TA, F ST TC and F ST AC windows) were analysed. The following over-representation tests were performed: “PANTHER GO-Slim Biological Process,” “PANTHER GO-Slim Molecular Function,” “PANTHER Protein Class,” “GO biological process complete,” “GO molecular function complete” and “GO cellular component complete”. Annotations from the human (all genes in the database) were used as a reference list. Only results of the over-representation test with P < 0.05 after Bonferroni correction were reported (Supplementary Table 9).

Brain-expressed genes

The genes found in the windows were checked for enrichment of genes expressed in the brain. Version 17 of Human Protein Atlas119(http://www.proteinatlas.org/) was used and downloaded from http://v17.proteinatlas.org/download/normal_tissue.tsv.zip. Brain tissues were considered to be caudate, cerebellum, cerebral cortex, hippocampus, hypothalamus and pituitary gland. All genes that have any expression level in any brain tissue except ‘none detected’ were included in the list of brain-expressed genes. Of the 12,976 genes in the version of the protein atlas with relevant data, there were 10,424 genes that showed expression in the brain. Among 15,694 annotated genes in 9,151 fox windows (15,826 high-confidence annotated genes total, but not all are in the windows used in the analysis), 10,991 have data in the Human Protein Atlas and 9,058 are brain-expressed (82.4%). There are 971 annotated genes in our significant windows, among which 698 have data in the Human Protein Atlas and 571 show brain expression (81.8%) (Supplementary Table 10). A hypergeometric test was conducted at https://www.geneprof.org/GeneProf/tools/hypergeometric.jsp and did not find enrichment for brain-expressed genes in significant windows (P = 0.69).

Genes from significant windows were also compared to genes involved in glutamatergic, serotonergic, dopaminergic, GABAergic and cholinergic synapses as listed in the KEGG database (KEGG last updated: 7 December 2017). The enrichment for synapse-related genes from KEGG database (Supplementary Table 11) was tested using a hypergeometric test (https://www.geneprof.org/GeneProf/tools/hypergeometric.jsp) and adjusted for multiple testing with Benjamini–Hochberg correction. No significant enrichment for genes in glutamatergic (adjusted P = 0.148), serotonergic (0.241), dopaminergic (0.381), GABAergic (0.148) and cholinergic (0.148) synapses was observed.

Comparison of fox significant windows with regions associated with domestication and positive selection in dogs

The positions of the 103 fox regions from Supplementary Table 7 were compared with the dog regions associated with domestication and positive selection from four publications42,43,44,45. In three of these studies the dog regions were reported according their location in CanFam242,44,45, the positions of these regions were identified in CanFam3.1 using the liftOver tool from the UCSC browser. The syntenic regions were then identified using an alignment between the fox and dog genomes. Fox windows located within 2 Mb of the fox syntenic positions of the dog regions were considered to be regions that overlap between fox and dog. To test whether this overlap occurred at a rate higher than expected by chance, the extent to which these regions would be expected to overlap was computed by permutation. We combined the four sets of reported dog regions42,43,44,45 into one set of regions for the permutation test. Our 103 fox regions were randomly permuted across the all 9,151 fox windows 10,000 times and the positions of the dog regions were held constant. The number of permuted fox regions that overlapped or were within 2 Mb of the dog regions was recorded for each permutation. The P-value for the actual number of overlap/close regions is the percentage of the 10,000 replications where the number of permuted regions marked as overlapping/close to the dog regions was at or higher than the actual number of overlapping/close regions.

Comparison of 103 fox regions from Supplementary Table 7 with fox behavioural QTL

The positions of fox regions from Supplementary Table 7 were compared with positions of nine fox behavioural QTL identified in previous studies26,46. Only QTL for behavioural phenotypes defined using PCA were included in this analysis. A QTL interval was defined as the genomic region extending 5–15 cM in both directions from the QTL peak, which is the cM position of the QTL with the most significant statistical support. The interval boundary on either side of the QTL peak was defined by the position of the mapped microsatellite marker46 located within the 5–15 cM interval from the QTL peak that was farthest from the QTL peak. For example, if there were three markers on the fox meiotic linkage map46 that fell on same side of the QTL peak at distances 7, 14 and 17 cM, respectively, the boundary of the QTL interval on this side would be placed at the position of the marker located 14 cM from the QTL peak. All microsatellite markers used for QTL mapping were dog-derived markers with known positions in the dog genome. Because the current QTL intervals are large and often correspond to several fox scaffolds, we used the locations of the microsatellite markers in the dog genome46 to define the length and positions of the dog genomic regions syntenic to the fox QTL intervals. These regions were then compared to the dog genomic coordinates of the 103 fox regions from Supplementary Table 7. This analysis identified 30 fox regions (positive regions) that overlap with five out of the nine fox behavioural QTL (Supplementary Table 14).

To test whether the observed overlap between the fox regions and fox QTL intervals is statistically significant, we compared the proportion of the dog genome represented in QTL intervals (that is, the length of all nine QTL intervals relative to the total length of dog autosomes and the X chromosome in CanFam3.1) to the proportion of the windows in the 103 regions from Supplementary Table 7 that overlap with the QTL intervals (that is, the number of windows that are located in the 30 positive regions and that overlap with QTL intervals relative to the total number of windows in 103 regions). The null hypothesis was that the proportion of windows that overlap with the QTL intervals would be similar to the proportion of the dog genome that is represented in the QTL intervals. Based on dog–fox synteny46, we estimated that the length of all nine QTL intervals corresponds to 474,130,369 bases in the dog genome; therefore, 20% of the dog genome is represented in QTL intervals (corresponding to 474,130,369 bases in the QTL intervals out of 2,327,633,984 bases in dog autosomes in CanFam3.1). Out of the 103 regions in Supplementary Table 7, 29 regions completely overlapped QTL intervals (that is, all windows in these regions overlap with QTL intervals) and one region (region 46) partly overlapped a QTL interval (61 out of 77 windows in that region overlapped the QTL interval). In total, the proportion of windows that overlapped QTL intervals was 40% (corresponding to 228 windows across the 30 positive regions overlapping the QTL intervals out of a total of 555 windows in the 103 regions). We performed a chi-square test (http://vassarstats.net/tab2x2.html) and found that that the proportion of the windows that overlap with the QTL intervals was significantly higher than would be expected by chance (χ2 = 82.84, d.f. = 1, P-value < 0.0001).

Functional analysis of intergenic SNPs in significant windows

We used the well-annotated dog genome for functional analysis of intergenic SNPs. As with variant calling in the fox de novo assembly, the reads obtained for the tame, aggressive and conventional populations were aligned to the dog genome (CanFam3.1) using Bowtie291,92, and SNPs were called using the UnifiedGenotyper tool from GATK106,107,108. Sequence variants that showed differences only between the dog and the fox (that is, positions where all foxes were identical and different from dog) were removed. The remaining SNPs were polymorphic in foxes and were filtered using VCFtools120 to include only those that had two alleles, a mean depth from 30–180 reads, and a quality of 100 or greater. This filtering step used the parameters: “–min-meanDP 10 –max-meanDP 60 –min-alleles 2 –max-alleles 2 –minQ 100”. The predicted effects of the SNPs that passed the filtering (Supplementary Table 22) were analysed with the program SNPeff121 using the CanFam3.1.82 database from SNPeff. To find the SNPs located in significant H p and F ST windows, we utilized the results of mapping the windows to the dog genome to extract the variants that were located in dog regions that mapped to our significant windows.

Fine mapping of the region on VVU15

Twenty-five short polymorphic indels (1–7 nucleotides) were identified by analysing the sequences of the re-sequenced foxes aligned to fox scaffold 1. Primers were designed with AmplifX v.1.7.0 (http://crn2m.univ-mrs.fr/pub/amplifx-dist) using the sequence of fox scaffold 1. Forward primers were tagged with fluorescent tags and markers were arranged into five multiplexes (Supplementary Table 19). PCR was performed at a volume of 15 µl using 20 ng of DNA, 1 × Promega GoTaq Colorless Master Mix (Promega), and 0.3 pMol each of the tagged forward and untagged reverse primer. The following conditions were used: 96 °C 2 m; 30 cycles of 96 °C (20 s), 58 °C (20 s), 72 °C (20 s); final extension of 72 °C 1 h. The PCR products were combined post-PCR and analysed on ABI3730 Genetic Analyzer (PE Biosystems). PCR products were sized relative to an internal size standard using ABI GeneMapper 3.5 software package (PE Biosystems). In total, 70 aggressive, 64 tame, 109 F 1 and 537 F 2 individuals were genotyped.

Haploview53 analysis of the tame and aggressive individuals was performed separately to determine the haplotypes in the two populations (Supplementary Fig. 9). Based on the Haploview data and the distances between the genotyped markers, three different sets of markers were chosen for haplotype analysis in the F 2 population. The three maker sets were: upstream (left) of SorCS1 (i13, i16, i17, i19, i20), over SorCS1 (i11, i10, i9, i7, i3, i4, i1, i12) and downstream (right) of SorCS1 (i34, i37, i45, i47, i49, i52) (Supplementary Tables 19, 20). The frequency of the haplotypes for these three marker sets in the tame and aggressive populations were calculated by Haploview, and the F 2 individuals were examined manually using the pedigree information to determine their haplotypes for each marker set.

The haplotype network for the middle haplotypes was calculated using Network 5122. The median-joining method was used to calculate the network, leaving all options at the default settings. All haplotypes that were found by Haploview were used in the calculation (Fig. 4).

The effect of haplotypes on behaviour was analysed in the F 2 population (see Supplementary Note 3 for details). F 2 individuals that were homozygous for any haplotype in any of the three regions (left of SorCS1, at SorCS1 (middle) and right of SorCS1) were identified. The haplotypes that were present in a homozygous state in more than 10 of the F 2 were selected for the analysis of their effect on DPC.1 phenotype46. The D.PC1 values of F 2 individuals from the groups homozygous for different haplotypes were compared using the Kruskal–Wallis test, and, for haplotypes found to be significant with Kruskal–Wallis, a post-hoc Dunn’s test was used to compare individual haplotypes to each other. This analysis used the kruskal.test and dunn.test functions in R.

Karyotype analysis

Chromosome preparation and banding techniques

A fibroblast cell line was established from an ear skin biopsy using conventional techniques123. Metaphase preparations were obtained as previously described29,124,125. Standard G- and C-bandings were made using the methods described in Seabright126 and Sumner127. Chromosomes were identified according to a previous study128.

Fluorescence in situ hybridization

Metaphase chromosomes from the fox primary fibroblast cell line were GTG-stained and captured. Slides were then washed in methanol–acetic acid fixative following xylol treatment. In situ hybridization was performed with a digoxigenin-11-dUTP-labelled (TTAGGG) n telomere repeats probe and a biotin-11-dUTP labelled 18s RNA plus 28 s RNA probe29,129. Hybridization signals were assigned to specific chromosomes or chromosome regions defined by G-banding patterns captured before hybridization.

Image capture

Digital images of the banded metaphase spreads and hybridization signals were captured as described29,125,130 using the VideoTest system with a CCD camera (Jenoptic) mounted on a Zeiss microscope Axioscope 2 (Carl Zeiss). Metaphase spreads images were edited by Corel Paint Shop Pro Photo X2.

Ethics statement

All animal procedures complied with standards for humane care and use of laboratory animals by foreign institutions.

Reporting Summary

Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.

Data availability

The red fox genome assembly and raw reads that were used to generate it are under NCBI project number PRJNA378561. The sequencing data for the tame, aggressive and conventional fox populations are under NCBI project number PRJNA376561. The RNA-seq data is under NCBI/GEO project number GSE76517. Scripts used for all analyses are available upon request.