Rabbits softly swept to domestication When people domesticate animals, they select for tameness and tolerance of humans. What else do they look for? To identify the selective pressures that led to rabbit domestication, Carneiro et al. sequenced a domestic rabbit genome and compared it to that of its wild brethren (see the Perspective by Lohmueller). Domestication did not involve a single gene changing, but rather many gene alleles changing in frequency between tame and domestic rabbits, known as a soft selective sweep. Many of these alleles have changes that may affect brain development, supporting the idea that tameness involves changes at multiple loci. Science, this issue p. 1074; see also p. 1000

Abstract The genetic changes underlying the initial steps of animal domestication are still poorly understood. We generated a high-quality reference genome for the rabbit and compared it to resequencing data from populations of wild and domestic rabbits. We identified more than 100 selective sweeps specific to domestic rabbits but only a relatively small number of fixed (or nearly fixed) single-nucleotide polymorphisms (SNPs) for derived alleles. SNPs with marked allele frequency differences between wild and domestic rabbits were enriched for conserved noncoding sites. Enrichment analyses suggest that genes affecting brain and neuronal development have often been targeted during domestication. We propose that because of a truly complex genetic background, tame behavior in rabbits and other domestic animals evolved by shifts in allele frequencies at many loci, rather than by critical changes at only a few domestication loci.

Domestication of animals (that is, the evolution of wild species into tame forms) has resulted in notable changes in behavior, morphology, physiology, and reproduction (1). The genetic underpinnings of the initial steps of animal domestication are poorly understood but probably involved changes in behavior that allowed the animals to survive and reproduce under conditions that might be too stressful for wild animals. Given the differences in behavior between wild and domesticated animals, we investigated to what extent this process involved fixation of new mutations with large phenotypic effects, as opposed to selection on standing variation. Such studies are hampered in most domestic animals due to ancient domestication events, extinct wild ancestors, or geographically widespread wild ancestors.

Rabbit domestication was initiated in monasteries in southern France as recently as ~1400 years ago (2). At this time, wild rabbits were mostly restricted to the Iberian Peninsula, where two subspecies occurred (Oryctolagus cuniculus cuniculus and O. c. algirus), and to France, colonized by O. c. cuniculus (Fig. 1B). Additionally, the area of origin of domestic rabbits is still populated with wild rabbits related to the ancestors of the domestic rabbit (3). This recent and well-defined origin provides a major advantage for inferring genetic changes underlying domestication.

Fig. 1 Experimental design and population data. (A) Images of the six rabbit breeds included in the study (sized to reflect differences in body weight) and of a wild rabbit. (B) Map of the Iberian Peninsula and southern France with sample locations marked (orange dots). Demographic history of this species is indicated, and a logarithmic time scale is shown at right. The hybrid zone between the two subspecies is marked with dashes. (C) Nucleotide diversities in domestic and wild populations. The French (FRW1 to FRW3) and Iberian (IW1 to IW11) wild rabbit populations are ordered according to a northeast-to-southwest transection. Sample locations are provided in table S4.

We performed Sanger sequencing and assembly of a female rabbit genome (4). The draft OryCun2.0 assembly size is 2.66 Gb, with a contig N50 size of 64.7 kb and a scaffold N50 size of 35.9 Mb (tables S1 and S2). The genome assembly was annotated using the Ensembl gene annotation pipeline (Ensembl release 73, September 2013) and with both rabbit RNA sequencing data and the annotation of human orthologs (4) (table S3). Our analysis of rabbit domestication used Ensembl annotations as well as a custom pipeline for annotation of untranslated regions (UTRs) (168,286 distinct features), noncoding RNA (n = 9666), and noncoding conserved elements (2,518,476 distinct features).

To identify genomic regions under selection during domestication, we performed whole-genome resequencing (10× coverage) of pooled samples (table S4) of six different breeds of domestic rabbits (Fig. 1A), 3 pools of wild rabbits from southern France, and 11 pools of wild rabbits from the Iberian Peninsula, representing both subspecies (Fig. 1B). We also sequenced a close relative, the snowshoe hare (Lepus americanus), to deduce the ancestral state at polymorphic sites. Short sequence reads were aligned to our assembly; single-nucleotide polymorphism (SNP) calling resulted in the identification of 50 million high-quality SNPs and 5.6 million insertion/deletion polymorphisms after filtering (table S5). The numbers of SNPs at noncoding conserved sites and in coding sequences were 719,911 and 154,489, respectively. The per-site nucleotide diversity (π) within populations of wild rabbits was in the range of 0.6 to 0.9% (Fig. 1C). Thus, the rabbit is one of the most polymorphic mammals sequenced so far, presumably due to a larger long-term effective population size relative to other sequenced mammals (5). Identity scores confirm that the domestic rabbit is most closely related to wild rabbits from southern France (fig. S1A), and we inferred a strong correlation (r = 0.94) in allele frequencies at most loci between these groups (fig. S1B). The average nucleotide diversity of each sequenced population is consistent with a bottleneck and reduction in genetic diversity when rabbits from the Iberian Peninsula colonized southern France and a second bottleneck during domestication (3) (Fig. 1, B and C).

