BACKGROUND
in dairy cattle, milk yield and composition are the most important economic traits.
compared with traditional dairy cattle breeding programs, dna-based marker-assisted selection has obvious advantages to shorten generation interval and enhance the accuracy of selection.
there are several main strategies of identifying the genes with large genetic effects on milk production traits, including marker-qtl linkage analysis , candidate gene approach, genome-wide association analysis  and next-generation sequencing .
many studies have been performed to investigate milk production traits in dairy cattle .
copy number variations  are dna segments that are 1 kb or larger and present at variable copy number in comparison with a reference genome  <cit> .
cnvs are widely distributed in the genome  <cit> .
as a complementary genetic variant to single nucleotide polymorphisms , cnvs have attracted increasing attention in recent years.
compared with the snp, cnv can affect a larger portion of the genome and cause effects, like changing gene structure and dosage, altering gene regulation and exposing recessive alleles  <cit> .
cnvs are also associated with various diseases  <cit>  and may contribute to a fraction of the missing heritability  <cit> .
along with the development of large-scale cnv studies in human  <cit> , substantial progress has been made in the cnv identification in domestic animals, including cattle , dog , sheep  <cit> , goat  <cit> , chicken  <cit>  and pig .
so far, there are four mechanisms known to allow formation of cnv, i.e., non-allelic homologous recombination , non-homologous end joining , fork stalling and template switching  and retrotransposition  <cit> .
in addition, previous studies also suggested that segmental duplications  are one of the catalysts and hotspots for cnv formation  <cit> .
traditionally, there were two array-based methods for cnv discovery, array comparative genomic hybridization  and snp arrays  <cit> .
although they promoted the progress of cnv studies, these two array-based methods have limitations  <cit> .
they cannot detect small cnvs  <cit> .
in addition, it is also a great challenge for microarrays to detect cnvs in the sd regions due to insufficient coverage  <cit>  although customized chips can be designed to cover sds.
recently, the advent of next-generation sequencing  technology has sped up the study of cnv  <cit> .
four basic strategies have been applied for detecting cnvs with ngs data in the  <dig> genomes project pilot studies  <cit> .
read pair  method  <cit>  analyzes discordant mapping pairs of clone reads or high-throughput sequencing fragments whose distances are different from the normal average insert size.
read-depth   analysis detects cnvs based on the read depth-of-coverage, i.e., the density of aligned reads along the chromosomes.
a random distribution  is assumed first in this method.
based on the depth-of-coverage, cnvs are detected with duplication regions showing high coverage, while deletion regions show low coverage.
split-read   <cit>  analysis can evaluate gapped sequence alignments for cnv detection.
this method first splits one read into multiple fragments randomly.
then the first and last fragment aligned along the reference genome respectively.
according to whether the fragments align or not, and the locations and directions if aligned, cnvs can be detected.
the mechanism of sr is similar to rp to some extent.
sequence assembly  method  <cit>  could discover all kinds of genetic variations theoretically because of its fine-scale working.
the direct assembly of short reads without reference genome is called de novo assembly and the general strategy is to reconstruct dna fragments, i.e., contigs, based on assembling overlapping reads firstly.
then by comparing the assembled fragments to the reference genome, the abnormal genomic regions with discordant copy number  can be identified.
additionally, as-based methods can also use a reference genome to improve the computational efficiency and contig quality.
rd methods applied in the  <dig> genomes project data have been shown to predict accurate copy number values due to its capability of high-resolution cnv calls  <cit> .
there have been several approaches based on rd, such as maq  <cit> , segseq  <cit> , mrfast  <cit>  and cnvnator  <cit> .
cnvnator can overcome some disadvantages, including unique regions of the genome  <cit> , poor breakpoint resolution  <cit> , and detect different sizes of cnvs, from a few hundred bases to megabases in the whole genome.
for cnvs accessible by rd described abyzov et al.
<cit> , cnvnator has high sensitivity , low false-discovery rate , high genotyping accuracy , and high resolution in breakpoint discovery.
in addition, they estimated that at least 11% of all cnv loci involve complex, multi-allelic events, a considerably higher estimate than reported earlier  <cit> .
for the cnv detection in the cattle genome, there have been several studies reported using such methods, including cgh  <cit> , bovinesnp <dig> beadchip  <cit> , bovinehd snp beadchip  <cit>  and ngs .
in this study, the objective was to identify candidate genes for milk protein and fat traits of dairy cattle through cnv detection based on ngs data of specific holstein bulls that have extremely high and low estimated breeding values  for milk protein and fat percentages.
methods
animals and re-sequencing
eight proven holstein bulls with high reliabilities  of estimated breeding values  for milk protein percentage  and fat percentage , born between  <dig> and  <dig>  were selected from the beijing dairy cattle center  according to their ebvs for pp and fp.
ebvs were calculated based on a multiple trait random regression test-day model using the software runge by the dairy data center of china .
the bulls were from two half sib families and two full sib families with two bulls in each family.
the two bulls in each group showed extremely high and low ebv for milk pp and fp, respectively.
the detailed information of the  <dig> bulls is present in table  <dig> table  <dig> the estimated breeding values and family information about  <dig> holstein bulls
re-sequencing, data filter and sequence alignment
genomic dna of each bull was extracted from frozen sperms by a standard phenol-chloroform method  <cit> .
dna degradation and contamination were monitored on 1% agarose gels and the concentration and purity were assessed on nanodrop  <dig> ; the high-quality dnas were then used for library construction.
two paired-end libraries were constructed for each individual, the read length was 2 × 100 bp, and whole genome sequencing was performed using illumina hiseq <dig> instruments .
all processes were performed according to the standard manufacturer’s protocols.
in order to get high-quality data, we removed low-quality reads and those containing primer/adaptor contamination which existing in the raw sequencing data by utilizing ngs qc toolkit with default parameters   <cit> .
after data filtering, we used the burrows-wheeler aligner  program  <cit>  with default parameters  to perform sequence alignment based on the umd <dig>  genome assembly which was retrieved from the ucsc website .
to save run time during the downstream analysis, we converted the sam files to bam files and then sorted and merged them by samtools  <cit> .
detection of cnv
cnvnator was run on merged bam files with a bin size of 200 bp following the authors’ recommendations  <cit> .
after calling, quality control was performed on the raw cnvs for each bull.
the filtering criteria included p-value < <dig>  , size >1 kb, and q0 <  <dig> .
p-value < <dig>  means that the region between two calls is not a same cnv and q <dig> means fraction of mapped reads with zero quality.
in addition, the cnvs that overlapped with gaps or unplaced chromosomes  were removed.
statistical analysis
according to the ebvs for pp and fp, the  <dig> holstein bulls were divided into  <dig> groups, high-group and low-group, and the differential cnvs between the high and low groups were obtained as the following steps.
here, a differential cnvr describes a cnvr that was segregating within the two populations.
firstly, as for any two cnvs from any two individuals of  <dig> bulls within each group, they were considered to be common cnvs if they have >30% reciprocal overlap, then we obtained the common cnv regions  by merging the common cnvs across the four individuals in either the high or low groups, respectively.
secondly, after getting the common cnvrs in each group, differential cnvrs were identified between the high and low groups of bulls with extremely high/low pp and fp.
to compare our results with previous studies, we used the ucsc liftover tool  <cit>  to convert the coordinates of cnvrs between umd <dig>  and btau <dig> .
quantitative pcr validation
quantitative pcr  was used to validate cnvrs detected by cnvnator.
a total of  <dig> cnvrs was randomly chosen.
for each cnvr, we firstly determined the best primers after designing multiple pairs of primers because of the uncertainty of the cnvr boundaries using primer <dig> webtool .
to ensure the amplification efficiencies of all pairs of primers, a serial diluted genomic dna sample from a common cattle was used as template for creating a standard curve of each pair of primer.
the basic transcription factor  <dig>  gene was chosen as the control with the assumption that there were two copies of dna segment in this region  <cit> .
with a total volume of 15 μl reagents in a 96-well plate, qpcr was conducted using sybr green chemistry in triplicate reactions on lightcycler®  <dig>  roche.
the condition for thermal cycle was as follows: 5 min at 95 °c followed by 45 cycles at 95 °c for 10 s, 60 °c for 10 s and 72 °c for 15 s. the 2-ΔΔct method was used to calculate the relative copy number for each test region.
first, we obtained the average ct value of three replications of each sample and normalized against the control gene.
then we calculated the Δct value between the test sample and reference sample detected with normal status  by cnvnator.
finally, a value around  <dig> or above was considered as gain and a value around  <dig> or below was considered as loss.
gene contents and functional annotation
using the biomart database, the genes within the detected cnvrs were retrieved based on umd <dig>  sequence assembly .
ensembl genes overlapping with cnvrs completely or partially were considered as copy number variable and selected for further analysis.
to provide insight into the functional enrichment of genes picked out above, we carried out annotation analysis, including go  and kegg , using kobas  <dig>   <cit> , which annotates an input set of genes with putative pathways and disease relationships based on mapping to genes with known annotation.
kobas  <dig>  accepts id and cross-species sequence similarity mapping and then performs statistical tests to identify statistically significantly enriched pathways and diseases.
kobas  <dig>  incorporates knowledge across  <dig> species from  <dig> pathway databases  and  <dig> human disease databases .
all annotated ensembl genes are used as background.
in addition, we compared cnvrs with the reported cattle qtls for milk pp and fp traits in the animal qtl database  <cit> .
RESULTS
sequencing data set statistics and cnv discovery
with illumina paired-end sequencing technology, we obtained ngs data from the  <dig> holstein bulls .
after we mapped them on the umd <dig>  bovine genome assembly and excluded potential pcr duplicates, the depth of coverage for each individual varied from  <dig> ×  to  <dig> × .
as shown previously, a 4x coverage is sufficient for cnv detection using a rd method  <cit> .
with cnvnator, cnvs were detected for  <dig> individuals.
after quality control, the number of duplication ranged from  <dig>  to  <dig> , and the number of deletion varied from  <dig>  to  <dig>  .
in order to determine how many cnvrs were detected from all  <dig> bulls, all the cnvs were merged if overlaps were 1 bp or greater, and a total of  <dig> cnvrs were obtained.
the detailed information about cnvs is shown in additional file 1: table s <dig>  from fig.
<dig>  we can see the cnv landscapes roughly.
although chromosome  <dig> was the longest, the number of cnvs it contained was not the most in any individual.
chromosome x occupied the most cnvs and simultaneously the largest cnvs.table  <dig> summary statistics of sequencing data and cnv of  <dig> holstein bulls
fig.
<dig> individualized cattle cnv map.
the tracks under every black bar represent the cnvs for sample  <dig> to sample  <dig> .
the colors for each bar in the animal data se tracks represent the average estimated cn for each cnv as shown in the legend
to confirm the cnvs detected by cnvnator, qpcr based on the relative comparative threshold cycle  method was performed to validate  <dig> randomly chosen predicted cnvrs representing different types of duplication, deletion and both, on the same  <dig> samples for whole genome sequencing .
it was shown that the validate rates of the  <dig> samples varied from  <dig> % to 90% with an average of  <dig> %.
identification of differential cnvrs between high and low groups
according to the experimental design and filtering standards, we first screened common cnvrs shared by the high and low groups.
then these common cnvrs were excluded from whole cnvrs of high and low group respectively, and  <dig> and  <dig> cnvrs as group-specific in high and low group were remained.
finally, a total of  <dig> differential cnvrs were obtained after merging two group-specific cnvrs if overlaps were 1 bp or greater, covering chromosomes 1-x , which amounted to  <dig>  mb of the cattle genome .
the length of cnvrs varied from  <dig>  kb to  <dig>  kb with an average of  <dig>  kb and a median of  <dig>  kb  and  <dig> % of all cnvrs had sizes ranging from 5 kb to 10 kb .
the cnvrs were divided into  <dig> categories, i.e.
<dig> deletions,  <dig> duplications and  <dig> both events .
in terms of count and length, deletion and duplication cnvrs were almost similar .fig.
<dig> different cnvrs between high and low group.
based on umd <dig> ,  <dig> cnvrs derived from difference between high and low group were shown in blue , red , and green
fig.
<dig> the length and frequency distribution of differential cnvrs
fig.
<dig> types of differential cnvrs
gene contents of differential and group-specific cnvrs
differential cnvrs
utilizing biomart in the ensembl database , we obtained the ids for the genes that are located within or overlapped with the detected cnvrs.
as a result, a total of  <dig> genes were identified, including  <dig> protein-coding genes,  <dig> mirna genes,  <dig> snrna,  <dig> pseudogenes,  <dig> rrna and  <dig> snorna  and  <dig> % of the cnvrs encompass  <dig> or more genes.
to know about the biological functions of these genes, go and kegg pathway analysis were performed with kobas.
we found that there were  <dig> significant go terms and  <dig> significant kegg pathways.
go terms related to protein and lipid metabolism were enriched , such as long-chain fatty acid binding, protein glycosylation, asymmetric protein localization, glycoprotein biosynthetic process, protein serine/threonine kinase activator activity and negative regulation of protein acetylation.
also, the enriched kegg pathways included several well-known protein and lipid metabolisms pathways  such as insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathway .
group-specific cnvrs
furthermore, we obtained  <dig> and  <dig> genes based on the  <dig> and  <dig> cnvrs across  <dig> individuals in the high and low groups, respectively.
in high group, there were  <dig> significant go terms, including lipid metabolism and glucose metabolic processes, and  <dig> significant kegg pathways, a well-known lipid metabolisms pathways   .
in low group,  <dig> significant go terms, i.e., dopamine biosynthetic process and insulin receptor binding, and  <dig> significant kegg pathway  were enriched  .
quantitative traits locus overlapped with differential and group-specific cnvrs
differential cnvrs
we compared the detected differential cnvrs between high and low groups with the previously reported cattle qtl regions for milk production traits  in order to further study the potential genetic effects of these cnvrs.
finally,  <dig> qtls for protein yield, protein percentage, fat yield and fat percentage were found to be overlapped with  <dig> cnvrs , implying the functional genes within these cnvrs are likely candidates for milk protein and fat traits.
group-specific cnvrs
in addition, we compared  <dig> cnvrs in the high group with cattle qtl regions for same traits as above.
totally,  <dig> qtls were obtained which overlap with  <dig> cnvrs .
similarly, we compared  <dig> cnvrs with corresponding traits in the low group and  <dig> qtls overlapped with  <dig> cnvrs were found .
discussion
in this study, we detected genome-wide cnvs of  <dig> holstein bulls with extremely high and low ebvs for pp and fp based on ngs using cnvnator.
we obtained 1834 ~  <dig> cnvs with an average of  <dig>  per bull.
compared with the previous methods based on snp chip and acgh of detecting cnv, ngs has many advantages in terms of both cnv numbers and sizes because the sequencing approach overcomes the sensitivity limits in the previous methods, and can more precisely identify cnv boundaries  <cit> .
with the ongoing developments and cost decreases in ngs, the sequencing approaches has become more and more popular for cnv detection.
due to the fact that it was not designed for cnv detection specifically and incomplete coverage of the whole genome, snp chip was restricted in the application of cnv detection.
based on the observation of cnv distribution, they were enriched in centromeric and the subtelomeric which is in agreement with the distribution of the cattle sd regions as reported before  <cit> .
the number of cnvs  identified in this study was more than the reports based on ngs data by bichkhart et al.
in angus, holstein, hereford and nelore cattle  <cit> , stothard et al.
in holstein and black angus  <cit> , zhan et al.
in holstein  <cit> , boussaha et al.
in holstein, montbéliarde and normande  <cit>  and ben et al.
in holstein  <cit> .
in addition, jiang et al.
detected cnvs based on illumina bovinesnp <dig>  and bovinehd chips  data in chinese holstein population  <cit> , which were less than what we detected in this study.
after converting  <dig> cnvrs to corresponding results based on btau <dig>  using the ucsc liftover tool with 50% of bases that must remap,  <dig> cnvrs of which were successfully converted amounting to about  <dig>  mb.
we found that ~80% of the  <dig> cnvrs overlapped with those reported by previous investigations  by 1 bp or greater , and the largest overlap was ~ <dig>  mb of the reported by bickhart et al.
<cit> .
as for the above inconsistencies, there are likely due to different detection methods and different samples.
bickhart et al.
used mrfast/mrsfast and wssd  <cit> , and both zhan et al.
<cit>  and stothard et al.
<cit>  performed cnv-seq, and boussaha et al.
<cit>  used gatk, pindel, and ben et al.
<cit>  performed control-free.
while in this study, we used cnvnator.
in addition, different cattle breed with specific genetic background may induce the inconsistencies of number of cnvs and cnvrs among various studies as well.
in this study, the qpcr validation rates of the detected cnvs was  <dig> % to 90%, which was similar to those reported by bickhart et al.
<cit> , zhan et al.
<cit> , stothard et al.
<cit>  and yi et al.
<cit> , showing the high accuracy of ngs-based cnv detection.
the relatively lower validation rate in this study may be due to the following reasons:  false positive in cnv calling even if cnvnator has a low false-discovery rate   <cit> ,  primers in qpcr experiment were not the best although we tried multiple primers.
as we know, there may be potential snps and small indels in the genome, and the negative impact of these potential variants could result in the reduced primer efficiency.fig.
<dig> comparison between  <dig> cnvrs in this study and the other cattle cnvr datasets on btau <dig>
genome wide cnv detection is also a strategy to identify the potential key genes for the traits of interest by mining the genes within the cnvrs in a specific experimental design.
hence, the different cnvrs between the high and low groups of holstein bulls with extremely high and low ebvs for pp and fp were used for candidate gene identification for milk protein and milk fat.
we determined a total of  <dig> differential cnvrs between the high and low groups, and further found that  <dig> functional genes were located within these cnvrs.
function analysis showed that the  <dig> genes were enriched in  <dig> significant go terms and  <dig> significant kegg pathways.
especially, the processes of long-chain fatty acid binding, protein glycosylation, asymmetric protein localization, glycoprotein biosynthetic process, protein serine/threonine kinase activator activity and negative regulation of protein acetylation were related to protein or lipid metabolism were enriched.
also among the kegg pathways we detected, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathway are well-known pathways for protein and lipid metabolism  <cit>  and  <dig> genes involved in the biological process such as cell apoptosis, protein modification, conversion of amino acid and metabolism of fatty acids were included.
these were ins, igf <dig>  foxo <dig>  th, scd <dig>  galnt <dig>  galnt <dig>  art <dig>  snca and wnt7a.
the bovine insulin  gene is close to the igf <dig> and th genes.
insulin binding to the insulin receptor  exerts biological function that maintaining the blood glucose concentration through multiple signaling pathways, such as ampk, insulin, mtor and pi3k-akt signaling pathways, which play critical roles in milk fat and protein synthesis in dairy cattle  <cit> .
the igf <dig> gene was related to breast epithelial and stromal cell proliferation in human  <cit> , and over-expression of igf <dig> increased breast cancer development  <cit> , thus, it was implied that igf <dig> may play an important role in maintaining bovine mammary gland epithelial cell function well.
forkhead box o <dig>  was also known as foxo3a, which was considered as a key downstream effector of the well-known signaling pathways for lipid and protein metabolism, i.e.
pi3k-akt, mapk and jak-stat  <cit> .
tyrosine hydroxylase  is the rate limiting enzyme for converting tyrosine to dopamine which was a crucial regulator of prolactin   <cit> .
prl is essential for mammary gland involution and lactation  <cit> .
the stearoyl-coa desaturase  <dig>  gene is located within a known qtl region for milk protein  <cit>  and fat yield  <cit> , but also near to the snps significantly associated with milk fat yield, protein yield, fat percentage and protein percentage identified by a previous gwa study  <cit> .
the adp-ribosyltransferase  <dig>  gene encodes the arginine-specific adp-ribosyltransferase that has impact on cell proliferation and apoptosis etc.
<cit> .
the wnt family member 7a  gene, as a member of wnt gene family, encodes secreted signaling proteins and is related to suppressing human lung cancer progression  <cit> .
the synuclein alpha  gene was found to be associated with parkinson’s and alzheimer’s diseases  <cit> .
among the candidate genes, ins, igf <dig>  foxo <dig>  th, scd <dig> were related with milk composition traits according to previous studies, and identification of them in current study confirmed their potential functions.
as for the remaining genes, there existed more or less indirect association with milk composition traits.
thereby, to gain further insights into the association of the  <dig> candidate genes with milk composition traits, we compared the chromosome positions of the  <dig> genes with the significant snps detected by previous gwas for milk production traits in dairy cattle  <cit> , and found that all genes were located near to multiple significant snps for milk protein and fat traits with  <dig>  mb to  <dig>  mb , suggesting their potential associations with milk compositions.
in the study of xu et al.
<cit> ,  <dig> cnvs were found significantly associated with milk production traits, of which  <dig> cnvs were included in the differential cnvrs identified in this study, i.e., cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig>  cnvr <dig> and cvnr <dig>  within cnvr <dig>  two candidate genes, ins and igf <dig> were enriched.
ben et al.
<cit>  identified two cnvrs associated with milk composition, including one  with milk fat yield and milk protein yield, and another  with milk fat yield, and  <dig> genes were enriched in these two regions, especially the mthfsd gene within the second cnvr belongs to the folate metabolism gene family and plays critical roles in regulating milk protein synthesis  <cit> .
although there was no overlap between these cnvrs and ours, two cnvrs in this study were located near to them with  <dig>  mb and  <dig>  mb, respectively.
the depdc <dig> gene overlapped with cnvr <dig> encodes a protein which was a component of the gap activity toward rags complex and is involved in mtorc <dig> pathway  <cit> .
in addition, we found that  <dig> differential cnvrs detected in this study were overlapped with  <dig> known qtls that have been shown to be associated with protein yield, protein percentage, fat yield and fat percentage in dairy cattle .
eight annotated genes were overlapped with these differential cnvrs .
CONCLUSIONS
in conclusion, based on ngs data of  <dig> holstein bulls with extremely high and low ebvs for milk pp and fp, we identified a total of  <dig> cnvs corresponding to  <dig> cnvrs.
of these,  <dig> differential cnvrs between the high and low groups were obtained.
of note, we further identified  <dig> annotated genes that were located in or overlapped with these differential cnvrs, including  <dig> genes significantly enriched for specific biological functions related to protein and lipid metabolism, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathways.
these genes included ins, igf <dig>  foxo <dig>  th, scd <dig>  galnt <dig>  galnt <dig>  art <dig>  snca and wnt7a, implying their potential association with milk protein and fat traits.
additional files
additional file 1: table s <dig>  summary of identified cnvs and cnvrs in the  <dig> bulls’ genomes.
additional file 2: table s <dig>  primers information and confirmation results of the  <dig> chosen cnvrs by qpcr analysis.
additional file 3: table s <dig>  general statistics of  <dig> differential cnvrs between high and low group based on umd <dig> .
additional file 4: table s <dig>  the detailed features of genes completely or partial overlapped with differential cnvrs.
additional file 5: table s <dig>  functional enrichment of go and kegg pathway analysis of genes covered by differential cnvrs.
additional file 6: table s <dig>  functional enrichment of go and kegg pathway analysis of genes covered by group-specific cnvrs.
additional file 7: table s <dig>  qtls covered by differential cnvrs.
additional file 8: table s <dig>  qtls covered by group-specific cnvrs.
additional file 9: table s <dig>  detailed information on the significant snps identified in previous gwas that are near to the  <dig> candidate genes identified in this study.
abbreviations
acgharray comparative genomic hybridization
asassembly
btf3basic transcription factor 3
bwaburrows-wheeler aligner
cncopy number
cnvcopy number variation
cnvrcnv region
ebvestimated breeding value
fostesfork stalling and template switching
fpfat percentage
gogene ontology
gwasgenome-wide association analysis
keggkyoto encyclopedia of genes and genomes
lalinkage analysis
nahrnon-allelic homologous recombination
ngsnext-generation sequencing
nhejnon-homologous end joining
ppprotein percentage
qpcrquantitative pcr
qtlquantitative traits locus
rdread-depth
rpread pair
sdsegmental duplications
snpsingle nucleotide polymorphisms
srsplit-read

 in dairy cattle, milk yield and composition are the most important economic traits. compared with traditional dairy cattle breeding programs, dna-based marker-assisted selection has obvious advantages to shorten generation interval and enhance the accuracy of selection. there are several main strategies of identifying the genes with large genetic effects on milk production traits, including marker-qtl linkage analysis, candidate gene approach, genome-wide association analysis and next-generation sequencing. many studies have been performed to investigate milk production traits in dairy cattle. copy number variations are dna segments that are 1 kb or larger and present at variable copy number in comparison with a reference genome. cnvs are widely distributed in the genome. as a complementary genetic variant to single nucleotide polymorphisms, cnvs have attracted increasing attention in recent years. compared with the snp, cnv can affect a larger portion of the genome and cause effects, like changing gene structure and dosage, altering gene regulation and exposing recessive alleles. cnvs are also associated with various diseases and may contribute to a fraction of the missing heritability. along with the development of large-scale cnv studies in human, substantial progress has been made in the cnv identification in domestic animals, including cattle, dog, sheep, goat, chicken and pig. so far, there are four mechanisms known to allow formation of cnv, i.e., non-allelic homologous recombination, non-homologous end joining, fork stalling and template switching and retrotransposition. in addition, previous studies also suggested that segmental duplications are one of the catalysts and hotspots for cnv formation. traditionally, there were two array-based methods for cnv discovery, array comparative genomic hybridization and snp arrays. although they promoted the progress of cnv studies, these two array-based methods have limitations. they cannot detect small cnvs. in addition, it is also a great challenge for microarrays to detect cnvs in the sd regions due to insufficient coverage although customized chips can be designed to cover sds. recently, the advent of next-generation sequencing technology has sped up the study of cnv. four basic strategies have been applied for detecting cnvs with ngs data in the genomes project pilot studies. read pair method analyzes discordant mapping pairs of clone reads or high-throughput sequencing fragments whose distances are different from the normal average insert size. read-depth analysis detects cnvs based on the read depth-of-coverage, i.e., the density of aligned reads along the chromosomes. a random distribution is assumed first in this method. based on the depth-of-coverage, cnvs are detected with duplication regions showing high coverage, while deletion regions show low coverage. split-read analysis can evaluate gapped sequence alignments for cnv detection. this method first splits one read into multiple fragments randomly. then the first and last fragment aligned along the reference genome respectively. according to whether the fragments align or not, and the locations and directions if aligned, cnvs can be detected. the mechanism of sr is similar to rp to some extent. sequence assembly method could discover all kinds of genetic variations theoretically because of its fine-scale working. the direct assembly of short reads without reference genome is called de novo assembly and the general strategy is to reconstruct dna fragments, i.e., contigs, based on assembling overlapping reads firstly. then by comparing the assembled fragments to the reference genome, the abnormal genomic regions with discordant copy number can be identified. additionally, as-based methods can also use a reference genome to improve the computational efficiency and contig quality. rd methods applied in the genomes project data have been shown to predict accurate copy number values due to its capability of high-resolution cnv calls. there have been several approaches based on rd, such as maq, segseq, mrfast and cnvnator. cnvnator can overcome some disadvantages, including unique regions of the genome, poor breakpoint resolution, and detect different sizes of cnvs, from a few hundred bases to megabases in the whole genome. for cnvs accessible by rd described abyzov et al., cnvnator has high sensitivity, low false-discovery rate, high genotyping accuracy, and high resolution in breakpoint discovery. in addition, they estimated that at least 11% of all cnv loci involve complex, multi-allelic events, a considerably higher estimate than reported earlier. for the cnv detection in the cattle genome, there have been several studies reported using such methods, including cgh, bovinesnp beadchip, bovinehd snp beadchip and ngs. in this study, the objective was to identify candidate genes for milk protein and fat traits of dairy cattle through cnv detection based on ngs data of specific holstein bulls that have extremely high and low estimated breeding values for milk protein and fat percentages. methods animals and re-sequencing eight proven holstein bulls with high reliabilities of estimated breeding values for milk protein percentage and fat percentage, born between and were selected from the beijing dairy cattle center according to their ebvs for pp and fp. ebvs were calculated based on a multiple trait random regression test-day model using the software runge by the dairy data center of china. the bulls were from two half sib families and two full sib families with two bulls in each family. the two bulls in each group showed extremely high and low ebv for milk pp and fp, respectively. the detailed information of the bulls is present in table  table the estimated breeding values and family information about holstein bulls re-sequencing, data filter and sequence alignment genomic dna of each bull was extracted from frozen sperms by a standard phenol-chloroform method. dna degradation and contamination were monitored on 1% agarose gels and the concentration and purity were assessed on nanodrop; the high-quality dnas were then used for library construction. two paired-end libraries were constructed for each individual, the read length was 2 × 100 bp, and whole genome sequencing was performed using illumina hiseq instruments. all processes were performed according to the standard manufacturer’s protocols. in order to get high-quality data, we removed low-quality reads and those containing primer/adaptor contamination which existing in the raw sequencing data by utilizing ngs qc toolkit with default parameters. after data filtering, we used the burrows-wheeler aligner program with default parameters to perform sequence alignment based on the umd genome assembly which was retrieved from the ucsc website. to save run time during the downstream analysis, we converted the sam files to bam files and then sorted and merged them by samtools. detection of cnv cnvnator was run on merged bam files with a bin size of 200 bp following the authors’ recommendations. after calling, quality control was performed on the raw cnvs for each bull. the filtering criteria included p-value <, size >1 kb, and q0 < . p-value < means that the region between two calls is not a same cnv and q means fraction of mapped reads with zero quality. in addition, the cnvs that overlapped with gaps or unplaced chromosomes were removed. statistical analysis according to the ebvs for pp and fp, the holstein bulls were divided into groups, high-group and low-group, and the differential cnvs between the high and low groups were obtained as the following steps. here, a differential cnvr describes a cnvr that was segregating within the two populations. firstly, as for any two cnvs from any two individuals of bulls within each group, they were considered to be common cnvs if they have >30% reciprocal overlap, then we obtained the common cnv regions by merging the common cnvs across the four individuals in either the high or low groups, respectively. secondly, after getting the common cnvrs in each group, differential cnvrs were identified between the high and low groups of bulls with extremely high/low pp and fp. to compare our results with previous studies, we used the ucsc liftover tool to convert the coordinates of cnvrs between umd and btau. quantitative pcr validation quantitative pcr was used to validate cnvrs detected by cnvnator. a total of cnvrs was randomly chosen. for each cnvr, we firstly determined the best primers after designing multiple pairs of primers because of the uncertainty of the cnvr boundaries using primer webtool. to ensure the amplification efficiencies of all pairs of primers, a serial diluted genomic dna sample from a common cattle was used as template for creating a standard curve of each pair of primer. the basic transcription factor gene was chosen as the control with the assumption that there were two copies of dna segment in this region. with a total volume of 15 μl reagents in a 96-well plate, qpcr was conducted using sybr green chemistry in triplicate reactions on lightcycler® roche. the condition for thermal cycle was as follows: 5 min at 95 °c followed by 45 cycles at 95 °c for 10 s, 60 °c for 10 s and 72 °c for 15 s. the 2-ΔΔct method was used to calculate the relative copy number for each test region. first, we obtained the average ct value of three replications of each sample and normalized against the control gene. then we calculated the Δct value between the test sample and reference sample detected with normal status by cnvnator. finally, a value around or above was considered as gain and a value around or below was considered as loss. gene contents and functional annotation using the biomart database, the genes within the detected cnvrs were retrieved based on umd sequence assembly. ensembl genes overlapping with cnvrs completely or partially were considered as copy number variable and selected for further analysis. to provide insight into the functional enrichment of genes picked out above, we carried out annotation analysis, including go and kegg, using kobas, which annotates an input set of genes with putative pathways and disease relationships based on mapping to genes with known annotation. kobas accepts id and cross-species sequence similarity mapping and then performs statistical tests to identify statistically significantly enriched pathways and diseases. kobas incorporates knowledge across species from pathway databases and human disease databases. all annotated ensembl genes are used as background. in addition, we compared cnvrs with the reported cattle qtls for milk pp and fp traits in the animal qtl database.  sequencing data set statistics and cnv discovery with illumina paired-end sequencing technology, we obtained ngs data from the holstein bulls. after we mapped them on the umd bovine genome assembly and excluded potential pcr duplicates, the depth of coverage for each individual varied from × to ×. as shown previously, a 4x coverage is sufficient for cnv detection using a rd method. with cnvnator, cnvs were detected for individuals. after quality control, the number of duplication ranged from to, and the number of deletion varied from to. in order to determine how many cnvrs were detected from all bulls, all the cnvs were merged if overlaps were 1 bp or greater, and a total of cnvrs were obtained. the detailed information about cnvs is shown in additional file 1: table s from fig.  we can see the cnv landscapes roughly. although chromosome was the longest, the number of cnvs it contained was not the most in any individual. chromosome x occupied the most cnvs and simultaneously the largest cnvs.table summary statistics of sequencing data and cnv of holstein bulls fig. individualized cattle cnv map. the tracks under every black bar represent the cnvs for sample to sample. the colors for each bar in the animal data se tracks represent the average estimated cn for each cnv as shown in the legend to confirm the cnvs detected by cnvnator, qpcr based on the relative comparative threshold cycle method was performed to validate randomly chosen predicted cnvrs representing different types of duplication, deletion and both, on the same samples for whole genome sequencing. it was shown that the validate rates of the samples varied from % to 90% with an average of %. identification of differential cnvrs between high and low groups according to the experimental design and filtering standards, we first screened common cnvrs shared by the high and low groups. then these common cnvrs were excluded from whole cnvrs of high and low group respectively, and and cnvrs as group-specific in high and low group were remained. finally, a total of differential cnvrs were obtained after merging two group-specific cnvrs if overlaps were 1 bp or greater, covering chromosomes 1-x, which amounted to  mb of the cattle genome. the length of cnvrs varied from  kb to  kb with an average of  kb and a median of  kb and % of all cnvrs had sizes ranging from 5 kb to 10 kb. the cnvrs were divided into categories, i.e. deletions, duplications and both events. in terms of count and length, deletion and duplication cnvrs were almost similar.fig. different cnvrs between high and low group. based on umd, cnvrs derived from difference between high and low group were shown in blue, red, and green fig. the length and frequency distribution of differential cnvrs fig. types of differential cnvrs gene contents of differential and group-specific cnvrs differential cnvrs utilizing biomart in the ensembl database, we obtained the ids for the genes that are located within or overlapped with the detected cnvrs. as a result, a total of genes were identified, including protein-coding genes, mirna genes, snrna, pseudogenes, rrna and snorna and % of the cnvrs encompass or more genes. to know about the biological functions of these genes, go and kegg pathway analysis were performed with kobas. we found that there were significant go terms and significant kegg pathways. go terms related to protein and lipid metabolism were enriched, such as long-chain fatty acid binding, protein glycosylation, asymmetric protein localization, glycoprotein biosynthetic process, protein serine/threonine kinase activator activity and negative regulation of protein acetylation. also, the enriched kegg pathways included several well-known protein and lipid metabolisms pathways such as insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathway. group-specific cnvrs furthermore, we obtained and genes based on the and cnvrs across individuals in the high and low groups, respectively. in high group, there were significant go terms, including lipid metabolism and glucose metabolic processes, and significant kegg pathways, a well-known lipid metabolisms pathways. in low group, significant go terms, i.e., dopamine biosynthetic process and insulin receptor binding, and significant kegg pathway were enriched. quantitative traits locus overlapped with differential and group-specific cnvrs differential cnvrs we compared the detected differential cnvrs between high and low groups with the previously reported cattle qtl regions for milk production traits in order to further study the potential genetic effects of these cnvrs. finally, qtls for protein yield, protein percentage, fat yield and fat percentage were found to be overlapped with cnvrs, implying the functional genes within these cnvrs are likely candidates for milk protein and fat traits. group-specific cnvrs in addition, we compared cnvrs in the high group with cattle qtl regions for same traits as above. totally, qtls were obtained which overlap with cnvrs. similarly, we compared cnvrs with corresponding traits in the low group and qtls overlapped with cnvrs were found. discussion in this study, we detected genome-wide cnvs of holstein bulls with extremely high and low ebvs for pp and fp based on ngs using cnvnator. we obtained 1834 ~  cnvs with an average of per bull. compared with the previous methods based on snp chip and acgh of detecting cnv, ngs has many advantages in terms of both cnv numbers and sizes because the sequencing approach overcomes the sensitivity limits in the previous methods, and can more precisely identify cnv boundaries. with the ongoing developments and cost decreases in ngs, the sequencing approaches has become more and more popular for cnv detection. due to the fact that it was not designed for cnv detection specifically and incomplete coverage of the whole genome, snp chip was restricted in the application of cnv detection. based on the observation of cnv distribution, they were enriched in centromeric and the subtelomeric which is in agreement with the distribution of the cattle sd regions as reported before. the number of cnvs identified in this study was more than the reports based on ngs data by bichkhart et al. in angus, holstein, hereford and nelore cattle, stothard et al. in holstein and black angus, zhan et al. in holstein, boussaha et al. in holstein, montbéliarde and normande and ben et al. in holstein. in addition, jiang et al. detected cnvs based on illumina bovinesnp and bovinehd chips data in chinese holstein population, which were less than what we detected in this study. after converting cnvrs to corresponding results based on btau using the ucsc liftover tool with 50% of bases that must remap, cnvrs of which were successfully converted amounting to about  mb. we found that ~80% of the cnvrs overlapped with those reported by previous investigations by 1 bp or greater, and the largest overlap was ~  mb of the reported by bickhart et al.. as for the above inconsistencies, there are likely due to different detection methods and different samples. bickhart et al. used mrfast/mrsfast and wssd, and both zhan et al. and stothard et al. performed cnv-seq, and boussaha et al. used gatk, pindel, and ben et al. performed control-free. while in this study, we used cnvnator. in addition, different cattle breed with specific genetic background may induce the inconsistencies of number of cnvs and cnvrs among various studies as well. in this study, the qpcr validation rates of the detected cnvs was % to 90%, which was similar to those reported by bickhart et al., zhan et al., stothard et al. and yi et al., showing the high accuracy of ngs-based cnv detection. the relatively lower validation rate in this study may be due to the following reasons: false positive in cnv calling even if cnvnator has a low false-discovery rate, primers in qpcr experiment were not the best although we tried multiple primers. as we know, there may be potential snps and small indels in the genome, and the negative impact of these potential variants could result in the reduced primer efficiency.fig. comparison between cnvrs in this study and the other cattle cnvr datasets on btau genome wide cnv detection is also a strategy to identify the potential key genes for the traits of interest by mining the genes within the cnvrs in a specific experimental design. hence, the different cnvrs between the high and low groups of holstein bulls with extremely high and low ebvs for pp and fp were used for candidate gene identification for milk protein and milk fat. we determined a total of differential cnvrs between the high and low groups, and further found that functional genes were located within these cnvrs. function analysis showed that the genes were enriched in significant go terms and significant kegg pathways. especially, the processes of long-chain fatty acid binding, protein glycosylation, asymmetric protein localization, glycoprotein biosynthetic process, protein serine/threonine kinase activator activity and negative regulation of protein acetylation were related to protein or lipid metabolism were enriched. also among the kegg pathways we detected, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathway are well-known pathways for protein and lipid metabolism and genes involved in the biological process such as cell apoptosis, protein modification, conversion of amino acid and metabolism of fatty acids were included. these were ins, igf foxo th, scd galnt galnt art snca and wnt7a. the bovine insulin gene is close to the igf and th genes. insulin binding to the insulin receptor exerts biological function that maintaining the blood glucose concentration through multiple signaling pathways, such as ampk, insulin, mtor and pi3k-akt signaling pathways, which play critical roles in milk fat and protein synthesis in dairy cattle. the igf gene was related to breast epithelial and stromal cell proliferation in human, and over-expression of igf increased breast cancer development, thus, it was implied that igf may play an important role in maintaining bovine mammary gland epithelial cell function well. forkhead box o was also known as foxo3a, which was considered as a key downstream effector of the well-known signaling pathways for lipid and protein metabolism, i.e. pi3k-akt, mapk and jak-stat. tyrosine hydroxylase is the rate limiting enzyme for converting tyrosine to dopamine which was a crucial regulator of prolactin. prl is essential for mammary gland involution and lactation. the stearoyl-coa desaturase gene is located within a known qtl region for milk protein and fat yield, but also near to the snps significantly associated with milk fat yield, protein yield, fat percentage and protein percentage identified by a previous gwa study. the adp-ribosyltransferase gene encodes the arginine-specific adp-ribosyltransferase that has impact on cell proliferation and apoptosis etc.. the wnt family member 7a gene, as a member of wnt gene family, encodes secreted signaling proteins and is related to suppressing human lung cancer progression. the synuclein alpha gene was found to be associated with parkinson’s and alzheimer’s diseases. among the candidate genes, ins, igf foxo th, scd were related with milk composition traits according to previous studies, and identification of them in current study confirmed their potential functions. as for the remaining genes, there existed more or less indirect association with milk composition traits. thereby, to gain further insights into the association of the candidate genes with milk composition traits, we compared the chromosome positions of the genes with the significant snps detected by previous gwas for milk production traits in dairy cattle, and found that all genes were located near to multiple significant snps for milk protein and fat traits with  mb to  mb, suggesting their potential associations with milk compositions. in the study of xu et al., cnvs were found significantly associated with milk production traits, of which cnvs were included in the differential cnvrs identified in this study, i.e., cnvr cnvr cnvr cnvr cnvr cnvr cnvr cnvr cnvr cnvr and cvnr within cnvr two candidate genes, ins and igf were enriched. ben et al. identified two cnvrs associated with milk composition, including one with milk fat yield and milk protein yield, and another with milk fat yield, and genes were enriched in these two regions, especially the mthfsd gene within the second cnvr belongs to the folate metabolism gene family and plays critical roles in regulating milk protein synthesis. although there was no overlap between these cnvrs and ours, two cnvrs in this study were located near to them with  mb and  mb, respectively. the depdc gene overlapped with cnvr encodes a protein which was a component of the gap activity toward rags complex and is involved in mtorc pathway. in addition, we found that differential cnvrs detected in this study were overlapped with known qtls that have been shown to be associated with protein yield, protein percentage, fat yield and fat percentage in dairy cattle. eight annotated genes were overlapped with these differential cnvrs.  in conclusion, based on ngs data of holstein bulls with extremely high and low ebvs for milk pp and fp, we identified a total of cnvs corresponding to cnvrs. of these, differential cnvrs between the high and low groups were obtained. of note, we further identified annotated genes that were located in or overlapped with these differential cnvrs, including genes significantly enriched for specific biological functions related to protein and lipid metabolism, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathways. these genes included ins, igf foxo th, scd galnt galnt art snca and wnt7a, implying their potential association with milk protein and fat traits. additional files additional file 1: table s summary of identified cnvs and cnvrs in the bulls’ genomes. additional file 2: table s primers information and confirmation results of the chosen cnvrs by qpcr analysis. additional file 3: table s general statistics of differential cnvrs between high and low group based on umd. additional file 4: table s the detailed features of genes completely or partial overlapped with differential cnvrs. additional file 5: table s functional enrichment of go and kegg pathway analysis of genes covered by differential cnvrs. additional file 6: table s functional enrichment of go and kegg pathway analysis of genes covered by group-specific cnvrs. additional file 7: table s qtls covered by differential cnvrs. additional file 8: table s qtls covered by group-specific cnvrs. additional file 9: table s detailed information on the significant snps identified in previous gwas that are near to the candidate genes identified in this study. abbreviations acgharray comparative genomic hybridization asassembly btf3basic transcription factor 3 bwaburrows-wheeler aligner cncopy number cnvcopy number variation cnvrcnv region ebvestimated breeding value fostesfork stalling and template switching fpfat percentage gogene ontology gwasgenome-wide association analysis keggkyoto encyclopedia of genes and genomes lalinkage analysis nahrnon-allelic homologous recombination ngsnext-generation sequencing nhejnon-homologous end joining ppprotein percentage qpcrquantitative pcr qtlquantitative traits locus rdread-depth rpread pair sdsegmental duplications snpsingle nucleotide polymorphisms srsplit-read 
