Guan et al. Genet Sel Evol (2020) 52:44 https://doi.org/10.1186/s12711-020-00564-4 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access A genome‑wide analysis of copy number variation in Murciano‑Granadina goats Dailu Guan1, Amparo Martínez2, Anna Castelló1,3, Vincenzo Landi2,4, María Gracia Luigi‑Sierra1, Javier Fernández‑Álvarez5, Betlem Cabrera1,3, Juan Vicente Delgado2, Xavier Such6, Jordi Jordana3 and Marcel Amills1,3* Abstract Background: In this work, our aim was to generate a map of the copy number variations (CNV) segregating in a population of Murciano-Granadina goats, the most important dairy breed in Spain, and to ascertain the main biologi‑ cal functions of the genes that map to copy number variable regions. Results: Using a dataset that comprised 1036 Murciano-Granadina goats genotyped with the Goat SNP50 Bead‑ Chip, we were able to detect 4617 and 7750 autosomal CNV with the PennCNV and QuantiSNP software, respec‑ tively. By applying the EnsembleCNV algorithm, these CNV were assembled into 1461 CNV regions (CNVR), of which 486 (33.3% of the total CNVR count) were consistently called by PennCNV and QuantiSNP and used in subsequent analyses. In this set of 486 CNVR, we identified 78 gain, 353 loss and 55 gain/loss events. The total length of all the CNVR (95.69 Mb) represented 3.9% of the goat autosomal genome (2466.19 Mb), whereas their size ranged from 2.0 kb to 11.1 Mb, with an average size of 196.89 kb. Functional annotation of the genes that overlapped with the CNVR revealed an enrichment of pathways related with olfactory transduction (fold-enrichment = 2.33, q-value = 1.61 × 10−10), ABC transporters (fold-enrichment = 5.27, q-value = 4.27 × 10−04) and bile secretion (fold- enrichment = 3.90, q-value = 5.70 × 10−03). Conclusions: A previous study reported that the average number of CNVR per goat breed was ~ 20 (978 CNVR/50 breeds), which is much smaller than the number we found here (486 CNVR). We attribute this difference to the fact that the previous study included multiple caprine breeds that were represented by small to moderate numbers of individuals. Given the low frequencies of CNV (in our study, the average frequency of CNV is 1.44%), such a design would probably underestimate the levels of the diversity of CNV at the within-breed level. We also observed that functions related with sensory perception, metabolism and embryo development are overrepresented in the set of genes that overlapped with CNV, and that these loci often belong to large multigene families with tens, hundreds or thousands of paralogous members, a feature that could favor the occurrence of duplications or deletions by non- allelic homologous recombination. Background Copy number variations (CNV) encompass genomic deletions or duplications, with sizes ranging from 50 base pairs (bp) to several megabases (Mb), and which display polymorphisms (in terms of copy number) among indi- *Correspondence: marcel.amills@uab.cat 1 viduals of a particular species [1–3]. In livestock, a broad Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain array of phenotypes related with, among others, mor- Full list of author information is available at the end of the article phology [4, 5], pigmentation [6–9], sexual development © The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat​iveco​ mmons​.org/licen​ses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creat​iveco​mmons​.org/publi​cdoma​in/ zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Guan et al. Genet Sel Evol (2020) 52:44 Page 2 of 10 [10] and susceptibility to disease [11] is caused by the magnitude of the CNV diversity at the within-breed level. segregation of CNV. Genome scans to detect structural Our goal was to fill this gap by analyzing a population of variations in cattle have revealed that CNV regions 1036 individuals from a single Spanish breed (Murciano- (CNVR) are often enriched in genes that are involved in Granadina), and to investigate the functional roles of immunity [12–15], metabolism [12, 13], embryo develop- genes that map to CNVR and compare these results with ment [12, 15] and sensory perception [13, 14]. There is data obtained in composite goat populations. evidence that the d ­ N/dS ratios of genes that map to tau- rine CNV are generally higher than those of genes that do Methods not overlap with CNV, which indicates that CNV genes Genomic DNA extraction and high‑throughput genotyping probably evolve under reduced selective constraint [13]. Blood samples from 1036 Murciano-Granadina female The analysis of gene networks has also shown that genes goats from 15 farms that are connected through the that co-localize with duplications tend to have fewer use of artificial insemination were collected in EDTA interactions with other genes than loci that do not over- K3 coated vacuum tubes and stored at − 20 °C before lap with CNV, reinforcing the idea that genes mapping processing. Genomic DNA was isolated by a modi- to duplicated regions have fewer essential housekeeping fied salting-out procedure [20]. Four volumes of red cell functions than non-CNV genes, and also have reduced lysis solution (Tris–HCl 10 mmol/L, pH = 6.5; EDTA pleiotropy [13]. 2 mmol/L; Tween 20 1%) were added to 3 mL of whole Although structural chromosomal variations can have blood, and this mixture was centrifuged at 2000×g. Pel- strong effects on gene expression and phenotypic vari- leted cells were resuspended in 3 mL lysis buffer (Tris– ability, technical limitations and the moderate quality HCl 200 mmol/L, pH = 8, EDTA 30 mmol/L, SDS 1%; of genome assemblies have hampered CNV mapping in NaCl 250 mmol/L) plus 100 µL proteinase K (20 mg/ livestock [1]. Until recently, this has been particularly mL). The resulting mixture was incubated at 55 °C for 3 true for goats. In 2010, Fontanesi et al. [16] published the h followed by centrifugation at 2000×g in the presence first caprine CNV map by identifying, with the Bovine of 1 mL of ammonium acetate (10 mol/L). The superna- 385 k aCGH array, 127 CNVR including 86 and 41 copy tant (~ 4 mL) was mixed with 3 mL of isopropanol 96%, loss and gain variants, respectively. Later on, resequenc- which was subsequently centrifuged at 2000×g for 3 min. ing the genome of individuals from several caprine The supernatant was removed and the DNA pellet was breeds made it possible to identify CNV that overlap washed with 3 mL of ethanol 70%. After centrifuging with 13 pigmentation genes and to detect an association at 2000×g for 1 min, the DNA precipitate was dried at between increased ASIP copy number and light pigmen- room temperature and resuspended in 1 mL of TE buffer tation [17]. The first worldwide survey of copy num- (10 mmol/L Tris, pH = 8.0; 1 mmol/L EDTA, pH = 8). ber variation in goats was performed within the Goat High-throughput genotyping of the 1036 Murciano- ADAPTmap Project (http://www.goata​daptm​ap.org), and Granadina DNA samples was carried out with the Goat involved the genome-wide genotyping of 1023 goats from SNP50 BeadChip [21] according to the manufacturer’s 50 breeds [18]. This study resulted in the identification of instructions (Illumina). Signal intensity ratios i.e. log R 978 CNVR among which several overlapped with genes Ratio or LRR (the total probe intensity of a SNP referred that are functionally related with local adaptation such to a canonical set of normal controls [22]), and B allele as coat color, muscle development, metabolic processes, frequencies or BAF (relative quantity of one allele com- and embryonic development [18]. Moreover, the patterns pared to the other one) [22], were exported for each of the diversity of CNV differed according to geographic single nucleotide polymorphism (SNP) with the Genom- origin, which indicates that they have been influenced eStudio software 2.0.4 (Illumina, https​://emea.illum​ina. by population history [18]. In another study on 433 indi- com). Then, SNP coordinates were converted to the latest viduals from 13 East African goat breeds, Nandolo et al. version of the goat reference genome (ARS1) [23]. After [19] detected 325 CNVR. More recently, Henkel et al. [8] filtering out unmapped and non-autosomal SNPs and demonstrated the existence of complex patterns of struc- those with a call rate lower than 98%, a set of 50,551 SNPs tural variation in the regions containing the caprine ASIP remained for CNV mapping. and KIT genes, with potential causal effects on pigmenta- tion. In spite of these efforts, the description of structural Copy number variant calling with PennCNV and QuantiSNP chromosomal variation in goats is still lagging behind Based on their excellent performance in compara- that of other domestic species. Most of the CNV surveys tive studies, we selected two software packages, Pen- in goats have analyzed large populations that represent a nCNV v1.0.5 [24] and QuantiSNP v2 [25], to call CNV mixture of different breeds each with a limited number of in the Murciano-Granadina population [26, 27]. The Pen- individuals [18, 19], thus making it difficult to assess the nCNV software [24] detects CNV by applying the default Guan et al. Genet Sel Evol (2020) 52:44 Page 3 of 10 parameters of the Hidden-Markov model. Population fre- maximum final volume of 15 µL. Each sample was ana- quencies of B alleles were compiled based on the BAF of lyzed in triplicate in order to obtain averaged copy num- each SNP in the population. We used the “–gcmodelfile” ber estimates. Reactions were loaded onto 384-well plates option to adjust “genomic waves” [28]. The number of and run in a QuantStudio 12 K Flex Real-Time PCR Sys- goat chromosomes was set with the “–lastchr 29” instruc- tem instrument (Applied Biosystems). The specificity of tion. The QuantiSNP analysis [25] assumes an objective the PCR reactions was evaluated with a melting curve Bayes hidden-Markov model to improve the accuracy of analysis procedure, and the efficiency (96.2-105.4%) was segmental aneuploidy identification and mapping. This assessed with standard curves. Thus, relative copy num- CNV calling software was run under default parameters ber was inferred with the qbase + software (Biogazelle, by modifying the “–chr 1:29” option. The CNV that were Ghent, Belgium) by using the 2 ­ −ΔΔCt approach [36]. Copy supported by less than three SNPs were removed from number values were calibrated by taking as a reference, the filtered set used here. four samples which, according to Goat SNP50 BeadChip data, had two copies of the investigated genomic loci. Definition and functional annotation of copy number variant regions Results We used the EnsembleCNV algorithm (beta version) [29] Detection of copy number variation to assemble CNVR. All CNV called by PennCNV and/ in Murciano‑Granadina goats or QuantiSNP were combined to generate a set of ini- The initial calling with PennCNV and QuantiSNP yielded tial CNVR by using the heuristic algorithm (threshold 4617 and 7750 autosomal CNV, respectively. By using of minimum overlap = 30%) described in [29]. Subse- the EnsembleCNV tool [29], we assigned these CNV quently, CNVR boundaries were refined by considering into 1461 CNVR with refined boundaries, of which 486 the local correlation structure of the LRR values of the (33.3% of the total CNVR count) were detected simul- SNPs mapping to CNVR [29]. Then, we reassigned the taneously by PennCNV and QuantiSNP. The resulting CNV calls that were initially obtained with PennCNV CNVR included 78 copy gain, 353 copy loss and 55 copy and QuantiSNP to each refined CNVR, so that the final gain/loss variants (Fig. 1, and Table 1) and (see Addi- set of CNVR comprised only those that were simulta- tional file 2: Table S2). The total length of the CNVR neously detected by both callers. The resulting CNVR covered 95.69 Mb (3.9%) of the goat autosomal genome were matched to gene features that are annotated in the (2466.19 Mb), whereas their individual size ranged from National Center for Biotechnology Information (NCBI, 2.0 kb to 11.1 Mb, with an average of 196.9 kb (Fig. 2a and https​://www.ncbi.nlm.nih.gov) by using BEDTools Table 1). Moreover, we found that 72.6% of the CNVR v2.25.0 [30]. In addition, we performed gene ontology showed minimum allele frequencies lower than 0.01, (GO) enrichment and pathway analyses using the DAVID with an average frequency of 1.44% (Fig. 2b). In addition, Bioinformatics Resources 6.8 [31, 32] based on human 10 CNVR with frequencies higher than 10% were distrib- and goat background gene sets. The statistical signifi- uted over seven caprine chromosomes. With a frequency cance was set to a q-value ≤ 0.05. of 41%, CNVR_1229_chr21 was the CNVR with the high- est frequency in the whole dataset (see Additional file 2: Confirmation of copy number variant regions Table S2). By using the BEDTools v2.25.0 program [30], by quantitative real‑time PCR 212 of the CNVR that we detected overlapped with 191 In order to evaluate the rate of false positives in our unique CNVR published by Liu et al. [18] (Fig. 1) and (see experiment, we conducted quantitative real-time PCR Additional file 2: Table S2). The CNVR that were detected (qPCR) experiments to obtain an independent esti- in both studies are referred to as “shared CNVR”, whereas mate of the copy number of four putative CNVR those that were identified in our study only are referred (CNVR_371_chr5, CNVR_506_chr6, CNVR_160_chr2 to as “non-shared CNVR” (Fig. 1). Six of the ten “shared and CNVR_1229_chr21). Primers were designed with the CNVR” with frequencies higher than 0.1 show positional Primer Express software (Applied Biosystems) to amplify concordance with six CNVR detected by Liu et al. [18] specific regions of the ADAMTS20, BST1, NCKAP5 (see Additional file 2: Table S2). and TNFAIP2 genes (see Additional file 1: Table S1). As reference genes, we used the melanocortin 1 recep- Functional annotation of the genes that are located in copy tor (MC1R) and glucagon (GCG​) genes (see Additional number variable regions file 1: Table S1) loci [18, 33–35]. Quantitative PCR reac- Within the CNVR defined in our study, we detected 779 tions contained 7.5 ng genomic DNA, 7.5 µL 2 × SybrSe- protein-coding genes according to the goat reference lect Master mix (Applied Biosystems), 4.5 pmol of each genome annotation (ARS1) [23] from the NCBI data- forward and reverse primer, and ultrapure water to a base (see Additional file 2: Table S2 and Additional file 3: Guan et al. Genet Sel Evol (2020) 52:44 Page 4 of 10 Fig. 1 Genomic distribution of 486 CNVR detected with the PennCNV and QuantiSNP software on the 29 caprine autosomes. Squares, triangles and circles represent copy number gain, loss and gain/loss events, respectively. Red and black colors represent shared and non-shared CNVR, respectively. Shared CNVR are those detected both in our study and in Liu et al. [18], while non-shared CNVR are those identified only in our study Table S3). In a survey of the diversity of CNV in goats with CNVR are olfactory receptors or olfactory receptor- with a worldwide distribution, Liu et al. [18] detected like genes (see Additional file 3: Table S3). Consistently, 1437 copy number variable genes, of which 116 were also the most significantly enriched pathway was “Olfactory identified in our study and are referred to as “shared copy transduction” (q-value = 1.61 × 10−10, Table 2), followed number variable genes” (see Additional file 3: Table S3). by “ABC transporter” (q-value = 4.27 × 10−4, Table 2). A Among the “shared copy number variable genes”, the significant pathway related with immunity (i.e. Fc epsilon ASIP and ADAMTS20 genes are particularly relevant: RI signaling, q-value = 0.02) was also identified based on they are involved in pigmentation [6, 8, 17, 35, 37–39] a human background gene set (Table 2). Several over- and co-localize with selection signals detected in a world- represented GO terms were related with embryonic skel- wide sample of goats [40]. In addition, we found that etal system morphogenesis (q-value = 1.22 × 10−3) and about 11.4% (89) of the annotated genes that co-localize G-protein coupled purinergic nucleotide receptor activity Guan et al. Genet Sel Evol (2020) 52:44 Page 5 of 10 Table 1 Main features of copy number variation regions (q-value = 6.22 × 10−3, Table 2). Interestingly, the copy (CNVR) detected in 1036 Murciano-Granadina goats number variable genes were also enriched in pathways Summary statistics Total Gain Loss Gain/loss with metabolic significance, such as prolactin signaling and insulin signaling, as well as GO terms related with Total length (Mb) 95.69 26.52 61.17 8 feeding behavior, but none of these pathways reached the Total number of CNVR 486 78 353 55 significance threshold (q-value ≤ 0.05) after correction Number of CNVR (< 10 kb) 1 1 0 0 for multiple testing (see Additional file 4: Table S4). Sev- Number of CNVR (10–50 kb) 4 2 1 1 eral of the pathways outlined in Additional file 4: Table S4 Number of CNVR (50–100 kb) 152 25 113 14 play important roles in immunity (e.g. chemokine sign- Number of CNVR (100–500 kb) 313 47 227 39 aling, B cell receptor signaling and T cell receptor sign- Number of CNVR (500 kb–1 Mb) 10 0 9 1 aling), cancer (e.g. endometrial cancer, proteoglycans in Number of CNVR (≥ 1 Mb) 6 3 3 0 cancer, thyroid cancer), as well as in oncogenic signal- Average number of SNPs per 5.59 9.01 5.03 4.35 ing (e.g. Ras and ErbB signaling) (see Additional file 4: CNVR Table S4), but most of them are not significant after cor- Minimum size of CNVR (kb) 2.04 2.04 23.2 43.1 rection for multiple testing. Maximum size of CNVR (kb) 11,124 11,124 1629.39 534.16 Average CNVR size (kb) 196.89 339.99 173.28 145.49 Standard deviation of CNVR 539.35 1299.49 156.89 91.51 size (kb) Fig. 2 Histograms displaying the distribution of CNVR according to their size (a) and frequency (b). CNVR that were longer than 1000 kb were included in the 1000-kb bin, whereas those with frequencies above 0.1 were grouped in the 0.1 bin. The histograms were drawn by using the ggplot2 package (http://ggplo​t2.tidyv​erse.org/) implemented in R (https​://www.r-proje​ct.org/) Guan et al. Genet Sel Evol (2020) 52:44 Page 6 of 10 Table 2 Functional enrichment of genes co-localizing with CNVR detected in 1036 Murciano-Granadina goats Background Category ID Term Number Fold enrichment P value q-value gene set of genes Goat KEGG chx04740 Olfactory transduction 69 2.33 1.26E−11 1.61E−10 Goat KEGG chx02010 ABC transporters 11 5.27 3.33E−05 4.27E−04 Goat KEGG chx04976 Bile secretion 11 3.90 4.46E−04 5.70E−03 Human KEGG hsa04664 Fc epsilon RI signaling pathway 8 4.71 1.40E−03 1.76E−02 Human GO/BP GO:0009952 Anterior/posterior pattern specification 12 5.56 9.36E−06 1.61E−04 Human GO/BP GO:0048704 Embryonic skeletal system morphogenesis 8 7.60 7.13E−05 1.22E−03 Human GO/BP GO:0035589 G-protein coupled purinergic nucleotide receptor 5 13.24 4.18E−04 7.16E−03 signaling pathway Human GO/CC GO:0016020 Membrane 81 1.40 1.45E−03 1.98E−02 Human GO/MF GO:0003677 DNA binding 67 1.48 1.10E−03 1.60E−02 Human GO/MF GO:0045028 G-protein coupled purinergic nucleotide receptor 5 13.19 4.24E−04 6.22E−03 activity KEGG: Kyoto Encyclopedia of Genes and Genomes pathway; GO/MF: gene ontology (GO) term related with molecular function; GO/BP: GO term related with biological process; GO/CC: GO term related with cellular component Validation of four copy number variants by real‑time to 2.38 copies relative to the calibrator (ADAMTS20), quantitative polymerase chain reaction 1.06 to 2.96 copies (BST1), 1.51 to 2.39 copies (NCKAP5) In order to confirm our results, we selected four CNVR and 1.83 to 3.28 copies (TNFAIP2). According to D’haene (i.e. CNVR_371_chr5, CNVR_506_chr6, CNVR_160_ et al. [36], copy number estimates between 1.414 and chr2 and CNVR_1229_chr21) that co-localized with 2.449 most likely correspond to a normal copy number of the ADAMTS20, BST1, NCKAP5 and TNFAIP2 genes, 2, whereas any number below or above these thresholds respectively (the primers used to amplify these CNVR could represent a deletion or a duplication, respectively. are listed in Additional file 1: Table S1). As shown in Thus, based on these values, evidence of copy number Fig. 3, the estimated copy numbers obtained by qPCR variation was inferred for three of the four genes ana- analysis of Murciano-Granadina goat samples were: 0.93 lyzed by qPCR. Fig. 3 Relative quantification of four copy number variation regions by real-time quantitative polymerase chain reaction analysis: a CNVR_371_chr5 (ADAMTS20), b CNVR_506_chr6 (BST1), c CNVR_160_chr2 (NCKAP5), d CNVR_1229_chr21 (TNFAIP2). The x and y axes represent sample ID and relative quantification of CNVR (mean ± standard error, with each sample analyzed in triplicate), respectively. As calibrator, we used the average of four samples estimated to have two copies (diploid status) based on the Goat SNP50 BeadChip analysis Guan et al. Genet Sel Evol (2020) 52:44 Page 7 of 10 Discussion lower CNVR frequencies ranging from 0.34% (Alpine In this work, our aim was to characterize copy number and Northern European goats) to 0.98% (Northwestern variation in Murciano-Granadina goats, a native Span- African goats). This decreased average CNVR frequency ish breed used for milk production. By genotyping 1036 is not very significant and probably reflects differences Murciano-Granadina goats with a SNP array, we were in sampling size and the use of composite populations able to identify 486 CNVR covering 3.9% of the goat with multiple breeds, each one with its specific CNVR genome, whereas Liu et al. [18] identified CNVR that frequencies. covered ~ 9% of the goat genome. The latter higher per- The CNVR detected in our study covered 779 protein- centage reported by Liu et al. [18] can be explained by coding genes. Pathway analyses reflected a substantial the fact that they analyzed 50 breeds with different geo- enrichment of genes that are involved in olfactory per- graphical origins, i.e. a composite population that is ception, which is consistent with previous reports in probably much more diverse than that used in our work. cattle [13, 14]. In this regard, there is an important dif- Besides, the pipeline that we used to identify CNVR ference between our results and those by Liu et al. [18]. is more stringent than that employed by Liu et al. [18], Whereas in the study of Liu et al. [18], the term “sensory removing CNVR that were not consistently detected by perception” was underrepresented among the CNV genes PennCNV and QuantiSNP. In the literature, estimates of (fold enrichment = 0.21), in our work the terms “olfactory 4.8 to 9.5% for CNVR coverage in the human genome are transduction” (fold enrichment = 2.33) and “G-protein reported [2]. Our results and those obtained by Liu et al. coupled purinergic nucleotide receptor activity” (fold [18] are consistent with these values. enrichment = 13.19) were overrepresented, and many Indeed, when Liu et al. [18] calculated the CNVR CNV genes were olfactory receptors. The two terms length for each breed normalized by the goat genome mentioned before are closely related because a broad size, their results agreed well with our estimate of 3.9%. array of purinergic receptors are differentially expressed For instance, this parameter reached values of 3.94% in the olfactory receptor neurons that modulate odor in goats from Southeastern Africa and 3.13% in goats responsiveness [42]. Moreover, purinergic nucleotides from Northwestern Africa and Eastern Mediterranean, are important neuromodulators of peripheral auditory whereas it was lowest (0.70%) for individuals from West and visual sensory systems [42]. In cattle, Keel et al. [13] Asia [18]. The number of CNV detected at the within- reported that “sensory perception of smell” and “G-pro- breed level by Liu et al. [18] was on average 126 CNV tein coupled receptor signaling pathway” were signifi- per breed and ranged from 6 to 714, whereas the aver- cantly overrepresented in the protein-coding genes that age number of CNVR was only ~ 20 per breed [18]. Since overlapped with CNVR. Similarly, Upadhyay et al. [14] the number of detected CNVR is proportional to popula- showed that “sensory perceptions of smell” and “chemi- tion size, for most of the breeds investigated in [18], the cal stimuli” are enriched in their set of CNV genes. A level of within-breed CNV variation is probably under- potential explanation for the underrepresentation of estimated. In summary, one important conclusion from the “sensory perception” functional category among the our study is that the magnitude of CNV diversity at the genes overlapping CNV reported by Liu et al. [18] could within-breed level is likely to be much larger than that be that in goats these genes are not well annotated yet, previously reported in studies that analyzed multiple so the majority of them are identified with a LOC pre- populations, each represented by a small or moderate fix and a number and, as a consequence of this, they are number of individuals. not correctly detected by PANTHER [43], thus biasing Most of the CNVR that we report here ranged in size the results obtained in the gene ontology enrichment from 50 to 500 kb, with a mean size of 196.89 kb. Simi- analysis. larly, the average CNVR size reported by Liu et al. [18] Loci belonging to large multigene families might be was 268 kb. Both estimates are quite large and reflect that more prone to co-localize with CNV because paralo- medium-density SNP arrays are not well suited to detect gous genes can act as templates in non-allelic homolo- small CNVR in spite of their high abundance. In cattle, gous recombination events, which promote increases the average sizes of CNVR detected with the Illumina or reductions in copy number [44]. It should be noted BovineHD Genotyping BeadChip (777 K SNPs) [14], Illu- that olfactory receptor genes constitute the largest gene mina whole-genome sequencing and PacBio sequencing superfamily, and in humans more than 900 genes and [41] were 66.15, 10 and 0.81 kb, respectively. Another pseudogenes have been identified [45]. In cattle, 1071 consistent feature of CNVR is that, in general, their fre- olfactory receptor genes and pseudogenes are distributed quencies are low or very low. In our study, approximately in 49 clusters across 26 bovine chromosomes [46], and 73% of the CNVR had frequencies lower than 1% and similar numbers have been reported for pigs [47]. More- the average frequency was 1.44%. Liu et al. [18] reported over, purifying selection against CNV is probably less Guan et al. Genet Sel Evol (2020) 52:44 Page 8 of 10 intense in regions that contain olfactory-receptor genes gene [6, 39]. Although some studies proposed that the than in genomic regions that contain genes with essential ASIP CNV might be associated with different pigmenta- functions [48]. Interestingly, copy number changes in the tion patterns in goats [8, 17, 37], no functional assay has olfactory receptor genes of wild and domestic mammals verified an association of ASIP copy number with ASIP might have consequences on food foraging as well as on mRNA levels. Another interesting shared copy number mate and predator recognition [49, 50]. variable gene is ADAMTS20, which was also identified in In the set of genes that co-localize with CNVR, we also two previous CNV surveys [17, 18]. This gene encodes a detected an enrichment of loci related with the multi- metalloproteinase with an important role in melanoblast gene family of ATP binding cassette (ABC) transport- survival by mediating Kit signaling [38] and in palatogen- ers, a result that agrees well with previous findings in esis [61]. Bertolini et al. [40] performed a selection scan humans [51–54] and cattle [14, 56]. In mammals, ABC in white vs. colored (black and red) goats and detected transporters fulfill the mission of carrying a broad array a selective sweep in the ADAMTS20 gene. In the light of of endogenous substrates, such as amino acids, peptides, these results, the potential involvement of a structural sugars, anions and hydrophobic compounds and metab- variation in ADAMTS20 in goat pigmentation should olites across lipid membranes. At least 49 ABC genes be explored further. Moreover, it is worthwhile to men- that belong to eight subfamilies have been identified in tion that several CNVR genes have functions related with the human genome [52]. Copy number variation in the production and reproduction traits. For instance, the human ABCC4 and ABCC6 genes is associated with sus- NCKAP5 gene, which co-localizes with CNVR_160_chr2 ceptibility to esophageal squamous cell carcinoma [51] (frequency = 0.1), is associated with milk fat percentage and to the rare autosomal recessive disease pseudoxan- in cattle [62]. Taking the above evidence into account, thoma elasticum [54], respectively. Moreover, large-scale the implication of structural chromosomal variations in deletions of the human ABCA1 gene are a causative fac- the genetic determinism of traits of economic interest tor for hypoalphalipoproteinemia [53], a disease that is with a complex inheritance deserves further exploration characterized by the complete absence of the apolipopro- by designing tools that allow inferring CNVR genotypes tein AI and extremely low levels of plasma high-density with high confidence. lipoprotein (HDL) cholesterol. We also found a highly significant enrichment of pathways related with embryo development (anterior/posterior pattern specification, Conclusions embryonic skeletal system morphogenesis), as previously With the PennCNV and QuantiSNP software, we reported [18]. These pathways are featured by genes that detected 486 CNVR in the genome of the Murciano- belong to the Hox multigene family of transcription fac- Granadina breed. In a previous study [18] that used a tors, possibly reflecting the genomic instability of certain less stringent pipeline (only PennCNV was used) and homeobox gene clusters as evidenced by the existence of included multiple populations with small to moder- many synteny/paralogy breakpoints and assembly gaps as ate sample sizes, the average number of CNVR events outlined in comparative studies [55]. per breed was ~ 20. One conclusion of our study is that Although not significant after correction for multiple CNV surveys, which are based on a broad array of breeds testing, we detected an enrichment of pathways with represented by only a few individuals, underestimate metabolic significance, such as prolactin and insulin the true levels of the CNV diversity at the within-breed signaling, which could have an impact on milk produc- level. The main reason for this outcome is that since the tion and growth [57–59]. Interestingly, the comparison majority of CNV have very low frequencies, they cannot of our work with that of Liu et al. [18] revealed 116 pro- be detected efficiently when sample size is small and, in tein-coding genes that co-localize with the set of shared consequence, much of the existing variation is missed. CNVR. One of the most relevant shared genes encodes We have also found that genes that overlap with CNV are ASIP, a protein that increases the ratio of pheomelanin functionally related with olfactory transduction, embryo to eumelanin by binding to the melanocortin 1 recep- development, ABC transporters and G-protein coupled tor and delivering an antagonist signal that blocks the purinergic nucleotide receptor activity. Most of these downstream expression of eumelanogenic enzymes [60]. genes belong to large multigene families encompassing Mutations in the ASIP gene play critical roles in animal tens, hundreds or thousands of paralogous genes that pigmentation [39]. For instance, the causal factor of the could act as substrates in non-allelic homologous recom- white color typical of many sheep breeds is the ubiqui- bination events, which is one of the main mechanisms tous expression of a duplicated copy of the ASIP coding generating duplications and deletions in humans and sequence, which is regulated by a duplicated promoter other species. Finally, we detected CNV that co-local- corresponding to the itchy E3 ubiquitin protein ligase ize with the ASIP and ADAMTS20 pigmentation genes, Guan et al. Genet Sel Evol (2020) 52:44 Page 9 of 10 which according to previous studies have been subjected Genética, Universidad de Córdoba, 14071 Córdoba, Spain. 3 Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma to positive selection for coat color in goats. de Barcelona, 08193 Bellaterra, Spain. 4 Department of Veterinary Medicine, University of Bari “Aldo Moro”, SP. 62 per Casamassima km. 3, 70010 Valenzano, Supplementary information BA, Italy. 5 Asociación Nacional de Criadores de Caprino de Raza Murciano-Gra‑ nadina (CAPRIGRAN), 18340 Granada, Spain. 6 Group of Research in Ruminants Supplementary information accompanies this paper at https​://doi. (G2R), Department of Animal and Food Science, Universitat Autònoma de org/10.1186/s1271​1-020-00564​-4. Barcelona (UAB), Bellaterra, Barcelona, Spain. Additional file 1: Table S1. List of primers used in the real-time quantita‑ Received: 2 March 2020 Accepted: 28 July 2020 tive PCR experiment to validate four putative copy number variable genes. Additional file 2: Table S2. List of copy number variation regions (CNVR) consistently detected with PennCNV and QuantiSNP in 1036 Murciano- Granadina goats. References Additional file 3: Table S3. List of copy number variable genes detected 1. Clop A, Vidal O, Amills M. Copy number variation in the genomes of in the current work and their concordance with those reported by Liu domestic animals. Anim Genet. 2012;43:503–17. et al. [18]. 2. Zarrei M, MacDonald JR, Merico D, Scherer SW. A copy number variation map of the human genome. Nat Rev Genet. 2015;16:172–83. Additional file 4: Table S4. Functional enrichment of genes co-localizing 3. Bickhart DM, Liu GE. The challenges and importance of structural varia‑ with copy number variation regions detected in 1036 Murciano-Granad‑ tion detection in livestock. Front Genet. 2014;5:37. ina goats. 4. Wright D, Boije H, Meadows JRS, Bed’hom B, Gourichon D, Vieaud A, et al. Copy number variation in intron 1 of SOX5 causes the pea-comb pheno‑ Acknowledgements type in chickens. PLoS Genet. 2009;5:e1000512. The authors are indebted to the Asociación Nacional de Criadores de Caprino 5. Chen C, Liu C, Xiong X, Fang S, Yang H, Zhang Z, et al. Copy number vari‑ de Raza Murciano-Granadina (CAPRIGRAN), and especially to Miguel García ation in the MSRB3 gene enlarges porcine ear size through a mechanism García and Teresa Novo Díaz, who collected all blood samples. We are also involving miR-584-5p. Genet Sel Evol. 2018;50:72. indebted to Dr. George Liu (USDA) for providing technical advice in CNV 6. Norris BJ, Whan VA. A gene duplication affecting expression of the calling. ovine ASIP gene is responsible for white and black sheep. Genome Res. 2008;18:1282–93. Authors’ contributions 7. Menzi F, Keller I, Reber I, Beck J, Brenig B, Schütz E, et al. Genomic amplifi‑ MA, JJ, VL and JVD conceived and designed the experiment; JF-A and XS cation of the caprine EDNRA locus might lead to a dose dependent loss contributed to the collection of materials and reagents, VL and AM performed of pigmentation. Sci Rep. 2016;6:28438. DNA extractions, AC and BC carried out DNA quality assessment and 8. Henkel J, Saif R, Jagannathan V, Schmocker C, Zeindler F, Bangerter E, et al. genotyping tasks, DG performed the bioinformatic analyses of the data, MGL Selection signatures in goats reveal copy number variants underlying contributed to the bioinformatic analyses of the data, AC and DG did the breed-defining coat color phenotypes. PLoS Genet. 2019;15:e1008536. quantitative PCR analyses, MA and DG wrote the manuscript. All authors read 9. Giuffra E, Törnsten A, Marklund S, Bongcam-Rudloff E, Chardon P, Kijas and approved the final manuscript. JMH, et al. A large duplication associated with dominant white color in pigs originated by homologous recombination between LINE elements Funding flanking KIT. Mamm Genome. 2002;13:569–77. This research was funded by the European Regional Development Fund 10. Pailhoux E, Vigier B, Chaffaux S, Servel N, Taourit S, Furet JP, et al. A 11.7- (FEDER)/Ministerio de Ciencia, Innovación y Universidades—Agencia Estatal kb deletion triggers intersexuality and polledness in goats. Nat Genet. de Investigación/Project Reference Grant: AGL2016-76108-R. We also acknowl‑ 2001;29:453–8. edge the support of the Spanish Ministry of Economy and Competitivity for 11. Sundström E, Imsland F, Mikko S, Wade C, Sigurdsson S, Pielberg GR, et al. the Center of Excellence Severo Ochoa 2016–2019 (SEV-2015-0533) grant Copy number expansion of the STX17 duplication in melanoma tissue awarded to the Centre for Research in Agricultural Genomics (CRAG, Bellaterra, from Grey horses. BMC Genomics. 2012;13:365. Spain) as well as to the CERCA Programme/Generalitat de Catalunya. Dailu 12. Letaief R, Rebours E, Grohs C, Meersseman C, Fritz S, Trouilh L, et al. Iden‑ Guan was funded by a Ph.D. fellowship from the China Scholarship Council tification of copy number variation in French dairy and beef breeds using (CSC). María Gracia Luigi-Sierra was funded with a FPI Ph.D. grant from the next-generation sequencing. Genet Sel Evol. 2017;49:77. Spanish Ministry of Economy and Competitivity (BES-2017-079709). 13. Keel BN, Lindholm-Perry AK, Snelling WM. Evolutionary and functional features of copy number variation in the cattle genome. Front Genet. Availability of data and materials 2016;7:207. The dataset supporting the conclusions of this article is accessible at Figshare 14. Upadhyay M, da Silva VH, Megens HJ, Visker MHPW, Ajmone-Marsan P, (https​://doi.org/10.6084/m9.figsh​are.12674​357). Bâlteanu VA, et al. Distribution and functionality of copy number variation across European cattle populations. Front Genet. 2017;8:108. Ethics approval and consent to participate 15. Hou Y, Bickhart DM, Hvinden ML, Li C, Song J, Boichard DA, et al. Fine The collection of blood is a routine procedure carried out by trained veterinar‑ mapping of copy number variations on two cattle genome assemblies ians working for the CAPRIGRAN association, so it does not require a permis‑ using high density SNP array. BMC Genomics. 2012;13:376. sion from the Committee on Ethics in Animal and Human Experimentation of 16. Fontanesi L, Martelli PL, Beretti F, Riggio V, Dall’Olio S, Colombo M, et al. the Universitat Autònoma de Barcelona. An initial comparative map of copy number variations in the goat (Capra hircus) genome. BMC Genomics. 2010;11:639. Consent for publication 17. Dong Y, Zhang X, Xie M, Arefnezhad B, Wang Z, Wang W, et al. Reference Not applicable. genome of wild goat (capra aegagrus) and sequencing of goat breeds provide insight into genic basis of goat domestication. BMC Genomics. Competing interests 2015;16:431. The authors declare that they have no competing interests. 18. Liu M, Zhou Y, Rosen BD, Van Tassell CP, Stella A, Tosser-Klopp G, et al. Diversity of copy number variation in the worldwide goat population. Author details Heredity. 2018;122:636–46. 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Uni‑ 19. Nandolo W, Lamuno D, Banda L, Gondwe T, Mulindwa H, Nakimbugwe H, versitat Autònoma de Barcelona, 08193 Bellaterra, Spain. 2 Departamento de et al. Distribution of copy number variants in the genomes of East African Guan et al. Genet Sel Evol (2020) 52:44 Page 10 of 10 goat breeds. In: Proceedings of the 11th world congress on genetics long-read and Illumina sequencing in New Zealand dairy cattle. J Dairy applied to livestock production, 11–16 Feb 2018, Auckland. 2018. Sci. 2017;100:5472–8. 20. Miller SA, Dykes DD, Polesky HF. A simple salting out procedure for 42. Hegg CC, Greenwood D, Huang W, Han P, Lucero MT. Activation of extracting DNA from human nucleated cells. Nucleic Acids Res. purinergic receptor subtypes modulates odor sensitivity. J Neurosci. 1988;16:1215. 2003;23:8291–301. 21. Tosser-Klopp G, Bardou P, Bouchez O, Cabau C, Crooijmans R, Dong Y, 43. Mi H, Muruganujan A, Huang X, Ebert D, Mills C, Guo X, Thomas PD. Proto‑ et al. Design and characterization of a 52 K SNP Chip for goats. PLoS One. col update for large-scale genome and gene function analysis with the 2014;9:e86227. PANTHER classification system (v.14.0). Nat Protoc. 2019;14:703–21. 22. Attiyeh EF, Diskin SJ, Attiyeh MA, Mossé YP, Hou C, Jackson EM, et al. 44. Hastings PJ, Lupski JR, Rosenberg SM, Ira G. Mechanisms of change in Genomic copy number determination in cancer cells from single gene copy number. Nat Rev Genet. 2009;10:551–64. nucleotide polymorphism microarrays based on quantitative genotyping 45. Glusman G, Yanai I, Rubin I, Lancet D. The complete human olfactory corrected for aneuploidy. Genome Res. 2009;19:276–83. subgenome. Genome Res. 2001;11:685–702. 23. Bickhart DM, Rosen BD, Koren S, Sayre BL, Hastie AR, Chan S, et al. Single- 46. Lee K, Nguyen DT, Choi M, Cha SY, Kim JH, Dadi H, et al. Analysis of cattle molecule sequencing and chromatin conformation capture enable de olfactory subgenome: the first detail study on the characteristics of the novo reference assembly of the domestic goat genome. Nat Genet. complete olfactory receptor repertoire of a ruminant. BMC Genomics. 2017;49:643–50. 2013;14:596. 24. Wang K, Li M, Hadley D, Liu R, Glessner J, Grant SF, et al. PennCNV: an inte‑ 47. Nguyen DT, Lee K, Choi H, Choi MK, Le MT, Song N, et al. The complete grated hidden Markov model designed for high-resolution copy number swine olfactory subgenome: expansion of the olfactory gene repertoire variation detection in whole-genome SNP genotyping data. Genome in the pig genome. BMC Genomics. 2012;13:584. Res. 2007;17:1665–74. 48. Young JM, Endicott RM, Parghi SS, Walker M, Kidd JM, Trask BJ. Extensive 25. Colella S, Yau C, Taylor JM, Mirza G, Butler H, Clouston P, et al. QuantiSNP: copy-number variation of the human olfactory receptor gene family. Am an Objective Bayes Hidden-Markov Model to detect and accurately map J Hum Genet. 2008;83:228. copy number variation using SNP genotyping data. Nucleic Acids Res. 49. Paudel Y, Madsen O, Megens HJ, Frantz LAF, Bosse M, Crooijmans RPMA, 2007;35:2013–25. et al. Copy number variation in the speciation of pigs: a possible promi‑ 26. Winchester L, Yau C, Ragoussis J. Comparing CNV detection methods for nent role for olfactory receptors. BMC Genomics. 2015;16:330. SNP arrays. Brief Funct Genomics. 2009;8:353–66. 50. Rinker DC, Specian NK, Zhao S, Gibbons JG. Polar bear evolution is 27. Pinto D, Darvishi K, Shi X, Rajan D, Rigler D, Fitzgerald T, et al. Compre‑ marked by rapid changes in gene copy number in response to dietary hensive assessment of array-based platforms and calling algorithms for shift. Proc Natl Acad Sci USA. 2019;116:13446–51. detection of copy number variants. Nat Biotechnol. 2011;29:512–20. 51. Sun Y, Shi N, Lu H, Zhang J, Ma Y, Qiao Y, et al. ABCC4 copy number 28. Diskin SJ, Li M, Hou C, Yang S, Glessner J, Hakonarson H, et al. Adjustment variation is associated with susceptibility to esophageal squamous cell of genomic waves in signal intensities from whole-genome SNP geno‑ carcinoma. Carcinogenesis. 2014;35:1941–50. typing platforms. Nucleic Acids Res. 2008;36:e126. 52. Vasiliou V, Vasiliou K, Nebert DW. Human ATP-binding cassette (ABC) 29. Zhang Z, Cheng H, Hong X, Di Narzo AF, Franzen O, Peng S, et al. transporter family. Hum Genomics. 2009;3:281–90. EnsembleCNV: an ensemble machine learning algorithm to identify and 53. Dron JS, Wang J, Berberich AJ, Iacocca MA, Cao H, Yang P, et al. Large- genotype copy number variation using SNP array data. Nucleic Acids Res. scale deletions of the ABCA1 gene in patients with hypoalphalipopro‑ 2019;47:e39. teinemia. J Lipid Res. 2018;59:1529–35. 30. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing 54. Kringen MK, Stormo C, Berg JP, Terry SF, Vocke CM, Rizvi S, et al. Copy genomic features. Bioinformatics. 2010;26:841–2. number variation in the ATP-binding cassette transporter ABCC6 gene 31. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis and ABCC6 pseudogenes in patients with pseudoxanthoma elasticum. of large gene lists using DAVID bioinformatics resources. Nat Protoc. Mol Genet Genomic Med. 2015;3:233–7. 2009;4:44–57. 55. Wilming LG, Boychenko V, Harrow JL. Comprehensive compara‑ 32. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: tive homeobox gene annotation in human and mouse. Database. paths toward the comprehensive functional analysis of large gene lists. 2015;2015:bav091. Nucleic Acids Res. 2009;37:1–13. 56. Liu GE, Hou Y, Zhu B, Cardone MF, Jiang L, Cellamare A, et al. Analysis 33. Ballester M, Castelló A, Ibáñez E, Sánchez A, Folch JM. Real-time quan‑ of copy number variations among diverse cattle breeds. Genome Res. titative PCR-based system for determining transgene copy number in 2010;20:693–703. transgenic animals. Biotechniques. 2004;37:610–3. 57. Fujita S, Rasmussen BB, Cadenas JG, Grady JJ, Volpi E. Effect of insulin on 34. Ramayo-Caldas Y, Castelló A, Pena RN, Alves E, Mercadé A, Souza CA, et al. human skeletal muscle protein synthesis is modulated by insulin-induced Copy number variation in the porcine genome inferred from a 60 K SNP changes in muscle blood flow and amino acid availability. Am J Physiol BeadChip. BMC Genomics. 2010;11:593. Endocrinol Metab. 2006;291:E745–54. 35. Fontanesi L, Beretti F, Riggio V, Gómez González E, Dall’Olio S, Davoli 58. Bequette BJ, Kyle CE, Crompton LA, Buchan V, Hanigan MD. Insulin R, et al. Copy number variation and missense mutations of the agouti regulates milk production and mammary gland and hind-leg amino acid signaling protein (ASIP) gene in goat breeds with different coat colors. fluxes and blood flow in lactating goats. J Dairy Sci. 2001;84:241–55. Cytogenet Genome Res. 2009;126:333–47. 59. Freeman ME, Kanyicska B, Lerant A, Nagy G. Prolactin: structure, function, 36. D’haene B, Vandesompele J, Hellemans J. Accurate and objective and regulation of secretion. Physiol Rev. 2000;80:1523–631. copy number profiling using real-time quantitative PCR. Methods. 60. Nasti TH, Timares L. MC1R, eumelanin and pheomelanin: their role in 2010;50:262–70. determining the susceptibility to skin cancer. Photochem Photobiol. 37. Zhang B, Chang L, Lan X, Asif N, Guan F, Fu D, et al. Genome-wide 2015;91:188–200. definition of selective sweeps reveals molecular evidence of trait-driven 61. Wolf ZT, Brand HA, Shaffer JR, Leslie EJ, Arzi B, Willet CE, et al. Genome- domestication among elite goat (Capra species) breeds for the produc‑ wide association studies in dogs and humans identify ADAMTS20 as a risk tion of dairy, cashmere, and meat. GigaScience. 2018;7:giy105. variant for cleft lip and palate. PLoS Genet. 2015;11:e1005059. 38. Silver DL, Hou L, Somerville R, Young ME, Apte SS, Pavan WJ. The secreted 62. Jiang J, Ma L, Prakapenka D, VanRaden PM, Cole JB, Da Y. A large-scale metalloprotease ADAMTS20 is required for melanoblast survival. PLoS genome-wide association study in U.S. Holstein cattle. Front Genet. Genet. 2008;4:e1000003. 2019;10:420. 39. Georges M, Charlier C, Hayes B. Harnessing genomic information for livestock improvement. Nat Rev Genet. 2019;20:135–56. 40. Bertolini F, Servin B, Talenti A, Rochat E, Kim ES, Oget C, et al. Signatures of Publisher’s Note selection and environmental adaptation across the goat genome post- Springer Nature remains neutral with regard to jurisdictional claims in pub‑ domestication. Genet Sel Evol. 2018;50:57. lished maps and institutional affiliations. 41. Couldrey C, Keehan M, Johnson T, Tiplady K, Winkelman A, Littlejohn MD, et al. Detection and assessment of copy number variation using PacBio