Selective sweeps occur when beneficial genetic variants increase in frequency due to positive selection together with linked neutral sequence variants (6). This results in genomic islands of reduced heterozygosity and increased differentiation between populations around the selected site. We compared genetic diversity between domestic rabbits as one group to wild rabbits representing 14 different locations in France and the Iberian Peninsula. We calculated fixation index (F ST ) values between wild and domestic rabbits and average pooled heterozygosity (H) in domestic rabbits in 50-kb sliding windows across the genome (hereafter referred to as the F ST -H outlier approach). We identified 78 outliers with F ST > 0.35 and H < 0.05 (Fig. 2A, fig. S2, database S1). We also used SweepFinder (7), which calculates maximum composite likelihoods for the presence of a selective sweep, taking into account the background pattern of genetic variation in the data and with a significance threshold set by coalescent simulations incorporating the recent demographic history of domestic and wild rabbits (figs. S3 and S4 and databases S1 and S2) (4). This analysis resulted in the identification of 78 significant sweeps (false discovery rate = 5%) (Fig. 2A, database S1). Thirty-one (40%) of these were also detected with the F ST -H approach (Fig. 2A). This incomplete overlap is probably explained by the fact that SweepFinder primarily assesses the distribution of genetic diversity within the selected population, whereas the F ST -H analysis identifies the most differentiated regions of the genome between wild and domestic rabbits. We carried out an additional screen using targeted sequence capture on an independent sample of individual French wild and domestic rabbits. We targeted more than 6 Mb of DNA sequence split into 5000 1.2-kb intronic fragments that were distributed across the genome and selected independently of the genome resequencing results above. Coalescent simulations, using the targeted resequencing data set and incorporating the recent demographic history of domestic rabbit as a null model (figs. S3 and S4 and databases S1 and S2) (4), revealed that the majority of the sweep regions detected by whole-genome resequencing showed levels and patterns of genetic variation that were observed less than 5% of the time in the simulated data set (76.0% with SweepFinder and 73.7% with F ST -H outlier regions, excluding regions without targeted fragments), a highly significant overlap (Fisher’s exact test, P < 1 × 10–5 for both tests). Furthermore, 26 of the 31 sweep regions detected with both SweepFinder and the F ST -H approach were targeted in the capture experiment, and an even greater proportion (88.5%) showed levels and patterns of genetic variation unlikely to be generated under the specified demographic model.

Fig. 2 Selective sweep and Δ allele frequency analyses. (A) Plot of F ST values between wild and domestic rabbits. Sweeps detected with the F ST -H outlier approach, SweepFinder, and their overlaps are marked on top. Unassigned scaffolds were not included in the analysis. (B and C) Selective sweeps at GRIK2 (B) and SOX2 (C). Heterozygosity plots for wild (red) and domestic (black) rabbits together with plots of F ST values and SNPs with ΔAF > 0.75 (HΔAF). The bottom panels show putative sweep regions, detected with the F ST -H outlier approach and SweepFinder, marked with horizontal bars. Gene annotations in sweep regions are indicated: * represents ENSOCUT000000; **SOX2-OT represents the manually annotated SOX2 overlapping transcript (4). (D) The majority of SNPs showed low ΔAF between wild and domestic rabbits. The black line indicates the number of SNPs in nonoverlapping ΔAF bins (left y axis). Colored lines denote M values (log2-fold changes) of the relative frequencies of SNPs at noncoding evolutionary conserved sites (blue), in UTRs (red), exons (yellow), and introns (green), according to ΔAF bins (right y axis). M values were calculated by comparing the frequency of SNPs in a given annotation category in a specific bin with the corresponding frequency across all bins. (E) Location of SNPs at conserved noncoding sites with ΔAF ≥ 0.8 SNPs (n = 1635) and ΔAF < 0.8 SNPs (n = 502,343) in relation to the TSS of the most closely linked gene. **P < 0.01.

