Analysis of copy number variations in the sheep genome using 50K SNP BeadChip array

In recent years, genome-wide association studies (GWAS) have successfully uncovered single-nucleotide polymorphisms (SNPs) associated with complex diseases or traits [1]. With the rapid development of chip array-based genotyping techniques, thousands of genomic submicroscopic structural variations have been found in the human genome [2]. As a main genetic form of submicroscopic structural variation copy number variations (CNVs) are widely distributed in the human genome and influence gene expression, phenotypic variation and adaptation by disrupting genes and altering gene dosage [2–5]. Numerous studies showed that CNVs contributed to both disease susceptibility and phenotypic diversity [2, 5]. Now, CNV is increasingly considered to be an important and abundant source of genetic variation and phenotypic diversity [5, 6].

Investigations on CNVs have been successively carried out in human and other species [7–13]. In the domestic animals, there are involving in cattle [14–20], dog [21], chicken [22], pig [23, 24], goat [25], sheep [26] and rabbit [27]. As for sheep, Fontanesi et al. [26] provided a first comparative map of CNVs of the sheep genome referred to the cattle genome using a cross-species array comparative genome hybridization(aCGH). However, the cross-species analysis based on heterologous hybridization couldn’t identify all detectable CNVRs due to low homology between cattle probes and sheep DNA for some regions and doesn’t show the CNVR distributions on the sheep genome. In addition to CGH, another major platform commonly used to identify CNVs is the SNP array. In SNP array, intensity values of SNPs derived from each sample are used to detect CNVs in each individual. Comparison with two panels, CGH array has excellent performance in signal-to noise ratios, while the SNP array based approach is more convenient for high-throughput analyses and follow-up association studies [28]. With the development of high density SNP arrays, higher resolution of genomic regions can be achieved [29]. Moreover, due to their low cost and high-density, commercial SNP arrays have been widely used for CNV detection in domestic animals, and CNV mapping and functional studies have made important progress. However, there are no reports on CNV detection of sheep based on SNP array data.

In this study, we will investigate genome-wide CNV in three sheep populations. To pursue convincing results, we firstly employ the PennCNV program to analyze Ovine SNP50 genotyping data, and then use other algorithm program, cnvPartition, to validate CNVRs detected by PennCNV. To our knowledge, we will construct the first sheep CNV map based on SNP array data. This research will provide useful addition to the sheep CNV maps, and will provide potential genetic markers for further investigation on the roles of CNV in sheep productive traits and evolutionary adaptation.

In recent years, CNVs have been increasingly recognized as an important source of genetic variation and phenotypic diversity [4, 5, 30]. A number of CNVs associated with diseases have been found in human genetic researches, involving autoimmune disorders, schizophrenia, autism, as well as infectious and cardiovascular diseases [32, 33]. In addition, CNVs have been shown to play a role in pharmacokinetics in terms of drug efficacy and toxicity [34]. Therefore, the systematic analysis and characterization of CNVs improves our understanding of genetic variation and is an important tool for deciphering the role of CNV in the heritability of complex traits. In the past few years, with the development of high-density genotyping arrays, detection of CNVs using SNP genotyping arrays has become a cost-effective and efficient approach. Fine-scale CNV maps for human and other species have been constructed and refined [18, 35, 36].

In this study, we used whole genome genotyping based on Ovine SNP50 BeadChip arrays to identify CNVs in the sheep genome. After genotyping 329 samples using the Illumina platform from three sheep breeds, the signal intensity (LRR) and allelic intensity (BAF) ratio of all samples were exported using the Illumina GenomeStudio software (Illumina Inc., San Diego, CA, USA). In order to obtain a higher accuracy, quality control and CNV detection optimization was conducted using the PennCNV software when generating CNV calls [16]. To identify reliable CNV data, other studies have used different criteria for filtering unreliable CNV calls. Jakobsson et al. restricted analyses to CNVs with a minimum of 10 markers per CNV to minimize the number of false positives [37]. In our study using PennCNV analysis, we tested the impact of genomic waves on CNV calls and the filtering parameters with or without the –gcmodel option set. Compared with the results from the PennCNV using the –gcmodel option, more CNVs were identified without the –gcmodel option (5418 CNVs) than with the –gcmodel option (5008 CNVs; data not shown). This result was similar to that reported by Hou et al.[16]. These discordant calls were likely due to false positives called from the differentiating signal intensities caused by ‘genomic waves’. This further demonstrated that genomic waves have a significant effect on CNV analysis. We also used the ‘filter_cnv.pl’ program to perform sample-level quality control (see Materials and Methods). After this step, 73 further samples were excluded because their intensity data failed to conform to the criteria. Finally, the remaining 256 samples were used for CNV detection. This study focused only on the CNVs in 26 autosomes and excluded chrX and chrUn from our analysis because the chrX sequence and annotations are still not well-known [38], furthermore, their sequences and SNPs were uncertain, including the SNP mapping.

