SNP Characteristics and Data Quality.

We genotyped ≈21,000 autosomal SNPs with Illumina’s Infinium CanineSNP20 BeadChip in a panel of 275 unrelated dogs from 10 phenotypically and genetically diverse breeds (Table 1). SNP markers were uniformly distributed throughout the genome, with a median SNP density of 103.5 ± 124.6 kb. Table 1 provides summary statistics of polymorphism for each breed. Note the average minor allele frequency was ≈25% across breeds, which reflects the ascertainment bias toward common alleles based on the SNP discovery strategy (12). Relationships among breeds were investigated by principal components analysis, which demonstrated that the German Shepherd, Shar-Pei, Beagle, and Greyhound were particularly genetically distinct (Fig. S1).

We performed several analyses to assess SNP data quality. First, four individuals were genotyped in duplicate, and the concordance among genotype calls was >99% across all replicates. Second, for each breed, arrays were performed on a trio of samples and non-Mendelian transmission, indicative of genotyping errors or copy number variants, was assessed. In total, ≈0.4% of markers exhibited Mendelian inconsistencies, consistent with the low genotyping error rate suggested by the replicate arrays. Finally, we assessed the genotyping call rate across all individuals and found uniformly high call rates (≥99%). Thus, these analyses suggest that the genotype data are of high quality.

Signatures of Selection in the Canine Genome.

A large number of statistical tests have been developed to detect deviations from neutrality (10). We developed a population-genomics strategy based on levels of population differentiation, as it is well suited to detect lineage-specific selective events and is robust to whether selection acts on newly arisen or preexisting variation (14). Specifically, for each SNP we defined a statistic, d_i, which is a function of pairwise F_ST (15) between breed i and the remaining breeds. A formal description of d_i is provided in Methods, but in words, d_i measures the standardized locus-specific deviation in levels of population structure for breed i relative to the genome-wide average, summed across all pairwise combinations involving breed i. Large positive values indicate loci, with high levels of population structure relative to the genome-at-large. Thus, it is particularly well suited for detecting selection specific to a particular breed, or subset of breeds, and isolating the direction of change, which is not possible when a single estimate of F_ST is calculated across all populations (16). To attenuate the stochastic variation inherent in single-locus estimates of population structure (17), we performed a sliding-window analysis in which d_i values were averaged in nonoverlapping 1-Mb windows throughout the genome.

The genome-wide distribution of d_i is shown in Fig. 1. We define candidate selection regions as outliers falling in the 99th percentile of the empirical distribution of d_i. In total, 155 out of the 1,933 windows met this criterion in one or more of the 10 breeds (Table S1). Several observations suggest that our set of outlier loci is enriched for targets of selection. First, all five genes that have been mapped to date through large-scale association studies of hallmark breed traits are among our list of most differentiated regions: IGF1 in breeds of small size (7), a locus on CFA 18 that is responsible for the characteristic short-limb phenotype in Daschshunds and other breeds (8), and three genes (RSPO2, FGF5, and KRT71) that influence coat phenotypes in many breeds (9).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Genomic distribution of population structure in 10 dog breeds. The distribution of d_i for each 1-Mb interval across all autosomes is shown for each breed. Alternating gray and black indicate values in d_i from adjacent chromosomes. The dashed red line denotes the 99th percentile for each breed. Breeds are abbreviated as described in Table 1.

Second, we performed extensive coalescent simulations that take into account SNP ascertainment and major demographic features, such as population structure and breed-specific bottlenecks. The neutral coalescent model closely recapitulates many characteristics of the observed data, such as average pairwise F_ST, average number of markers per 1-Mb window and distribution of minor allele frequencies (Fig. S2). The observed data contains significantly more highly differentiated loci (P = 1.3 × 10⁻⁷) compared with the simulated data.

Third, we observed a significant enrichment of signatures of selection in or around genes relative to putatively neutrally evolving regions (P = 2.95 × 10⁻³), which is expected if adaptive variation is overrepresented in genic regions. Similar observations have also been made in genome-wide scans of selection in humans (18, 19). Collectively, these observations support the hypothesis that the most differentiated regions of the canine genome are enriched for targets of selection.

Shared Versus Unique Selective Events.