An example of a selective sweep overlapping the 3′-part of GRIK2 (glutamate receptor, ionotropic, kainate 2) is shown in Fig. 2B. Parts of this region have low heterozygosity in domestic rabbits, and at position chr12:90,153,383 base pairs, domestic rabbits carry a nearly fixed derived allele at a site with 100% sequence conservation among 29 mammals except for the allele present in domestic rabbits (8), suggesting functional importance. GRIK2 encodes a subunit of a glutamate receptor that is highly expressed in the brain and has been associated with recessive mental retardation in humans (9). Both SweepFinder and the F ST -H outlier analysis identified two sweeps near SOX2 (SRY-BOX 2), separated by a region of high heterozygosity (Fig. 2C). SOX2 encodes a transcription factor that is required for stem cell maintenance (10).

Given the comprehensive sampling in our study and the correlation in allele frequencies between domestic and French wild rabbits (fig. S1B), highly differentiated individual SNPs are likely either to have been directly targeted by selection or to occur in the vicinity of loci under selection. For each SNP, we calculated the absolute allele frequency difference between wild and domestic rabbits (ΔAF) and sorted these into 5% bins (ΔAF = 0 to 0.05, etc.). The majority of SNPs showed low ΔAF between wild and domestic rabbits (Fig. 2D). We examined exons, introns, UTRs, and evolutionarily conserved sites for enrichment of SNPs with high ΔAF, as would be expected under directional selection on many independent mutations (Fig. 2D and table S6). We observed no consistent enrichment for high ΔAF SNPs in introns, but we found significant enrichments in exons, UTRs, and conserved noncoding sites (χ2 test, P < 0.05). We detected a significant excess of SNPs at conserved noncoding sites for each bin ΔAF > 0.45 (χ2 test, P = 1.8 × 10–3 to 7.3 × 10–17), whereas in coding sequence, a significant excess was found only at ΔAF > 0.80 (χ2 test, P = 3.0 × 10–2 to 1.0 × 10–3). Compared to the relative proportions in the entire data set, there was an excess of 3000 SNPs at conserved noncoding sites with ΔAF > 0.45, whereas for exonic SNPs with ΔAF > 0.80, the excess was only 83 SNPs (table S6). Thus, changes at regulatory sites have played a much more prominent role in rabbit domestication, at least numerically, than changes in coding sequences.

We selected the 1635 SNPs at conserved noncoding sites with ΔAF > 0.80, which represent 681 nonoverlapping 1-Mb blocks of the rabbit genome. So as not to inflate significances due to inclusion of SNPs in strong linkage disequilibrium, we selected only one SNP per 50 kb, leaving 1071 SNPs. More than 60% of the SNPs were located 50 kb or more from the closest transcriptional start site (TSS) (Fig. 2E), suggesting that many differentiated SNPs are located in long-range regulatory elements. A gene ontology (GO) overrepresentation analysis (11) examining all genes located within 1 Mb from high-ΔAF SNPs showed that the most enriched categories of biological processes involved “cell fate commitment” (Bonferroni P = 3.1 × 10–3 to 5.4 × 10–5) (Table 1 and database S3), whereas the statistically most significant categories involved brain and nervous system cell development (Bonferroni P = 2.9 × 10–3 to 3.7 × 10–10). Many of the mouse orthologs of genes associated with noncoding high-ΔAF SNPs were expressed in the brain or sensory organs, and this enrichment was highly significant (Table 1). We also examined phenotypes observed in mouse mutants (www.informatics.jax.org) for these genes, revealing a significant (Bonferroni P = 3.7 × 10–2 to 7.5 × 10–17) enrichment of genes associated with defects in brain and neuronal development, development of sensory organs (hearing and vision), ectoderm development, and respiratory system phenotypes (fig. S5). These highly significant overrepresentations were obtained because there were many genes in the overrepresented categories (Table 1). For example, we observed high-ΔAF SNPs associated with 191 genes (113 expected by chance) from the nervous system–development GO category (Bonferroni P = 3.7 × 10–10). Thus, rabbit domestication must have a highly polygenic basis with many loci responding to selection and where genes affecting brain and neuronal development have been particularly targeted.

Table 1 Summary of results from enrichment analysis of ΔAF > 0.8 SNPs located in conserved noncoding elements. One significantly enriched term was chosen from each group of significantly enriched intercorrelated terms. Full lists of enriched terms and intercorrelations are presented in database S3, and the most enriched intercorrelated terms are presented in fig. S5. P values are Bonferroni-corrected. O/R, number of distinct nonoverlapping 1-Mb windows observed (O) and the average number of 1-Mb windows observed in 1000 random (R) samplings of the same number of genes (rounded to the nearest integer). TS, Thieler stage. View this table:

