Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Genome-wide patterns of copy number variation in the diversified chicken genomes using next-generation sequencing

Genome-wide patterns of copy number variation in the diversified chicken genomes using... Background: Copy number variation (CNV) is important and widespread in the genome, and is a major cause of disease and phenotypic diversity. Herein, we performed a genome-wide CNV analysis in 12 diversified chicken genomes based on whole genome sequencing. Results: A total of 8,840 CNV regions (CNVRs) covering 98.2 Mb and representing 9.4% of the chicken genome were identified, ranging in size from 1.1 to 268.8 kb with an average of 11.1 kb. Sequencing-based predictions were confirmed at a high validation rate by two independent approaches, including array comparative genomic hybridization (aCGH) and quantitative PCR (qPCR). The Pearson’s correlation coefficients between sequencing and aCGH results ranged from 0.435 to 0.755, and qPCR experiments revealed a positive validation rate of 91.71% and a false negative rate of 22.43%. In total, 2,214 (25.0%) predicted CNVRs span 2,216 (36.4%) RefSeq genes associated with specific biological functions. Besides two previously reported copy number variable genes EDN3 and PRLR,we also found some promising genes with potential in phenotypic variation. Two genes, FZD6 and LIMS1, related to disease susceptibility/resistance are covered by CNVRs. The highly duplicated SOCS2 may lead to higher bone mineral density. Entire or partial duplication of some genes like POPDC3 may have great economic importance in poultry breeding. Conclusions: Our results based on extensive genetic diversity provide a more refined chicken CNV map and genome-wide gene copy number estimates, and warrant future CNV association studies for important traits in chickens. Keywords: Copy number variation, Whole genome sequencing, aCGH, Genetic diversity, Chicken Background sequences (insertions) of ≥1 kb in length and ≥90% se- Copy number variations (CNVs) are defined as gains or quence identity are also suggested to be one of the major losses of DNA fragments of 50 bp or longer in length in catalysts and hotspots for CNV formation [11,12], mainly comparison with reference genome [1,2]. CNVs contrib- because the genomic regions flanked by SDs are suscep- ute significantly to both disease susceptibility/resistance tible to recurrent rearrangement by NAHR [11,13]. In and normal phenotypic variability in humans [3-5] and terms of total bases involved, the percentage of the genome animals [6-9]. Four major mechanisms have been found affected by CNVs is higher than that of single nucleotide to be related to CNV formation including non-allelic hom- polymorphism (SNP) markers. Although SNPs are gener- ologous recombination (NAHR), non-homologous end ally considered as more suitable markers in the genome- joining (NHEJ), Fork Stalling and Template Switching wide association studies (GWASs), most reported SNP (FoSTeS) and LINE1 Retrotransposition [4,10]. Addition- variants have relatively limited effects and explain only a ally, segmental duplications (SDs) which are duplicated small proportion of phenotypic variance [14]. Further, CNVs encompassing genes or regulatory elements are be- * Correspondence: [email protected] lieved to exert potentially larger effects on gene expression Equal contributors through changing gene structure and dosage, altering gene Department of Animal Genetics and Breeding, College of Animal Science regulation, exposing recessive alleles and other mechanisms and Technology, China Agricultural University, Beijing, China © 2014 Yi et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Yi et al. BMC Genomics 2014, 15:962 Page 2 of 16 http://www.biomedcentral.com/1471-2164/15/962 [1,4,15,16]. CNVs are also found to alone capture 18 to 30% 12 chickens from multiple breeds with extensive genetic of the total detected genetic variation in gene expression in diversity, including seven Chinese indigenous breeds humans and animals, and may contribute to a fraction of [37], four commercial breeds and one Red Jungle Fowl. the missing heritability [17,18].Therefore,identificationof Then we applied NGS-based method to construct a CNVs is essential in whole genome fine-mapping of CNVs more refined and individualized chicken CNV map, in- and association studies for important phenotypes. vestigate genome-wide CNV characteristics and estimate Originally, two cost-effective and high-throughput genome-wide gene copy number. The results will enable methods including array comparative genomic hybridization us to better understand the patterns of CNVs in the (aCGH) and commercial SNP microarrays are used for chicken genome and future CNV association studies. CNV screening [19,20]. However, due to the limitation in resolution and sensitivity, it is difficult with the two ap- Results proaches to detect small CNVs shorter than 10 kb in Mapping statistics and CNV detection length and identify the precise boundaries of CNVs We performed whole genome sequencing in 12 different [21,22]. The two analytical platforms also reveal inconsist- breeds of female chickens using Illumina paired-end li- ent results with poor overlaps owing to different designs braries and obtained a total of 12.9 Gb high quality se- and probe densities and coordinates [18,20]. Furthermore, quence data per individual after quality filtering. After the presence of SD regions is also a common challenge for sequence alignment and removing potential PCR dupli- microarrays, because a considerable proportion of CNVs cates, the sequencing depth varied from 8.2× (CS) to fall into SD regions not well-covered by microarrays 12.4× (WR), which was sufficient for CNV detection, and [2,23]. Recently, a variety of CNV detection approaches the average coverage with respect to the chicken reference based on next-generation sequencing (NGS) were pro- genome sequence was 97.2% (Table 1). We calculated the posed and offer promising alternatives as they have higher average read depth (RD) of 5 kb non-overlapping windows effective resolution to discover CNVs with more types and for all autosomes and performed GC correction. The GC- wider size ranges [24]. One leading method is read depth adjusted RD mean and standard deviation (STDEV) for (RD) (also known as depth of coverage (DOC)) with the each individual are listed in Table 1. We applied the pro- capability of inferring gains or losses of DNA segments gram CNVnator to 12 individuals and the average number and determining absolute copy number values, which de- of CNVs per individual was 1,328, ranging from 644 in tects CNVs by analyzing the number of reads that fall into WL to 1,921 CNVs in BY. A detailed description of CNV each pre-specified window with a certain size [25,26]. calls can be found in Additional file 1: Table S1. For all the Hence, the advent of NGS technologies and suitable ana- autosomal CNVs classified as duplications, the average lytical methods can promote to systematically identify copy number value of all the individuals was 3.78 and the CNVs at higher resolution and sensitivity. maximum copy number estimate was 40.8 on chromo- At present, the three aforementioned high-throughput some 2 (chr2) in RJF. platforms have been applied to livestock genomics for A total of 8,840 CNV regions (CNVRs) allowing for CNV detection, such as sheep [27], horse [28] and cattle CNV overlaps of 1 bp or greater were obtained, covering [2], and uncover several CNVs associated with important chromosomes 1–28, two linkage groups and sex chromo- phenotypes. Some CNVs are also found to be the genetic somes, which amounted to 98.2 Mb of the chicken gen- foundation of phenotypic variation in chickens. A dupli- ome and corresponded to 9.4% of the genome sequence cated sequence close to the first intron of SOX5 is asso- (Additional file 1: Table S1). The individualized chicken ciated with the chicken pea-comb phenotype [29] and CNV map across the genome is shown in Additional file 2: an inverted duplication containing EDN3 causes dermal Figure S1. The length of CNVRs ranged from 1.1 to hyperpigmentation [30]. Partial duplication of the PRLR 268.8 kb with an average of 11.1 kb and a median of 6.6 kb. is related to the late feathering [31]. In total, 6,137 (69.4%) out of all CNVRs had sizes varying A genome-wide chicken CNV analysis is desired since from 1.1 to 10 kb (Figure 1A). Although chr1 had a max- chicken is not only an economically important farm ani- imum of 1,928 CNVRs, the two largest CNVR density mal but also a valuable biomedical model [9,32]. How- values, defined as the average distance between CNVRs, ever, some previous CNV studies in chickens based on were 35.7 kb and 32.0 kb on the chr16 and chrLGE64, re- aCGH and SNP platforms mainly suffered from low spectively (Additional file 3: Table S2). The number of resolution and sensitivity [9,32-35]. A latest study exhib- CNVRs in different individuals varied greatly, ranging from ited the detection of four main types of genetic variation 629 in WL to 1,890 in BY. Among all CNVRs, 6,083 from whole genome sequencing data using two chickens (68.8%) were present in a single individual, 1,423 (16.1%) [36], suggesting the efficiency of CNV detection via deep shared in two individuals, and 1,334 (15.1%) shared in at sequencing. Considering that a great number of CNVs least three individuals (Figure 1B). Further, the mean and appears to be segregating in distinct breeds, we selected median lengths of the unique CNVRs were 8.9 kb and Yi et al. BMC Genomics 2014, 15:962 Page 3 of 16 http://www.biomedcentral.com/1471-2164/15/962 Table 1 Summary statistics for sequencing and CNVs of 12 individuals Chicken Numbers of Depth Coverage (%) Autosome reads Autosome Duplications Deletions Sequence a b abbreviation mapped reads per 5 kb window reads STDEV covered (Mb) BY 102,002,937 9.7 97.0 489.29 110.73 1,319 602 34.1 CS 85,383,494 8.2 96.9 409.93 101.42 1,132 663 26.8 DX 129,847,015 12.4 97.4 623.50 130.46 552 820 8.2 LX 105,152,881 10.0 97.3 503.82 112.74 898 821 11.7 RIR 102,464,756 9.8 97.3 490.96 108.21 578 669 8.3 RJF 105,517,587 10.1 97.2 504.23 113.52 702 620 9.8 SG 85,987,827 8.2 96.6 412.27 87.66 470 553 7.2 SK 95,322,371 9.1 97.1 457.21 100.61 773 657 12.3 TB 107,535,104 10.3 97.3 515.68 108.07 607 679 8.5 WC 119,116,969 11.4 97.4 572.35 121.39 710 768 9.8 WL 118,689,980 11.3 97.5 567.18 118.63 203 441 3.3 WR 130,307,416 12.4 97.6 625.01 132.32 224 477 3.3 BY Beijing You, CS Cornish, DX Dongxiang, LX Luxi Game, RIR Red Island Rhode, RJF Red Jungle Fowl, SG Shouguang, SK Silkie, TB Tibetan, WC Wenchang, WL White Leghorn, WR White Plymouth Rock. The number of reads per 5 kb windows after GC correction. 5.8 kb, whereas the shared CNVRs sizes were 15.9 kb on In total, 7,530 out of 8,487 (88.7%) autosomal CNVRs were average and 9.5 kb as the median. According to the type of converted successfully. The detailed comparison results are CNVRs, they were divided into three categories, including presented in Table 2 and Additional file 4: Table S3. In our 4,761 gain, 3,773 loss and 306 both (gain and loss) events. results, 1,052 (14.0%) CNVRs with the total length of Gain events possessed longer genomic sequences than 19.7 Mb were reported by eight previous studies, and the losses both on average (14.2 kb vs. 5.4 kb) and in total remaining 6,478 (86.0%) were regarded as novel CNVRs. It (67.6 Mb vs. 20.3 Mb). In addition, the count of CNVRs on should be noted that the novel CNVRs had slightly smaller each chromosome was directly proportional to the chromo- sizes (10.6 kb) on average than those reported CNVRs some length, and five macrochromosomes (chr1-5) pos- (18.8 kb). As a special and important chromosome in the sessed a large proportion (61.8%) of all putative CNVRs. chicken genome, chr16 encompassed some CNVRs which could be confirmed by different platforms. Comparison with previous chicken CNV studies Considering that most of the previous studies excluded CNV quality assessment by CNVnator, aCGH and qPCR the CNVRs on sex chromosomes and unassigned linkage The copy number values of diploid regions on auto- groups, we migrated our autosomal CNVR coordinates somes theoretically equal to two, so we could inspect the from galGal4 to galGal3 using the UCSC liftOver tool. potential of CNVnator to generate false positive results by Figure 1 The length and frequency distribution of CNVRs. (A) 6,137 (69.4%) CNVR events are shorter than 10 kb, and the number of CNVRs longer than 50 kb is only 291 (3.3%). (B) 6,083 (68.8%) CNVR occur in only one individual and 2,757 (31.2%) CNVRs are shared in at least two individuals. Yi et al. BMC Genomics 2014, 15:962 Page 4 of 16 http://www.biomedcentral.com/1471-2164/15/962 Table 2 Comparison between autosomal CNVRs identified in this study and other chicken studies Platforms Results from different studies Overlapped CNVRs in this study Study Breed Samples Number Total length (Mb) Number Pct. of number (%) Total length (Mb) Pct. of length (%) Sequencing-based study This study 12 12 7,530 88.12 CGH-based studies Wang et al. 2010 [40] 3 10 91 15.72 162 2.15 2.66 3.02 Wang et al. 2012 [9] 3 6 130 3.34 83 1.10 0.92 1.04 Crooijmans et al. 2013 [34] 7 64 1,504 57.44 721 9.58 8.30 9.42 Luo et al. 2013 [39] 4 6 29 1.46 21 0.28 0.42 0.47 Tian et al. 2013 [33] 11 22 308 10.81 166 2.20 2.00 2.27 Abernathy et al. 2014 [41] 2 12 147 4.18 68 0.90 0.63 0.71 SNP-based study (60 K) Jia et al. 2013 [32] 2 746 209 13.55 141 1.87 1.75 1.99 Sequencing-based study Fan et al. 2013 [36] 2 2 415 3.17 96 1.27 0.80 0.90 Total 1,052 13.97 12.45 14.13 Yi et al. BMC Genomics 2014, 15:962 Page 5 of 16 http://www.biomedcentral.com/1471-2164/15/962 evaluating these two copies regions. For all 12 individuals, average of 0.647. BY (0.502), SK (0.435) and WR (0.491) we selected all 5 kb non-overlapping windows on auto- showed lower correlation close to 0.500, and we found the somes and excluded the windows intersecting with pre- mean of all probes log ratio values in the three aforemen- dicted CNVs and gaps, and then estimated their average tioned individuals were 1.05, 0.85 and 1.05 respectively, CN. The average CN and STDEV per individual was which were larger than the values of others that were close 2.077 ± 0.291, varied from 2.041 ± 0.226 in WR to 2.104 ± to zero. 0.299 in RJF, showing low variability within the predicted In addition, we randomly chose 15 predicted CNVRs neutral regions. Further, we validated sequencing-based representing different types and frequencies for qPCR as- CNV predictions by two independent experiment ap- says, and tested all 12 samples for each CNVR. Two dis- proaches including aCGH and qPCR. We performed 11 tinct pairs of primers were designed for each predicted pairwise aCGH experiments using RJF as the reference CNVR (Additional file 6: Table S4). The proportion of and all others as the test samples. Considering that we es- confirmed positive samples (positive predictive value) var- timated the CN of selected individuals with respect to ref- ied from 50 to 100%, with an average of 91.71%. However, erence genome which cannot be used for the aCGH some negative samples were also confirmed to contain reference sample, we calculated the predicted log CN ra- CNVs, and the false negative rate varied from 0 to 60%, tios for the 11 test samples against RJF to make the CN with an average of 22.43%. We illustrated the qPCR results values comparable with the aCGH results, which was des- for three confirmed CNVRs of different types (gain, loss ignated as digital aCGH approach [12,38]. We focused on and both) (Additional file 7: Figure S3). the autosomal CNVs to avoid the impact of gender and unassigned linkage groups. For pairwise samples (each of Copy number polymorphic genes the 11 test samples and RJF), there were two types of We obtained 6,086 non-redundant RefSeq gene tran- CNV events, i.e., overlapping and unique segments. For scripts retrieved from the UCSC Genome Browser and the overlapping CNV segments, we first split them into estimated the copy number values of all genes in differ- non-overlapping subsegments. And then we estimated the ent individuals by CNVnator. A total of 2,216 (36.4%) CN of each subsegment and unique segment longer than genes overlapped with 2,214 (25.0%) predicted CNVRs. 1,000 bp for each of the two pairwise samples, and divided Among them, 537 genes were found to be completely the copy number estimates of the test sample by that of covered by CNVRs. The overlapping genes were found RJF and calculated log CN ratios as digital aCGH values. not to be highly duplicated sequences, and the max- Then we compared the digital values with aCGH probe imum copy number estimates was only 12.0. We exam- log ratios which were defined as the average of all probes ined the 25 most variable genes according to the STDEV log ratio values in the corresponding segments. We per- of their copy number estimates in different individuals formed a simple linear regression analysis to assess the (Additional file 8: Table S5), and found that these genes correlation between the two values. The Pearson’scorrel- were mainly involved in immune response and keratin ation coefficient (r) ranged from 0.435 in SK to 0.755 in formation. It should be noted that the keratin gene fam- DX (Figure 2 and Additional file 5: Figure S2), with an ilies were detected to have large CN values and variance. Figure 2 Correlation between digital aCGH and whole genome aCGH among Luxi Game and White Leghorn compared with Red Jungle Fowl (RJF). RJF is selected as the reference sample in each aCGH experiment. Digital aCGH values are estimated using calculated log CN ratios in which CN are estimated for identified CNV segments of two individuals and divided by the corresponding CN of RJF. Whole genome aCGH values are defined as the average of all probes log ratio values in the same segments as the digital aCGH. 2 Yi et al. BMC Genomics 2014, 15:962 Page 6 of 16 http://www.biomedcentral.com/1471-2164/15/962 Two significant CNVRs associated with dermal hyperpig- to perform qPCR experiments using the same two pairs of mentation were located on chr20 at positions 11,217,001 primers listed in Additional file 6: Table S4. Two qPCR re- to 11,272,200 (CNVR7962) and 11,651,801 to 11,822,900 sults demonstrated the copy number estimate of almost (CNVR7968), respectively, which had already been de- each LX was larger than the others (Figure 4), and the scribed in detail in a previous study [30], and the distance average copy number estimates (5.0 and 5.2 for two pairs between the two loci was 379.6 kb. SLMO2 and TUBB1 of primers, respectively) of all LX were significant larger were completely covered by the first region which was than those (2.6 and 2.6) in other individuals using the predicted to be about twice as many copies in DX and SK two-sample t-test (P-value =0.003 and 0.001). Additionally, as in other individuals (Figure 3A and Additional file 9: other identified CNV-gene overlaps could be potentially Figure S4a). The functional gene EDN3 (endothelin 3) is responsible for certain economic traits, as these genes were not archived because the predicted gene is not available involved in lipid metabolism, muscle development and pro- for UCSC RefSeq database. We found that only BY had tein secretion process. For example, our results suggested this CNVR while SK and DX as two typical breeds with higher copy number for the POPDC3 gene in WL (n =4.2) dermal hyperpigmentation did not. So we further checked than in the other 11 genomes (n =2.3) (Figure 3C and the raw results before removing CNVs overlapping with Additional file 9: Figure S4c). Similarly, the WL genome gaps. Two nearly identical CNVs comprising two gaps showed the greatest number of AVR2 copies (n =2.0) on (>100 bp) were found, one at positions 11,111,501 to chrZ compared with the others (n =1.1). Two promising 11,238,600 in DX and the other at positions 11,111,401 genes involving in lipid metabolism, AP2M1 and LBFABP, to 11,238,900 in SK, which were also confirmed by our were identified as the largest copy number (n =3.0 and whole genome aCGH experiments (Figure 3A and 3.2) in meat-type chicken (CS) compared with those in the Additional file 9: Figure S4a). The distance between the others. raw CNVR and the second region (CNVR7968) was 412.9 kb, which perfectly supported the reported results Heatmap visualization [30]. Conversely, the first CNVR in BY (11,217,001 to We performed a hierarchical clustering analysis and gen- 11,272,200) showing normal skin color does not contain erated heatmaps based on Pearson’s correlation coeffi- the EDN3 gene (11,148,025 to 11,160,484), which also pro- cient using the CN values of selected gene loci, in order vides evidence that the EDN3 with copy number poly- to infer the potential relationship of selected genes morphism is the causal mutation resulting in dermal among 12 individuals. The loci encompassing SLMO2 hyperpigmentation. Another previously identified CNVR and TUBB1 in DX and SK were found to be highly du- involving the PRLR (prolactin receptor) gene on chrZ [31] plicated regions and the two individuals were clustered was also detected in our study, and the copy number es- into one group (Figure 5A). Another promising gene, timates of PRLR in WC and WL were twice as many as SOCS2, was also confirmed for the difference in copy in other individuals. The sex-linked K allele containing number between LX and the others (Figure 5B). Mean- two copies of PRLR in females is associated with the late while, WL showed specific expansion in the POPDC3 feathering and used widely for sexing hatchlings. Our locus and was split into a separate clade (Figure 5C). sequencing-based and qPCR results showed that WC and WL should exhibit the late feathering phenotype, Gene content and QTL analysis of CNVRs which were supported by the actual phenotype record. A total of 2,216 RefSeq genes overlapped with putative In addition, we found that some genes related to the CNVRs. Then, we performed gene ontology (GO) and host immune and inflammatory response had CNVR Kyoto Encyclopedia of Genes and Genomes (KEGG) path- overlaps, like CD8A, FZD6, LIMS1, TNFSF13B and some way analysis for these genes. The GO analysis revealed MHC-related genes associated with Marek’s disease 646 GO terms, of which 175 were statistically significant (MD). The SOCS2 involving in the regulation of bone after Benjamini correction (Additional file 10: Table S6). growth and density was predicted to have the largest CN The significant GO terms were mainly involved in positive value in LX (n =6.4), while DX (n =3.0) and TB (n =3.6) regulation of macromolecule metabolic process and gene also possessed the duplicated sequences in this locus expression, plasma membrane, protein localization, en- compared with the neutral regions in other individuals zyme binding, response to oxidative stress and immune (Figure 3B and Additional file 9: Figure S4b). LX repre- system development. The KEGG pathway analysis indi- sents a characteristic breed for cockfighting in which bone cated that these genes were overrepresented in nine path- strength is an essential feature for selection. To validate ways, but none of which was significant after Benjamini the highly duplicated sequence (CNVR410) found only in correction. According to our artificial QTL filtering cri- LX, we selected another 16 individuals, i.e., eight LX (four teria, we identified 595 high-confidence QTLs in total, of males and four females) and other eight females consisting which 300 (50.4%) were found to overlap with 560 (6.3%) of one CS, one DX, one SG, one SK, two TB and two WL, CNVRs (Additional file 11: Table S7). These QTLs were Yi et al. BMC Genomics 2014, 15:962 Page 7 of 16 http://www.biomedcentral.com/1471-2164/15/962 Figure 3 (See legend on next page.) Yi et al. BMC Genomics 2014, 15:962 Page 8 of 16 http://www.biomedcentral.com/1471-2164/15/962 (See figure on previous page.) Figure 3 Visual examination by read depth, whole-genome aCGH and digital aCGH around three loci for five representative chicken genomes. The uppermost gene image is generated with the UCSC Genome Browser (http://genome.ucsc.edu/) using the galGal4 assembly. The track below the gene region is depth of coverage for all five individual genomes. Red indicates regions of excess read depth (> mean +3 × STDEV), whereas gray indicates intermediate read depth (mean +2 × STDEV < × < mean +3 × STDEV), and green indicates normal read depth (mean ± 2 × STDEV). All read depth values based on 1 kb non-overlapping windows are corrected by GC content. Whole-genome aCGH and digital aCGH values are depicted as the red-green histograms and correspond to a gain colored in green (>0.5), a loss colored in red (<−0.5) and normal status colored in gray (−0.5 < × <0.5). (A) Two previously reported CNVs (chr20: 11,111,401-11,238,900 and chr20: 11,651,801-11,822,900) associated with dermal hyperpigmentation. The DX and SK genomes show two additional copies of the two regions compared with RJF, and are also validated by whole-genome aCGH. (B) A higher copy number increase for the SOCS2 locus (chr1: 44,764,280-44,765,955) is predicted in LX than in other individuals. (C) The POPDC3 gene (chr3: 68,255,196-68,259,535) is predicted to be duplicated status only in WL. mainly involved in production and health traits, such as RD method employed in our work has advantages in both growth, body weight, abdominal fat weight, egg number technology platform and genetic diversity compared with and Marek’s disease-related traits. the eight previous reports [9,32-34,36,39-41]. Because a significant fraction of CNVs falls into genomic regions not Discussion well-covered by microarrays, especially for SD regions This study performed genome-wide CNV detection, deter- lacking sufficient probes [16,23], CNV as a major source mined absolute copy number and constructed the first in- of genetic variation is complementary to SNP and could dividualized chicken CNV map. The NGS technology and account for a substantial part of missing heritability [14]. Figure 4 Validation of CNVR410 by qPCR in another 16 chickens. X-axis represents all 16 samples and Y-axis represents normalized ratios (NR) estimated by qPCR. NR around 2 indicates normal status (2 copies), NR around 0 or 1 indicates loss status (0 copies or 1 copy), and NR around 3 or more indicates gain status (3 or more copies). (A) qPCR results confirmed by primer410_2. (B) qPCR results confirmed by primer410_3. Yi et al. BMC Genomics 2014, 15:962 Page 9 of 16 http://www.biomedcentral.com/1471-2164/15/962 Figure 5 (See legend on next page.) Yi et al. BMC Genomics 2014, 15:962 Page 10 of 16 http://www.biomedcentral.com/1471-2164/15/962 (See figure on previous page.) Figure 5 Hierarchical clustered heatmaps of preselected genetic loci for 12 chicken genomes. Every block in the heatmap indicates estimated CN values of 1 kb non-overlapping windows in the preselected region. These heatmaps are generated from hierarchical cluster analysis using Pearson’s correlation coefficient of the CN values. The colors for each bar denote different copy number (CN). (A) DX and SK which are predicted to be doubled within dermal hyperpigmentation loci are clustered together. (B) Upstream and downstream of the SOCS2 locus reveal higher CN values in DX, TB and WC especially LX. (C) WL shows specific expansion in the POPDC3 locus and is split into a separate clade. Most CNV studies to date have been discovery studies ra- larger than losses because chromosomal deletion can ther than association studies, mainly due to the limitations lead to a variety of serious malformations and disorders of CNV resolution and genotyping in each individual [3]. and is subjected to purifying selection [13,47]. In gen- The high-resolution individualized chicken CNV map eral, the length of chromosome is positively correlated based on extensive genetic diversity not only enriches the with the number of CNVRs. The chr16 (a microchromo- current genetic variation database but also encourages the some) was found to have the second densest CNVRs, future development of assays for accurately genotyping possibly owing to the highly variable major histocom- CNVs, enabling systematic exploration about CNV associ- patibility complex (MHC) regions and higher recombin- ation studies similar to SNPs. In future, integration of ation rate [48], which also results in repeatedly finding CNVs with SNPs maybeaneffectiveand promisingway to the same CNVRs on chr16 among different studies. elucidate the causes of complex diseases and traits [16,17]. Quality assessment and experimental validation Genome-wide CNV landscape in the chicken genome It is generally believed that the CN of neutral regions is The number of CNVs and CNVRs in each individual var- between 1.5 and 2.5 [25] and the mean ± 2 × STDEV in ied greatly, and all individuals shared a small number of our results corresponded closely to the hypothesis, them, suggesting that a considerable proportion of CNVs which demonstrates that CNVnator has efficient per- likely segregated among 12 distant breeds [2,34], therefore formance on CNV detection and CN estimation and can a larger population and multiple breeds are crucial to con- generate most reliable results. For CNV quality assess- struct a more complete chicken CNV map. The high per- ment by aCGH, the positive correlation values between centage of unique CNVRs could also be partly explained computational and experimental log CN ratios in our by the high recombination rate in the chicken genome study were higher than the previous results [2], mainly (2.5-21 cM/Mb), because recombination-based mecha- owing to the aCGH platform with higher resolution in nisms such as non-allelic homologous recombination our analysis. The slightly low correlation coefficients in (NAHR) are the major causes leading to CNVs [42]. Simi- BY, SK and WR might disclose certain experimental larly, the high recombination rate may induce more denser noises and biases resulting in misgenotyping in corre- CNVRs in microchromosomes [43]. These unique CNVRs sponding aCGH experiments [16], and particularly may be recent events in evolution and contribute to highly duplicated regions and rare deletions [15,25]. In breed-specific phenotype and performance [44]. Com- addition, the average positive predicted value of the 15 pared with the eight previous chicken CNV studies chosen CNVRs was 91.71%, similar to some previous [9,32-34,36,39-41], far more CNVRs both on average and results in animals [7,33,45], suggesting that most of the in total were found. A total of 6,478 (86.0%) autosomal positive samples detected by sequencing-based method CNVRs with slightly smaller average size (10.6 kb) were arehighlyconsistentwiththe qPCR experiments. We novel, likely due to the higher resolution and sensitivity of also estimated the false negative error rates as it is a NGS method than aCGH and SNP array. These novel common problem in CNV detection [7,49], and the CNVRs enrich significantly the published chicken CNV average percentage of false negative results was 22.43%. database. The low concordance between different studies The discrepancies between NGS results and qPCR val- results from the differences in technical issues, CNV call- idation may be due to the negative impact of potential ing algorithms as well as study populations [45], and can SNPs and small INDELs, which result in the reduced also indicate that numerous CNVs in the chicken genome primer efficiency. are still expected to be discovered. We found both maximum and mean copy number es- Promising candidate genes covered by CNVRs timates of autosomal duplicated sequences in chickens CNV is a significant source of genetic variation account- were less than those in mammalians [2,12], which might ing for disease and phenotypic diversity, owing to the be related to the relatively smaller genome size (only duplication or deletion of covered genes or regulation el- one third of a typical mammalian genome) and the lower ements [4]. Our results showed that 36.4% RefSeq genes repetitive DNA content in the chicken genome [46]. In intersected with 25.0% predicted CNVRs. It is probable addition, both the counts and sizes of gain events were that CNVs, especially deletions, are located preferably in Yi et al. BMC Genomics 2014, 15:962 Page 11 of 16 http://www.biomedcentral.com/1471-2164/15/962 gene-poor regions [13,47], because gene-rich CNVRs are in LX was larger than others. We suspect that the copy more likely to be pathogenic than gene-poor CNVRs number polymorphic locus is ubiquitous in the chicken and these deleterious CNVRs would be removed by genome, but the particularly high gene duplication in LX purifying selection [47,50]. Meanwhile, the maximum may be the result of the genetic effect of long-term artifi- CN of all genes covered by CNVRs was 12.0, suggesting cial selection such as crossing between the individuals again that the chicken genome has lower repetitive DNA with stronger bone. content [46]. It is noted that nine out of the 25 most Additionally, the copy number estimates of POPDC3 variable genes belong to four keratin subfamilies (claw, (popeye domain containing 3) in WL were found to be feather, feather-like and scale). In birds, skin appendages about twice as many as other individuals. The POPDC3 such as claws, scales, beaks and feathers are composed of gene belongs to the Popeye family encoding proteins beta (β) keratins and can prevent water loss and provide a with three potential transmembrane domains with a high barrier between the organism and external environment degree of sequence conservation, and is preferentially [51]. The avian keratin genes are significantly over- expressed in heart and skeletal muscle cells as well as represented with respect to mammals [34,48]. These highly smooth muscle cells [56]. It has been reported that the variable keratin genes suggest the scenario for the evolu- expression of two Popeye family members was upregu- tion of the β-keratin gene family through gene duplication lated in the uterus of pregnant mice [57]. The uterus has and divergence for their adaptive benefits [4,51]. Addition- been thought to be an important organ composed of ally, the four subfamilies of β-keratin genes form a cluster smooth muscle and containing the shell gland in favor on chr25, one of the more GC-rich chromosomes and con- of depositing eggshell [58]. Considering that WL is the taining a relatively larger number of minisatellites [51], most prolific egg laying chicken due to the fact that it which also result in high copy number of genes. has been extensively bred for egg production, the dupli- We validated two well-known causative genes with copy cation of the POPDC3 gene may reveal the important number polymorphism, EDN3 [30] and PRLR [31], in- differences in abilities like myometrium maturation and volved in dermal hyperpigmentation and late feathering, labor, protein secretion and eggshell formation between respectively. In our study, we used hierarchical clustering WL and other breeds. analysis based on CN content to visualize the potential re- Moreover, these enriched GO terms were mainly in- lationship among 12 breeds. For example, the heatmap for volved in cellular regulation and structure, various binding dermal hyperpigmentation grouped DX and SK together, functions as well as stress and immune responses, which and both of which are distributed in the Jiangxi province of are consistent with several previous studies [9,32-34], sug- China, suggesting that DX and SK may have a close evolu- gesting that the copy number variable genes may influence tionary relationship or purposely bred dermal hyperpig- the responses to external stimuli and provide the muta- mentation into different strains. In addition, two reported tional flexibility to adapt rapidly to changing selective copy number variable genes associated with Marek’sdis- pressures due to evolutionary adaption [59]. Most CNVRs ease, namely FZD6 (frizzled family receptor 6) and LIMS1 also spanned some QTL regions, which indicated that (LIM and senescent cell antigen-like domains 1) [39,52,53] these CNVRs may exert significant effects on traits of eco- were also found in our results. nomic interest in chickens. Furthermore, we also found some novel CNV-gene overlaps as potential candidates linked to some important Conclusions traits. For example, the SOCS2 (suppressor of cytokine sig- In this study, we performed genome-wide CNV detection naling 2) is a member of the suppressor of cytokine signal- and estimated the absolute copy number of the corre- ing family, and the related proteins are implicated in the sponding genetic locus based on whole genome sequen- negative regulation of cytokine action through inhibition cing data of 12 chickens abundant in genetic diversity, and of the JAK/STAT pathway (Janus kinase/signal transduc- constructed the highest-resolution individualized chicken ers and activators of transcription) [54]. Dual x-ray ab- CNV map so far. A total of 8,840 CNVRs were identified, sorptiometry (DXA) analysis demonstrated that SOCS2 and most of them were novel variants which could enrich inactivation resulted in reduced trabecular and cortical the current CNV database. The high CNVR confirmation volumetric bone mineral density (BMD) in SOCS2-defi- rates by aCGH and qPCR suggested that sequencing- cient mice [55]. We found that the SOCS2 had higher CN based method was more sensitive and efficient for CNV (n =6.4) in LX than in other individuals, which is particu- discovery. We detected 2,216 RefSeq genes overlapping larly interesting as the LX is known for cockfighting in with CNVRs, including genes involved in well-known phe- which the chickens with higher BMD have advantage over notypes such as dermal hyperpigmentation and late feath- others. The gene expansions were also supported by the ering. In addition, some novel genes like POPDC3 and heatmap. Additional qPCR experiments in 16 other indi- LBFABP covered by CNVRs may play an important role in viduals revealed that the increased copy number of SOCS2 production traits, and the highly duplicated SOCS2 may Yi et al. BMC Genomics 2014, 15:962 Page 12 of 16 http://www.biomedcentral.com/1471-2164/15/962 serve as an excellent candidate for bone mineral density. processing was performed using CNVnator software Our study lays the foundation for comprehensive under- based on RD method as previously described [25]. standing of copy number variation in the chicken genome CNVnator firstly calculated the counts of mapped reads and is beneficial to future association studies between within user specified non-overlapping bins of equal size CNV and important traits of chickens. as the RD signal, and then adjusted the signal in consid- eration of the potential correlation between RD signal Methods and GC content of the underlying genomic sequence. Sample collection and sequencing The mean-shift algorithm was employed to segment the We selected a total of 12 female chickens from different signal with presumably different underlying CN. Then types and genetic sources representing modern chicken CNVs were predicted by applying statistical significance populations, i.e., a Red Jungle Fowl (RJF, the ancestor of tests to the segments. A more detailed description about domestic chickens), seven Chinese indigenous chickens this method could be found in the CNVnator paper [25]. including Beijing You (BY), Dongxiang (DX), Luxi Game We ran CNVnator with a bin size of 100 bp for our data. (LX), Shouguang (SG), Silkie (SK), Tibetan (TB) and CNV calls were filtered using stringent criteria including Wenchang (WC), and four commercial breeds including P-value <0.01 and size >1 kb, and calls with >50% of q0 Cornish (CS), Rhode Island Red (RIR), White Leghorn (zero mapping quality) reads within the CNV regions (WL) and White Plymouth Rock (WR). The whole blood were removed (q0 filter). All CNV calls overlapping with samples were collected from brachial veins by standard gaps in the reference genome were excluded from con- venepuncture along with regular quarantine inspection sideration. For unlocalized and unplaced chromosomes of the experimental station of China Agricultural Uni- (chrN_random and chrun_random in UCSC, chrUn), we versity, and genomic DNA was isolated using the stand- removed them for further analysis due to the shorter ard phenol/chloroform extraction method. Whole length of the chrUn contigs and mapping ambiguity of genome sequencing for all 12 individuals was performed chrUn sequence reads. Meanwhile, we performed geno- on the HiSeq 2000 system (Illumina Inc., San Diego, typing of all 5 kb non-overlapping windows which did CA, USA). Two genomic DNA libraries of 500 bp insert not overlap with putative CNVs and gaps on autosomes. size per individual were constructed and sequenced with In order to compare our results with previous studies, 100 bp paired-end reads, and each library dataset was we converted all autosomal CNVRs from galGal4 to gal- generated with five-fold coverage depth. Library prepar- Gal3 using the UCSC liftOver tool [63]. ation and all Illumina runs were performed as the stand- ard manufacturer’s protocols. Array CGH for assessing genome-wide CNVs We conducted CNV consistency evaluation using two Quality control and sequence alignment similar whole genome tiling arrays based on galGal4 2011 For ensuring high-quality data, we used NGS QC Toolkit build. One of them is the NimbleGen aCGH (Madison, with default parameters to perform quality control of raw WI, USA), a custom-designed 3*1.4 M array containing a sequencing data, mainly by removing low-quality reads total of 1,425,178 50-75mer probes with the mean and and reads containing primer/adaptor contamination [60]. median interval of 734 bp and 700 bp. The other is the All high-quality Illumina sequence reads were aligned Agilent custom-designed 1*1.0 M array (Agilent Technol- against the galGal4 assembly by using the Burrows- ogy Inc., CA, USA), with the mean and median probe spa- Wheeler Aligner (BWA) program [61] with default param- cing of 1,056 bp and 1,050 bp. It should be noted that the eters. The draft genome sequence was retrieved from the average physical distance of the closest SNP probes be- UCSC website (http://hgdownload.soe.ucsc.edu/golden- tween two arrays was 262.6 bp and 95.2% distance inter- Path/galGal4/bigZips/). During the construction of a gen- vals were shorter than 500 bp. Meanwhile we only omic library, Illumina platform was likely to generate analyzed raw aCGH log ratio values instead of processed/ some duplicate reads named ‘PCR and optical duplicates’ normalized data. These cases could ensure reasonable ex- which imposed negative impact on the downstream ana- planation for our results although using different arrays. All lysis. So we first used SAMtools [62] to convert the .sam processing steps like DNA labeling (Cy3 for samples and files of different libraries belonging to the same individual Cy5 for references), array hybridization, data normalization to .bam files and sort and merge them, followed by re- and scanning analysis were performed following standard moval of potential PCR duplicates using Picard (http:// procedure. In each aCGH experiment, we chose the RJF as broadinstitute.github.io/picard/). the same reference sample. CNV detection Quantitative PCR confirmation Following the above filtering steps, the resulting .bam We also performed qPCR confirmation of 15 CNVRs files were utilized for CNV calling and genotyping, post- chosen from the CNVRs detected by CNVnator. Most Yi et al. BMC Genomics 2014, 15:962 Page 13 of 16 http://www.biomedcentral.com/1471-2164/15/962 chosen CNVRs have not been reported in the previous Hierarchical cluster analysis studies and are also adjacent to annotated genes. Two We used the heatmap.2() function of the gplots package distinct pairs of PCR primers were designed to target (http://cran.r-project.org/web/packages/gplots/index.html) each CNVR using Primer5.0 software for the uncer- to generate heatmap figures. We first specified the regions tainty of the CNVR boundaries. Furthermore, the UCSC extending 30 kb on each side of interested genes and used In-Silico PCR tool was used for in silico analysis of the estimated CN values of 1 kb non-overlapping windows primers specificity and sensitivity [63]. The PCCA gene in each individual for post-analysis, mainly considering which was previously identified as a non-CNV locus that the regulatory elements may be included in the up- was chosen as the control region [40]. Quality control stream or downstream of a gene. No reordering of those of all primer sets was evaluated using an 8-point stand- windows representing corresponding chromosome loca- ard curve in duplicate to ensure the similar amplifica- tions in the heatmap was made for the sake of clarity. The tion efficiencies between target and control primers. All Pearson’s correlation coefficient of the CN values was used qPCR experiments were conducted on the ABI Prism as the distance measure of the agglomerative hierarchical 7500 sequence detection system (Applied Biosystems clustering with average linkage, and to generate hierarch- group) using SYBR green chemistry in triplicate reac- ical cluster dendrograms. tions, each with a reaction volume of 15 μl in a 96-well plate. The condition for thermal cycle was as follows: Availability of supporting data 1 cycle of pre-incubation at 50°C for 2 min and 95°C for All raw sequence data has been deposited in NCBI Se- 10 min, 40 cycles of amplification (95°C for 15 s and 60°C quence Read Achieve (SRA) under the Bioproject number (1 - ΔΔCt) for 1 min). We used the formula 2 method to cal- PRJNA232548. The experiment numbers for the 12 chick- culate the relative copy number for each test region. The ens are SRX408161-SRX408172. All aCGH data have been cycle threshold (Ct) value of each test sample was first submitted to the NCBI Gene Expression Omnibus (GEO) normalized against the control region, and then the ΔCt (http://www.ncbi.nlm.nih.gov/geo/) under accession num- value was calculated between the test sample and a prese- ber GSE54119. lected reference sample predicted with normal copy num- ber status by CNVnator. The golden standard of each Additional files diploid CNV was generally considered to have two copies for autosomes or one copy when the locus was on Z Additional file 1: Table S1. Summary of identified CNVs and CNVRs in the 12 chicken genomes. chromosome (chrZ) of a female in chickens. Additional file 2: Figure S1. Individualized chicken CNV map in the chicken genome. The horizontal black lines represent the draft chicken genome (UCSC version galGal4). Tracks under the chromosomes indicate Gene contents and functional annotation corresponding CNV status of all individuals kept in the alphabetical order from top to bottom, for BY, CS, DX, LX, RIR, RJF, SG, SK, TB, WC, WL and WR. The RefSeq gene list was retrieved from the UCSC Merged CNVRs from all individuals are depicted above chromosomes. The RefSeq database [63]. All miRNA genes were excluded colors for each bar denote different copy number (CN) in CNV legend and because the nucleotide sequences were too short to es- different types of CNVRs. The downmost axis shows the chromosome, CNV and CNVR coordinates. Left-hand chromosomes are ordered from left to timate reliable copy number. We analyzed the propor- right, and the right-hands are just reversed. tion of the RefSeq genes overlapping with putative Additional file 3: Table S2. General statistics of the CNVRs on CNVRs and performed CN estimates for all 6,086 non- each chromosome. redundant RefSeq gene transcripts. In addition, to provide Additional file 4: Table S3. Summary of novel or reported CNVRs insight into the functional enrichment of the RefSeq genes on autosomes. covered by CNVRs, we performed Gene Ontology (GO) Additional file 5: Figure S2. Correlation between digital aCGH and whole-genome aCGH among nine individuals compared with Red Jungle functional annotation and Kyoto Encyclopedia of Genes Fowl (RJF). RJF is selected as the reference sample in each aCGH experiment. and Genomes (KEGG) pathway analysis employing the Digital aCGH values are estimated using calculated log CN ratios in which web-accessible program DAVID [64]. Statistical signifi- CN are estimated for identified CNV segments of nine individuals and divided by the corresponding CN of RJF. Whole genome aCGH values are cance was accessed by using a modified Fisher’sexact defined as the average of all probes log ratio values in the same segments test and Benjamini correction for multiple testing (P- as the digital aCGH. value <0.05). We also compared the CNVRs identified Additional file 6: Table S4. Primers information and confirmation in this study with the reported QTLs obtained from the results of the 15 chosen CNVRs by qPCR analysis. chicken QTL database [65]. We focused on the QTLs Additional file 7: Figure S3. Illustrating of qPCR confirmation results for three selected CNVRs of different types. X-axis represents all 12 samples and with confidence interval less than 10 Mb and consid- Y-axis represents normalized ratios (NR) estimated by qPCR. NR around 2 ered those QTLs with overlapped confidence intervals indicates normal status (2 copies), NR around 0 or 1 indicates loss status greater than 50% as the same QTL [45], because the (0 copies or 1 copy), and NR around 3 or more indicates gain status (3 or more copies). (A) Results for a gain status of CNVR3588. (B) Results for a loss QTL confidence intervals were too large to be used effi- status of CNVR6695. (C) Results for a both status of CNVR410. ciently in the post-processing. Yi et al. BMC Genomics 2014, 15:962 Page 14 of 16 http://www.biomedcentral.com/1471-2164/15/962 cattle genomes using next-generation sequencing. Genome Res 2012, Additional file 8: Table S5. The detailed features of RefSeq genes 22(4):778–790. completely or partial overlapped with CNVRs. 3. McCarroll SA, Altshuler DM: Copy-number variation and association Additional file 9: Figure S4. Visual examination by read depth, studies of human disease. Nat Genet 2007, 39(7 Suppl):S37–S42. whole-genome aCGH and digital aCGH around three loci for 12 4. Zhang F, Gu W, Hurles ME, Lupski JR: Copy number variation in human chicken genomes. The uppermost gene image is generated with the health, disease, and evolution. Annu Rev Genomics Hum Genet 2009, UCSC Genome Browser (http://genome.ucsc.edu/) using the galGal4 10:451–481. assembly. The track below the gene region is depth of coverage for all 5. Altshuler DM, Gibbs RA, Peltonen L, Dermitzakis E, Schaffner SF, Yu F, 12 individual genomes. Red indicates regions of excess read depth Bonnen PE, de Bakker PI, Deloukas P, Gabriel SB, Gwilliam R, Hunt S, Inouye M, (> mean +3 × STDEV), whereas gray indicates intermediate read depth Jia X, Palotie A, Parkin M, Whittaker P, Chang K, Hawes A, Lewis LR, Ren Y, (mean +2 × STDEV < × < mean +3 × STDEV), and green indicates Wheeler D, Muzny DM, Barnes C, Darvishi K, Hurles M, Korn JM, Kristiansson K, normal read depth (mean ± 2 × STDEV). All read depth values based Lee C, McCarrol SA, et al: Integrating common and rare genetic variation in on 1 kb non-overlapping windows are corrected by GC content. diverse human populations. Nature 2010, 467(7311):52–58. Whole-genome aCGH and digital aCGH values are depicted as the 6. Yalcin B, Wong K, Agam A, Goodson M, Keane TM, Gan X, Nellaker C, red-green histograms and correspond to a gain colored in green Goodstadt L, Nicod J, Bhomra A, Hernandez-Pliego P, Whitley H, Cleak J, (>0.5), a loss colored in red (<−0.5) and normal status colored in gray Dutton R, Janowitz D, Mott R, Adams DJ, Flint J: Sequence-based (−0.5 < x <0.5). (A) Two previously reported CNVs (chr20: 11,111,401-11,238,900 characterization of structural variation in the mouse genome. and chr20: 11,651,801-11,822,900) associated with dermal hyperpigmentation. Nature 2011, 477(7364):326–329. The DX and SK genomes show two additional copies of the two regions 7. Wang J, Jiang J, Fu W, Jiang L, Ding X, Liu JF, Zhang Q: A genome-wide compared with RJF, and are also validated by whole-genome aCGH. (B) A detection of copy number variations using SNP genotyping arrays in higher copy number increase for the SOCS2 locus (chr1: 44,764,280-44,765,955) swine. BMC Genomics 2012, 13:273. is predictedinLXthaninother individuals. (C)The POPDC3 gene (chr3: 8. Liu GE, Hou Y, Zhu B, Cardone MF, Jiang L, Cellamare A, Mitra A, Alexander 68,255,196-68,259,535) is predicted to be duplicated status only in WL. LJ, Coutinho LL, Dell'Aquila ME, Gasbarre LC, Lacalandra G, Li RW, Additional file 10: Table S6. Functional enrichment of GO and KEGG Matukumalli LK, Nonneman D, Regitano LC, Smith TP, Song J, Sonstegard pathway analysis of RefSeq genes covered by CNVRs. TS, Van Tassell CP, Ventura M, Eichler EE, McDaneld TG, Keele JW: Analysis of copy number variations among diverse cattle breeds. Genome Res Additional file 11: Table S7. The overlap information of QTLs and 2010, 20(5):693–703. CNVRs across the chicken genome. 9. Wang Y, Gu X, Feng C, Song C, Hu X, Li N: A genome-wide survey of copy number variation regions in various chicken breeds by array comparative Abbreviations genomic hybridization method. Anim Genet 2012, 43(3):282–289. CNV: Copy number variation; CNVR: Copy number variation region; 10. Hastings PJ, Ira G, Lupski JR: A microhomology-mediated break-induced aCGH: array comparative genomic hybridization; qPCR: quantitative replication model for the origin of human copy number variation. polymerase chain reaction; NAHR: Non-allelic homologous recombination; PLoS Genet 2009, 5(1):e1000327. NHEJ: Non-homologous end joining; FoSTes: Fork stalling and template 11. Sharp AJ, Locke DP, McGrath SD, Cheng Z, Bailey JA, Vallente RU, Pertz LM, switching; SD: Segmental duplication; GWAS: Genome-wide association Clark RA, Schwartz S, Segraves R, Oseroff VV, Albertson DG, Pinkel D, Eichler studies; RD: Read depth; SNP: Single nucleotide polymorphism; NGS: EE: Segmental duplications and copy-number variation in the human Next-generation sequencing; GO: Gene ontology; KEGG: Kyoto encyclopedia genome. Am J Hum Genet 2005, 77(1):78–88. of genes and genomes; QTL: Quantitative trait loci; MHC: Major 12. Alkan C, Kidd JM, Marques-Bonet T, Aksay G, Antonacci F, Hormozdiari F, histocompatibility complex; MD: Marek’s disease; BMD: Bone mineral density. Kitzman JO, Baker C, Malig M, Mutlu O, Sahinalp SC, Gibbs RA, Eichler EE: Personalized copy number and segmental duplication maps using Competing interests next-generation sequencing. Nat Genet 2009, 41(10):1061–1067. The authors declare that they have no competing interests. 13. Freeman JL, Perry GH, Feuk L, Redon R, McCarroll SA, Altshuler DM, Aburatani H, Jones KW, Tyler-Smith C, Hurles ME, Carter NP, Scherer SW, Lee Authors’ contributions C: Copy number variation: new insights in genome diversity. Genome Res NY and LQ conceived and designed all experiments. GY, LQ and YY 2006, 16(8):949–961. performed bioinformatics and statistical analysis with the help of JL, and 14. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, carried out aCGH and qPCR experiments. GX provided samples. GY and LQ McCarthy MI, Ramos EM, Cardon LR, Chakravarti A, Cho JH, Guttmacher AE, drafted the manuscript. NY revised the paper. All authors read and approved Kong A, Kruglyak L, Mardis E, Rotimi CN, Slatkin M, Valle D, Whittemore AS, the final manuscript. Boehnke M, Clark AG, Eichler EE, Gibson G, Haines JL, Mackay TF, McCarroll SA, Visscher PM: Finding the missing heritability of complex diseases. Nature 2009, 461(7265):747–753. Acknowledgements 15. Conrad DF, Pinto D, Redon R, Feuk L, Gokcumen O, Zhang Y, Aerts J, We thank Dr. Xiquan Zhang for providing some samples. This work was Andrews TD, Barnes C, Campbell P, Fitzgerald T, Hu M, Ihm CH, Kristiansson funded in part by grants of National High Technology Development Plan of K, Macarthur DG, Macdonald JR, Onyiah I, Pang AW, Robson S, Stirrups K, China (2013AA102501), Natural Science Foundation of China (31320103905), Valsesia A, Walter K, Wei J, Tyler-Smith C, Carter NP, Lee C, Scherer SW, Programs for Changjiang Scholars and Innovative Research in University Hurles ME: Origins and functional impact of copy number variation in (IRT1191), and China Agriculture Research Systems (CARS-41). the human genome. Nature 2010, 464(7289):704–712. Received: 4 March 2014 Accepted: 13 October 2014 16. Liu GE, Bickhart DM: Copy number variation in the cattle genome. Published: 7 November 2014 Funct Integr Genomics 2012, 12(4):609–624. 17. Stranger BE, Forrest MS, Dunning M, Ingle CE, Beazley C, Thorne N, Redon R, Bird CP, de Grassi A, Lee C, Tyler-Smith C, Carter N, Scherer SW, Tavare S, References Deloukas P, Hurles ME, Dermitzakis ET: Relative impact of nucleotide and 1. Redon R, Ishikawa S, Fitch KR, Feuk L, Perry GH, Andrews TD, Fiegler H, copy number variation on gene expression phenotypes. Science 2007, Shapero MH, Carson AR, Chen W, Cho EK, Dallaire S, Freeman JL, Gonzalez 315(5813):848–853. JR, Gratacos M, Huang J, Kalaitzopoulos D, Komura D, MacDonald JR, 18. Henrichsen CN, Chaignat E, Reymond A: Copy number variants, diseases Marshall CR, Mei R, Montgomery L, Nishimura K, Okamura K, Shen F, and gene expression. Hum Mol Genet 2009, 18(R1):R1–R8. Somerville MJ, Tchinda J, Valsesia A, Woodwark C, Yang F, et al: Global 19. LaFramboise T: Single nucleotide polymorphism arrays: a decade of variation in copy number in the human genome. Nature 2006, biological, computational and technological advances. Nucleic Acids Res 444(7118):444–454. 2009, 37(13):4181–4193. 2. Bickhart DM, Hou Y, Schroeder SG, Alkan C, Cardone MF, Matukumalli LK, 20. Pinto D, Darvishi K, Shi X, Rajan D, Rigler D, Fitzgerald T, Lionel AC, Song J, Schnabel RD, Ventura M, Taylor JF, Garcia JF, Van Tassell CP, Thiruvahindrapuram B, Macdonald JR, Mills R, Prasad A, Noonan K, Gribble S, Sonstegard TS, Eichler EE, Liu GE: Copy number variation of individual Yi et al. BMC Genomics 2014, 15:962 Page 15 of 16 http://www.biomedcentral.com/1471-2164/15/962 Prigmore E, Donahoe PK, Smith RS, Park JH, Hurles ME, Carter NP, Lee C, 39. Luo J, Yu Y, Mitra A, Chang S, Zhang H, Liu G, Yang N, Song J: Genome- Scherer SW, Feuk L: Comprehensive assessment of array-based platforms wide copy number variant analysis in inbred chickens lines with and calling algorithms for detection of copy number variants. different susceptibility to Marek’s disease. G3 (Bethesda) 2013, Nat Biotechnol 2011, 29(6):512–520. 3(2):217–223. 21. Bentley DR, Balasubramanian S, Swerdlow HP, Smith GP, Milton J, Brown 40. Wang X, Nahashon S, Feaster TK, Bohannon-Stewart A, Adefope N: An initial CG, Hall KP, Evers DJ, Barnes CL, Bignell HR, Boutell JM, Bryant J, Carter RJ, map of chromosomal segmental copy number variations in the chicken. Keira Cheetham R, Cox AJ, Ellis DJ, Flatbush MR, Gormley NA, Humphray SJ, BMC Genomics 2010, 11:351. Irving LJ, Karbelashvili MS, Kirk SM, Li H, Liu X, Maisinger KS, Murray LJ, 41. Abernathy J, Li X, Jia X, Chou W, Lamont SJ, Crooijmans R, Zhou H: Copy Obradovic B, Ost T, Parkinson ML, Pratt MR, et al: Accurate whole human number variation in fayoumi and leghorn chickens analyzed using array genome sequencing using reversible terminator chemistry. Nature 2008, comparative genomic hybridization. Anim Genet 2014, 45(3):400–411. 456(7218):53–59. 42. Munoz-Amatriain M, Eichten SR, Wicker T, Richmond TA, Mascher M, 22. Yoon S, Xuan Z, Makarov V, Ye K, Sebat J: Sensitive and accurate detection Steuernagel B, Scholz U, Ariyadasa R, Spannagl M, Nussbaumer T, Mayer KF, of copy number variants using read depth of coverage. Genome Res Taudien S, Platzer M, Jeddeloh JA, Springer NM, Muehlbauer GJ, Stein N: 2009, 19(9):1586–1592. Distribution, functional impact, and origin mechanisms of copy number 23. Campbell CD, Sampas N, Tsalenko A, Sudmant PH, Kidd JM, Malig M, Vu TH, variation in the barley genome. Genome Biol 2013, 14(6):R58. Vives L, Tsang P, Bruhn L, Eichler EE: Population-genetic properties of 43. Skinner BM, Al Mutery A, Smith D, Volker M, Hojjat N, Raja S, Trim S, Houde differentiated human copy-number polymorphisms. Am J Hum Genet P, Boecklen WJ, Griffin DK: Global patterns of apparent copy number 2011, 88(3):317–332. variation in birds revealed by cross-species comparative genomic 24. Teo SM, Pawitan Y, Ku CS, Chia KS, Salim A: Statistical challenges hybridization. Chromosome Res 2014, 22(1):59–70. associated with detecting copy number variations with next-generation 44. Clop A, Vidal O, Amills M: Copy number variation in the genomes of sequencing. Bioinformatics 2012, 28(21):2711–2718. domestic animals. Anim Genet 2012, 43(5):503–517. 25. Abyzov A, Urban AE, Snyder M, Gerstein M: CNVnator: an approach to 45. Jiang L, Jiang J, Yang J, Liu X, Wang J, Wang H, Ding X, Liu J, Zhang Q: discover, genotype, and characterize typical and atypical CNVs from family Genome-wide detection of copy number variations using high-density and population genome sequencing. Genome Res 2011, 21(6):974–984. SNP genotyping platforms in Holsteins. BMC Genomics 2013, 14:131. 26. Szatkiewicz JP, Wang W, Sullivan PF, Sun W: Improving detection of 46. Burt DW: Chicken genome: current status and future opportunities. copy-number variation by simultaneous bias correction and read-depth Genome Res 2005, 15(12):1692–1698. segmentation. Nucleic Acids Res 2013, 41(3):1519–1532. 47. Conrad DF, Andrews TD, Carter NP, Hurles ME, Pritchard JK: A 27. Norris BJ, Whan VA: A gene duplication affecting expression of the ovine high-resolution survey of deletion polymorphism in the human genome. ASIP gene is responsible for white and black sheep. Genome Res 2008, Nat Genet 2006, 38(1):75–81. 18(8):1282–1293. 48. International Chicken Genome Sequencing Consortium: Sequence and 28. Rosengren Pielberg G, Golovko A, Sundstrom E, Curik I, Lennartsson J, comparative analysis of the chicken genome provide unique Seltenhammer MH, Druml T, Binns M, Fitzsimmons C, Lindgren G, Sandberg K, perspectives on vertebrate evolution. Nature 2004, 432(7018):695–716. Baumung R, Vetterlein M, Stromberg S, Grabherr M, Wade C, Lindblad-Toh K, 49. Nicholas TJ, Cheng Z, Ventura M, Mealey K, Eichler EE, Akey JM: The Ponten F, Heldin CH, Solkner J, Andersson L: A cis-acting regulatory mutation genomic architecture of segmental duplications and associated copy causes premature hair graying and susceptibility to melanoma in the horse. number variants in dogs. Genome Res 2009, 19(3):491–499. Nat Genet 2008, 40(8):1004–1009. 50. Lee C, Iafrate AJ, Brothman AR: Copy number variations and clinical 29. Wright D, Boije H, Meadows JR, Bed'hom B, Gourichon D, Vieaud A, Tixier- cytogenetic diagnosis of constitutional disorders. Nat Genet 2007, Boichard M, Rubin CJ, Imsland F, Hallbook F, Andersson L: Copy number 39(7 Suppl):S48–S54. variation in intron 1 of SOX5 causes the Pea-comb phenotype in 51. Greenwold MJ, Sawyer RH: Genomic organization and molecular chickens. PLoS Genet 2009, 5(6):e1000512. phylogenies of the beta (beta) keratin multigene family in the chicken 30. Dorshorst B, Molin AM, Rubin CJ, Johansson AM, Stromstedt L, Pham MH, (Gallus gallus) and zebra finch (Taeniopygia guttata): implications for Chen CF, Hallbook F, Ashwell C, Andersson L: A complex genomic feather evolution. BMC Evol Biol 2010, 10:148. rearrangement involving the endothelin 3 locus causes dermal 52. Wang-Rodriguez J, Dreilinger AD, Alsharabi GM, Rearden A: The signaling hyperpigmentation in the chicken. PLoS Genet 2011, 7(12):e1002412. adapter protein PINCH is up-regulated in the stroma of common 31. Elferink MG, Vallee AA, Jungerius AP, Crooijmans RP, Groenen MA: Partial cancers, notably at invasive edges. Cancer 2002, 95(6):1387–1395. duplication of the PRLR and SPEF2 genes at the late feathering locus in 53. Chen K, Tu Y, Zhang Y, Blair HC, Zhang L, Wu C: PINCH-1 regulates the chicken. BMC Genomics 2008, 9:391. ERK-Bim pathway and contributes to apoptosis resistance in cancer cells. 32. Jia X, Chen S, Zhou H, Li D, Liu W, Yang N: Copy number variations J Biol Chem 2008, 283(5):2508–2517. identified in the chicken using a 60 K SNP BeadChip. Anim Genet 2013, 54. Metcalf D, Greenhalgh CJ, Viney E, Willson TA, Starr R, Nicola NA, Hilton DJ, 44(3):276–284. Alexander WS: Gigantism in mice lacking suppressor of cytokine 33. Tian M, Wang Y, Gu X, Feng C, Fang S, Hu X, Li N: Copy number variants in signalling-2. Nature 2000, 405(6790):1069–1073. locally raised Chinese chicken genomes determined using array 55. Lorentzon M, Greenhalgh CJ, Mohan S, Alexander WS, Ohlsson C: Reduced comparative genomic hybridization. BMC Genomics 2013, 14(1):262. bone mineral density in SOCS-2-deficient mice. Pediatr Res 2005, 34. Crooijmans RP, Fife MS, Fitzgerald TW, Strickland S, Cheng HH, Kaiser P, 57(2):223–226. Redon R, Groenen MA: Large scale variation in DNA copy number in 56. Brand T: The Popeye domain-containing gene family. Cell Biochem Biophys chicken breeds. BMC Genomics 2013, 14:398. 2005, 43(1):95–103. 35. Griffin DK, Robertson LB, Tempest HG, Vignal A, Fillon V, Crooijmans RP, 57. Andree B, Hillemann T, Kessler-Icekson G, Schmitt-John T, Jockusch H, Groenen MA, Deryusheva S, Gaginskaya E, Carre W, Waddington D, Talbot R, Arnold HH, Brand T: Isolation and characterization of the novel popeye Volker M, Masabanda JS, Burt DW: Whole genome comparative studies gene family expressed in skeletal muscle and heart. Dev Biol 2000, between chicken and turkey and their implications for avian genome 223(2):371–382. evolution. BMC Genomics 2008, 9:168. 58. Hincke MT, Nys Y, Gautron J, Mann K, Rodriguez-Navarro AB, McKee MD: 36. Fan WL, Ng CS, Chen CF, Lu MY, Chen YH, Liu CJ, Wu SM, Chen CK, Chen JJ, The eggshell: structure, composition and mineralization. Front Biosci Mao CT, Lai YT, Lo WS, Chang WH, Li WH: Genome-wide patterns of (Landmark Ed) 2012, 17:1266–1280. genetic variation in two domestic chickens. Genome Biol Evol 2013, 59. Gokcumen O, Babb PL, Iskow RC, Zhu Q, Shi X, Mills RE, Ionita-Laza I, Vallender EJ, 5(7):1376–1392. Clark AG, Johnson WE, Lee C: Refinement of primate copy number variation 37. Qu L, Li X, Xu G, Chen K, Yang H, Zhang L, Wu G, Hou Z, Yang N: hotspots identifies candidate genomic regions evolving under positive Evaluation of genetic diversity in Chinese indigenous chicken breeds selection. Genome Biol 2011, 12(5):R52. using microsatellite markers. Sci China C Life Sci 2006, 49(4):332–341. 60. Patel RK, Jain M: NGS QC Toolkit: a toolkit for quality control of next 38. Sudmant PH, Kitzman JO, Antonacci F, Alkan C, Malig M, Tsalenko A, generation sequencing data. PLoS One 2012, 7(2):e30619. Sampas N, Bruhn L, Shendure J, Eichler EE: Diversity of human copy 61. Li H, Durbin R: Fast and accurate short read alignment with burrows- number variation and multicopy genes. Science 2010, 330(6004):641–646. wheeler transform. Bioinformatics 2009, 25(14):1754–1760. Yi et al. BMC Genomics 2014, 15:962 Page 16 of 16 http://www.biomedcentral.com/1471-2164/15/962 62. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The sequence alignment/map format and SAMtools. Bioinformatics 2009, 25(16):2078–2079. 63. Karolchik D, Barber GP, Casper J, Clawson H, Cline MS, Diekhans M, Dreszer TR, Fujita PA, Guruvadoo L, Haeussler M, Harte RA, Heitner S, Hinrichs AS, Learned K, Lee BT, Li CH, Raney BJ, Rhead B, Rosenbloom KR, Sloan CA, Speir ML, Zweig AS, Haussler D, Kuhn RM, Kent WJ: The UCSC Genome Browser database: 2014 update. Nucleic Acids Res 2014, 42(Database issue):D764-770. 64. Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009, 4(1):44–57. 65. Hu ZL, Park CA, Wu XL, Reecy JM: Animal QTLdb: an improved database tool for livestock animal QTL/association data dissemination in the post-genome era. Nucleic Acids Res 2013, 41(Database issue):D871–D879. doi:10.1186/1471-2164-15-962 Cite this article as: Yi et al.: Genome-wide patterns of copy number variation in the diversified chicken genomes using next-generation sequencing. BMC Genomics 2014 15:962. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Genomics Springer Journals

