DNA methylation in canine brains is related to domestication and dog-breed formation

Epigenetic factors such as DNA methylation act as mediators in the interaction between genome and environment. Variation in the epigenome can both affect phenotype and be inherited, and epigenetics has been suggested to be an important factor in the evolutionary process. During domestication, dogs have evolved an unprecedented between-breed variation in morphology and behavior in an evolutionary short period. In the present study, we explore DNA methylation differences in brain, the most relevant tissue with respect to behavior, between wolf and dog breeds. We optimized a combined method of genotype-by-sequencing (GBS) and methylated DNA immunoprecipitation (MeDIP) for its application in canines. Genomic DNA from the frontal cortex of 38 dogs of 8 breeds and three wolves was used. GBS and GBS-MeDIP libraries were prepared and sequenced on Illuma HiSeq2500 platform. The reduced sample represented 1.18 ± 0.4% of the total dog genome (2,4 billion BP), while the GBS-MeDIP covered 11,250,788 ± 4,042,106 unique base pairs. We find substantial DNA methylation differences between wolf and dog and between the dog breeds. The methylation profiles of the different groups imply that epigenetic factors may have been important in the speciation from dog to wolf, but also in the divergence of different dog breeds. Specifically, we highlight methylation differences in genes related to behavior and morphology. We hypothesize that these differences are involved in the phenotypic variation found among dogs, whereas future studies will have to find the specific mechanisms. Our results not only add an intriguing new dimension to dog breeding but are also useful to further understanding of epigenetic involvement.


Introduction
During animal domestication, wild species are adapted to a life in human proximity. Through correlated selection responses to tameness selection, differences in phenotype can be acquired in only a few generations [1][2][3]. Even though behavior is profoundly modified during domestication, differences in DNA sequence between domesticates and their wild ancestor species are relatively few, while differences in brain gene expression are all the larger [4]. For

Ethical statement
Canine brains were donated by their owners (dogs) and by the zoo (wolves), and all samples were collected in connection with veterinary motivated procedures. No ethical license was therefore required.

Subjects and sampling
Brain tissue from 41 canine individuals was used in this study; three female wolves (Canis lupus) and 38 dogs (Canis familiaris) of eight different breeds, which represents the samples available at the time. . Dog samples were donated by owners after the dogs had been euthanized by veterinarians after decision by the owner. Samples were collected by Dog Genetics Project, Swedish University of Agricultural Sciences and Uppsala University, in connection to a previous study [25]. The ages of the dogs varied from one to ten years, with the exception of one rottweiler (3 months), and three German shepherds (6, 7 and 8 months). None of the dogs had been euthanized due to any brain disease, and the most common reasons were tumors outside of the nervous system and cardiovascular diseases. None of the dogs had metastases in the brains or circulating tumor emboli according the autopsy protocols. The brains were stored in a -80 freezer at Uppsala University until our sample dissection. Three wolf brain samples were donated from Borås Zoo in connection with routine euthanasia of two 2 years old and one 9 year old surplus animals.
For the present study, a piece of tissue was dissected from the medial prefrontal cortex of the left cerebral hemisphere. This represents a part of the brain generally believed to be involved in cognitive aspects of behaviour. The brain tissue was kept frozen either on dry ice, in liquid nitrogen or in a -80 freezer until DNA extraction.

DNA isolation
For DNA extraction, we used Qiagen All-Prep RNA/DNA kit and followed the instruction from the manufacturer. DNA quality was measured using ThermoFischer Scientific Nano-Drop1 ND-2000c and concentration was quantified using a fluorometer (Qubit1 Fluorometric Quantitation).

GBS and MeDIP
DNA from the prefrontal cortex was used for the genomic and epigenomic analyses. Because pyramidal cells are the most common neurons in the cerebral cortex, being its building block [34], it is expected that most of the epigenetic signal obtained belongs to these cells. To identify differences in methylation patterns, we have used a novel approach that combines genotyping by sequencing (GBS) [35] and methylated DNA immunoprecipitation (MeDIP) [36]. This method has previously been used in chickens (Gallus gallus) [37] and is here further optimized for application in canines based on in silico and in vitro digestion tests. Reducing the genome and pooling barcoded individual samples makes this a cost-efficient genotyping and epigenotyping method that allows comparisons of methylation profiles. Importantly, the GBS reduced representation method cleaves the genome at the same recognition sites in all individuals, which results in a non-random small fraction that constitutes a broadly representative genome sample. This sample can then be used for phylogenetic comparisons and also as input to represent the genetic background of each individual in the methylation analysis. The restriction sites used by the restriction enzyme PstI are unrelated to CpGs and, consequently, the positions are unbiased towards CpG islands.
To give an overview of the GBS method, the DNA is cleaved using PstI (Thermo Scientific) and enriched for fragments sized 200-500 basepairs (bp) long, suitable for Illumina sequencing. A DNA barcode, unique for each individual sample, together with a common adapter for Illumina sequencing barcoding system, is ligated to the fragments (Poland and Rife, 2012). The output of the barcoding and pooling is used as input for MeDIP. Two libraries are then prepared, one GBS library and one GBS-MeDIP library, and both libraries are sequenced. In our samples, sequencing was paired-end sequenced (read length: 125 bp) on the Illumina HiSeq2500 platform and performed at SNP&SEQ Technology Platform, SciLifeLab, Uppsala, Sweden. A more detailed description of the procedure can be found in S1 Table.