We identified a total of 238 CNVRs in three sheep breeds using the PennCNV detection algorithm. Of these, 219 CNVRs were losses, higher than the proportion of gains (Figure 2). This result might be due to the fact that losses are easier to detect than gains, because the exponential intensity data is linearly correlated with the copy number [39]. To exclude the potential false positive signals, we used another CNV detection software, CNVpartition, which detected 52 CNVRs. Forty-seven regions were consistently identified by both CNV detection software packages. Compared to CNVpartition based on other algorithms, such as QuantiSNP [40], Birdsuite [41], and GADA [42], PennCNV combines multiple sources of information, including the total signal intensity and allelic intensity ratios at each SNP marker to generate a hidden state for copy neutral loss of heterozygosity, the distance between neighboring SNPs, and the allele frequency of SNPs. PennCNV also integrates a computational approach by fitting regression models with GC content to overcome ‘genomic waves’ [43].

We compared our CNV length dataset with another sheep CNV study that was based on a cross species cattle-sheep aCGH experiment using a tiling oligonucleotide array with approximately 385,000 probes designed using the bovine genome. Because our dataset excluded CNV calls in the chrX and chrUn, only autosomal CNVR calls (238 CNVRs) were compared to the autosomal CNV calls of that study (135 CNVRs). As in other CNV studies using aCGH experiments [15, 20, 25, 26], the number of loss events in our dataset was larger than the number of gain events. However, further validation is required as it is also possible that purifying selection occurred [5]. The CNV sizes in our study ranged from 100 to 500 kb (average, 144.6 kb), which differed from those in the work of Fontanesi et al. (average, 73.9 kb) [26]. This CNV size difference was likely due to sampling differences or to resolution and genome coverage differences between the two techniques. The SNP genotyping resolution was 50.9 kb (mean probe spacing), whereas that of the aCGH platform was 6.3 kb. A SNP array lacked non-polymorphic probes designed specifically for CNV identification. Thus, only the large CNVRs were identified with the Ovine SNP50 assay. In future experiments, high-density SNP arrays combined with improved CNV calling algorithms could remedy these differences. In addition, we investigated the distribution pattern of these CNVRs in 26 autosomes and determined whether the CNV region was common in three sheep breeds. Interestingly, 75 of the 238 CNVRs had frequency rates of > 3%, whereas 121 of the 238 CNVRs had a frequency rate of < 1% in 256 individuals. This result might have been influenced by the innate limitation of Ovine SNP 50 arrays, and therefore, many common CNVs could have been missed.