Genome-wide patterns of copy number variation in the diversified chicken genomes using next-generation sequencing

Download PDF
16 pages

Loading next page...
 
/lp/springer-journals/genome-wide-patterns-of-copy-number-variation-in-the-diversified-WPC84d886r

References (86)

Publisher
Springer Journals
Copyright
2014 Yi et al.; licensee BioMed Central Ltd.
eISSN
1471-2164
DOI
10.1186/1471-2164-15-962
Publisher site
See Article on Publisher Site

Abstract

Background: Copy number variation (CNV) is important and widespread in the genome, and is a major cause of disease and phenotypic diversity. Herein, we performed a genome-wide CNV analysis in 12 diversified chicken genomes based on whole genome sequencing. Results: A total of 8,840 CNV regions (CNVRs) covering 98.2 Mb and representing 9.4% of the chicken genome were identified, ranging in size from 1.1 to 268.8 kb with an average of 11.1 kb. Sequencing-based predictions were confirmed at a high validation rate by two independent approaches, including array comparative genomic hybridization (aCGH) and quantitative PCR (qPCR). The Pearson’s correlation coefficients between sequencing and aCGH results ranged from 0.435 to 0.755, and qPCR experiments revealed a positive validation rate of 91.71% and a false negative rate of 22.43%. In total, 2,214 (25.0%) predicted CNVRs span 2,216 (36.4%) RefSeq genes associated with specific biological functions. Besides two previously reported copy number variable genes EDN3 and PRLR,we also found some promising genes with potential in phenotypic variation. Two genes, FZD6 and LIMS1, related to disease susceptibility/resistance are covered by CNVRs. The highly duplicated SOCS2 may lead to higher bone mineral density. Entire or partial duplication of some genes like POPDC3 may have great economic importance in poultry breeding. Conclusions: Our results based on extensive genetic diversity provide a more refined chicken CNV map and genome-wide gene copy number estimates, and warrant future CNV association studies for important traits in chickens. Keywords: Copy number variation, Whole genome sequencing, aCGH, Genetic diversity, Chicken Background sequences (insertions) of ≥1 kb in length and ≥90% se- Copy number variations (CNVs) are defined as gains or quence identity are also suggested to be one of the major losses of DNA fragments of 50 bp or longer in length in catalysts and hotspots for CNV formation [11,12], mainly comparison with reference genome [1,2]. CNVs contrib- because the genomic regions flanked by SDs are suscep- ute significantly to both disease susceptibility/resistance tible to recurrent rearrangement by NAHR [11,13]. In and normal phenotypic variability in humans [3-5] and terms of total bases involved, the percentage of the genome animals [6-9]. Four major mechanisms have been found affected by CNVs is higher than that of single nucleotide to be related to CNV formation including non-allelic hom- polymorphism (SNP) markers. Although SNPs are gener- ologous recombination (NAHR), non-homologous end ally considered as more suitable markers in the genome- joining (NHEJ), Fork Stalling and Template Switching wide association studies (GWASs), most reported SNP (FoSTeS) and LINE1 Retrotransposition [4,10]. Addition- variants have relatively limited effects and explain only a ally, segmental duplications (SDs) which are duplicated small proportion of phenotypic variance [14]. Further, CNVs encompassing genes or regulatory elements are be- * Correspondence: [email protected] lieved to exert potentially larger effects on gene expression Equal contributors through changing gene structure and dosage, altering gene Department of Animal Genetics and Breeding, College of Animal Science regulation, exposing recessive alleles and other mechanisms and Technology, China Agricultural University, Beijing, China © 2014 Yi et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Yi et al. BMC Genomics 2014, 15:962 Page 2 of 16 http://www.biomedcentral.com/1471-2164/15/962 [1,4,15,16]. CNVs are also found to alone capture 18 to 30% 12 chickens from multiple breeds with extensive genetic of the total detected genetic variation in gene expression in diversity, including seven Chinese indigenous breeds humans and animals, and may contribute to a fraction of [37], four commercial breeds and one Red Jungle Fowl. the missing heritability [17,18].Therefore,identificationof Then we applied NGS-based method to construct a CNVs is essential in whole genome fine-mapping of CNVs more refined and individualized chicken CNV map, in- and association studies for important phenotypes. vestigate genome-wide CNV characteristics and estimate Originally, two cost-effective and high-throughput genome-wide gene copy number. The results will enable methods including array comparative genomic hybridization us to better understand the patterns of CNVs in the (aCGH) and commercial SNP microarrays are used for chicken genome and future CNV association studies. CNV screening [19,20]. However, due to the limitation in resolution and sensitivity, it is difficult with the two ap- Results proaches to detect small CNVs shorter than 10 kb in Mapping statistics and CNV detection length and identify the precise boundaries of CNVs We performed whole genome sequencing in 12 different [21,22]. The two analytical platforms also reveal inconsist- breeds of female chickens using Illumina paired-end li- ent results with poor overlaps owing to different designs braries and obtained a total of 12.9 Gb high quality se- and probe densities and coordinates [18,20]. Furthermore, quence data per individual after quality filtering. After the presence of SD regions is also a common challenge for sequence alignment and removing potential PCR dupli- microarrays, because a considerable proportion of CNVs cates, the sequencing depth varied from 8.2× (CS) to fall into SD regions not well-covered by microarrays 12.4× (WR), which was sufficient for CNV detection, and [2,23]. Recently, a variety of CNV detection approaches the average coverage with respect to the chicken reference based on next-generation sequencing (NGS) were pro- genome sequence was 97.2% (Table 1). We calculated the posed and offer promising alternatives as they have higher average read depth (RD) of 5 kb non-overlapping windows effective resolution to discover CNVs with more types and for all autosomes and performed GC correction. The GC- wider size ranges [24]. One leading method is read depth adjusted RD mean and standard deviation (STDEV) for (RD) (also known as depth of coverage (DOC)) with the each individual are listed in Table 1. We applied the pro- capability of inferring gains or losses of DNA segments gram CNVnator to 12 individuals and the average number and determining absolute copy number values, which de- of CNVs per individual was 1,328, ranging from 644 in tects CNVs by analyzing the number of reads that fall into WL to 1,921 CNVs in BY. A detailed description of CNV each pre-specified window with a certain size [25,26]. calls can be found in Additional file 1: Table S1. For all the Hence, the advent of NGS technologies and suitable ana- autosomal CNVs classified as duplications, the average lytical methods can promote to systematically identify copy number value of all the individuals was 3.78 and the CNVs at higher resolution and sensitivity. maximum copy number estimate was 40.8 on chromo- At present, the three aforementioned high-throughput some 2 (chr2) in RJF. platforms have been applied to livestock genomics for A total of 8,840 CNV regions (CNVRs) allowing for CNV detection, such as sheep [27], horse [28] and cattle CNV overlaps of 1 bp or greater were obtained, covering [2], and uncover several CNVs associated with important chromosomes 1–28, two linkage groups and sex chromo- phenotypes. Some CNVs are also found to be the genetic somes, which amounted to 98.2 Mb of the chicken gen- foundation of phenotypic variation in chickens. A dupli- ome and corresponded to 9.4% of the genome sequence cated sequence close to the first intron of SOX5 is asso- (Additional file 1: Table S1). The individualized chicken ciated with the chicken pea-comb phenotype [29] and CNV map across the genome is shown in Additional file 2: an inverted duplication containing EDN3 causes dermal Figure S1. The length of CNVRs ranged from 1.1 to hyperpigmentation [30]. Partial duplication of the PRLR 268.8 kb with an average of 11.1 kb and a median of 6.6 kb. is related to the late feathering [31]. In total, 6,137 (69.4%) out of all CNVRs had sizes varying A genome-wide chicken CNV analysis is desired since from 1.1 to 10 kb (Figure 1A). Although chr1 had a max- chicken is not only an economically important farm ani- imum of 1,928 CNVRs, the two largest CNVR density mal but also a valuable biomedical model [9,32]. How- values, defined as the average distance between CNVRs, ever, some previous CNV studies in chickens based on were 35.7 kb and 32.0 kb on the chr16 and chrLGE64, re- aCGH and SNP platforms mainly suffered from low spectively (Additional file 3: Table S2). The number of resolution and sensitivity [9,32-35]. A latest study exhib- CNVRs in different individuals varied greatly, ranging from ited the detection of four main types of genetic variation 629 in WL to 1,890 in BY. Among all CNVRs, 6,083 from whole genome sequencing data using two chickens (68.8%) were present in a single individual, 1,423 (16.1%) [36], suggesting the efficiency of CNV detection via deep shared in two individuals, and 1,334 (15.1%) shared in at sequencing. Considering that a great number of CNVs least three individuals (Figure 1B). Further, the mean and appears to be segregating in distinct breeds, we selected median lengths of the unique CNVRs were 8.9 kb and Yi et al. BMC Genomics 2014, 15:962 Page 3 of 16 http://www.biomedcentral.com/1471-2164/15/962 Table 1 Summary statistics for sequencing and CNVs of 12 individuals Chicken Numbers of Depth Coverage (%) Autosome reads Autosome Duplications Deletions Sequence a b abbreviation mapped reads per 5 kb window reads STDEV covered (Mb) BY 102,002,937 9.7 97.0 489.29 110.73 1,319 602 34.1 CS 85,383,494 8.2 96.9 409.93 101.42 1,132 663 26.8 DX 129,847,015 12.4 97.4 623.50 130.46 552 820 8.2 LX 105,152,881 10.0 97.3 503.82 112.74 898 821 11.7 RIR 102,464,756 9.8 97.3 490.96 108.21 578 669 8.3 RJF 105,517,587 10.1 97.2 504.23 113.52 702 620 9.8 SG 85,987,827 8.2 96.6 412.27 87.66 470 553 7.2 SK 95,322,371 9.1 97.1 457.21 100.61 773 657 12.3 TB 107,535,104 10.3 97.3 515.68 108.07 607 679 8.5 WC 119,116,969 11.4 97.4 572.35 121.39 710 768 9.8 WL 118,689,980 11.3 97.5 567.18 118.63 203 441 3.3 WR 130,307,416 12.4 97.6 625.01 132.32 224 477 3.3 BY Beijing You, CS Cornish, DX Dongxiang, LX Luxi Game, RIR Red Island Rhode, RJF Red Jungle Fowl, SG Shouguang, SK Silkie, TB Tibetan, WC Wenchang, WL White Leghorn, WR White Plymouth Rock. The number of reads per 5 kb windows after GC correction. 5.8 kb, whereas the shared CNVRs sizes were 15.9 kb on In total, 7,530 out of 8,487 (88.7%) autosomal CNVRs were average and 9.5 kb as the median. According to the type of converted successfully. The detailed comparison results are CNVRs, they were divided into three categories, including presented in Table 2 and Additional file 4: Table S3. In our 4,761 gain, 3,773 loss and 306 both (gain and loss) events. results, 1,052 (14.0%) CNVRs with the total length of Gain events possessed longer genomic sequences than 19.7 Mb were reported by eight previous studies, and the losses both on average (14.2 kb vs. 5.4 kb) and in total remaining 6,478 (86.0%) were regarded as novel CNVRs. It (67.6 Mb vs. 20.3 Mb). In addition, the count of CNVRs on should be noted that the novel CNVRs had slightly smaller each chromosome was directly proportional to the chromo- sizes (10.6 kb) on average than those reported CNVRs some length, and five macrochromosomes (chr1-5) pos- (18.8 kb). As a special and important chromosome in the sessed a large proportion (61.8%) of all putative CNVRs. chicken genome, chr16 encompassed some CNVRs which could be confirmed by different platforms. Comparison with previous chicken CNV studies Considering that most of the previous studies excluded CNV quality assessment by CNVnator, aCGH and qPCR the CNVRs on sex chromosomes and unassigned linkage The copy number values of diploid regions on auto- groups, we migrated our autosomal CNVR coordinates somes theoretically equal to two, so we could inspect the from galGal4 to galGal3 using the UCSC liftOver tool. potential of CNVnator to generate false positive results by Figure 1 The length and frequency distribution of CNVRs. (A) 6,137 (69.4%) CNVR events are shorter than 10 kb, and the number of CNVRs longer than 50 kb is only 291 (3.3%). (B) 6,083 (68.8%) CNVR occur in only one individual and 2,757 (31.2%) CNVRs are shared in at least two individuals. Yi et al. BMC Genomics 2014, 15:962 Page 4 of 16 http://www.biomedcentral.com/1471-2164/15/962 Table 2 Comparison between autosomal CNVRs identified in this study and other chicken studies Platforms Results from different studies Overlapped CNVRs in this study Study Breed Samples Number Total length (Mb) Number Pct. of number (%) Total length (Mb) Pct. of length (%) Sequencing-based study This study 12 12 7,530 88.12 CGH-based studies Wang et al. 2010 [40] 3 10 91 15.72 162 2.15 2.66 3.02 Wang et al. 2012 [9] 3 6 130 3.34 83 1.10 0.92 1.04 Crooijmans et al. 2013 [34] 7 64 1,504 57.44 721 9.58 8.30 9.42 Luo et al. 2013 [39] 4 6 29 1.46 21 0.28 0.42 0.47 Tian et al. 2013 [33] 11 22 308 10.81 166 2.20 2.00 2.27 Abernathy et al. 2014 [41] 2 12 147 4.18 68 0.90 0.63 0.71 SNP-based study (60 K) Jia et al. 2013 [32] 2 746 209 13.55 141 1.87 1.75 1.99 Sequencing-based study Fan et al. 2013 [36] 2 2 415 3.17 96 1.27 0.80 0.90 Total 1,052 13.97 12.45 14.13 Yi et al. BMC Genomics 2014, 15:962 Page 5 of 16 http://www.biomedcentral.com/1471-2164/15/962 evaluating these two copies regions. For all 12 individuals, average of 0.647. BY (0.502), SK (0.435) and WR (0.491) we selected all 5 kb non-overlapping windows on auto- showed lower correlation close to 0.500, and we found the somes and excluded the windows intersecting with pre- mean of all probes log ratio values in the three aforemen- dicted CNVs and gaps, and then estimated their average tioned individuals were 1.05, 0.85 and 1.05 respectively, CN. The average CN and STDEV per individual was which were larger than the values of others that were close 2.077 ± 0.291, varied from 2.041 ± 0.226 in WR to 2.104 ± to zero. 0.299 in RJF, showing low variability within the predicted In addition, we randomly chose 15 predicted CNVRs neutral regions. Further, we validated sequencing-based representing different types and frequencies for qPCR as- CNV predictions by two independent experiment ap- says, and tested all 12 samples for each CNVR. Two dis- proaches including aCGH and qPCR. We performed 11 tinct pairs of primers were designed for each predicted pairwise aCGH experiments using RJF as the reference CNVR (Additional file 6: Table S4). The proportion of and all others as the test samples. Considering that we es- confirmed positive samples (positive predictive value) var- timated the CN of selected individuals with respect to ref- ied from 50 to 100%, with an average of 91.71%. However, erence genome which cannot be used for the aCGH some negative samples were also confirmed to contain reference sample, we calculated the predicted log CN ra- CNVs, and the false negative rate varied from 0 to 60%, tios for the 11 test samples against RJF to make the CN with an average of 22.43%. We illustrated the qPCR results values comparable with the aCGH results, which was des- for three confirmed CNVRs of different types (gain, loss ignated as digital aCGH approach [12,38]. We focused on and both) (Additional file 7: Figure S3). the autosomal CNVs to avoid the impact of gender and unassigned linkage groups. For pairwise samples (each of Copy number polymorphic genes the 11 test samples and RJF), there were two types of We obtained 6,086 non-redundant RefSeq gene tran- CNV events, i.e., overlapping and unique segments. For scripts retrieved from the UCSC Genome Browser and the overlapping CNV segments, we first split them into estimated the copy number values of all genes in differ- non-overlapping subsegments. And then we estimated the ent individuals by CNVnator. A total of 2,216 (36.4%) CN of each subsegment and unique segment longer than genes overlapped with 2,214 (25.0%) predicted CNVRs. 1,000 bp for each of the two pairwise samples, and divided Among them, 537 genes were found to be completely the copy number estimates of the test sample by that of covered by CNVRs. The overlapping genes were found RJF and calculated log CN ratios as digital aCGH values. not to be highly duplicated sequences, and the max- Then we compared the digital values with aCGH probe imum copy number estimates was only 12.0. We exam- log ratios which were defined as the average of all probes ined the 25 most variable genes according to the STDEV log ratio values in the corresponding segments. We per- of their copy number estimates in different individuals formed a simple linear regression analysis to assess the (Additional file 8: Table S5), and found that these genes correlation between the two values. The Pearson’scorrel- were mainly involved in immune response and keratin ation coefficient (r) ranged from 0.435 in SK to 0.755 in formation. It should be noted that the keratin gene fam- DX (Figure 2 and Additional file 5: Figure S2), with an ilies were detected to have large CN values and variance. Figure 2 Correlation between digital aCGH and whole genome aCGH among Luxi Game and White Leghorn compared with Red Jungle Fowl (RJF). RJF is selected as the reference sample in each aCGH experiment. Digital aCGH values are estimated using calculated log CN ratios in which CN are estimated for identified CNV segments of two individuals and divided by the corresponding CN of RJF. Whole genome aCGH values are defined as the average of all probes log ratio values in the same segments as the digital aCGH. 2 Yi et al. BMC Genomics 2014, 15:962 Page 6 of 16 http://www.biomedcentral.com/1471-2164/15/962 Two significant CNVRs associated with dermal hyperpig- to perform qPCR experiments using the same two pairs of mentation were located on chr20 at positions 11,217,001 primers listed in Additional file 6: Table S4. Two qPCR re- to 11,272,200 (CNVR7962) and 11,651,801 to 11,822,900 sults demonstrated the copy number estimate of almost (CNVR7968), respectively, which had already been de- each LX was larger than the others (Figure 4), and the scribed in detail in a previous study [30], and the distance average copy number estimates (5.0 and 5.2 for two pairs between the two loci was 379.6 kb. SLMO2 and TUBB1 of primers, respectively) of all LX were significant larger were completely covered by the first region which was than those (2.6 and 2.6) in other individuals using the predicted to be about twice as many copies in DX and SK two-sample t-test (P-value =0.003 and 0.001). Additionally, as in other individuals (Figure 3A and Additional file 9: other identified CNV-gene overlaps could be potentially Figure S4a). The functional gene EDN3 (endothelin 3) is responsible for certain economic traits, as these genes were not archived because the predicted gene is not available involved in lipid metabolism, muscle development and pro- for UCSC RefSeq database. We found that only BY had tein secretion process. For example, our results suggested this CNVR while SK and DX as two typical breeds with higher copy number for the POPDC3 gene in WL (n =4.2) dermal hyperpigmentation did not. So we further checked than in the other 11 genomes (n =2.3) (Figure 3C and the raw results before removing CNVs overlapping with Additional file 9: Figure S4c). Similarly, the WL genome gaps. Two nearly identical CNVs comprising two gaps showed the greatest number of AVR2 copies (n =2.0) on (>100 bp) were found, one at positions 11,111,501 to chrZ compared with the others (n =1.1). Two promising 11,238,600 in DX and the other at positions 11,111,401 genes involving in lipid metabolism, AP2M1 and LBFABP, to 11,238,900 in SK, which were also confirmed by our were identified as the largest copy number (n =3.0 and whole genome aCGH experiments (Figure 3A and 3.2) in meat-type chicken (CS) compared with those in the Additional file 9: Figure S4a). The distance between the others. raw CNVR and the second region (CNVR7968) was 412.9 kb, which perfectly supported the reported results Heatmap visualization [30]. Conversely, the first CNVR in BY (11,217,001 to We performed a hierarchical clustering analysis and gen- 11,272,200) showing normal skin color does not contain erated heatmaps based on Pearson’s correlation coeffi- the EDN3 gene (11,148,025 to 11,160,484), which also pro- cient using the CN values of selected gene loci, in order vides evidence that the EDN3 with copy number poly- to infer the potential relationship of selected genes morphism is the causal mutation resulting in dermal among 12 individuals. The loci encompassing SLMO2 hyperpigmentation. Another previously identified CNVR and TUBB1 in DX and SK were found to be highly du- involving the PRLR (prolactin receptor) gene on chrZ [31] plicated regions and the two individuals were clustered was also detected in our study, and the copy number es- into one group (Figure 5A). Another promising gene, timates of PRLR in WC and WL were twice as many as SOCS2, was also confirmed for the difference in copy in other individuals. The sex-linked K allele containing number between LX and the others (Figure 5B). Mean- two copies of PRLR in females is associated with the late while, WL showed specific expansion in the POPDC3 feathering and used widely for sexing hatchlings. Our locus and was split into a separate clade (Figure 5C). sequencing-based and qPCR results showed that WC and WL should exhibit the late feathering phenotype, Gene content and QTL analysis of CNVRs which were supported by the actual phenotype record. A total of 2,216 RefSeq genes overlapped with putative In addition, we found that some genes related to the CNVRs. Then, we performed gene ontology (GO) and host immune and inflammatory response had CNVR Kyoto Encyclopedia of Genes and Genomes (KEGG) path- overlaps, like CD8A, FZD6, LIMS1, TNFSF13B and some way analysis for these genes. The GO analysis revealed MHC-related genes associated with Marek’s disease 646 GO terms, of which 175 were statistically significant (MD). The SOCS2 involving in the regulation of bone after Benjamini correction (Additional file 10: Table S6). growth and density was predicted to have the largest CN The significant GO terms were mainly involved in positive value in LX (n =6.4), while DX (n =3.0) and TB (n =3.6) regulation of macromolecule metabolic process and gene also possessed the duplicated sequences in this locus expression, plasma membrane, protein localization, en- compared with the neutral regions in other individuals zyme binding, response to oxidative stress and immune (Figure 3B and Additional file 9: Figure S4b). LX repre- system development. The KEGG pathway analysis indi- sents a characteristic breed for cockfighting in which bone cated that these genes were overrepresented in nine path- strength is an essential feature for selection. To validate ways, but none of which was significant after Benjamini the highly duplicated sequence (CNVR410) found only in correction. According to our artificial QTL filtering cri- LX, we selected another 16 individuals, i.e., eight LX (four teria, we identified 595 high-confidence QTLs in total, of males and four females) and other eight females consisting which 300 (50.4%) were found to overlap with 560 (6.3%) of one CS, one DX, one SG, one SK, two TB and two WL, CNVRs (Additional file 11: Table S7). These QTLs were Yi et al. BMC Genomics 2014, 15:962 Page 7 of 16 http://www.biomedcentral.com/1471-2164/15/962 Figure 3 (See legend on next page.) Yi et al. BMC Genomics 2014, 15:962 Page 8 of 16 http://www.biomedcentral.com/1471-2164/15/962 (See figure on previous page.) Figure 3 Visual examination by read depth, whole-genome aCGH and digital aCGH around three loci for five representative chicken genomes. The uppermost gene image is generated with the UCSC Genome Browser (http://genome.ucsc.edu/) using the galGal4 assembly. The track below the gene region is depth of coverage for all five individual genomes. Red indicates regions of excess read depth (> mean +3 × STDEV), whereas gray indicates intermediate read depth (mean +2 × STDEV < × < mean +3 × STDEV), and green indicates normal read depth (mean ± 2 × STDEV). All read depth values based on 1 kb non-overlapping windows are corrected by GC content. Whole-genome aCGH and digital aCGH values are depicted as the red-green histograms and correspond to a gain colored in green (>0.5), a loss colored in red (<−0.5) and normal status colored in gray (−0.5 < × <0.5). (A) Two previously reported CNVs (chr20: 11,111,401-11,238,900 and chr20: 11,651,801-11,822,900) associated with dermal hyperpigmentation. The DX and SK genomes show two additional copies of the two regions compared with RJF, and are also validated by whole-genome aCGH. (B) A higher copy number increase for the SOCS2 locus (chr1: 44,764,280-44,765,955) is predicted in LX than in other individuals. (C) The POPDC3 gene (chr3: 68,255,196-68,259,535) is predicted to be duplicated status only in WL. mainly involved in production and health traits, such as RD method employed in our work has advantages in both growth, body weight, abdominal fat weight, egg number technology platform and genetic diversity compared with and Marek’s disease-related traits. the eight previous reports [9,32-34,36,39-41]. Because a significant fraction of CNVs falls into genomic regions not Discussion well-covered by microarrays, especially for SD regions This study performed genome-wide CNV detection, deter- lacking sufficient probes [16,23], CNV as a major source mined absolute copy number and constructed the first in- of genetic variation is complementary to SNP and could dividualized chicken CNV map. The NGS technology and account for a substantial part of missing heritability [14]. Figure 4 Validation of CNVR410 by qPCR in another 16 chickens. X-axis represents all 16 samples and Y-axis represents normalized ratios (NR) estimated by qPCR. NR around 2 indicates normal status (2 copies), NR around 0 or 1 indicates loss status (0 copies or 1 copy), and NR around 3 or more indicates gain status (3 or more copies). (A) qPCR results confirmed by primer410_2. (B) qPCR results confirmed by primer410_3. Yi et al. BMC Genomics 2014, 15:962 Page 9 of 16 http://www.biomedcentral.com/1471-2164/15/962 Figure 5 (See legend on next page.) Yi et al. BMC Genomics 2014, 15:962 Page 10 of 16 http://www.biomedcentral.com/1471-2164/15/962 (See figure on previous page.) Figure 5 Hierarchical clustered heatmaps of preselected genetic loci for 12 chicken genomes. Every block in the heatmap indicates estimated CN values of 1 kb non-overlapping windows in the preselected region. These heatmaps are generated from hierarchical cluster analysis using Pearson’s correlation coefficient of the CN values. The colors for each bar denote different copy number (CN). (A) DX and SK which are predicted to be doubled within dermal hyperpigmentation loci are clustered together. (B) Upstream and downstream of the SOCS2 locus reveal higher CN values in DX, TB and WC especially LX. (C) WL shows specific expansion in the POPDC3 locus and is split into a separate clade. Most CNV studies to date have been discovery studies ra- larger than losses because chromosomal deletion can ther than association studies, mainly due to the limitations lead to a variety of serious malformations and disorders of CNV resolution and genotyping in each individual [3]. and is subjected to purifying selection [13,47]. In gen- The high-resolution individualized chicken CNV map eral, the length of chromosome is positively correlated based on extensive genetic diversity not only enriches the with the number of CNVRs. The chr16 (a microchromo- current genetic variation database but also encourages the some) was found to have the second densest CNVRs, future development of assays for accurately genotyping possibly owing to the highly variable major histocom- CNVs, enabling systematic exploration about CNV associ- patibility complex (MHC) regions and higher recombin- ation studies similar to SNPs. In future, integration of ation rate [48], which also results in repeatedly finding CNVs with SNPs maybeaneffectiveand promisingway to the same CNVRs on chr16 among different studies. elucidate the causes of complex diseases and traits [16,17]. Quality assessment and experimental validation Genome-wide CNV landscape in the chicken genome It is generally believed that the CN of neutral regions is The number of CNVs and CNVRs in each individual var- between 1.5 and 2.5 [25] and the mean ± 2 × STDEV in ied greatly, and all individuals shared a small number of our results corresponded closely to the hypothesis, them, suggesting that a considerable proportion of CNVs which demonstrates that CNVnator has efficient per- likely segregated among 12 distant breeds [2,34], therefore formance on CNV detection and CN estimation and can a larger population and multiple breeds are crucial to con- generate most reliable results. For CNV quality assess- struct a more complete chicken CNV map. The high per- ment by aCGH, the positive correlation values between centage of unique CNVRs could also be partly explained computational and experimental log CN ratios in our by the high recombination rate in the chicken genome study were higher than the previous results [2], mainly (2.5-21 cM/Mb), because recombination-based mecha- owing to the aCGH platform with higher resolution in nisms such as non-allelic homologous recombination our analysis. The slightly low correlation coefficients in (NAHR) are the major causes leading to CNVs [42]. Simi- BY, SK and WR might disclose certain experimental larly, the high recombination rate may induce more denser noises and biases resulting in misgenotyping in corre- CNVRs in microchromosomes [43]. These unique CNVRs sponding aCGH experiments [16], and particularly may be recent events in evolution and contribute to highly duplicated regions and rare deletions [15,25]. In breed-specific phenotype and performance [44]. Com- addition, the average positive predicted value of the 15 pared with the eight previous chicken CNV studies chosen CNVRs was 91.71%, similar to some previous [9,32-34,36,39-41], far more CNVRs both on average and results in animals [7,33,45], suggesting that most of the in total were found. A total of 6,478 (86.0%) autosomal positive samples detected by sequencing-based method CNVRs with slightly smaller average size (10.6 kb) were arehighlyconsistentwiththe qPCR experiments. We novel, likely due to the higher resolution and sensitivity of also estimated the false negative error rates as it is a NGS method than aCGH and SNP array. These novel common problem in CNV detection [7,49], and the CNVRs enrich significantly the published chicken CNV average percentage of false negative results was 22.43%. database. The low concordance between different studies The discrepancies between NGS results and qPCR val- results from the differences in technical issues, CNV call- idation may be due to the negative impact of potential ing algorithms as well as study populations [45], and can SNPs and small INDELs, which result in the reduced also indicate that numerous CNVs in the chicken genome primer efficiency. are still expected to be discovered. We found both maximum and mean copy number es- Promising candidate genes covered by CNVRs timates of autosomal duplicated sequences in chickens CNV is a significant source of genetic variation account- were less than those in mammalians [2,12], which might ing for disease and phenotypic diversity, owing to the be related to the relatively smaller genome size (only duplication or deletion of covered genes or regulation el- one third of a typical mammalian genome) and the lower ements [4]. Our results showed that 36.4% RefSeq genes repetitive DNA content in the chicken genome [46]. In intersected with 25.0% predicted CNVRs. It is probable addition, both the counts and sizes of gain events were that CNVs, especially deletions, are located preferably in Yi et al. BMC Genomics 2014, 15:962 Page 11 of 16 http://www.biomedcentral.com/1471-2164/15/962 gene-poor regions [13,47], because gene-rich CNVRs are in LX was larger than others. We suspect that the copy more likely to be pathogenic than gene-poor CNVRs number polymorphic locus is ubiquitous in the chicken and these deleterious CNVRs would be removed by genome, but the particularly high gene duplication in LX purifying selection [47,50]. Meanwhile, the maximum may be the result of the genetic effect of long-term artifi- CN of all genes covered by CNVRs was 12.0, suggesting cial selection such as crossing between the individuals again that the chicken genome has lower repetitive DNA with stronger bone. content [46]. It is noted that nine out of the 25 most Additionally, the copy number estimates of POPDC3 variable genes belong to four keratin subfamilies (claw, (popeye domain containing 3) in WL were found to be feather, feather-like and scale). In birds, skin appendages about twice as many as other individuals. The POPDC3 such as claws, scales, beaks and feathers are composed of gene belongs to the Popeye family encoding proteins beta (β) keratins and can prevent water loss and provide a with three potential transmembrane domains with a high barrier between the organism and external environment degree of sequence conservation, and is preferentially [51]. The avian keratin genes are significantly over- expressed in heart and skeletal muscle cells as well as represented with respect to mammals [34,48]. These highly smooth muscle cells [56]. It has been reported that the variable keratin genes suggest the scenario for the evolu- expression of two Popeye family members was upregu- tion of the β-keratin gene family through gene duplication lated in the uterus of pregnant mice [57]. The uterus has and divergence for their adaptive benefits [4,51]. Addition- been thought to be an important organ composed of ally, the four subfamilies of β-keratin genes form a cluster smooth muscle and containing the shell gland in favor on chr25, one of the more GC-rich chromosomes and con- of depositing eggshell [58]. Considering that WL is the taining a relatively larger number of minisatellites [51], most prolific egg laying chicken due to the fact that it which also result in high copy number of genes. has been extensively bred for egg production, the dupli- We validated two well-known causative genes with copy cation of the POPDC3 gene may reveal the important number polymorphism, EDN3 [30] and PRLR [31], in- differences in abilities like myometrium maturation and volved in dermal hyperpigmentation and late feathering, labor, protein secretion and eggshell formation between respectively. In our study, we used hierarchical clustering WL and other breeds. analysis based on CN content to visualize the potential re- Moreover, these enriched GO terms were mainly in- lationship among 12 breeds. For example, the heatmap for volved in cellular regulation and structure, various binding dermal hyperpigmentation grouped DX and SK together, functions as well as stress and immune responses, which and both of which are distributed in the Jiangxi province of are consistent with several previous studies [9,32-34], sug- China, suggesting that DX and SK may have a close evolu- gesting that the copy number variable genes may influence tionary relationship or purposely bred dermal hyperpig- the responses to external stimuli and provide the muta- mentation into different strains. In addition, two reported tional flexibility to adapt rapidly to changing selective copy number variable genes associated with Marek’sdis- pressures due to evolutionary adaption [59]. Most CNVRs ease, namely FZD6 (frizzled family receptor 6) and LIMS1 also spanned some QTL regions, which indicated that (LIM and senescent cell antigen-like domains 1) [39,52,53] these CNVRs may exert significant effects on traits of eco- were also found in our results. nomic interest in chickens. Furthermore, we also found some novel CNV-gene overlaps as potential candidates linked to some important Conclusions traits. For example, the SOCS2 (suppressor of cytokine sig- In this study, we performed genome-wide CNV detection naling 2) is a member of the suppressor of cytokine signal- and estimated the absolute copy number of the corre- ing family, and the related proteins are implicated in the sponding genetic locus based on whole genome sequen- negative regulation of cytokine action through inhibition cing data of 12 chickens abundant in genetic diversity, and of the JAK/STAT pathway (Janus kinase/signal transduc- constructed the highest-resolution individualized chicken ers and activators of transcription) [54]. Dual x-ray ab- CNV map so far. A total of 8,840 CNVRs were identified, sorptiometry (DXA) analysis demonstrated that SOCS2 and most of them were novel variants which could enrich inactivation resulted in reduced trabecular and cortical the current CNV database. The high CNVR confirmation volumetric bone mineral density (BMD) in SOCS2-defi- rates by aCGH and qPCR suggested that sequencing- cient mice [55]. We found that the SOCS2 had higher CN based method was more sensitive and efficient for CNV (n =6.4) in LX than in other individuals, which is particu- discovery. We detected 2,216 RefSeq genes overlapping larly interesting as the LX is known for cockfighting in with CNVRs, including genes involved in well-known phe- which the chickens with higher BMD have advantage over notypes such as dermal hyperpigmentation and late feath- others. The gene expansions were also supported by the ering. In addition, some novel genes like POPDC3 and heatmap. Additional qPCR experiments in 16 other indi- LBFABP covered by CNVRs may play an important role in viduals revealed that the increased copy number of SOCS2 production traits, and the highly duplicated SOCS2 may Yi et al. BMC Genomics 2014, 15:962 Page 12 of 16 http://www.biomedcentral.com/1471-2164/15/962 serve as an excellent candidate for bone mineral density. processing was performed using CNVnator software Our study lays the foundation for comprehensive under- based on RD method as previously described [25]. standing of copy number variation in the chicken genome CNVnator firstly calculated the counts of mapped reads and is beneficial to future association studies between within user specified non-overlapping bins of equal size CNV and important traits of chickens. as the RD signal, and then adjusted the signal in consid- eration of the potential correlation between RD signal Methods and GC content of the underlying genomic sequence. Sample collection and sequencing The mean-shift algorithm was employed to segment the We selected a total of 12 female chickens from different signal with presumably different underlying CN. Then types and genetic sources representing modern chicken CNVs were predicted by applying statistical significance populations, i.e., a Red Jungle Fowl (RJF, the ancestor of tests to the segments. A more detailed description about domestic chickens), seven Chinese indigenous chickens this method could be found in the CNVnator paper [25]. including Beijing You (BY), Dongxiang (DX), Luxi Game We ran CNVnator with a bin size of 100 bp for our data. (LX), Shouguang (SG), Silkie (SK), Tibetan (TB) and CNV calls were filtered using stringent criteria including Wenchang (WC), and four commercial breeds including P-value <0.01 and size >1 kb, and calls with >50% of q0 Cornish (CS), Rhode Island Red (RIR), White Leghorn (zero mapping quality) reads within the CNV regions (WL) and White Plymouth Rock (WR). The whole blood were removed (q0 filter). All CNV calls overlapping with samples were collected from brachial veins by standard gaps in the reference genome were excluded from con- venepuncture along with regular quarantine inspection sideration. For unlocalized and unplaced chromosomes of the experimental station of China Agricultural Uni- (chrN_random and chrun_random in UCSC, chrUn), we versity, and genomic DNA was isolated using the stand- removed them for further analysis due to the shorter ard phenol/chloroform extraction method. Whole length of the chrUn contigs and mapping ambiguity of genome sequencing for all 12 individuals was performed chrUn sequence reads. Meanwhile, we performed geno- on the HiSeq 2000 system (Illumina Inc., San Diego, typing of all 5 kb non-overlapping windows which did CA, USA). Two genomic DNA libraries of 500 bp insert not overlap with putative CNVs and gaps on autosomes. size per individual were constructed and sequenced with In order to compare our results with previous studies, 100 bp paired-end reads, and each library dataset was we converted all autosomal CNVRs from galGal4 to gal- generated with five-fold coverage depth. Library prepar- Gal3 using the UCSC liftOver tool [63]. ation and all Illumina runs were performed as the stand- ard manufacturer’s protocols. Array CGH for assessing genome-wide CNVs We conducted CNV consistency evaluation using two Quality control and sequence alignment similar whole genome tiling arrays based on galGal4 2011 For ensuring high-quality data, we used NGS QC Toolkit build. One of them is the NimbleGen aCGH (Madison, with default parameters to perform quality control of raw WI, USA), a custom-designed 3*1.4 M array containing a sequencing data, mainly by removing low-quality reads total of 1,425,178 50-75mer probes with the mean and and reads containing primer/adaptor contamination [60]. median interval of 734 bp and 700 bp. The other is the All high-quality Illumina sequence reads were aligned Agilent custom-designed 1*1.0 M array (Agilent Technol- against the galGal4 assembly by using the Burrows- ogy Inc., CA, USA), with the mean and median probe spa- Wheeler Aligner (BWA) program [61] with default param- cing of 1,056 bp and 1,050 bp. It should be noted that the eters. The draft genome sequence was retrieved from the average physical distance of the closest SNP probes be- UCSC website (http://hgdownload.soe.ucsc.edu/golden- tween two arrays was 262.6 bp and 95.2% distance inter- Path/galGal4/bigZips/). During the construction of a gen- vals were shorter than 500 bp. Meanwhile we only omic library, Illumina platform was likely to generate analyzed raw aCGH log ratio values instead of processed/ some duplicate reads named ‘PCR and optical duplicates’ normalized data. These cases could ensure reasonable ex- which imposed negative impact on the downstream ana- planation for our results although using different arrays. All lysis. So we first used SAMtools [62] to convert the .sam processing steps like DNA labeling (Cy3 for samples and files of different libraries belonging to the same individual Cy5 for references), array hybridization, data normalization to .bam files and sort and merge them, followed by re- and scanning analysis were performed following standard moval of potential PCR duplicates using Picard (http:// procedure. In each aCGH experiment, we chose the RJF as broadinstitute.github.io/picard/). the same reference sample. CNV detection Quantitative PCR confirmation Following the above filtering steps, the resulting .bam We also performed qPCR confirmation of 15 CNVRs files were utilized for CNV calling and genotyping, post- chosen from the CNVRs detected by CNVnator. Most Yi et al. BMC Genomics 2014, 15:962 Page 13 of 16 http://www.biomedcentral.com/1471-2164/15/962 chosen CNVRs have not been reported in the previous Hierarchical cluster analysis studies and are also adjacent to annotated genes. Two We used the heatmap.2() function of the gplots package distinct pairs of PCR primers were designed to target (http://cran.r-project.org/web/packages/gplots/index.html) each CNVR using Primer5.0 software for the uncer- to generate heatmap figures. We first specified the regions tainty of the CNVR boundaries. Furthermore, the UCSC extending 30 kb on each side of interested genes and used In-Silico PCR tool was used for in silico analysis of the estimated CN values of 1 kb non-overlapping windows primers specificity and sensitivity [63]. The PCCA gene in each individual for post-analysis, mainly considering which was previously identified as a non-CNV locus that the regulatory elements may be included in the up- was chosen as the control region [40]. Quality control stream or downstream of a gene. No reordering of those of all primer sets was evaluated using an 8-point stand- windows representing corresponding chromosome loca- ard curve in duplicate to ensure the similar amplifica- tions in the heatmap was made for the sake of clarity. The tion efficiencies between target and control primers. All Pearson’s correlation coefficient of the CN values was used qPCR experiments were conducted on the ABI Prism as the distance measure of the agglomerative hierarchical 7500 sequence detection system (Applied Biosystems clustering with average linkage, and to generate hierarch- group) using SYBR green chemistry in triplicate reac- ical cluster dendrograms. tions, each with a reaction volume of 15 μl in a 96-well plate. The condition for thermal cycle was as follows: Availability of supporting data 1 cycle of pre-incubation at 50°C for 2 min and 95°C for All raw sequence data has been deposited in NCBI Se- 10 min, 40 cycles of amplification (95°C for 15 s and 60°C quence Read Achieve (SRA) under the Bioproject number (1 - ΔΔCt) for 1 min). We used the formula 2 method to cal- PRJNA232548. The experiment numbers for the 12 chick- culate the relative copy number for each test region. The ens are SRX408161-SRX408172. All aCGH data have been cycle threshold (Ct) value of each test sample was first submitted to the NCBI Gene Expression Omnibus (GEO) normalized against the control region, and then the ΔCt (http://www.ncbi.nlm.nih.gov/geo/) under accession num- value was calculated between the test sample and a prese- ber GSE54119. lected reference sample predicted with normal copy num- ber status by CNVnator. The golden standard of each Additional files diploid CNV was generally considered to have two copies for autosomes or one copy when the locus was on Z Additional file 1: Table S1. Summary of identified CNVs and CNVRs in the 12 chicken genomes. chromosome (chrZ) of a female in chickens. Additional file 2: Figure S1. Individualized chicken CNV map in the chicken genome. The horizontal black lines represent the draft chicken genome (UCSC version galGal4). Tracks under the chromosomes indicate Gene contents and functional annotation corresponding CNV status of all individuals kept in the alphabetical order from top to bottom, for BY, CS, DX, LX, RIR, RJF, SG, SK, TB, WC, WL and WR. The RefSeq gene list was retrieved from the UCSC Merged CNVRs from all individuals are depicted above chromosomes. The RefSeq database [63]. All miRNA genes were excluded colors for each bar denote different copy number (CN) in CNV legend and because the nucleotide sequences were too short to es- different types of CNVRs. The downmost axis shows the chromosome, CNV and CNVR coordinates. Left-hand chromosomes are ordered from left to timate reliable copy number. We analyzed the propor- right, and the right-hands are just reversed. tion of the RefSeq genes overlapping with putative Additional file 3: Table S2. General statistics of the CNVRs on CNVRs and performed CN estimates for all 6,086 non- each chromosome. redundant RefSeq gene transcripts. In addition, to provide Additional file 4: Table S3. Summary of novel or reported CNVRs insight into the functional enrichment of the RefSeq genes on autosomes. covered by CNVRs, we performed Gene Ontology (GO) Additional file 5: Figure S2. Correlation between digital aCGH and whole-genome aCGH among nine individuals compared with Red Jungle functional annotation and Kyoto Encyclopedia of Genes Fowl (RJF). RJF is selected as the reference sample in each aCGH experiment. and Genomes (KEGG) pathway analysis employing the Digital aCGH values are estimated using calculated log CN ratios in which web-accessible program DAVID [64]. Statistical signifi- CN are estimated for identified CNV segments of nine individuals and divided by the corresponding CN of RJF. Whole genome aCGH values are cance was accessed by using a modified Fisher’sexact defined as the average of all probes log ratio values in the same segments test and Benjamini correction for multiple testing (P- as the digital aCGH. value <0.05). We also compared the CNVRs identified Additional file 6: Table S4. Primers information and confirmation in this study with the reported QTLs obtained from the results of the 15 chosen CNVRs by qPCR analysis. chicken QTL database [65]. We focused on the QTLs Additional file 7: Figure S3. Illustrating of qPCR confirmation results for three selected CNVRs of different types. X-axis represents all 12 samples and with confidence interval less than 10 Mb and consid- Y-axis represents normalized ratios (NR) estimated by qPCR. NR around 2 ered those QTLs with overlapped confidence intervals indicates normal status (2 copies), NR around 0 or 1 indicates loss status greater than 50% as the same QTL [45], because the (0 copies or 1 copy), and NR around 3 or more indicates gain status (3 or more copies). (A) Results for a gain status of CNVR3588. (B) Results for a loss QTL confidence intervals were too large to be used effi- status of CNVR6695. (C) Results for a both status of CNVR410. ciently in the post-processing. Yi et al. BMC Genomics 2014, 15:962 Page 14 of 16 http://www.biomedcentral.com/1471-2164/15/962 cattle genomes using next-generation sequencing. Genome Res 2012, Additional file 8: Table S5. The detailed features of RefSeq genes 22(4):778–790. completely or partial overlapped with CNVRs. 3. McCarroll SA, Altshuler DM: Copy-number variation and association Additional file 9: Figure S4. Visual examination by read depth, studies of human disease. Nat Genet 2007, 39(7 Suppl):S37–S42. whole-genome aCGH and digital aCGH around three loci for 12 4. Zhang F, Gu W, Hurles ME, Lupski JR: Copy number variation in human chicken genomes. The uppermost gene image is generated with the health, disease, and evolution. Annu Rev Genomics Hum Genet 2009, UCSC Genome Browser (http://genome.ucsc.edu/) using the galGal4 10:451–481. assembly. The track below the gene region is depth of coverage for all 5. Altshuler DM, Gibbs RA, Peltonen L, Dermitzakis E, Schaffner SF, Yu F, 12 individual genomes. Red indicates regions of excess read depth Bonnen PE, de Bakker PI, Deloukas P, Gabriel SB, Gwilliam R, Hunt S, Inouye M, (> mean +3 × STDEV), whereas gray indicates intermediate read depth Jia X, Palotie A, Parkin M, Whittaker P, Chang K, Hawes A, Lewis LR, Ren Y, (mean +2 × STDEV < × < mean +3 × STDEV), and green indicates Wheeler D, Muzny DM, Barnes C, Darvishi K, Hurles M, Korn JM, Kristiansson K, normal read depth (mean ± 2 × STDEV). All read depth values based Lee C, McCarrol SA, et al: Integrating common and rare genetic variation in on 1 kb non-overlapping windows are corrected by GC content. diverse human populations. Nature 2010, 467(7311):52–58. Whole-genome aCGH and digital aCGH values are depicted as the 6. Yalcin B, Wong K, Agam A, Goodson M, Keane TM, Gan X, Nellaker C, red-green histograms and correspond to a gain colored in green Goodstadt L, Nicod J, Bhomra A, Hernandez-Pliego P, Whitley H, Cleak J, (>0.5), a loss colored in red (<−0.5) and normal status colored in gray Dutton R, Janowitz D, Mott R, Adams DJ, Flint J: Sequence-based (−0.5 < x <0.5). (A) Two previously reported CNVs (chr20: 11,111,401-11,238,900 characterization of structural variation in the mouse genome. and chr20: 11,651,801-11,822,900) associated with dermal hyperpigmentation. Nature 2011, 477(7364):326–329. The DX and SK genomes show two additional copies of the two regions 7. Wang J, Jiang J, Fu W, Jiang L, Ding X, Liu JF, Zhang Q: A genome-wide compared with RJF, and are also validated by whole-genome aCGH. (B) A detection of copy number variations using SNP genotyping arrays in higher copy number increase for the SOCS2 locus (chr1: 44,764,280-44,765,955) swine. BMC Genomics 2012, 13:273. is predictedinLXthaninother individuals. (C)The POPDC3 gene (chr3: 8. Liu GE, Hou Y, Zhu B, Cardone MF, Jiang L, Cellamare A, Mitra A, Alexander 68,255,196-68,259,535) is predicted to be duplicated status only in WL. LJ, Coutinho LL, Dell'Aquila ME, Gasbarre LC, Lacalandra G, Li RW, Additional file 10: Table S6. Functional enrichment of GO and KEGG Matukumalli LK, Nonneman D, Regitano LC, Smith TP, Song J, Sonstegard pathway analysis of RefSeq genes covered by CNVRs. TS, Van Tassell CP, Ventura M, Eichler EE, McDaneld TG, Keele JW: Analysis of copy number variations among diverse cattle breeds. Genome Res Additional file 11: Table S7. The overlap information of QTLs and 2010, 20(5):693–703. CNVRs across the chicken genome. 9. Wang Y, Gu X, Feng C, Song C, Hu X, Li N: A genome-wide survey of copy number variation regions in various chicken breeds by array comparative Abbreviations genomic hybridization method. Anim Genet 2012, 43(3):282–289. CNV: Copy number variation; CNVR: Copy number variation region; 10. Hastings PJ, Ira G, Lupski JR: A microhomology-mediated break-induced aCGH: array comparative genomic hybridization; qPCR: quantitative replication model for the origin of human copy number variation. polymerase chain reaction; NAHR: Non-allelic homologous recombination; PLoS Genet 2009, 5(1):e1000327. NHEJ: Non-homologous end joining; FoSTes: Fork stalling and template 11. Sharp AJ, Locke DP, McGrath SD, Cheng Z, Bailey JA, Vallente RU, Pertz LM, switching; SD: Segmental duplication; GWAS: Genome-wide association Clark RA, Schwartz S, Segraves R, Oseroff VV, Albertson DG, Pinkel D, Eichler studies; RD: Read depth; SNP: Single nucleotide polymorphism; NGS: EE: Segmental duplications and copy-number variation in the human Next-generation sequencing; GO: Gene ontology; KEGG: Kyoto encyclopedia genome. Am J Hum Genet 2005, 77(1):78–88. of genes and genomes; QTL: Quantitative trait loci; MHC: Major 12. Alkan C, Kidd JM, Marques-Bonet T, Aksay G, Antonacci F, Hormozdiari F, histocompatibility complex; MD: Marek’s disease; BMD: Bone mineral density. Kitzman JO, Baker C, Malig M, Mutlu O, Sahinalp SC, Gibbs RA, Eichler EE: Personalized copy number and segmental duplication maps using Competing interests next-generation sequencing. Nat Genet 2009, 41(10):1061–1067. The authors declare that they have no competing interests. 13. Freeman JL, Perry GH, Feuk L, Redon R, McCarroll SA, Altshuler DM, Aburatani H, Jones KW, Tyler-Smith C, Hurles ME, Carter NP, Scherer SW, Lee Authors’ contributions C: Copy number variation: new insights in genome diversity. Genome Res NY and LQ conceived and designed all experiments. GY, LQ and YY 2006, 16(8):949–961. performed bioinformatics and statistical analysis with the help of JL, and 14. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, carried out aCGH and qPCR experiments. GX provided samples. GY and LQ McCarthy MI, Ramos EM, Cardon LR, Chakravarti A, Cho JH, Guttmacher AE, drafted the manuscript. NY revised the paper. All authors read and approved Kong A, Kruglyak L, Mardis E, Rotimi CN, Slatkin M, Valle D, Whittemore AS, the final manuscript. Boehnke M, Clark AG, Eichler EE, Gibson G, Haines JL, Mackay TF, McCarroll SA, Visscher PM: Finding the missing heritability of complex diseases. Nature 2009, 461(7265):747–753. Acknowledgements 15. Conrad DF, Pinto D, Redon R, Feuk L, Gokcumen O, Zhang Y, Aerts J, We thank Dr. Xiquan Zhang for providing some samples. This work was Andrews TD, Barnes C, Campbell P, Fitzgerald T, Hu M, Ihm CH, Kristiansson funded in part by grants of National High Technology Development Plan of K, Macarthur DG, Macdonald JR, Onyiah I, Pang AW, Robson S, Stirrups K, China (2013AA102501), Natural Science Foundation of China (31320103905), Valsesia A, Walter K, Wei J, Tyler-Smith C, Carter NP, Lee C, Scherer SW, Programs for Changjiang Scholars and Innovative Research in University Hurles ME: Origins and functional impact of copy number variation in (IRT1191), and China Agriculture Research Systems (CARS-41). the human genome. Nature 2010, 464(7289):704–712. Received: 4 March 2014 Accepted: 13 October 2014 16. Liu GE, Bickhart DM: Copy number variation in the cattle genome. Published: 7 November 2014 Funct Integr Genomics 2012, 12(4):609–624. 17. Stranger BE, Forrest MS, Dunning M, Ingle CE, Beazley C, Thorne N, Redon R, Bird CP, de Grassi A, Lee C, Tyler-Smith C, Carter N, Scherer SW, Tavare S, References Deloukas P, Hurles ME, Dermitzakis ET: Relative impact of nucleotide and 1. Redon R, Ishikawa S, Fitch KR, Feuk L, Perry GH, Andrews TD, Fiegler H, copy number variation on gene expression phenotypes. Science 2007, Shapero MH, Carson AR, Chen W, Cho EK, Dallaire S, Freeman JL, Gonzalez 315(5813):848–853. JR, Gratacos M, Huang J, Kalaitzopoulos D, Komura D, MacDonald JR, 18. Henrichsen CN, Chaignat E, Reymond A: Copy number variants, diseases Marshall CR, Mei R, Montgomery L, Nishimura K, Okamura K, Shen F, and gene expression. Hum Mol Genet 2009, 18(R1):R1–R8. Somerville MJ, Tchinda J, Valsesia A, Woodwark C, Yang F, et al: Global 19. LaFramboise T: Single nucleotide polymorphism arrays: a decade of variation in copy number in the human genome. Nature 2006, biological, computational and technological advances. Nucleic Acids Res 444(7118):444–454. 2009, 37(13):4181–4193. 2. Bickhart DM, Hou Y, Schroeder SG, Alkan C, Cardone MF, Matukumalli LK, 20. Pinto D, Darvishi K, Shi X, Rajan D, Rigler D, Fitzgerald T, Lionel AC, Song J, Schnabel RD, Ventura M, Taylor JF, Garcia JF, Van Tassell CP, Thiruvahindrapuram B, Macdonald JR, Mills R, Prasad A, Noonan K, Gribble S, Sonstegard TS, Eichler EE, Liu GE: Copy number variation of individual Yi et al. BMC Genomics 2014, 15:962 Page 15 of 16 http://www.biomedcentral.com/1471-2164/15/962 Prigmore E, Donahoe PK, Smith RS, Park JH, Hurles ME, Carter NP, Lee C, 39. Luo J, Yu Y, Mitra A, Chang S, Zhang H, Liu G, Yang N, Song J: Genome- Scherer SW, Feuk L: Comprehensive assessment of array-based platforms wide copy number variant analysis in inbred chickens lines with and calling algorithms for detection of copy number variants. different susceptibility to Marek’s disease. G3 (Bethesda) 2013, Nat Biotechnol 2011, 29(6):512–520. 3(2):217–223. 21. Bentley DR, Balasubramanian S, Swerdlow HP, Smith GP, Milton J, Brown 40. Wang X, Nahashon S, Feaster TK, Bohannon-Stewart A, Adefope N: An initial CG, Hall KP, Evers DJ, Barnes CL, Bignell HR, Boutell JM, Bryant J, Carter RJ, map of chromosomal segmental copy number variations in the chicken. Keira Cheetham R, Cox AJ, Ellis DJ, Flatbush MR, Gormley NA, Humphray SJ, BMC Genomics 2010, 11:351. Irving LJ, Karbelashvili MS, Kirk SM, Li H, Liu X, Maisinger KS, Murray LJ, 41. Abernathy J, Li X, Jia X, Chou W, Lamont SJ, Crooijmans R, Zhou H: Copy Obradovic B, Ost T, Parkinson ML, Pratt MR, et al: Accurate whole human number variation in fayoumi and leghorn chickens analyzed using array genome sequencing using reversible terminator chemistry. Nature 2008, comparative genomic hybridization. Anim Genet 2014, 45(3):400–411. 456(7218):53–59. 42. Munoz-Amatriain M, Eichten SR, Wicker T, Richmond TA, Mascher M, 22. Yoon S, Xuan Z, Makarov V, Ye K, Sebat J: Sensitive and accurate detection Steuernagel B, Scholz U, Ariyadasa R, Spannagl M, Nussbaumer T, Mayer KF, of copy number variants using read depth of coverage. Genome Res Taudien S, Platzer M, Jeddeloh JA, Springer NM, Muehlbauer GJ, Stein N: 2009, 19(9):1586–1592. Distribution, functional impact, and origin mechanisms of copy number 23. Campbell CD, Sampas N, Tsalenko A, Sudmant PH, Kidd JM, Malig M, Vu TH, variation in the barley genome. Genome Biol 2013, 14(6):R58. Vives L, Tsang P, Bruhn L, Eichler EE: Population-genetic properties of 43. Skinner BM, Al Mutery A, Smith D, Volker M, Hojjat N, Raja S, Trim S, Houde differentiated human copy-number polymorphisms. Am J Hum Genet P, Boecklen WJ, Griffin DK: Global patterns of apparent copy number 2011, 88(3):317–332. variation in birds revealed by cross-species comparative genomic 24. Teo SM, Pawitan Y, Ku CS, Chia KS, Salim A: Statistical challenges hybridization. Chromosome Res 2014, 22(1):59–70. associated with detecting copy number variations with next-generation 44. Clop A, Vidal O, Amills M: Copy number variation in the genomes of sequencing. Bioinformatics 2012, 28(21):2711–2718. domestic animals. Anim Genet 2012, 43(5):503–517. 25. Abyzov A, Urban AE, Snyder M, Gerstein M: CNVnator: an approach to 45. Jiang L, Jiang J, Yang J, Liu X, Wang J, Wang H, Ding X, Liu J, Zhang Q: discover, genotype, and characterize typical and atypical CNVs from family Genome-wide detection of copy number variations using high-density and population genome sequencing. Genome Res 2011, 21(6):974–984. SNP genotyping platforms in Holsteins. BMC Genomics 2013, 14:131. 26. Szatkiewicz JP, Wang W, Sullivan PF, Sun W: Improving detection of 46. Burt DW: Chicken genome: current status and future opportunities. copy-number variation by simultaneous bias correction and read-depth Genome Res 2005, 15(12):1692–1698. segmentation. Nucleic Acids Res 2013, 41(3):1519–1532. 47. Conrad DF, Andrews TD, Carter NP, Hurles ME, Pritchard JK: A 27. Norris BJ, Whan VA: A gene duplication affecting expression of the ovine high-resolution survey of deletion polymorphism in the human genome. ASIP gene is responsible for white and black sheep. Genome Res 2008, Nat Genet 2006, 38(1):75–81. 18(8):1282–1293. 48. International Chicken Genome Sequencing Consortium: Sequence and 28. Rosengren Pielberg G, Golovko A, Sundstrom E, Curik I, Lennartsson J, comparative analysis of the chicken genome provide unique Seltenhammer MH, Druml T, Binns M, Fitzsimmons C, Lindgren G, Sandberg K, perspectives on vertebrate evolution. Nature 2004, 432(7018):695–716. Baumung R, Vetterlein M, Stromberg S, Grabherr M, Wade C, Lindblad-Toh K, 49. Nicholas TJ, Cheng Z, Ventura M, Mealey K, Eichler EE, Akey JM: The Ponten F, Heldin CH, Solkner J, Andersson L: A cis-acting regulatory mutation genomic architecture of segmental duplications and associated copy causes premature hair graying and susceptibility to melanoma in the horse. number variants in dogs. Genome Res 2009, 19(3):491–499. Nat Genet 2008, 40(8):1004–1009. 50. Lee C, Iafrate AJ, Brothman AR: Copy number variations and clinical 29. Wright D, Boije H, Meadows JR, Bed'hom B, Gourichon D, Vieaud A, Tixier- cytogenetic diagnosis of constitutional disorders. Nat Genet 2007, Boichard M, Rubin CJ, Imsland F, Hallbook F, Andersson L: Copy number 39(7 Suppl):S48–S54. variation in intron 1 of SOX5 causes the Pea-comb phenotype in 51. Greenwold MJ, Sawyer RH: Genomic organization and molecular chickens. PLoS Genet 2009, 5(6):e1000512. phylogenies of the beta (beta) keratin multigene family in the chicken 30. Dorshorst B, Molin AM, Rubin CJ, Johansson AM, Stromstedt L, Pham MH, (Gallus gallus) and zebra finch (Taeniopygia guttata): implications for Chen CF, Hallbook F, Ashwell C, Andersson L: A complex genomic feather evolution. BMC Evol Biol 2010, 10:148. rearrangement involving the endothelin 3 locus causes dermal 52. Wang-Rodriguez J, Dreilinger AD, Alsharabi GM, Rearden A: The signaling hyperpigmentation in the chicken. PLoS Genet 2011, 7(12):e1002412. adapter protein PINCH is up-regulated in the stroma of common 31. Elferink MG, Vallee AA, Jungerius AP, Crooijmans RP, Groenen MA: Partial cancers, notably at invasive edges. Cancer 2002, 95(6):1387–1395. duplication of the PRLR and SPEF2 genes at the late feathering locus in 53. Chen K, Tu Y, Zhang Y, Blair HC, Zhang L, Wu C: PINCH-1 regulates the chicken. BMC Genomics 2008, 9:391. ERK-Bim pathway and contributes to apoptosis resistance in cancer cells. 32. Jia X, Chen S, Zhou H, Li D, Liu W, Yang N: Copy number variations J Biol Chem 2008, 283(5):2508–2517. identified in the chicken using a 60 K SNP BeadChip. Anim Genet 2013, 54. Metcalf D, Greenhalgh CJ, Viney E, Willson TA, Starr R, Nicola NA, Hilton DJ, 44(3):276–284. Alexander WS: Gigantism in mice lacking suppressor of cytokine 33. Tian M, Wang Y, Gu X, Feng C, Fang S, Hu X, Li N: Copy number variants in signalling-2. Nature 2000, 405(6790):1069–1073. locally raised Chinese chicken genomes determined using array 55. Lorentzon M, Greenhalgh CJ, Mohan S, Alexander WS, Ohlsson C: Reduced comparative genomic hybridization. BMC Genomics 2013, 14(1):262. bone mineral density in SOCS-2-deficient mice. Pediatr Res 2005, 34. Crooijmans RP, Fife MS, Fitzgerald TW, Strickland S, Cheng HH, Kaiser P, 57(2):223–226. Redon R, Groenen MA: Large scale variation in DNA copy number in 56. Brand T: The Popeye domain-containing gene family. Cell Biochem Biophys chicken breeds. BMC Genomics 2013, 14:398. 2005, 43(1):95–103. 35. Griffin DK, Robertson LB, Tempest HG, Vignal A, Fillon V, Crooijmans RP, 57. Andree B, Hillemann T, Kessler-Icekson G, Schmitt-John T, Jockusch H, Groenen MA, Deryusheva S, Gaginskaya E, Carre W, Waddington D, Talbot R, Arnold HH, Brand T: Isolation and characterization of the novel popeye Volker M, Masabanda JS, Burt DW: Whole genome comparative studies gene family expressed in skeletal muscle and heart. Dev Biol 2000, between chicken and turkey and their implications for avian genome 223(2):371–382. evolution. BMC Genomics 2008, 9:168. 58. Hincke MT, Nys Y, Gautron J, Mann K, Rodriguez-Navarro AB, McKee MD: 36. Fan WL, Ng CS, Chen CF, Lu MY, Chen YH, Liu CJ, Wu SM, Chen CK, Chen JJ, The eggshell: structure, composition and mineralization. Front Biosci Mao CT, Lai YT, Lo WS, Chang WH, Li WH: Genome-wide patterns of (Landmark Ed) 2012, 17:1266–1280. genetic variation in two domestic chickens. Genome Biol Evol 2013, 59. Gokcumen O, Babb PL, Iskow RC, Zhu Q, Shi X, Mills RE, Ionita-Laza I, Vallender EJ, 5(7):1376–1392. Clark AG, Johnson WE, Lee C: Refinement of primate copy number variation 37. Qu L, Li X, Xu G, Chen K, Yang H, Zhang L, Wu G, Hou Z, Yang N: hotspots identifies candidate genomic regions evolving under positive Evaluation of genetic diversity in Chinese indigenous chicken breeds selection. Genome Biol 2011, 12(5):R52. using microsatellite markers. Sci China C Life Sci 2006, 49(4):332–341. 60. Patel RK, Jain M: NGS QC Toolkit: a toolkit for quality control of next 38. Sudmant PH, Kitzman JO, Antonacci F, Alkan C, Malig M, Tsalenko A, generation sequencing data. PLoS One 2012, 7(2):e30619. Sampas N, Bruhn L, Shendure J, Eichler EE: Diversity of human copy 61. Li H, Durbin R: Fast and accurate short read alignment with burrows- number variation and multicopy genes. Science 2010, 330(6004):641–646. wheeler transform. Bioinformatics 2009, 25(14):1754–1760. Yi et al. BMC Genomics 2014, 15:962 Page 16 of 16 http://www.biomedcentral.com/1471-2164/15/962 62. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The sequence alignment/map format and SAMtools. Bioinformatics 2009, 25(16):2078–2079. 63. Karolchik D, Barber GP, Casper J, Clawson H, Cline MS, Diekhans M, Dreszer TR, Fujita PA, Guruvadoo L, Haeussler M, Harte RA, Heitner S, Hinrichs AS, Learned K, Lee BT, Li CH, Raney BJ, Rhead B, Rosenbloom KR, Sloan CA, Speir ML, Zweig AS, Haussler D, Kuhn RM, Kent WJ: The UCSC Genome Browser database: 2014 update. Nucleic Acids Res 2014, 42(Database issue):D764-770. 64. Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009, 4(1):44–57. 65. Hu ZL, Park CA, Wu XL, Reecy JM: Animal QTLdb: an improved database tool for livestock animal QTL/association data dissemination in the post-genome era. Nucleic Acids Res 2013, 41(Database issue):D871–D879. doi:10.1186/1471-2164-15-962 Cite this article as: Yi et al.: Genome-wide patterns of copy number variation in the diversified chicken genomes using next-generation sequencing. BMC Genomics 2014 15:962. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit

Journal

BMC Genomics – Springer Journals

Published: Dec 1, 2014

Keywords: life sciences, general; microarrays; proteomics; animal genetics and genomics; microbial genetics and genomics; plant genetics and genomics

There are no references for this article.