To investigate how frequently selective events were unique or shared among breeds, we calculated the number of overlapping signatures of selection for each of the 155 significant 1-Mb windows (Fig. 2A). Approximately 103 of the 155 significant windows (∼66%) were observed in just one or two breeds (Fig. 2A). These loci likely contain genes that confer breed-restricted phenotypes, such as skin wrinkling in the Shar-Pei (see below). Conversely, 16 of the 155 significant windows (∼10%) exhibited signatures of selection in five or more breeds. Such pervasive differentiation at a single locus is consistent with the action of a gene that generally sorts individuals into phenotypic classes and breed groups. For example, one window with strong evidence of selection in multiple breeds is located on CFA15 (43.6–44.6 Mb) and contains the IGF1 gene, which governs the miniature size of breeds in the “toy” group (7). Interestingly, a region on CFA 3 (44.6–45.6 Mb) that includes the IGF1R gene also shows a strong signature of selection in the Dachshund and Brittany, suggesting that multiple steps in the insulin growth-factor signaling pathway have been substrates of artificial selection in dogs.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Shared versus unique signatures of selection. (A) The number of overlapping signatures of selection in each 1-Mb window is shown. We define an overlapping signature of selection for each window if the empirical P value is ≤ 0.01 in one breed and ≤ 0.05 in another breed. Alternating white and vertical light yellow rectangles indicate adjacent chromosomes. The red arrow indicates the chromosomal region shown in B. (B) (Upper) Sliding-window analyses of pairwise F_ST among German Shepherds, Jack Russell Terriers, and Beagles. Gray boxes indicate two distinct peaks of differentiation (at ≈10.3–10.4 Mb and 11.3–11.4 Mb). Note the different patterns of pairwise F_ST in the two peaks of differentiation. Additional breed comparisons have been omitted for clarity. (Lower) Unrooted neighbor joining trees are shown for all breeds inferred from markers in each of the peaks of differentiation above. Note the distinct topology between the two trees (BGL, JRT, DSH, and BRT lineages are shown in blue).

One of the most differentiated regions of the canine genome that shows evidence of selection in multiple breeds occurs in three contiguous windows on CFA 10 (Fig. 2B). Sliding-window analyses of pairwise F_ST across the 3-Mb interval suggests two or more independent selective events, reflected by two peaks of differentiation with distinct patterns of allele frequency divergence among breeds (Fig. 2B). The peak of differentiation observed from 11.2 to 11.3 Mb coincides with the HMGA2 gene, whose protein product is an integral component of enhanceosomes and regulates gene expression (20). In mice, mutations in HMGA2 result in the pygmy phenotype (21), characterized by aberrations in adiposity and disrupted growth leading to dwarfism. In our data, the small-sized breeds (Dachshund, Beagle, Jack Russell Terrier, and Brittany) show high levels of differentiation at HMGA2 compared to the larger-sized breeds (Fig. 2B). At the most differentiated SNP near HMGA2, allele frequency is significantly correlated with body weight (Pearson r² = 0.68, P = 0.003). Thus, HMGA2 is a strong candidate for mediating variation in size among dogs. The second peak of differentiation in the CFA 10 region (10.35–10.45 Mb), in which the German Shepherd, Jack Russell Terrier, Border Collie, and Greyhound are strongly differentiated from the Dachshund, Beagle, Brittany, and Shar-Pei (Fig. 2B), overlaps two genes, GNS and RASSF3. GNS is a particularly interesting candidate as its avian ortholog (QSulf1), regulates WNT signaling during embryogenesis in myogenic somite progenitors (22).

Overview of Candidate Selection Genes.

The 155 candidate selection loci contain 1,630 known or predicted protein-coding genes. To obtain a broad overview into the molecular functions of these genes and to test the hypothesis that particular functional classes are enriched in the most differentiated regions of the canine genome, we performed a gene ontology (GO) analysis. Table S2 summarizes GO molecular function and biological process terms that are significantly enriched among genes in candidate-selection regions. Similar to analyses of selection in natural populations (23), we find that genes involved in immunity and defense are also significantly overrepresented in the 155 candidate selection regions. This is somewhat surprising, as natural and artificial selection would not necessarily be expected a priori to act on similar classes of genes, and suggests that immune related genes are pervasive targets of selection because of their critical role in pathogen defense or propensity for pleiotropic effects (24).

The average number of genes in each of the 155 candidate selection regions was ≈11. Thus, it is difficult to precisely identify the specific gene that has been influenced by selection. Nonetheless, the 155 most differentiated loci possess many strong candidate genes that influence phenotypes that vary conspicuously among breeds, such as size (HMGA2 and IGF1R), coat color and texture (SILV and MITF), behavior (CDH9, DRD5, and HTR2A), skeletal morphology (SOX9), and physiology (FTO, SLC2A9, and SLC5A2).

However, more definitive inferences can be made for eight regions, in which there is only a single protein-coding gene located within the interval (Table S3). Possible phenotypes that each gene may influence are listed in Table S3, and more detailed information is provided in the SI Text. Interestingly, three of the eight genes are transcription factors (ZFHX3, SOX9, and SATB1). There has been considerable debate about the relative contribution of changes in gene regulation versus protein structure as mechanisms of evolutionary change (25, 26). Similar to analyses of artificial selection in other domesticated species (27–29), our data suggest that tinkering (30) with gene-expression networks may have played a prominent role in the rapid phenotypic diversification of modern dog breeds.