Bioinformatic analyses
We used CASAVA (Illumina) for the initial processing of the samples and converted ".bcl" (base call files) to ".fastq" extensions. These are compatible with programs used for alignment reading. The quality of the reads was checked using FastQC v.0.11.3 34 and we performed quality trimming in short read sequences. Quality-trimmed reads were aligned against the canine reference genome (CanFam3, NCBI) using Bowtie2 tool v.2-2.2.9 [38] with default parameters for very-sensitive-local alignment. Coverage depth was determined by using Samtools v.1.3.1 [39].
TASSEL-GBS Discovery Pipeline [40] was used to process the GBS data. For SNP calling, default filtering parameters were used, except for 5% for minimum minor allele frequency (mnMAF), 70% of minimum taxon coverage (mnTCov), and 70% for minimum site coverage (mnScov). With Tassel we checked for allele changes between wolves and dogs and among dog breeds by comparing the number of fixed SNPs within each group. Tassel was additionally used to create a cladogram. It was generated by Neighbor Joining distance matrix and was plotted using Archaeopteryx tree.
For epigenetics analysis, reads from fastq files were demultiplexed using Stacks v.1.46 [41]. Uncalled and low-quality score bases were eliminated using process radtags function from Stacks. MEDIPS R-package from Bioconductor [42] was used for basic data processing, quality controls, normalization, and identification of differentially methylated regions (DMRs). We followed the same specific parameters from MEDIP package as described previously in Pértille, Brantsaeter [37]. The genome was divided into adjacent windows of pre-defined length size of 100 bp and differential methylation analysis used a weighted trimmed mean of the log expression ratios (trimmed mean of M values) [43].
We compared DNA methylation of the adjacent 100 bp windows between several groups: 1) wolves were compared to all female dogs, 2) wolves were compared to females of each breed, 3) female dogs were compared to male dogs, and 4) breeds were compared to breeds, including both male and female individuals. The reason for only comparing wolves to female dogs was that our samples only included female wolves, and we wanted to avoid finding sexspecific effects instead of species effects. DMRs with p<0.0005 were considered significant. Manhattan plots were created using the R package qqman [44].
The DMRs were annotated against the dog reference genome (CanFam3, NCBI) using the Variant Effect Predictor (VEP) tool [45]. The annotations were analyzed through different online bioinformatic tools. Reactome pathway browser [ [46]; reactome.org] was used to find overrepresented pathways related to genes in the wolf-dog comparison. An overrepresentation enrichment gene ontology analysis was performed using WEB-based GEne SeT AnaLysis

Sequencing results
On average, the GBS resulted in a coverage of 28,226,049 ± 9,895,466 (SD) unique base pairs (breadth) which represents 1.18 ± 0.4% of the total dog genome (2,4 billion BP). The GBS-Me-DIP covered 11,250,788 ± 4,042,106 unique base pairs which represents 0.47 ± 0.17% of the total genome. Table 1 presents average coverage for each dog breed and the wolf for both GBS and GBS-MeDIP.
The regions obtained via GBS include 1,001,135 CpG regions which represents 3.98 ± 1.48% of the total number of CpGs (25 millions) in the canine genome. GBS-MeDIP include 149,275 ± 93,255 CpGs, 0.59 ± 0.29% of the total. The ratios (enrichment scores) between the number of observed CpGs and the expected from the reference genome were 1.37 and 2.16 for GBS and GBS-MeDIP, respectively, indicating an enrichment in CpGs for the GBS-MeDIP.
A neighbor-joining tree was plotted based on 62,452 single nucleotide polymorphisms (SNPs) across all the 41 individuals from the GBS sequencing results. Although we only used a reduced genome, individuals from the same species/breed grouped together (except for two outliers), indicating that the GBS produced results representative of the genetic variability in each breed (Fig 1).
The number of SNPs in the GBS reduced genome that was fixed within species and breed groups was compared across groups. A total of 2479 SNPs with allele changes were found between the wolf and all dogs combined. The results from additional comparisons are presented in Table 2. For example, in wolf versus breed comparisons, boxer was the breed with greatest number of allele differences (2940 fixed SNPs) and walker hound the breed with least (745 SNPs) compared to wolves. When comparing the breeds, walker hound compared to  Although only a reduced representation sample was used, individuals from the same species/breed group together, except two outliers (one Labrador (turquoise) and one walker hound (grey)).