We used UCSC gene annotation (http://genome.ucsc.edu) to identify genes that were located within or partially overlapped with CNVRs. Initially, we used a database of known sheep genes to annotate the gene content of CNVRs. However, the annotation results showed only four CNVRs overlapping with five known RefSeq genes in sheep CNVRs (see Additional file 5: Table S6). These results are likely attributable to the fact that the sheep genome is not well-annotated compared to the human genome or that of other domesticated animals [44–48]. The current sheep genome version (ISGC Ovis aries 3.1) is incomplete and known related genes are rarely identified in the sheep RefSeq gene database [38]. Therefore, we performed a BLASTN search for homologous human and cattle sequences using the UCSC table browser tool, and identified 563 genes homologous to human IDs identified in DAVID. We conducted a GO analysis to determine the biological effects of the 563 copy number variant genes (homologous to human genes). Similar to other results for GO analyses, the enriched genes were related to those involved in olfactory receptor (OR) activity, sensory perception of smell, sensory perception of chemical stimulus, sensory perception, cognition, G-protein coupled receptor protein signaling pathway, neurological system processes, cell surface receptor-linked signal transduction, plasma membrane and integral membrane components. Table 4 showed that the genes involved in environmental responses were over-represented in the CNVRs. In this study, functional analysis results were similar to those previously reported in CNV studies in human and other mammals [11, 15, 19, 26, 49]. This significant (p < 3.51 × 10^-158) enrichment for OR activity might represent the frequent occurrence of CNVs in gene clusters for OR [50], as previously observed in cattle [19]. In addition, we observed 563 human orthologous genes in sheep CNVRs, of which 341 genes were previously reported in the DGV database (Additional file 6: Table S7). Moreover, we identified CNV regions that may overlap with Online Mendelian Inheritance in Man (OMIM) genes influencing disease susceptibility (Table 4). Notwithstanding, in the great majority of cases, these CNVs encompass genes and thereby directly affect gene dosage by loss or gain in the level of gene expression [2, 3, 51]. However, some other reports have suggested that CNVs are located preferably in gene-poor regions [52, 53] or in noncoding regions, and some studies provide evidence that the CNVs are considered nonpathogenic. Recent studies have described genome-wide distribution of CNVs in regions that encompass noncoding sequences, thereby affecting the regulation of distant target genes [54–57]. This suggested impact of CNVs in noncoding regions requires further elucidation in future studies. As discussed above, functional analysis studies indicated that the CNV genes possess a wide spectrum of molecular functions. However, since the sheep genome is less well-defined, Sheep CNVRs warrant further investigation for their roles in complex traits.

To verify CNV obtained by PennCNV, we conducted quantitative PCR (qPCR) on ten selected CNVRs and compared them with a sample control region known to have no CNVs (the DGAT1 gene fragment). We found that 26.7% of our qPCR assays (16 confirmed/60 total) agree with our CNVRs predictions in these regions. In total, six regions of these CNVRs were validated in at least one qPCR assay. However, it should be noted that another four CNVRs were not confirmed by qPCR (Additional file 4: Table S5). Reasons for this discrepancy between the CNVR analysis of Ovine SNP50 BeadChip data and the qPCR experiment may be due to the fact that the low probe density of the Ovine SNP50 BeadChip made it difficult to identify the true boundaries of CNVRs, resulting in an overestimation of their actual sizes. In addition, the primer pairs might have been designed outside the copy number polymorphic region.

In this research, only a small proportion of CNVRs identified overlapped with other studies. A similar situation was also encountered in human and other mammal CNV studies [10, 16, 58]. This indicated that a vast amount of CNVs existing had not been detected in the Ovine genome. Furthermore, a comparative analysis of CNVRs identified in sheep, goat and cattle could reveal the evolutionary mechanisms determining CNV formation during the mammalian evolution.

We firstly constructed a sheep CNV map based on the Ovine SNP50 array. 238 CNVRs were totally identified in the sheep autosomal genome, these CNVRs were non-randomly distributed on chromosomes, and 70% of which (166/238 CNVRs) were shared in > two animals. We also confirmed 6 CNVRs in total of 10 selected CNVRs by qPCR method. Compared with the results of CNV studies based on aCGH, the currently available genome wide SNP assays can infer Ovine CNV efficiently. However, it should be noticed that smaller size CNVs (< 10 kb) were not identified using this SNP panel and only 16 out of the 60 qPCR assays confirmed in this study was likely to be a overestimation of the true number and true boundary of CNVRs in the sheep genome. On the other hand, in this study, we using two detection programs: PennCNV and CNVpartition, to identify the sheep genome CNVs, results indicated that 47 CNVRs identified by the PennCNV (19.7%) overlapped with the CNVpartition data (90.4%), which highlighted the differences and commonalities of the two detection methods. Follow-up studies for high-resolution CNV mapping require improved SNP assay and next-generation sequencing to improve the accuracy of CNV calling. In addition, integration of different algorithms can enhance the detection of genomic structure variations. Furthermore, our findings would be of help for understanding the ovine genome and provide preliminary foundation for carrying out the CNVs association studies with economically important phenotypes of sheep in the future.

The experiments were performed on trait records and DNA samples that had been collected previously from animals born in 2010. Because no new animals were handled in these experiments, there are no needs of Animal Care and Use Committee approval for this study.