We note that the stringent threshold used to define candidate selection regions has likely excluded genuine substrates of selection. For example, regions on CFA 9 and 27 lie just beyond our threshold of significance in Poodles (empirical P-values = 0.021 and 0.014, respectively). These regions contain numerous keratin gene family members, which are important structural proteins of the skin, nails, and hair (Fig. S3). Of particular interest are members of the type-I hair keratins on CFA 9 (KRT25, KRT27, KRT28, KRT32, KRT35, and KRT36) and type-II hair keratins on CFA 27 (KRT71, KRT72, KRT73, KRT74, KRT82, KRT84, and KRT85). Recently, variation in KRT71 has been associated with curly coat phenotypes in several breeds (9), which validates our CFA 27 results. Our data suggest that additional keratin genes on CFA 9 are also strong candidates for contributing to the curly coat phenotype.

Regulatory Variation in HAS2 Is Associated with Skin Wrinkling of Shar-Peis.

To characterize candidate selection genes in more detail, we focused on a region on CFA 13 with evidence of selection in the Shar-Pei (Fig. 3A) that contains three genes (SNTB1, FTSJ1, and HAS2). A distinguishing characteristic of the Shar-Pei is cutaneous mucinosis, or excessive skin wrinkling. The degree of skin folds correlates with high mucin content histologically and elevated levels of hyaluronic acid biochemically (31). HAS2, which is a hyaluronic acid synthase, was thus a strong candidate gene. In addition, rare mutations in human HAS2 have been described that result in severe cutaneous mucinosis (32).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Genetic variation in HAS2 is associated with skin wrinkling in Shar-Pei. (A) Single locus estimates of F_ST between Shar-Pei and Dachshund across a 1-Mb window. Similar patterns were observed for Shar-Pei compared to other breeds, but have been omitted for clarity. The location of all protein-coding genes are shown as rectangular boxes. (B) An example of smooth (Left) and wrinkled (Right) Shar-Pei dogs. (C) Exon structure of HAS2. Conservation values obtained from the University of California Santa Cruz genome browser are shown below. Black horizontal lines indicate sequenced regions. (D) Genotype frequencies of the intron 2 indel (site 13805 in Table S4) in smooth and wrinkled Shar-Pei, which are significantly different (P = 6.28 × 10⁻⁵). Deletion and insertion alleles are denoted as “D” and “d,” respectively.

To test the hypothesis that genetic variation in HAS2 contributes to skin wrinkling, we exploited the intrabreed phenotypic variation that exists in the degree of wrinkling within Shar-Pei (Fig. 3B). Specifically, we sequenced ≈3.7 kb of HAS2 (including all exons, intron/exon boundaries, and untranslated regions) (Fig. 3C) from 32 wrinkled and 18 smooth-coated purebred Shar-Pei (Fig. 3B). In total, we discovered five polymorphisms, none of which are located in coding regions (Table S4). One of the upstream polymorphisms was nearly fixed in both wrinkled and smooth dogs and was not considered further (Table S4). Association mapping was performed with a permutation-based Cochran-Armitage trend test (33) on the four remaining SNPs, all of which demonstrated significant differences in genotype frequencies between wrinkled versus smooth Shar-Pei (Table S4).

We next sequenced all of the HAS2 amplicons in a diverse panel of 94 dogs derived from 20 breeds (Table S5). The most differentiated SNP between the Shar-Pei and other breeds is a 2-bp indel ≈86 bp 3′ of exon 2 (Table S5). The deletion allele is significantly associated with the wrinkling phenotype (P = 6.28 × 10⁻⁵; see site 13805 in Table S4 and Fig. 3D), where the frequency of the deletion allele is ≈0.91 and 0.53 in wrinkled and smooth Shar-Pei, respectively. The deletion allele is rare outside of the Shar-Pei (∼1.6%) (Table S5) and no homozygous deletions were found in any of the 94 dogs.

Although experimental studies will ultimately be necessary to determine whether the polymorphisms described in Table S4 are functionally important, it seems unlikely that any of them are causally related to skin wrinkling in the Shar-Pei. The most strongly associated polymorphism (site −424) is common across breeds (Table S5). Even though the intron 2 polymorphism does possess patterns of variation between the Shar-Pei and other breeds expected for a causal polymorphism, it is not located in a region of high sequence conservation and is not embedded in any obvious regulatory elements. Therefore, we hypothesize that the polymorphisms in Table S4, and in particular the 2-bp intron 2 indel, are in linkage disequilibrium (LD) with the unidentified causative allele. As no variation was found in the HAS2 coding region, the causal allele is likely a regulatory polymorphism. Consistent with this hypothesis, of the 50 Shar-Pei dogs in the resequencing panel, 22 overlap with the set of individuals used for Illumina SNP genotyping. As shown in Fig. S4, the strongest associations for this subset of individuals in the SNP data occur upstream of the HAS2 gene, suggesting the causal polymorphism lies 5′ to HAS2.