None of the coding SNPs that differed between wild and domestic rabbits was a nonsense or frame-shift mutation, consistent with data from chicken (12) and pigs (13), suggesting that gene loss has not played a major role during animal domestication. This is an important finding, as it has been suggested that gene inactivation could be an important mechanism for rapid evolutionary change, for instance, during domestication (14). Of 69,985 autosomal missense mutations, there were no fixed differences, and only 14 showed a ΔAF above 90%. On the basis of poor sequence conservation, similar chemical properties of the substituted amino acids, and/or the derived state of the domestic allele, we assume that most of these result from hitchhiking rather than being functionally important (database S4). However, two missense mutations stand out; these may be direct targets of selection because at these two positions the domestic rabbit differs from all other sequenced mammals (>40 species). The first is a Gln813→Arg813 substitution in TTC21B (tetratricopeptide repeat domain 21B protein), which is part of the ciliome and modulates sonic hedgehog signaling during embryonic development (15). The other is an Arg1627→Trp1627 substitution in KDM6B (lysine-specific demethylase 6B), which encodes a histone H3K27 demethylase involved in HOX gene regulation during development (16).

Deletions distinct to domestic rabbits were difficult to identify because the genome assembly is based on a domestic rabbit, but some convincing duplications were detected with marked frequency differences between wild and domestic rabbits (database S5). We observed a one–base pair insertion/deletion polymorphism located within an intron of IMMP2L (inner mitochondrial membrane peptidase-2 like protein), where domestic and wild rabbits were fixed for different alleles. The polymorphism occurs in a sweep region and is the sequence polymorphism with highest ΔAF in the region (fig. S6). Mutations in IMMP2L have been associated with Tourette syndrome and autism in humans (17).

Cell fate determination was a strongly enriched GO category (enrichment factor = 4.9) (database S3) for genes near variants with high ΔAF. We examined the functional importance of 12 SOX2, 4 KLF4, and 1 PAX2 high ΔAF SNPs associated with this GO category and where all 17 SNPs were distinct to domestic rabbits compared with other sequenced mammals. Electrophoretic mobility shift assay (EMSA) with nuclear extracts from mouse embryonic stem cell–derived neural stem cells revealed specific DNA-protein interactions (Fig. 3, fig. S7, table S7). Four probes, all from the SOX2 region, showed a gel shift difference between wild and domestic alleles. Nuclear extracts from a mouse P19 embryonic carcinoma cell line before and after neuronal differentiation recapitulated these four gel shifts and revealed three additional probes, one in PAX2 and two more in SOX2, that showed gel shift differences between wild-type and mutant probes only after neuronal differentiation. Thus, altered DNA-protein interactions were identified for 7 of the 17 high ΔAF SNPs that we tested, qualifying them as candidate causal SNPs that may have contributed to rabbit domestication.

Fig. 3 Bioinformatic and functional analysis of candidate causal mutations. Three examples of SNPs near SOX2 and PAX2 where the domestic allele differs from other mammals. The locations of these three SNPs assessed with EMSA are indicated by red crosses on top. EMSA with nuclear extracts from embryonic stem cell–derived neural stem cells or from mouse P19 embryonic carcinoma cells before (un-diff) or after neuronal differentiation (diff) are shown for three SNPs. Exact nucleotide positions of polymorphic sites are indicated. Allele-specific gel shifts are indicated by arrows. WT, wild-type allele; Dom, domestic, the most common allele in domestic rabbits. Cold probes at 100-fold excess were used to verify specific DNA-protein interactions.

Our results show that very few loci have gone to complete fixation in domestic rabbits and none at coding sites or at noncoding conserved sites. However, allele frequency shifts were detected at many loci spread across the genome, and the great majority of domestic alleles were also found in wild rabbits, implying that directional selection events associated with rabbit domestication are consistent with polygenic and soft sweep modes of selection (18) that primarily acted on standing genetic variation in regulatory regions of the genome. This stands in contrast with breed-specific traits in many domesticated animals that often show a simple genetic basis with complete fixation of causative alleles (19). Our finding that many genes affecting brain and neuronal development have been targeted during rabbit domestication is fully consistent with the view that the most critical phenotypic changes during the initial steps of animal domestication probably involved behavioral traits that allowed animals to tolerate humans and the environment humans offered. On the basis of these observations, we propose that the reason for the paucity of specific fixed domestication genes in animals is that no single genetic change is either necessary or sufficient for domestication. Because of the complex genetic background for tame behavior, we propose that domestic animals evolved by means of many mutations of small effect, rather than by critical changes at only a few domestication loci.