these genes included ins, igf foxo th, scd galnt galnt art snca and wnt7a, implying their potential association with milk protein and fat traits. in addition, cnvrs were overlapped with known qtls that are associated with milk protein and fat traits of dairy cattle. in conclusion, based on ngs of holstein bulls with extremely high and low ebvs for milk pp and fp, we identified a total of cnvs, differential cnvrs between groups, and genes, which were suggested as promising candidate genes for milk protein and fat traits. by gene ontology and kegg pathway analyses, we found that genes were significantly enriched for specific biological functions related to protein and lipid metabolism, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathways. herein, we present a genome-wide assessment of cnvs that are potentially associated with milk composition traits in dairy cattle. in this study, cnvs were detected based on whole genome re-sequencing data of eight holstein bulls from four half- and/or full-sib families, with extremely high and low estimated breeding values of milk protein percentage and fat percentage. copy number variations are important and widely distributed in the genome. copy number variationchinese holsteinwhole genome re-sequencinghttp://dx.doi.org/ /501100002855ministry of science and technology of the people's republic of china2013aa 2011bad28b 2014zx08009-053bsun dongxiao http://dx.doi.org/ /501100001809national natural science foundation of china 31472065sun dongxiao http://dx.doi.org/ /501100004826natural science foundation of beijing municipality6152013sun dongxiao beijing municipal bureau of agriculturebaic06-2016sun dongxiao http://dx.doi.org/ /501100002338ministry of education of the people's republic of chinairt1191sun dongxiao issue-copyright-statement© the author 2017 annotation of these differential cnvrs were performed based on the cattle genome reference assembly and totally functional genes were found within the cnvrs. cnv detection opens a new avenue for exploring genes associated with complex traits in humans, animals and plants. among them, differential cnv regions comprising ~  mb of the cattle genome were observed between the high and low groups.
 copy number variations are important and widely distributed in the genome. cnv detection opens a new avenue for exploring genes associated with complex traits in humans, animals and plants. herein, we present a genome-wide assessment of cnvs that are potentially associated with milk composition traits in dairy cattle.  in this study, cnvs were detected based on whole genome re-sequencing data of eight holstein bulls from four half- and/or full-sib families, with extremely high and low estimated breeding values of milk protein percentage and fat percentage. the range of coverage depth per individual was – ×. using cnvnator, we identified a total of cnvs, including duplications and deletions. among them, differential cnv regions comprising ~  mb of the cattle genome were observed between the high and low groups. annotation of these differential cnvrs were performed based on the cattle genome reference assembly and totally functional genes were found within the cnvrs. by gene ontology and kegg pathway analyses, we found that genes were significantly enriched for specific biological functions related to protein and lipid metabolism, insulin/igf pathway-protein kinase b signaling cascade, prolactin signaling pathway and ampk signaling pathways. these genes included ins, igf foxo th, scd galnt galnt art snca and wnt7a, implying their potential association with milk protein and fat traits. in addition, cnvrs were overlapped with known qtls that are associated with milk protein and fat traits of dairy cattle.  in conclusion, based on ngs of holstein bulls with extremely high and low ebvs for milk pp and fp, we identified a total of cnvs, differential cnvrs between groups, and genes, which were suggested as promising candidate genes for milk protein and fat traits. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords copy number variationchinese holsteinwhole genome re-sequencinghttp://dx.doi.org/ /501100002855ministry of science and technology of the people's republic of china2013aa 2011bad28b 2014zx08009-053bsun dongxiao http://dx.doi.org/ /501100001809national natural science foundation of china 31472065sun dongxiao http://dx.doi.org/ /501100004826natural science foundation of beijing municipality6152013sun dongxiao beijing municipal bureau of agriculturebaic06-2016sun dongxiao http://dx.doi.org/ /501100002338ministry of education of the people's republic of chinairt1191sun dongxiao issue-copyright-statement© the author 2017 