González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 DOI 10.1186/s12711-016-0258-1 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 On the performance of tests for the detection of signatures of selection: a case study with the Spanish autochthonous beef cattle populations Aldemar González‑Rodríguez1, Sebastián Munilla1,2, Elena F. Mouresan1, Jhon J. Cañas‑Álvarez3, Clara Díaz4, Jesús Piedrafita3, Juan Altarriba1,5, Jesús Á. Baro6, Antonio Molina7 and Luis Varona1,5* Abstract Background: Procedures for the detection of signatures of selection can be classified according to the source of information they use to reject the null hypothesis of absence of selection. Three main groups of tests can be identified that are based on: (1) the analysis of the site frequency spectrum, (2) the study of the extension of the linkage disequi‑ librium across the length of the haplotypes that surround the polymorphism, and (3) the differentiation among popu‑ lations. The aim of this study was to compare the performance of a subset of these procedures by using a dataset on seven Spanish autochthonous beef cattle populations. Results: Analysis of the correlations between the logarithms of the statistics that were obtained by 11 tests for detecting signatures of selection at each single nucleotide polymorphism confirmed that they can be clustered into the three main groups mentioned above. A factor analysis summarized the results of the 11 tests into three canoni‑ cal axes that were each associated with one of the three groups. Moreover, the signatures of selection identified with the first and second groups of tests were shared across populations, whereas those with the third group were more breed-specific. Nevertheless, an enrichment analysis identified the metabolic pathways that were associated with each group; they coincided with canonical axes and were related to immune response, muscle development, protein biosynthesis, skin and pigmentation, glucose metabolism, fat metabolism, embryogenesis and morphology, heart and uterine metabolism, regulation of the hypothalamic–pituitary–thyroid axis, hormonal, cellular cycle, cell signaling and extracellular receptors. Conclusions: We show that the results of the procedures used to identify signals of selection differed substantially between the three groups of tests. However, they can be classified using a factor analysis. Moreover, each canoni‑ cal factor that coincided with a group of tests identified different signals of selection, which could be attributed to processes of selection that occurred at different evolutionary times. Nevertheless, the metabolic pathways that were associated with each group of tests were similar, which suggests that the selection events that occurred during the evolutionary history of the populations probably affected the same group of traits. *Correspondence: lvarona@unizar.es 1 Departamento de Anatomía, Embriología y Genética, Universidad de Zaragoza, 50013 Saragossa, Spain Full list of author information is available at the end of the article © The Author(s) 2016. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/ publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 2 of 12 Background SNP genotyping and phasing The evolutionary history of animal populations involves Genomic DNA was extracted by standard protocols. both natural and artificial selection. These processes High-density SNP genotyping was performed at a com- not only affect the allelic frequencies at causal poly- mercial laboratory (Xenética Fontao, Lugo, Spain) by morphisms, but also the surrounding genomic regions using the BovineHD BeadChip (Illumina Inc, USA) due to the so-called “hitchhiking” effect. Thus, they may according to the manufacturer’s protocol; this HD chip is leave detectable signals on the structure of the genome designed to genotype 777,962 SNPs. The SNPs that were that can be identified by using appropriate procedures retained for our study were located on autosomal chro- [1, 2]. mosomes at a single position. Additional requirements The vast majority of the procedures used to detect were a Mendelian error rate lower than 0.05, and SNP signatures of selection [2] is based on the null hypoth- and individual call rates higher than 0.95. Quality control esis of absence of selection, which relies on the neutral was performed by using PLINK software [17] and finally, model of evolution [3]. In fact, these procedures can be 703,707 SNPs that covered 2,510,606 kb were available classified according to the source of information they for the analyses with on average one SNP per 3.567 kb. use to reject the null hypothesis. Based on the literature Haplotypes for the parental chromosomes were derived [2], three main groups of tests can be identified: the first with Beagle software [18] using the “TRIO” option. group is based on the analysis of the site frequency spec- trum [4–6], the second group focuses on the study of the Detection of signatures of selection extension of the linkage disequilibrium across the length The data were analysed using the following procedures of the haplotypes that surround a polymorphism [7, 8], for the detection of signatures of selection. and the third group is based on several measures of dif- ferentiation among populations [9–11]. In addition, the Tajima results of all these tests can be affected to some degree The procedure that was developed by Tajima [4] com- by demographic events and by the ascertainment bias pares two statistics to estimate the scaled mutation rate. caused by the procedure used to select the single nucleo- The first statistic (θπ ) is based on the number of segre- tide polymorphisms (SNPs) for the genotyping chip [12]. gating sites within a genomic region and the second (θκ ) Thus, the results of each test may not be fully consistent is the average heterozygosity at segregating sites in the with each other [13], which has led to propose strategies sample. The standardized difference between these two for summarizing results into a single statistic that either values, D = θπ − θκ , is used to infer departures from does [13] or does not [14, 15] account for the correlations neutrality. Theoretically, if D < 0 either the popula- between the results from different methods. tion has suffered expansion after a recent bottleneck or The aim of our study was to compare the performance a recent selective sweep has taken place; on the contrary of a subset of these procedures by using a dataset on if D > 0, the population has either experienced a sudden seven autochthonous beef cattle populations (Asturiana population contraction or is under balancing selection. de los Valles, Avileña-Negra Ibérica, Bruna dels Pirineus, The analysis was performed over sliding windows of 100 Morucha, Pirenaica, Retinta and Rubia Gallega) which SNPs by using own software. share close genetic relationships between them [16]. A second objective was to identify candidate genes and/or Fay and Wu metabolic processes that are associated with the regions This procedure [6] calculates the following statistic involved in the selection processes that occurred during D = θπ − θH , where θH depends on the number of sites at the evolution of these populations. which a derived allele is present within a genomic region. In the analysis, the ancestral alleles were extracted from Methods the study of Rocha et al. [19]. This test was computed Animals and sample size over sliding windows of 100 SNPs by using own software. A total of 171 sire/dam/offspring triplets were collected from seven Spanish beef cattle populations, includ- Fu and Li ing Asturiana de los Valles (AV, n = 25), Avileña-Negra This procedure [5] is based on counting the number Ibérica (ANI, n = 24), Bruna dels Pirineus (BP, n = 25), of singletons or alleles present in only one phase. The Morucha (Mo, n = 24), Pirenaica (Pi, n = 24), Retinta rationale is that a selection process will extend time to (Re, n = 24) and Rubia Gallega (RG, n = 24) breeds. The coalescence so that a larger number of mutations may selected parents were chosen as unrelated as possible to take place in new or external branches of the tree and fully represent the diversity of the populations. thus appear only once in the observed sample. As before, González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 3 of 12 the analysis was performed over sliding windows of 100 Selestim SNPs by using own software. This procedure [10] assumes a hierarchical Bayesian model to distinguish selected polymorphisms from the iHS background of neutral (or almost neutral) polymor- This procedure [8] calculates the ratio of the integrated phisms and also to estimate the intensity of selection in haplotype score (iHH) for the ancestral allele and the each population. The model assumes a binomial distribu- derived allele at a given SNP. The iHH is the integral tion of the allele counts at each locus and for each popu- (area) of the observed decay of the EHH (extended hap- lation, and the prior distribution of allelic frequencies is lotype homozygosity) as defined by Sabeti et al. [7]. As modeled under the assumption of a stationary density in the previous test, ancestral alleles were extracted from of the diffusion process [10]. The model is implemented Rocha et al. [19]. The iHS was calculated with the selscan by using a Markov chain Monte Carlo method. SelEstim software [20] using the parameters recommended by the software (http://www1.montpellier.inra.fr/CBGP/soft- authors. For further calculations, we used the |iHS|. ware/selestim/) was used for this purpose with the stand- ard parameters that are recommended by the authors. nSL Among the outputs provided by the Selestim approach, This procedure was recently presented by Ferrer- we extracted the σij parameter [10], which represents the Admetlla et al. [21]. The procedure of calculation is simi- coefficient of selection for the ith subpopulation and the lar to iHS, but replaces IHH by an alternative statistic jth locus. (SL) that measures the length of a segment of haplotype homozygosity in terms of segregating sites. The main XP‑CLR advantage of nSL over iHS is that it uses segregating sites This approach [23] assumes that the allele frequen- as a measure of distance, while iHS needs the recombina- cies of two populations that diverge from an ancestral tion distance. Thus, the iHS is more sensitive to recombi- population follow a Gaussian distribution for which nation rate [21]. The analysis used the same parameters the variance contains information on the history of as in the iHS test and own software. As before, we used the populations since they split. Under the assump- |nSL|. tion that the evolutionary process is reversible, the procedure defines the distribution of allelic frequen- H12 cies in the first population (reference) given the allele This method was recently proposed [22] with the H12 frequencies in the second population (objective). We statistic being defined as: calculated XP-CLR by taking each pair of populations  as objective and reference with the software XPCLR H 12 = (p1 + p2 )2 + pj2 , (http://genetics.med.harvard.edu/reich/Reich_Lab/ j>2 Software.html). Then, we averaged the six available where pj is the frequency of the jth most common haplo- tests for each population that was treated as an objec- type in the population. Here, the frequencies of the first tive population, to infer the signatures of selection for and second most common haplotypes were combined each breed. into a single frequency. The calculation was performed over sliding windows of 100 SNPs by using own software. XP‑EHH This approach [7] is also computed for each pair of pop- Fixation index (FST) ulations. For each population, as in the iHS test, it cal- This procedure was described by Wright [9] and is the culates the EHH between a core SNP and a set of SNPs most classical approach to study the pattern of differen- within a predefined genomic interval and integrates it tiation between populations. The fixation index FST is cal- with respect to genetic distance to calculate the integrated culated for each SNP and for each pair of populations as haplotype score (IHH) for populations A and B. Then, the FST = (HO − HE )/HE , where HO and HE are the observed statistic is computed as XPEHHAB = ln (IHHA /IHHB ). and expected heterozygosities, respectively. Estimates for As previously, we computed this statistic for each SNP FST were averaged over sliding windows of 100 SNPs and and each pair of populations and the results were aver- assigned to the central SNP in each window. The proce- aged over the six comparisons to generate a unique dure was computed with own software. Finally, the results result for each population. We used the software selscan for each population were computed by averaging the [20] with the parameters that are recommended by the paired FST estimates with the other six populations. authors. González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 4 of 12 VarLD the relaxation of the empirical threshold was to capture This procedure [24] evaluates the magnitude of the dif- softer signals of selection. With these selected genomic ferences in linkage disequilibrium between a pair of regions for each canonical axis, we used the software populations. It calculates the linkage disequilibrium as WebGestalt [28] (http://bioinfo.vanderbilt.edu/webge- the correlation coefficient between pairs of SNPs within stalt/) by setting the Homo sapiens genome as the refer- a genomic region and creates a matrix of those correla- ence genome. In addition, we used a hypergeometric p tions for each population. Then, it evaluates the differ- value to correct for multiple-testing. The results included ences between the matrices of both populations as the the top 10 pathways (WikiPathways). difference between its eigenvalues. The procedure was computed using the software VarLD [25] over sliding Results and discussion windows of 100 SNPs. Summarizing footprints of selection detected by 11 For all the above-described methods, we used the procedures empirical distribution of the results generated along the A large set of procedures is available for the identification genome as the null distribution of the test, in order to of footprints of selection across the genome [2]. Most of reduce the possible effects of the demographic history these procedures are based on the rejection of the null or the ascertainment bias. The underlying hypothesis is hypothesis of absence of selection based on the neutral that, on average, both demographic events and ascertain- theory of evolution [3]. However, each of these methods ment bias affect all the genome in a similar way, and thus, calculates a different statistic to test this hypothesis. In deviations or extreme values of the empirical distribution addition, they are influenced to varying degrees by demo- could be understood as signals of selection events. graphic history and ascertainment bias caused by the selection of SNPs [12]. Thus, it is expected that each test Summary of signals of selection provides a different output as confirmed by the correla- In order to detect communalities and summarize the tions between the results obtained by the 11 procedures results of the 11 procedures for ease of interpretation, we used in this study (Fig. 1) and by the Manhattan plots normalized these results for each SNP using a logarithm generated with the results for each test and population transformation to make the scale of the different results (see Additional file 1: Figure S1, Additional file 2: Figure comparable and, then, we calculated the correlation S2, Additional file 3: Figure S3, Additional file 4: Figure between the logarithms (or the negative of logarithms for S4, Additional file 5: Figure S5, Additional file 6: Figure the Tajima, Fu and Li and Fay and Wu procedures) for the S6, Additional file 7: Figure S7, Additional file 8: Figure 703,707 SNPs. In a confirmatory analysis, provided that S8, Additional file 9: Figure S9, Additional file 10: Figure the methods used to detect signatures of selection were S10, Additional file 11: Figure S11). classified into three groups, we performed a factor analy- In order to summarize the signals of selection that were sis restricted to a subspace of three axes using a varimax detected by the 11 tests, there are procedures to con- rotation [26]. The analysis was done with R [27] by using dense such results into a single statistic by using Bayes the function principal() included in the package psych. factors [14] or a combination of p values [13, 15]. How- ever, these strategies imply that the signals of selection Selection of candidate genes that are captured by the different methods are compara- First, we identified candidate genes based on the empiri- ble. Nevertheless, as the definition of the null hypothesis cal distribution of the output of the three canonical varies between tests, the signals of selection identified axes of the factor analysis. Thus, we defined a very strict by each procedure may correspond to different types of threshold by selecting the genomic 1-Mb regions with at selection events. In fact, some authors [29] pointed out least 25 SNPs that were in the top 0.1% of the results for that within-population haplotype length methods [8, 21] each axis. Then, we used the Ensembl-Biomart database can detect only very recent selection processes, because to identify the genes that were present in those genomic they become ineffective when the selected alleles reach regions and compared our results with those in the litera- fixation or are very close to fixation. The same authors ture to identify potential candidate genes for selection in [29] indicated that signals of selection that are based on the bovine populations. a reduction of genetic diversity [4, 5] persist for a longer period of time and these methods can detect older sig- Enrichment analysis nals of selection, while tests that are based on population Finally, in order to obtain a clearer picture of the meta- differentiation [9, 10] occupy an intermediate position. bolic pathways that were affected by the selection pro- In this study, the correlations of the absolute loga- cesses, we identified the genomic regions that were rithm of the results between the 11 methods used were above the top 5% of each canonical axis. The objective of low or even negative (Fig. 1). However, there are some González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 5 of 12 Fig. 1 Heatmap of the correlations between logarithms of the results of the 11 tests applied for the detection of signatures of selection and their clustering remarkable exceptions such as the correlations between Such a structure of the correlations between the results of Tajima and Fu and Li, and iHS and nSL tests, that were these methods indicates that the implementation of a fac- remarkably high. It should be noted that the Tajima and tor analysis, as suggested by Simianer et al. [30], could be Fu and Li tests are both based on the analysis of the appropriate to summarize the results into a few canonical site frequency spectrum, and that the nSL test is just a axes. In addition, as described below, each axis was asso- modification of iHS where the map distance is replaced ciated with signatures of selection of a different kind. In by the number of segregating sites [21]. A more detailed particular, we applied a factor analysis restricted to three analysis of the structure of these correlations allows to canonical axes using a varimax approximation [26] that identify three main groups of tests, one group based on explains up to the 56% of variation. the site frequency spectrum (Tajima, Fu and Li, Fay and Table 1 shows the loadings for the canonical axes that Wu); a second group based on the haplotype length (iHS resulted from the factor analysis for each of the 11 tests and nSL), and a third group that focuses on the differen- used to detect selection signatures. Moreover, Table 1 tiation between populations (FST, SelEstim and XP-CLR). presents the correlations between the canonical axes The remaining tests (VarLD, H12 and XP-EHH) are in an and each specific test, which are fully consistent with intermediate position between the latter two, although the results of the correlations presented in Fig. 1. The slightly closer to tests based on population differentiation. first axis explains 20% of the variation and shows a high González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 6 of 12 Table 1 Weights in the factor analysis with, between and explains 17% of the total variation. For each test, the parentheses, the correlation between the results of each three axes explain between 32 (XP-CLR) and 90% (iHS) test and the canonical axis, and percentage of variance of the variation. The Manhattan plots of the results that explained by the three axes relate to the three canonical axes are in Figs. 2, 3 and 4. Method First axis Second axis Third axis % variance The first two axes presented a higher level of shared sig- nals between populations (see Figs. 2, 3) whereas the Tajima −0.07 (0.07) 0.42 (0.85) −0.01 (0.06) 73 results of the third axis were, in general, breed-specific. Fu-Li −0.08 (−0.02) 0.35 (0.68) −0.05 (−0.04) 47 This statement is supported by the results in Fig. 5, which Fay-Wu −0.05 (0.09) 0.38 (0.77) −0.00 (0.07) 61 shows the correlations of the results obtained for the Selestim 0.28 (0.58) −0.10 (−0.05) 0.02 (0.13) 36 first, second and third canonical axes between popula- XPCLR 0.22 (0.51) 0.06 (0.24) −0.02 (0.06) 32 tions. An average correlation of 0.50 was found for the H12 0.29 (0.67) 0.08 (0.32) −0.04 (0.06) 55 first axis [ranging from 0.39 (BP and Re) to 0.71 (AV and IHS −0.06 (0.07) −0.04 (0.03) 0.53 (0.95) 90 RG)]. Furthermore, the second axis also showed high NSL −0.04 (0.11) −0.02 (0.07) 0.52 (0.94) 89 correlations between populations that ranged from 0.37 FST 0.38 (0.77) −0.10 (−0.02) −0.03 (0.08) 60 (Re and Pi) to 0.60 (AV and RG) with an average of 0.49. XP-EHH 0.17 (0.47) 0.14 (0.40) 0.01 (0.13) 40 On the contrary, the correlations for the third axis were VarLD 0.31 (0.59) −0.11 (−0.09) −0.10 (−0.08) 36 lower with an average value of 0.08 and ranged from 0.05 (Pi and Re) to 0.16 (AV and BP). In addition, the structure of the correlations (Fig. 5) confirmed the classification of correlation with the procedures based on the analysis the populations into two main clusters, one composed by of population differentiation (FST, SelEstim, XP-EHH, the ANI, Mo and Re populations and the other by the Pi, XP-CLR and VarLD) and H12; the second axis explains BP, RG and ANI populations, as previously reported by 19% of the variation and is correlated with the methods Cañas-Álvarez et al. [16] based on distance measures and based on the site frequency spectrum (Tajima, Fu and Li, admixture analysis. and Fay and Wu); and, finally, the third axis is strongly As in the study of Sabeti et al. [29], our results may correlated with methods based on the extension of link- indicate that old selection or adaptation processes age disequilibrium or haplotype length (iHS and nSL) that occurred before breed differentiation or during Fig. 2 Manhattan plots for the results of the first axis (a) and genomic regions identified with at least 25 SNPs within the top 0.1% of the results (b) González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 7 of 12 Fig. 3 Manhattan plots for the results of the second axis (a) and genomic regions identified with at least 25 SNPs within the top 0.1% of the results (b) Fig. 4 Manhattan plots for the results of the third axis (a) and genomic regions identified with at least 25 SNPs within the top 0.1% of the results (b) González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 8 of 12 Fig. 5 Correlations of the results from the first, second and third canonical axes between populations speciation were only detected by the site frequency spec- In fact, the most remarkable signal of selection from trum methods that are associated with the second canon- this third axis was identified on chromosome 2 around ical axis. The signals of selection that were generated the myostatin (MTSN) gene (between 6,213,566 and by later isolation and recent selection events within the 6,220,196 bp) in the AV population, where double-mus- populations are identified by the differentiation meth- cling is included as a criterion of selection in its breed- ods, which are linked to the first canonical axis. Finally, ing program. This specific genomic region can be used the haplotype length methods, summarized in the third to illustrate the timing of the signatures of selection that canonical axis, identified more recent and, in general, less were detected by each group of methods. In the AV pop- intense selection events that are mostly specific to each ulation, two large signatures of selection were detected population. The absence of regions with strong recent with the first and third canonical axes, respectively. The signals of selection agrees with the postulate that artifi- first axis is related to processes that were involved in the cial selection processes do not leave relevant signatures creation of the breed and the third axis to recent selec- of selection [31]. The main reason of this absence can be tion. In addition, a large signature of selection associ- due to the polygenic nature of most of the traits associ- ated with the first canonical axis in the RG population ated with current selection processes [32–34] or to the was observed. However, in this population, there is no effect of epistasis [35]. relevant signature in the results of the third axis. This González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 9 of 12 result is consistent with several previous studies [36] factor 7) genes, which are close to a strong signature of that reported the presence of haplotypes associated with selection that was reported by Gautier [44] and is associ- double-muscling also in the RG population. This may ated with the VDAC1 (voltage-dependent anion-selective indicate that some degree of selection around this gene channel protein 1) gene. A strong signature of selection may have occurred during the process of breed forma- was observed for the genomic region on BTA13, in the tion, but the current breeding program does no longer Re population, where is located the END3 (endothelin 3) put any selection pressure on double-muscling. Finally, gene that plays a role in melanocyte development [45] the results from the second canonical axis are less rele- and was recently associated with piebald pattern [44]. vant (AV) or even absent (RG), which indicates that, for Finally, there is a very strong signature of selection on this group of tests, the selection effects may be diluted BTA11 for the BP population, where the closest gene to because a larger number of generations without selection the maximum signal is BMP10 (bone morphogenetic pro- is considered. tein 10). Within this genomic region on BTA11, some authors [46, 47] identified signatures of selection and sev- Candidate genes and metabolic paths eral genes that could be associated with fertility: PROKR1 The results of the first axis (Fig. 2) allowed us to highlight (prokineticin receptor 1), GFPT1 (glutamine-fructose- seven relevant genomic regions on Bos taurus chromo- 6-phosphate transaminase 1), GMCL1 (germ cell-less some (BTA) 2 (between 1,047,347 and 11,899,039 bp), spermatogenesis associated 1), PCBP1 (poly(rC) binding BTA5 (between 15,920,995 and 20,321,882 bp), BTA6 protein 1) and EHD3 (EH-domain containing 3). (between 37,853,912 and 41,160,000 bp), BTA7 (between The results of the second canonical axis (Fig. 3) 47,276,124 and 47,745,164 bp), BTA11 (between confirmed some of the signals of selection that were 65,077,840 and 72,203,248 bp), BTA13 (between detected in the first axis, but also revealed several new 57,430,392 and 57,754,760 bp) and BTA18 (between ones that are shared by several populations and located 12,675,262 and 16,202,289 bp). In some cases, these on five chromosomes: BTA2 (between 61,684,232 and genomic regions were extremely large, because of strong 62,199,344 and between 72,158,144 and 73,356,296 bp), signatures of selection such as those on BTA2 for the BTA7 (between 20,612,988 and 21,163,812 bp), BTA13 AV and RG populations or on BTA5, BTA6 and BTA11 (between 11,860,881 and 12,062,522 bp), BTA16 for the BP population. Such huge signatures of selection (between 44,612,592 and 45,846,144 bp) and BTA21 imply that large genomic regions included SNPs that (between 32,207,264 and 32,414,316 bp). Previously, two were associated with results above the top 0.1% of the meta-analyses [38, 43] showed that these regions were empirical distribution along the genome. However, the also associated with signatures of selection in other popu- localization of the strongest signals within each genomic lations. Among the genes included in these regions, some region and for each population allowed us to narrow of them may be good candidates for being affected by down the genomic regions (Fig. 2b), which are similar to selection i.e.: (1) genes that are related to energy balance those reported in a previous study on the differentiation and homeostasis: R3HGM1 (R3H domain containing 1) between populations [37]. These regions included well- on BTA2 [48]; CAMK1D (calcium/calmodulin-dependent known genes that were previously reported as potential protein kinase ID) on BTA13 [49]; and SLC25A33 (solute candidates of selection signatures in cattle [38], such as carrier family 25 (pyrimidine nucleotide carrier), member MTSN (myostatin) on BTA2, suggested in several beef 33) and SLC2A5 (solute carrier family 2 (facilitated glu- cattle populations [39–41], KIT-LG (kit-ligand) on BTA5 cose/fructose transporter), member 5) on BTA16 [49, 50]; with a very large peak in the BP population, MC1R (mel- (2) PLIN5 (perilipin 5) on BTA7, which is involved in the anocortin 1 receptor) on BTA18, which controls the regulation of lipid metabolism in rat and pigs [51, 52]; (3) production of eumelanin (black) or pheomelanin (red) SCAPER (s-phase cyclin A-associated protein in the endo- pigments [42] and appears to be relevant in populations plasmic reticulum) on BTA21, which regulates cell cycle with black (AV and Mo) or red (Re) coat color. Moreo- progression [53]; and (4) there are two other signatures ver, it should be also highlighted that the region on BTA6 of selection that are worth noting i.e. one detected for that includes LAP3 (leucine aminopeptidase3), LCORL the RG population on BTA19 (between 27,941,270 and (ligand dependent nuclear receptor corepressor-like) and 28,571,032 bp) where ALOX15B (arachidonate 15-lipoxy- NCAPG (non-SMC condensing I complex, subunit G) genase, type B) and ALOX12B (arachidonate 12-lipoxyge- and was identified in two meta-analyses [38, 43] as one nase, 12R type) are located and are involved in immune the genomic regions that is most frequently identified response [54], and one for the ANI population on BTA20 with signatures of selection in the bovine genome. The (between 40,854,136 and 40,996,384 bp) where the NPR3 genomic region identified on BTA7 includes the CAMLG (natriuretic peptide receptor 3) gene is located, which is (calcium modulating ligand) and TCF (transcription related with cattle stature [55]. González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 10 of 12 Although less relevant, the results of the third canoni- Table 2 Top 10 enriched pathways for the three axes cal axis (Fig. 4) confirm the signatures of selection Pathway Ngenesa Totalb around the MTSN gene for the AV, RG and Pi popula- tions, and around the complex LAP–LCORL–NCAG for First axis the BP population. There are several other interesting Focal adhesion 48 185 signatures of selection such as those located on BTA27 Integrated pancreatic cancer pathway 46 181 (between 36,466,580 and 40,862,444 bp) and BTA28 MAPK signalling pathway 43 163 (between 41,643,416 and 45,215,488 bp) for the ANI and Lymphocite TarBase 96 533 Mo populations, respectively. The first genomic region Epithelium TarBase 69 340 includes the IKBKB (inhibitor of kappa light polypep- TSH signalling pathway 25 70 tide gene enhancer in B-cells, kinase beta), DKK4 (dick- Adipogenesis 36 130 kopf WNT signaling pathway inhibitor 4) and VDAC3 Cytoplasmic ribosomal proteins 27 88 (voltage-dependent anion-selective channel protein 3) Muscle cell TarBase 77 424 genes that are associated with immune response to GPCRs, class A rhodopsin-like 53 259 trypanosoma infection in African populations [56], and Second axis the second region contains the ALOX5 (arachidonate Epithelium TarBase 51 340 5-lipoxygenase) and RASSF4 (ras association (RalGDS/ Lymphocyte TarBase 62 533 AF-6) domain family member 4) genes, that are related Translation factors 15 51 with growth [57] and feed conversion [58], respectively. Focal adhesion 27 185 Finally, we performed a pathway enrichment analysis Adipogenesis 21 130 [59] for the identified genomic regions by applying a less Muscle cell TarBase 44 424 restrictive criterion (top 5%). The objective of this anal- Notch signalling pathway 11 45 ysis was to identify a larger number of genomic regions Integrated pancreatic cancer pathway 23 181 associated with signatures of selection although with less G1 to S cell cycle control 14 77 strong signals. The results of the top 10 pathways that Eukaryotic transcription initiation 10 41 were identified by the enrichment analysis are in Table 2. Third axis In general, the pathways associated with each canonical Lymphocyte TarBase 304 533 axis are coincident, which indicates that the metabolic MAPK signalling pathway 124 165 pathways that were involved in old and recent selection Insulin signalling 123 163 events are similar, although probably with variable inten- Muscle cell TarBase 240 424 sities and directions [60]. The enrichment analysis iden- Calcium regulation in the cardiac cell 116 151 tified pathways that are related with immune response Focal adhesion 132 185 (lymphocyte TarBase), muscle development (muscle cell Integrated pancreatic cancer pathway 130 181 TarBase), protein biosynthesis (translation factors, cyto- Myometrial relaxation and contraction pathways 116 162 plasmic ribosomal proteins), skin and pigmentation (epi- Adipogenesis 99 130 thelium TarBase), glucose metabolism (insulin signaling, Epithelium TarBase 191 340 integrated pancreatic cancer pathway), fat metabolism a Ngenes: number of genes present in the genomic regions (adipogenesis), embryogenesis and morphology (focal b Total: number genes in the pathway adhesion), heart (calcium regulation in the cardiac cell) and uterine metabolism (myometrial relaxation and contraction pathways), regulation of the hypothalamic– that the metabolic pathways involved in old and recent pituitary–thyroid axis (TSH signalling pathway), hor- processes of selection are similar to those detected by monal, cellular cycle (MAPK-signaling pathway, G1 to using equivalent approaches in other cattle populations. S cell cycle control, eukaryotic transcription initiation), cell signaling (notch signaling pathway) and extracellular Conclusions receptors (GPCR, class A rhodopsin-like). Among these, In this study, we confirm that the results of various proce- 10 pathways (focal adhesion, integrated pancreatic cancer dures used to identify signatures of selection varied largely pathway, adipogenesis, myometrial relaxation and con- among groups of tests depending on the source of infor- traction pathways, adipogenesis, lymphocyte TarBase, mation they use to reject the null hypothesis of absence insulin signaling, MAPK signaling pathway, focal adhe- of selection. However, we observed some correlations sion, epithelium TarBase) had been previously identified between the results of each test. Accordingly, these tests in a meta-analysis based on a very large number of stud- could be clustered into three groups that matched with ies on selection signatures in cattle [38], which confirmed the three canonical axes of a factor analysis. Moreover, González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 11 of 12 each canonical factor (or group of tests) identified differ- 289592—Gene2Farm, from FEDER funds and from the regional government of Aragon (DGA) through its research group A51. The collaboration of Breed Soci‑ ent signals of selection, which were assigned to selection eties in collecting samples and the support of FEAGAS are also acknowledged. events that occurred at different evolutionary times. In The authors wish to thank JJ Arranz (Universidad de Leon) for facilitating the fact, older selection events generated signatures of selec- sampling of the Morucha population. J. J. Cañas-Álvarez acknowledges the financial support provided by COLCIENCIAS through the Francisco José de tion that presented communalities between populations, Caldas fellowship 497/2009 and A González-Rodríguez acknowledges the whereas more recent selection events were detected spe- financial support from the Spanish Government BES-2011-045434. cifically for each population. Nevertheless, the enriched Competing interests metabolic pathways associated to each group of tests The authors declare that they have no competing interests. showed an important degree of agreement which suggests that the traits involved in the selection events were similar Received: 22 March 2016 Accepted: 18 October 2016 during the evolutionary history of the populations. Additional files References 1. Nielsen R. Molecular signatures of natural selection. Annu Rev Genet. Additional file 1: Figure S1. Manhattan plots of the results along the 2005;39:197–218. autosomal genome obtained with the Tajima procedure. 2. Qanbari S, Simianer H. Mapping signatures of positive selection in the Additional file 2: Figure S2. Manhattan plots of the results along the genome of livestock. Livest Sci. 2014;166:133–43. autosomal genome obtained with the Fu and Li procedure. 3. Kimura M. The neutral theory of molecular evolution. Cambridge: Cam‑ Additional file 3: Figure S3. Manhattan plots of the results along the bridge University Press; 1983. autosomal genome obtained with the Fay and Wu procedure. 4. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95. Additional file 4: Figure S4. Manhattan plots of the results along the 5. Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. autosomal genome obtained with the SelEstim procedure. 1993;133:693–709. Additional file 5: Figure S5. Manhattan plots of the results along the 6. Fay JC, Wu CI. Hitchhiking under positive Darwinian selection. Genetics. autosomal genome obtained with the XP-CLR procedure. 2000;155:1405–13. 7. Sabeti PC, Reich DE, Higgins JM, Levine HZP, Richter DJ, Schaffner SF, Additional file 6: Figure S6 Manhattan plots of the results along the et al. Detecting recent positive selection in the human genome from autosomal genome obtained with the H12 procedure. haplotype structure. Nature. 2002;419:832–7. Additional file 7: Figure S7. Manhattan plots of the results along the 8. Voight BF, Kudaravalli S, Wen X, Pritchard JK. A map of recent positive autosomal genome obtained with the iHS procedure. selection in the human genome. PLoS Biol. 2006;4:e72. 9. Wright S. Isolation by distance. Genetics. 1943;28:114–38. Additional file 8: Figure S8. Manhattan plots of the results along the 10. Vitalis R, Gautier M, Dawson KJ, Beaumont MA. Detecting and measuring autosomal genome obtained with the nSL procedure. selection from gene frequency data. Genetics. 2014;196:799–817. Additional file 9: Figure S9. Manhattan plots of the results along the 11. Beaumont MA, Balding DJ. Identifying adaptive genetic divergence autosomal genome obtained with the FST procedure. among populations from genome scans. Mol Ecol. 2004;13:969–80. 12. Nielsen R. Population genetic analysis of ascertained SNP data. Hum Additional file 10: Figure S10. Manhattan plots of the results along the Genomics. 2004;1:218–24. autosomal genome obtained with the XP-EHH procedure. 13. Ma Y, Ding X, Qanbari S, Weigend S, Zhang Q, Simianer H. Properties of Additional file 11: Figure S11. Manhattan plots of the results along the different selection signature statistics and a new strategy for combining autosomal genome obtained with the VarLD procedure. them. Heredity (Edinb). 2015;115:426–36. 14. Grossman SR, Shlyakhter I, Karlsson EK, Byrne EH, Morales S, Frieden G, et al. A composite of multiple signals distinguishes causal variants in regions of positive selection. Science. 2006;327:883–6. Authors’ contributions 15. Utsunomiya YT, Perez-O’Brien AM, Sonstegard TS, Van Tassell CP, do LV conceived the study; CD, JP and JA contributed to the study design; AG, Carmo AS, Mészáros G, et al. Detecting loci under recent positive selec‑ SM, EFM, JJC and LV developed and applied the methods; AG, CD, JP and LV tion in dairy and beef cattle by combining different genome-wide scan drafted the manuscript; CD, JP, JA, AM and JAB contributed with the genera‑ methods. PLoS One. 2013;8:e64280. tion of experimental data. All authors read and approved the final manuscript. 16. Cañas-Álvarez JJ, González-Rodríguez A, Munilla S, Varona L, Díaz C, et al. Genetic diversity and divergence among Spanish beef cattle breeds Author details assessed by a bovine high-density SNP chip. J Anim Sci. 2015;93:5164–74. 1 17. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. Departamento de Anatomía, Embriología y Genética, Universidad de Zaragoza, 50013 Saragossa, Spain. 2 Departamento de Producción Animal, PLINK: a tool set for whole-genome association and population-based Facultad de Agronomía, Universidad de Buenos Aires, 1417 Buenos Aires, linkage analyses. Am J Hum Genet. 2007;81:559–75. Argentina. 3 Grup de Recerca en Remugants, Departament de Ciència Animal i 18. Browning BL, Browning SR. A unified approach to genotype imputation dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, and haplotype-phase inference for large data sets of trios and unrelated Spain. 4 Departamento de Mejora Genética Animal, INIA, 28040 Madrid, Spain. individuals. Am J Hum Genet. 2009;84:210–23. 5 Instituto Agroalimentario de Aragón (IA2), 50013 Saragossa, Spain. 6 Depar‑ 19. Rocha D, Billerey C, Samson F, Boichard D, Boussaha M. Identification of tamento de Ciencias Agroforestales, Universidad de Valladolid, 34004 Palencia, the putative ancestral allele of bovine single-nucleotide polymorphism. J Spain. 7 MERAGEM, Universidad de Córdoba, 14071 Córdoba, Spain. Anim Breed Genet. 2014;131:483–6. 20. Szpiech ZA, Hernández RD. Selscan: an efficient multi-threaded program Acknowledgements to calculate EHH-based scans for positive selection. Mol Biol Evol. The research leading to these results received funding from a Ministerio de 2014;31:2824–7. Ciencia e Innovación AGL 2010-15903 Grant from the Spanish Government, 21. Ferrer-Admetlla A, Liang M, Korneliussen T, Nielsen R. On detecting from the European Union’s Seventh Framework Programme for research, incomplete soft or hard selective sweeps using haplotype structure. Mol technological development and demonstration under Grant Agreement No. Biol Evol. 2014;31:1275–91. González‑Rodríguez et al. Genet Sel Evol (2016) 48:81 Page 12 of 12 22. Garud NR, Messer PW, Buzbas EO, Petrov DA. Recent selective sweeps in 46. Stella A, Ajmone-Marsan P, Lazzari B, Boettcher P. Identification of selec‑ North American Drosophila melanogaster show signatures of soft sweeps. tion signatures in cattle breeds selected for dairy production. Genetics. PLoS Genet. 2015;11:e1005004. 2010;185:1451–61. 23. Chen H, Patterson N, Reich D. Population differentiation as a test for 47. Rothammer S, Seichter D, Förster M, Medugorac I. A genomic scan for selective sweeps. Genome Res. 2010;20:393–402. signatures of differential artificial selection in ten cattle breeds. BMC 24. Teo YY, Fry AE, Bhattacharya K, Small KS, Kwiatkowski DP. Genome- Genomics. 2013;14:908. wide comparisons of variation in linkage disequilibrium. Genome Res. 48. Burt DW. The cattle genome reveals its secrets. J Biol. 2009;8:36. 2009;19:1849–60. 49. Randhawa IAS, Khatkar MS, Thomson PC, Raadsma HW. Composite selec‑ 25. Ong R, Teo YY. VarLD: a program for quantifying variation in link‑ tion signals can localize the trait specific genomic regions in multi-breed age disequilibrium patterns between populations. Bioinformatics. populations of cattle and sheep. BMC Genet. 2014;15:34. 2010;26:1269–70. 50. Bomba L, Nicolazzi EL, Milanesi M, Negrini R, Mancini G, Biscarini F, et al. 26. Kaiser HF. The varimax criterion for analytic rotation in factor analysis. Relative extended haplotype homozygosity signals across breeds reveal Psychometrika. 1958;23:187–200. dairy and beef specific signatures of selection. Genet Sel Evol. 2015;47:25. 27. R Development Core Team. R: a language and environment for statistical 51. Bosma M, Sparks LM, Hooiveld GJ, Jorgensen JA, Houten SM, Schrauwen computing. Vienna: R Foundation for Statistical Computing; 2008. ISBN P, et al. Overexpression of PLIN5 in skeletal muscle promotes oxidative 3-900051-07-0. http://www.R-project.org. gene expression and intramyocellular lipid content without compromis‑ 28. Wang J, Duncan D, Shi Z, Zhang B. WEB-based GEne SeT AnaLysis Toolkit ing insulin sensitivity. Biochim Biophys Acta. 2013;1831:844–52. (WebGestalt). Nucleic Acids Res. 2013;41:W77–83. 52. Puig-Oliveras A, Ramayo-Caldas Y, Corominas J, Estellé J, Pérez-Montarelo 29. Sabeti PC, Schaffner SF, Fry B, Lohmueller J, Varilly P, Shamovsky O, et al. D, Hudson NJ, et al. Differences in muscle composition transcriptome Positive natural selection in human lineage. Science. 2006;312:1614–20. among pigs phenotypically extreme for fatty acid composition. PLoS 30. Simianer H, Ma Y, Qanbari S. Statistical problems in livestock popula‑ One. 2014;9:e99720. tion genomics. In: Proceedings, 10th World Congress of genetics 53. Tsang WY, Wang L, Chen Z, Sanchez I, Dynlacht BD. SCAPER, a novel cyclin applied to livestock production. Vancouver; 2014. https://asas.org/docs/ A-interacting protein that regulates cell cycle progression. J Cell Biol. default-source/wcgalp-proceedings-oral/202_paper_10373_manu‑ 2007;178:621–33. script_1346_0.pdf?sfvrsn=2. 54. Mashima R, Okuyama T. The role of lipoxygenasas in pathophysiology; 31. Kemper KE, Saxton SJ, Bolormaa S, Hayes BJ, Goddard ME. Selection for new insights and future perspectives. Redox Biol. 2015;6:297–310. complex traits leaves little or no classic signatures of selection. BMC 55. Randhawa IAS, Khatkar MS, Thomson PC, Raadsma HW. Composite selec‑ Genomics. 2014;15:246. tion signals for complex traits exemplified through bovine stature using 32. Cole JB, VanRaden PM, O’Connell JR, Van Tassell CP, Sonstegard TS, Schna‑ multibreed cohorts of European and African Bos taurus. G3 (Bethesda). bel RD, et al. Distribution and location of genetic effects for dairy traits. J 2015;5:1391–401. Dairy Sci. 2009;92:2931–46. 56. Noyes H, Brass A, Obara I, Anderson S, Archibald AL, Bradley DG, et al. 33. Hayes BJ, Pryce J, Chamberlain AJ, Bowman PJ, Goddard ME. Genetic Genetic and expression analysis of cattle identifies candidate genes in architecture of complex traits and accuracy of genomic prediction: coat pathways responding to Trypanosoma congolense infection. Proc Natl colour, milk-fat percentage and type in Holstein cattle as contrasting Acad Sci USA. 2011;108:9304–9. model traits. PLoS Genet. 2010;23:e1001139. 57. Pareek CS, Michno J, Smoczynski R, Tyburski J, Golebiewski M, Piechocki 34. Pimentel EG, Erbe M, König S, Simianer H. Genome partitioning of K, et al. Identification of predicted genes expressed differentially in genetic variation for milk production and composition traits in Holstein pituitary gland tissue of young growing bulls revealed by the cDNA-AFLP cattle. Front Genet. 2011;2:19. technique. Czech J Anim Sci. 2013;58:147–58. 35. Shao H, Burrage LC, Sinasac DS, Hill AE, Ernest SR, O’Brien W, et al. Genetic 58. Yao C, Spurlock DM, Armentano LE, Page CD Jr, VandeHaar MJ, Bickhart architecture of complex traits: large phenotypic effects and pervasive DM, et al. Random Forests approach for identifying additive and epistatic epistasis. Proc Natl Acad Sci USA. 2008;105:19910–4. single nucleotide polymorphisms associated with residual feed intake in 36. Dunner S, Miranda ME, Amigues Y, Cañon J, Georges M, Hanset R, et al. dairy cattle. J Dairy Sci. 2013;96:6716–29. Haplotype diversity of the myostatin gene among beef cattle breeds. 59. Flori L, Fritz S, Jaffrézic F, Boussaha M, Gut I, Heath S, et al. The genome Genet Sel Evol. 2003;35:103–18. response to artificial selection: a case study in dairy cattle. PLoS One. 37. González-Rodríguez A, Toro MA, Varona L, Carabaño MJ, Cañas-Álvarez JJ, 2009;4:e6595. Altarriba J, et al. Genome-wide analysis of genetic diversity in Autoch‑ 60. Fellius M, Beerling ML, Buchanan DS, Theunissen B, Kollmees PA, Lenstra thonous Spanish populations of beef cattle. In: Proceedings of the 10th JA. On the history of cattle genetic resources. Diversity. 2014;6:705–50. World Congress of genetics applied to livestock production. Vancouver; 2014. https://asas.org/docs/default-source/wcgalp-proceedings- oral/255_paper_9607_manuscript_1524_0.pdf?sfvrsn=2. 38. Gutierrez-Gil B, Arranz JJ, Wiener P. An interpretive review of selective sweep studies in Bos taurus cattle populations: identification of unique and shared selection signals across breeds. Front Genet. 2015;6:167. 39. Wiener P, Gutierrez-Gil B. Assessment of selection mapping near the myostatin gene (GDF-8) in cattle. Anim Genet. 2009;40:598–608. 40. Boitard S, Rocha D. Detection of signatures of selective sweeps in the Submit your next manuscript to BioMed Central Blonde d’Aquitaine cattle breed. Anim Genet. 2013;44:579–83. 41. Druet T, Pérez-Pardal L, Charlier C, Gautier M. Identification of large and we will help you at every step: selective sweeps associated with major genes in cattle. Anim Genet. • We accept pre-submission inquiries 2013;44:758–62. 42. Werth LA, Hawkins GA, Eggen A, Petit E, Elduque C, Kreigesmann B, et al. • Our selector tool helps you to find the most relevant journal Rapid communication: melanocyte stimulating hormone receptor (MC1R) • We provide round the clock customer support maps to bovine chromosome 18. J Anim Sci. 1996;74:262. • Convenient online submission 43. Randhawa IAS, Khatkar MS, Thomson PC, Raadsma HW. A meta-assembly of selection signatures in cCattle. PLoS One. 2016;11:e0152013. • Thorough peer review 44. Gautier M. Genome-wide scan for adaptive divergence and association • Inclusion in PubMed and all major indexing services with population-specific covariates. Genetics. 2015;201:1555–79. • Maximum visibility for your research 45. Saldana-Caboverde A, Kos L. Roles of endothelin signaling in mel‑ anocyte development and melanoma. Pigment Cell Melanoma Res. Submit your manuscript at 2010;23:160–70. www.biomedcentral.com/submit