Wolf-dog methylation differences
For the comparisons between wolves and dogs, only female individuals were included, since we only had brains from female wolves. The differentially methylated regions (DMRs) comparing wolf and dog are visualized in Manhattan and volcano plots in Fig 2A and 2B, respectively. In this reduced representation sample, there were 64 significant DMRs (p<0.0005) across 15 chromosomes between wolf and dog, located in, or close to, 11 unique genes ( Table 3). All of these were hypermethylated in the wolf compared to dogs. The position of the DMRs relative to their closest genes is shown in Fig 2C. As can be seen, most DMRs are found in introns. DNA methylation was also compared between wolf and each of the eight breeds separately (Fig 2D). In contrast to the DMRs in the wolf-dog comparison, in the breed comparisons, a majority of the significant DMRs were hypomethylated in the wolf (represented by a negative log fold change). Across all breeds, there were in total 116 significant DMRs in 17 genes (Table 3). Interestingly, few DMRs overlap across the different comparisons between wolves and dogs.
The gene list from the first comparison, wolf versus all dogs, was further analyzed in GO analysis and pathway analyses. 183 genes (DMRs p<0.05) were used as input, of which 40 were annotated to a specific GO function and 65 assigned specific pathway. The GO analysis identified an enrichment for, among other things, regulation of nervous system development, synapse structure, and cell development. GO terms related to anatomical structures were also found (p<0.05; Table 4). However, none reached adjusted significance (false discovery rate, FDR). The pathway analysis identified two pathways that reached adjusted significance (FDR<0.05): 1) TNF signaling and 2) SHC-related events triggered by IGF1R. Both of these are signal transduction pathways. TNF signaling pathway is involved in, for example, cell growth and death and immune and inflammatory responses [50], and IGF1R signaling promotes cell growth and differentiation. Interestingly, IGF1 is a determinant of small size in dogs [51]. Other pathways that we found were related to, among others, developmental biology (axon guidance), immune response and metabolism (p<0.05) ( Table 5). High level GO terms for genes from both wolf versus dog and wolf versus breed comparisons are presented in S2 Table. The genes are related to a wide array of processes, but noteworthy are genes involved in neurological processes, stress response, and reproduction. It is important to mention that due to the nature of these analyses, the pathways showing enrichment are biased towards genes that have been extensively investigated in the literature.

Methylation differences across dog groups
Each breed was also compared to each of the other breeds. Volcano plots for each comparison are presented in Fig 3A. The boxer was the most divergent and had the largest number of DMRs that were hypermethylated compared to the other breeds. German shepherd dog had the least DMRs, and Beagle, Great Dane, and pitbull terrier were hypomethylated ( Fig 3B). Together, the DMRs were located in 272 genes and a list of these is available in the S3 Table. High-level GO terms are presented in S4 Table. Noteworthy are genes involved in anatomical structures, behavior, neurotransmitter secretion, and stress response.
Additionally, female and male dogs were compared, and we found several differences. In the reduced representation sample, 75 DMRs differed between the sexes, situated across 17 chromosomes and located in, or close to, 21 unique genes (S5 Table). Females were more methylated than males at a majority of the DMRs (64%). The differences are visualized in a Manhattan plot (Fig 4A) and a volcano plot (Fig 4B). Most DMRs were found in introns followed by intergenic and upstream regions (Fig 4C). High level GO terms for groups of genes are presented in S6 Table. Although breeds were not balanced for sex, no overlaps of DMRs (p<0.005) were found between breed and sex analyses.

Discussion
The role of epigenetic factors in evolution, domestication, and selection is receiving increasing attention and it has been suggested that epigenetics may play a greater role than previously assumed. In the present study, we explored differences in DNA methylation in the brain of domesticated dogs and their ancestor species, the grey wolf, and between breeds of dogs, which reflects a more recent selection. By utilizing a reduced fraction of the canine genome, we found distinct DNA methylation profiles in the brains of the wolf and dog as well as for the different dog breeds, suggesting that epigenetics has been important in the divergent selection during dog domestication and breed formation. Our results, however, cannot discriminate between true methylation alterations and potential confounding factors such as copy number variations in the regions showing methylation differences [52].
DNA methylation differences between ancestral and domesticated populations have been shown in previous studies comparing DNA methylation in blood and buccal samples of wolves and dogs [28,29], and comparing domesticated chickens to their ancestor [18,19]. Methylation differences have also been demonstrated between populations of Red Junglefowl selected for either high or low tameness over only five generations [53]. Although our present results do not allow any conclusions about specific causative effects, they strongly suggest that DNAmethylation is affecting or being affected by selection during dog domestication.
Because species-specific DNA methylation patterns have been shown to be stable across generations [18,28], we can, based on our limited but very valuable biological material, compare species and patterns of DNA methylation related to the breeds. Our results show that the approach of combining genotype-by-sequencing and methylated DNA immunoprecipitation can be used in canines as a genotyping and epigenotyping method. GBS has been used previously to test domestication scenarios in chickens (Gallus gallus) [19] and in Lima beans (Phaseolus lunatus L.) [54], and the combined GBS-MeDIP has been used to study effects of different rearing conditions on chicken DNA methylation profiles in erythrocytes [37]. The reduced genome from the GBS procedure represented approx. 1% of the canine genome and combined with MeDIP it represented approx. 0.5%. Although only a small fraction of the genome was sequenced, the cladogram groups the individuals correctly according to group (wolf and breeds), similarly to previous whole-genome studies [22,23]. Hence, in spite of the fact that occasional recognitions sites for the enzyme may differ between individuals due to sequence differences, the reduced fractions are sufficiently similar across the populations in order to produce a correct cladogram.
The methylated regions of the female wolves were compared to female dogs in two ways: 1) to the complete set of dogs and 2) to each of the eight breeds. In both types of analyses, we found several DMRs between wolves and dogs. However, different patterns emerged from the different analyses. In the comparison with all dogs, DMRs were hypermethylated in the wolves  Table 3. List of genes with significant DMRs between wolf and dogs and wolf and breeds. Location is presented as chromosome:start-stop and adjacent regions have been merged. Strand is indicated by 1 or -1. Significant DMRs in intergenic regions are not included. P<0.0005 is considered significant. whereas the opposite was found when comparing the wolf to each of the breeds. Also, there were few overlaps in significant DMRs across the wolf-breed analyses. When we compare wolves to a pooled sample of different dog populations, we explore mainly epigenetic effects related to domestication, whereas when we compare wolves to each breed, we rather explore breed diversification. This implies that DNA methylation may be an important factor in the recent intense dog selection process, and this may perhaps be valid for speciation in general. This is further supported by our finding of substantial differences in methylation patterns between dog breeds. Previous studies have identified somatic DNA methylation differences as being important in speciation. Skinner, Gurerrero-Bosagna [55] found that epigenome differences between different species of Darwin's finches were a better match than DNA sequence variation with respect to the evolutionary relationship. In line with this, Smith, Martin [56] found that the epigenome across and within darter populations, a stream fish of genus Etheostomais sp., were more diverse and changed faster than the genome, and that behavioral isolation of a population increased with DNA methylation differences. In chicken domestication and breed formation, CpG-related SNPs are reported to increase with genetic distance, which suggests an important role for DNA methylation in speciation [19]. Selection for behavior has been proposed as being the main driver in domestication (Trut et [30] and the brain is the organ of interest for behavioral effects. Therefore, the tissue we have used to measure methylation differences is the brain. It is decidedly interesting that our results have highlighted genes and processes that may have been important for domestication and breed divergence. These include genes related to neurological processes, behavior, and synapse activity. Interestingly, Janowitz Koch, Clark [28] found genes related to the neurotransmitters GABA and glutamate to be differentially methylated between wolf and dog. This was, however, from blood samples. In the present study, using brain samples, we find several DMRs in genes related to neurotransmission between the domestic dog and their ancestor, giving further evidence that domestication has targeted neurological functions. For example, TEPSIN protein is part of the adaptor protein complex-4 that interacts with a glutamate receptor [57]; the expression of CAMKV is regulated by a glutamate receptor (AMPA) [58]; PLXNB2 regulates GABAergic and glutamatergic synapse development [59]; GRIN2B is an effector in a pathway activated by, among others, neurotransmitters, and that leads to neurite outgrowth [60]; JAK2 and ROCK are parts of pathways that regulate GABA expression [61]. We also identified differentially methylated genes between dog breeds that are interesting from a behavioral viewpoint, indicating that DNA methylation of these regions is important also for breed differences. For example, SLC17A5 (sialin) transports neurotransmitters glutamate and aspartate [62], and PTPRZ1, that has been associated to schizophrenia. Mice overexpressing this gene show hyperactive behavior and have an altered glutamatergic, GABAergic and dopaminergic activity [63]. It has been shown that DNA methylation can affect behavior. In rats, maternal care affects methylation of the glucocorticoid receptor gene and alters stress reactivity [8]; in great tits, explorative behavior is affected by methylation level at a dopamine receptor gene [64], and in dogs, DNA methylation in the promoter region of the oxytocin receptor gene affects dog's social behavior [10].

Comparison
Wolf and dog breeds do not only differ in behavior, but also in morphological traits. Therefore, it is also of interest that we have identified DMRs in genes related to morphology.  From the pathway analysis, an IGF1 signaling pathway was significant. The IGF1 gene is involved in several prenatal and postnatal processes and it regulates imprinted genes, many of which control metabolism and growth [65]. The gene is a known determinant of small size in dogs [51] and has also been associated with dog anxiety [66]. In the GO analysis, the most relevant was 'pattern specific process'. Hox genes, for example, are involved in forming the appendicular skeleton and mutant mice show dramatic phenotypic differences [67]. Interestingly, the dog breed pug shows typical vertebral anomalies consistent with altered hox genes [68]. Furthermore, we found DMRs in PER2 gene in the wolf-boxer comparison. In humans, DNA methylation in this gene has been associated to obesity [69].
Most DMRs were found in intronic regions. It has been reported that DNA methylation of introns regulates gene expression. Anastasiadi, Esteve-Codina [70] found, across species, that the methylation level of the first intron was inversely linked with expression levels. Specifically, it has been suggested that neurological-related processes may be regulated by intronic DNA methylation [71]. DNA methylation can also affect splicing outcomes [72], and mutations in introns may affect splicing with consequences for RNA processing [73]. Interestingly, mutations in splicing sites have been associated with, for instance, human neurological diseases [74,75]. The relevance of genetic and epigenetic changes of intronic and intergenic regions merits further investigation.
As expected, and in accordance with previous studies on other species [32,33,53], we also found sexually dimorphic DNA methylation profiles in dogs. DMRs are found both in autosomes and sex chromosomes. Female DMRs were more methylated than males which has previously been shown in, for example, human and rat brains [32]. There are behavioral differences between female and male dogs [e.g. [76][77][78], where epigenetic factors seem to be important. For example, the level of methylation of the oxytocin receptor gene has been reported to be linked to differences in the social behavior in male and female dogs [29]. Moreover, it has been shown that DNA methylation of specific regions are strongly associated with sex-specific gene expression [33], and environmental challenges may have sexually dimorphic effects on DNA methylation and, consequently, on gene expression and behavior [79].
Our study has a number of important limitations. Firstly, the samples were all obtained from ad hoc donated brains, so we can not be certain that they are truly representative of the different breeds. However, the phylogenetic analysis clearly showed that the samples clustered according to breed with respect to DNA-sequence variation, so we feel relatively confident in assuming that they should grossly represent the breeds. Secondly, the dogs were of different age, and it is well known that methylation patterns may change over life, as has been demonstrated in leucocytes [80]. However, the dynamic part of the methylome represents a relatively minor fraction of the methylation differences, so we believe that the majority of the DMR's detected here are permanent, so called obligatory epigenetic variation that are considerably more stable over time. It is also to be expected that methylation in brain neurons are more stable than in leucocytes. Furthermore, there were no systematic age differences between the dogs from the different breeds. This also pertains to the third important caveat, and that is the fact that some epigenetic variation may occur due to environmental impact. Again, this will cause dynamic, or facultative, stochastic methylation differences and should logically not affect the present results to any major degree.
In conclusion, we show that the epigenome of the brain, in the form of DNA methylation pattern, differs substantially between the wolf and the domestic dog, and between different dog breeds. This tentatively suggests that epigenetic factors may have been, and still are, important mechanisms involved in dog domestication and dog-breed formation. Specifically, we highlight methylation differences in genes related to behavior and morphology. Although the causal effects of methylation differences between and within canines remains to be studied, these may alter gene expression and, thereby, phenotype. We hypothesize that these differences are involved in the great phenotypic variation found among dogs and future studies should further explore the impact of epigenetic factors in relation to genetic factors. We also show that the new method of combining genotype-by-sequencing and methylated DNA immunoprecipitation can be used as an efficient epigenotyping method in canines. Our results are useful to further our understanding of epigenetic involvement in evolution and speciation, and they add an intriguing dimension to dog breeding.