A whole genome Bayesian scan for adaptive genetic divergence in West African cattle
© Gautier et al. 2009
Received: 29 May 2009
Accepted: 21 November 2009
Published: 21 November 2009
Skip to main content
© Gautier et al. 2009
Received: 29 May 2009
Accepted: 21 November 2009
Published: 21 November 2009
The recent settlement of cattle in West Africa after several waves of migration from remote centres of domestication has imposed dramatic changes in their environmental conditions, in particular through exposure to new pathogens. West African cattle populations thus represent an appealing model to unravel the genome response to adaptation to tropical conditions. The purpose of this study was to identify footprints of adaptive selection at the whole genome level in a newly collected data set comprising 36,320 SNPs genotyped in 9 West African cattle populations.
After a detailed analysis of population structure, we performed a scan for SNP differentiation via a previously proposed Bayesian procedure including extensions to improve the detection of loci under selection. Based on these results we identified 53 genomic regions and 42 strong candidate genes. Their physiological functions were mainly related to immune response (MHC region which was found under strong balancing selection, CD79A, CXCR4, DLK1, RFX3, SEMA4A, TICAM1 and TRIM21), nervous system (NEUROD6, OLFM2, MAGI1, SEMA4A and HTR4) and skin and hair properties (EDNRB, TRSP1 and KRTAP8-1).
The main possible underlying selective pressures may be related to climatic conditions but also to the host response to pathogens such as Trypanosoma(sp). Overall, these results might open the way towards the identification of important variants involved in adaptation to tropical conditions and in particular to resistance to tropical infectious diseases.
Cattle are still playing a major role in Africa for food supply, to generate income and draught power or for ceremonial purposes. Archaeological, historical and anthropological evidence combined with recent genetic data  have provided insights into the complex origins of present day West-African cattle diversity. Indeed, although their wild ancestor Bos primigenius was not native to sub-Saharan Africa, West African cattle populations are representative of both shorthorn (Bos taurus brachyceros) and longhorn (Bos taurus longifrons) humpless taurines, humped zebus (Bos indicus) and zebu/taurine hybrid cattle. This early suggested that West African cattle has originated from several successive and recent colonization events [2, 3]. Briefly, shorthorn taurines were introduced from the Middle-East and possibly North Africa around 4,000 years BP [3, 4] while longhorn taurine probably arrived at an earlier period (5,000 years BP) following different migration routes . Although, zebu cattle first penetrated through the Horn of Africa in the late 2nd millennium BC, the major wave of indicine introgression really started with the Arab settlements along the East Coast of Africa from the end of the 7th century AD. Zebu cattle spread even more recently over West Africa with movements of pastoralist people such as the Fulani .
West African cattle populations thus represent an appealing model to unravel the genome response to adaptation to tropical conditions. The purpose of this study was to perform a whole genome scan for footprints of adaptive selection based on a newly collected genotyping data set containing 36,320 SNPs genotyped on 9 West African cattle populations from different bovine sub-species and agro-ecological areas (Figure 1). In particular, we sampled populations on both side of the tsetse infested zone. Based on this large data set, we first carried out a detailed analysis of the genetic structure of these populations. We next performed a scan for differentiation among SNPs under a full Bayesian framework [7, 8] for which we proposed additional extensions to improve the detection of loci under selection. We then annotated regions containing SNPs subjected to selection, adopting a systems biology approach to highlight the main underlying physiological functions.
In total, 437 individuals belonging to 12 different cattle populations, nine of which originating from West Africa (Figure 1), were genotyped for the BovineSNP50 chip assay containing 54,001 SNPs mainly derived from sequences available in European cattle breeds . Among the autosomal SNPs that passed Quality Control analyses (see Methods), we retained the 36,320 SNPs polymorphic (Minor Allele Frequency or MAF above 0.01) in at least one West-African taurine and one West African zebu populations. As expected and shown in Figure S1A (additional file 1), this SNP selection procedure, leading mainly to the elimination of SNPs from European origin, resulted in a relative increase in calculated heterozygosity for all but outgroup taurine populations (Aubrac or AUB and Oulmès Zaer or OUL). In addition, as revealed by the distribution of the MAF for the different populations (Figure S1B in additional file 1), the remaining differences in heterozygosity among populations were mainly explained by the proportion of SNPs with rare variants (MAF < 0.05). In particular, marker polymorphism remained clearly lower in West African taurine populations (ND1, ND2, LAG, SOM and BAO), the proportion of SNPs with a MAF < 0.05 being the highest (36%) in LAG which also displayed the lowest average heterozygosity (0.118). These observations are in agreement with previous studies based on microsatellite markers on the same populations  and might essentially be explained by the recent demographic history of West African taurine cattle characterized by a marked isolation associated with a strong effective population size decrease (e.g. ). Although West African zebus (ZFU and ZCH) and hybrids (BOR and KUR) displayed more genetic diversity than West African taurines in our study (Figure S1 in additional file 1), ZFU, ZCH and ZMA levels of polymorphism were clearly lower than, and BOR and KUR ones similar to the taurine outgroup (AUB and OUL) ones (Figure S1A in additional file 1). Nevertheless, studies based on microsatellite data (e.g. ) revealed that heterozygosities of ZCH, ZFU, BOR and KUR were higher than European breeds' ones (AUB) which were themselves similar to ZMA. The lack of zebu specific sequences in the SNP discovery process might explain these apparent discrepancies.
As reported in Table S1 (additional file 2), within population F IS were almost null (< 0.006) for ZFU and BOR or very close to zero (from -0.0121 to 0.0198) for ND1, BAO, KUR and ZCH. The slightly positive F IS (0.0421) observed in SOM might result from the sampling area being apart the border between Togo and Benin. The positive F IS (0.0414) of similar magnitude observed in LAG might rather originate from a higher level of inbreeding as suggested by the extent of Linkage Disequilibrium (LD) within this population (see below). Conversely, for ND2, the negative F IS (-0.108) might be explained by recent outcrossing events. Indeed ND2 individuals derived from founders recently collected in different and distant Ivorian villages (see Methods). This hypothesis is also in agreement with their higher dispersion compared to other West African taurines in the PCA (Figure 2) and a higher level of overall extent of LD (see below). Nevertheless, the global F IS coefficient remained close to zero (0.0103) while the global F ST and F IT were respectively equal to 0.132 and 0.141. Such a level of differentiation among the West African populations was in agreement with the one computed based on microsatellite markers  suggesting a low impact of the SNP ascertainment procedure on these estimates.
We further studied the differentiation among SNPs and populations based on a Bayesian hierarchical model derived from  and presented in more details in the Methods section. Among the factors influencing population differences in allele counts (differentiation), the model aims at distinguishing locus-specific from population-specific factors, such as migration or drift . Because most SNPs originated prior to breed differentiation, the effect of different mutation rates across loci might be negligible and the SNP-specific effect in the model is mainly related to selection. Notice that the model assumed the F IS to be null within each population which appeared sound given the characteristics of our data (see Figure S3 in additional file 1 and Methods). We also ignored ascertainment biases originating both from the SNP discovery process and our selection of SNPs "informative" in West African populations. These two different sources are expected to favour polymorphic SNPs of ancient origin (co-segregating in European and African cattle breeds) and thus modifies to some extent the distribution of allele frequencies in the gene pool (x following our model notations) . Recently, efficient algorithms were proposed to account for ascertainment biases , providing it can be modelled (which is difficult in our case at least for the first identified source). Nevertheless, for the purpose of this study we were mainly concerned with possible biases introduced in the locus-specific effect estimation which might not be of importance . At least for our second source of ascertainment bias, simulated data allowed to confirm this statement (Figure S4 in additional file 1) since estimation of x was found robust to departure from the prior distribution assumed. Similarly, although more variable, locus-specific effect estimates were highly correlated with their corresponding simulated values (r > 0.8).
FDR (FNR) for different thresholds on BF.
BF threshold (in dB units)
P = 0.1%
P = 1%
P = 10%
P = 18.5%
Relationship between selection coefficient (s) and locus effect α i
Number of loci simulated
Average FST (sd)
FST range (median)
Regions under selection identified at the 5% local FDR (q-value) thresholds.
Position (size) in Mb
Peak Position in Mb
BF value of the SNP at the peak position (q-value of the smoothed signal)
Gene at or close to the peak
The 36,320 SNPs were representative of about 7,177 different genes corresponding to approximately one third of the total number expected in the genome. For each of these, Table S5 (additional file 2) summarizes the number of SNPs they contained and results from the differentiation analysis. Individual SNP information made it nevertheless difficult to propose a list of genes subjected to selection and this strategy might be more sensitive to possible false positives introduced by hierarchical structure among the populations under study . Alternatively, a great proportion of such genes could include SNPs appearing as significantly differentiated as a result of hitchhiking with a favourable variant located nearby (see above). Hence, among the 191 (853) genes containing or close to at least one SNP with a BF > 30 (BF > 15), 103 (237) mapped within one of the 53 regions identified above. For functional and network analyses, we thus decided to mainly concentrate on the 42 genes which were located for 25 of them (represented in Table S5 in additional file 2), under the peak of each region, or for the remaining 17 less than 50 kb from the peak of each region (Table 3) (and not represented by any of the 36,320 SNPs surveyed). These 42 genes were viewed as the strongest candidates underlying the observed footprints of selection. Among these, 37 genes were eligible for network analysis, 29 being also eligible for functional analysis. The five remaining ones (CCDC46, CTXN3, KRTAP8-1, RNF180, TEKT3) were not eligible for any network or functional analyses.
Several candidate genes (such as CD79A, CXCR4, DLK1, RFX3, SEMA4A, TICAM1 and TRIM21) belonging to N directly participate in antigen recognition, a key process underlying the development of immune response. For instance, TRIM21, a newly identified type of IgG Fc receptors  participates to both the host anti-viral response through the modulation of IRF3 functions  and mediation of autoimmune diseases. Footprints of selection were recently reported in primates within another member of the same tripartite motif (TRIM) family of proteins which has antiretroviral properties [29, 30]. TICAM1 is a Toll like receptor (TLR) adapter and thus plays a role in innate immunity. It possesses, in particular, a high IFN type I inducing activity during viral infection . DLK1, a member of the epidermal growth factor-like gene family, has been shown in mouse to be involved in B lymphocytes differentiation and function . More specifically, other genes (such as CD79A, CXCR4, RFX3 and SEMA4A) are related to antigen presentation to T cells in the context of MHC. Although no MHC genes were present in the list analyzed, the corresponding region (#49 in Table 3) and a SNP (with BF = 17 and F ST = 0.005) within the bovine HLA-DMB ortholog (Table S4 in additional file 2) were found under strong balancing selection. Among the molecules participating to antigen presentation, the chemokine receptor CXCR4 is co-recruited with CCR5 in the immunological synapse, and it participates to modulation of T cells response . These two molecules were shown to act as co-receptors for HIV entry in human cells . A strong signature of balancing selection was reported in the 5' cis-regulatory region of CCR5 in human , this region displaying as well a selective sweep in chimpanzee . We also recently observed a significant signal of selection on a microsatellite located less than 1 Mb from CXCR4 in West African taurine  and defining the boundary of a QTL underlying the trypanotolerance trait in cattle . Similarly, RFX3 (43.235-43.563 Mb on BTA08) is located within another QTL identified in this latter study . This member of the regulatory factor X (RFX) gene family participates with CIITA in the regulation of MHC genes [38–40]. CD79A via CD79α-CD79β heterodimer located on the surface of B cells is required for antigen presentation through MHC class II . SEMA4A a semaphorin expressed on dendritic cells and B cells is involved in T cell priming and Th1 differentiation through its interaction with Tim2 expressed on T cells .
Although heat stress response may have imposed immunological parameter modifications , infectious and parasitical diseases are likely to have represented selective pressures acting on genes involved in such functions . As previously mentioned, trypanosomiasis represents one of the well known examples of such selective pressure. Nevertheless, the importance of other diseases such as tick-borne diseases (anaplasmosis, babesiosis, cowdriosis) or anthrax should not be overlooked. For instance, anthrax is hyperendemic (especially in Chad, Togo and Ivory Coast) or endemic in West Africa http://www.vetmed.lsu.edu/whocc/mp_world.htm. Notice that the ANTRXR2, a receptor of the anthrax toxin was among our candidate genes (see network N, Figure 6).
In addition to the exposure to new pathogens, the recent settlement of cattle to Sub Saharan Africa might have imposed a dramatic change in their environmental conditions . Hence, cattle (in particular taurine) were exposed to warmer temperature, to a thermal amplitude decrease, to a day length change, to differences in solar radiation and also to new feeding conditions. As a result, adaptations needed modifications in global appearance (such as morphology or coat color), in body temperature regulation, in circadian clock and also in reproductive behaviour and abilities. Such selective pressures might partly explain the physiological functions related to nervous system (neurogenesis and eyesight). Alternatively, neurotropism of several parasites such as Trypanosoma(sp) have been demonstrated at least in human. Some of the candidate genes related to nervous system development and function belonged to N3 (e.g. NEUROD6 and OLFM2) or N (e.g. MAGI1, SEMA4A, HTR4 and EDNRB). NEUROD6 is a member of the neurogenic differentiation transcription factor family  while OLFM2 is a secreted glycoprotein belonging to the noelins family which modulate the timing of neuronal differentiation during development . Interestingly, MAGI1 is a scaffolding protein present in tight junction of epithelial cells , some transcripts of the corresponding gene being only expressed in brain. Moreover it is located within a QTL underlying the trypanotolerance trait . Besides their role in immune response (see above), semaphorins represented by SEMA4A were initially characterized for their role in the guidance of axonal migration during neuronal development. Mutations in the conserved semaphorin domain of SEMA4A are associated with two retinal degenerative diseases, retinis pigmentosa and cone-rod dystrophy . Of more particular interest is HTR4 which is a serotonin receptor participating to the serotonergic system. Recently, a mutation within HTR2A, a receptor of the same family was reported associated with the chronic fatigue syndrome corresponding to a dysregulation of hypothalamic-pituitary-adrenal (HPA) axis and serotonergic system . A reduced responsiveness of the HPA axis after corticotropin-releasing hormone (CRH) challenge was observed in Boran cattle infected with Trypanosoma congolense . Taken together, a possible involvement of HTR4 in trypanotolerance/susceptibility in cattle might be hypothesized though we cannot exclude an implication of this serotonin receptor in more general aspects of circadian rythmicity .
Although EDNRB plays an essential role in the development of enteric neurons, it is also involved in the development of epidermal melanocytes, both cell lineages being derived from the neural crest . As an illustration, several EDNRB mutations are associated with an auditory pigmentary syndrome caused by the absence of melanocytes . However, given the variety of cattle populations surveyed (Figure 1), the footprint of selection observed in the EDNRB region is more likely related to differences in coat color and more particularly to the spotting pattern (SOM, BAO and LAG) or the white color (KUR and ZFU). Indeed, EDNRB is also referred to as piebald or S locus in mouse  and a null mutation induce a white coat color in rat . In the Hereford cattle breed, although not yet fully characterized, locus S seems responsible for white spotting pattern . In addition to EDNRB involved in skin and hair pigmentation, other genes such as TRSP1 which belonged to N and KRTAP8-1 which was not eligible for network analysis play a role in hair development in human and mouse. Indeed, TRPS1 is a zinc finger transcription factor implicated in growth and trichosis, some of its variants being associated with trichorhinophalangeal syndromes [54, 55]. It is down-regulated in patients affected by hypertrichosis and in a mouse model of hypertrichosis . KRTAP8-1, a keratin associated protein, plays a role in the formation of hair shafts .
The skin color and thickness, the hair size and sleekness and the number of sweat glands directly influence thermo-resistance of cattle living in the tropics . Compared to European taurine breeds, zebus which have a higher density of sweat glands and a smoother and shinier hair coat were reported to better regulate body temperature and more efficiently maintain cellular function during heat stress . Similarly, slick-haired Holstein cattle are more able to regulate their body temperature than wild-type . Some skin and hair properties might also confer higher resistance to tick infestation [60, 61]. Interestingly, some genes involved in hair development also participate in cornification. Hence the observed footprint of selection within the cluster of keratin associated proteins (represented in the candidate genes list by KRTPA8-1, the closest to the peak) might thus also be related to horn morphology differences . A striking example is the floating horns of KUR which could represent original adaptation facilitating swimming in their swamp living area . Note that the identified region is also localized less than 5 Mb from the bovine polled locus 
West African cattle provide a valuable resource to better understand the genomic response to selective pressure arising in tropical conditions. This large-scale whole genome Bayesian scan for adaptive differentiation in populations representative of the current breed diversity allowed us to identify candidate genes involved in several key physiological functions. These results might open the way towards the identification of variants underlying these footprints of selection, in particular those involved in the resistance to tropical infectious diseases in cattle but also in other mammals such as human populations subjected to similar pressures.
Figure 1 provides information on the origin of the nine different West African populations sampled in our study and representatives of the three main types of cattle found in this area. Longhorn taurines were represented by N'Dama individuals from two distinct origins (ND1 and ND2). ND1 individuals (n = 14) originated from the Gaoua herd constituted by acquisitions coming from different villages in the Pays Lobi (Burkina Faso). ND2 (n = 17) samples were collected in the Samandeni herd (Burkina Faso) which was constituted in the early eighties by acquisitions from the Marahoue ranch (Ivory Coast) and originated from northwestern villages from Ivory Coast where zebu had not been introduced . Shorthorn taurines were represented by three different breeds: Baoulé (BAO) samples (n = 29) originated, as ND1, from the Gaoua herd , Somba (SOM) individuals (n = 44) were sampled in the breed birthplace across the border near Nadoba (Togo) and Boukoumbe (Benin), and Lagune (LAG) individuals (n = 44) were sampled in the Porto Novo district (Benin) . West African zebus were represented by two populations: Sudanese Fulani or ZFU (n = 43) which was sampled in the Malanville region in Benin  and Choah zebus or ZCH (n = 59) which was sampled in the Bol district (Chad) . This latter population was initially sub-divided into two breeds (M'Bororo and Choah Zebus), however, as previously shown  and confirmed in our analysis with a dense marker set, no clear differentiation appeared among these. Finally, two hybrid populations were considered: Borgou or BOR (n = 45) which is a stabilized crossbred between shorthorn taurines (LAG or BAO) and was sampled in its region of origin around Parakou in Benin  and Kuri or KUR (n = 47) which was sampled from different islands of Lake Chad around Bol . This latter breed was sometimes referred to as a particular longhorn taurine , however, several molecular analyses showed a high level of Zebu admixture in it [4, 10]. Moreover, unlike West African taurines, KUR is known to be trypano-susceptible. Three other breeds were considered as outgroups for the detailed genetic structure analysis: ZMA (n = 35) sampled in the Madagascar Island  and representing pure zebu , AUB (n = 20) sampled in the birthplace of the breed  and representing European taurine, and OUL (n = 40) sampled in the North of Morocco to represent North African taurine cattle.
Individuals were genotyped on the Illumina BovineSNP50 chip assay  at the Centre National de Génotypage (CNG) platform (Evry, France) using standard procedures http://www.illumina.com. Among the 54,001 SNPs included in the chip, only the 51,581 mapping to a bovine autosome on the latest bovine genome assembly Btau_4.0  were retained for further analysis. To limit ascertainment bias favouring SNPs from European origin, we subsequently discarded 13,786 SNPs (~25%) which were not polymorphic (MAF > 0.01) in at least one West African taurine (ND1, ND2, LAG, SOM or BAO) and one West African zebu (ZCH or ZFU). Among the remaining ones, 1,422 SNPs which were genotyped on less than 85% of the individuals from at least one of the nine West African breeds, were also eliminated. An exact test for Hardy-Weinberg Equilibrium (HWE)  was subsequently carried out within each breed separately. Based on the obtained p-values, q-values  were estimated for each SNP using the R package qvalue http://cran.r-project.org/web/packages/qvalue/index.html. Fifty three SNPs with a q-value < 0.01 in at least one breed were then discarded from further analysis. In total, 36,320 SNPs were finally considered for the study leading to an average marker density of 1 SNP every 70 kb over the genome (Table S8 in additional file 2). Moreover, as shown in Figure S6 (additional file 1) and detailed in Table S8 (additional file 2), the genome coverage was homogeneous with a median distance between consecutive SNPs equal to 47.54 kb. Few large gaps between SNPs were present since the 95th (99th) percentile of this distance was 189 kb (385 kb), the largest gap localized on BTA10 being 2 Mb long. Conversely, less than 0.5% of the distances between successive SNPs were shorter than 20 kb. Genotyping data are given in additional file 3.
ASD were computed for each pair of individuals using all available SNP information by a simple counting algorithm: for a given pair of individuals i and j, ASD was defined as 1-x ij where x ij represents the proportion of alleles alike in state averaged over all genotyped SNPs. A neighbor-joining tree was computed based on the resulting distance matrix using the R package APE . We subsequently performed a PCA based on all available SNP information using the SMARTPCA software package . As suggested by Patterson et al. , we performed LD correction by replacing individual SNP values with the residuals from a multivariate regression without intercept on the two preceding SNPs on the map, providing they were less than 200 kb apart. Results were further visualized using functionalities available in the R package ade4 . The global F-statistics F IT , F ST and F IS were estimated respectively in the form of F, θ and f  using the program GENEPOP 4.0 . GENEPOP 4.0 was also used to estimate diversity for each locus and population both within individuals and among individuals within a population. The within breed F IS was derived from the average of these two quantities over all the SNPs. In order to evaluate the reliability of the F IS estimates we computed the mean and standard deviation over 10,000 samples of 5,000 randomly chosen SNPs. F ST statistics between populations were estimated using both SMARTPCA  and GENEPOP  which are based on two different models of population divergence.
Individual genotyping data were modelled according to the reparameterized extension, recently proposed by Riebler et al. , of the initial Bayesian hierarchical model developed by Beaumont and Balding . Briefly and in the bi-allelic SNP case, let a ij be the observed reference allele (defined arbitrarily) count in population j = 1, ..., J at locus i = 1, ..., I. The conditional distribution given the true (unobserved) allele frequency p ij in that population at that locus is assumed to be binomial with parameters n ij (twice the number of genotyped individuals in population j at locus i) and p ij : a ij | p ij ~Bin(n ij , p ij ) (1). Note that 1) implicitly assumes populations are in HWE or their respective inbreeding coefficients (F IS ) are null. Non null F IS could be taken into account in the model by considering instead that the three possible genotypes are drawn from a multinomial distribution with parameters corresponding to the number of individuals genotyped and genotype probabilities: (1-F IS j )p ij 2 + F IS j p ij ; 2(1-F IS j )p ij (1-p ij ) and (1-F IS j )(1-p ij )2 + F IS j (1-p ij ) . Nevertheless, for co-dominant markers such as SNP and given the range of F IS values in our study (see Results), the binomial distribution seems to be reasonable (Figure S3 in additional file 1). The second step assumes that the true allele frequencies p ij are themselves sampled from a Beta distribution: p ij | λ ij , x i ~Beta(λ ij x i , λ ij (1-x i )) (2). This distribution relies on an infinite Wright island model involving mutation, drift and migration at its equilibrium state [76, 77]. Under this model, x i might be interpreted as the frequency of the chosen reference allele at locus i in the gene pool from which each allele frequency p ij is sampled. The scaling parameter λ ij reflects the gene flow from the gene pool to population j. Under Wright's model, this parameter is specific to each population j but remains homogeneous over loci so that any departure from that property may indicate that locus i is no longer neutral. Note that under the Beta distribution, the first two moments can be simply expressed as E(p ij ) = x i and Var(p ij ) = x i (1-x i )(1+λ ij )-1. Actually (1+λ ij )-1 can be identified as a F ST ij coefficient . The next level of the model consists of specifying the distributions of λ ij (or F ST ij ) and of x i . These frequencies are nuisance parameters and uncertainty about them will be taken into account by assigning to them non informative prior such as a Beta(1,1) distribution (e.g. ). Following Beaumont and Balding , the λ ij 's are modelled via a linear model on the logistic transformation of the F ST ij . Since F ST ij /(1-F ST ij ) = 1/λ ij , we can write this model in terms of η ij = log(F ST ij /(1-F ST ij )) = -log(λ ij ) as: η ij = α i + β j + γ ij where α i is a locus effect, β j a population effect and γ ij an error term corresponding to a departure of the logit of F ST ij from the additive decomposition. For theoretical and computational reasons, it is more efficient to consider the previous decomposition under a hierarchical Bayesian structure than under its basic form, and implement it as follows: η ij | α i , β j , σ γ 2~ iid N(α i + β j , σ γ 2) with β j | μ β , σ β 2~ iid N(μ β , σ β 2) and α i | σ α 2~ iid N(0, σ α 2). Introduction of additional levels by defining priors on the variance components σ α 2, σ β 2 and σ γ 2 is straightforward but adds marginal gain in the estimates and requires high supplementary computational costs (Gautier and Foulley, unpublished data). We thus gave for these components the known values proposed initially by Beaumont and Balding (2004): σ α 2 = 1, σ β 2 = 3.24 and σ γ 2 = 0.25. The prior mean μ for the β j was also taken as -2.0. Recently, Riebler et al.  introduced in the above logistic model an auxiliary indicator variable δ i attached to each locus specifying whether it can be regarded as selected (δ i = 1) or neutral (δ i = 0). They demonstrated the efficiency of this approach for improving the power of the statistical procedure, particularly for the detection of loci under balancing selection. Under this reparameterized model, the previous parameters α i are written as: α i = δ i α i * with α i * | σ α *2~ iid N(0, σ α *2). The model further assumes a Bernoulli distribution for the indicator δ i variable with parameter P: δ i | P ~Bin(1, P). P is itself assumed to be Beta distributed to take into account uncertainty on this crucial parameter. Here we took P ~Beta(0.2,1.8) , thus assuming that a priori 10% of the loci are on average under selection, but within a very large range of possible values as the 95% credibility interval goes from almost 0 to 65%. Note that the value σ α *2 can be derived from σ α 2 since by construction σ α 2 = E(P) σ α *2 (hence σ α *2 = 10).
The posterior distributions of the different parameters of interest were estimated via MCMC procedures as previously described  from 10,000 post burn-in samples (with a burn-in period of 3,000 iterations) and a thinning interval of k = 25. Convergence was checked using standard criterion. Parameter estimates for the F ST ij were taken as the median of the posterior distribution.
where π(P) is the density of the distribution of P defined above. To make interpretation of the BF i easier, we expressed them in deciban units (dB) ie dB i = 10log 10(BF i ) so that 10 dB corresponds to an odds ratio of 10, 20 dB to an odds ratio of 100, 30 dB to an odds ratio of 1,000 and so forth. We then applied the Jeffreys' rule  which quantifies the strength of evidence (here to consider a SNP as being under selection) based on BF using the following scale: "strong evidence" when 10<dB i <15, "very strong evidence" when 15<dB i <20 and "decisive evidence" when dB i >20. Notice that the BF captures the change in the odds in favor of δ i = 1 as we move from the prior to the posterior. In that way, it makes this criterion independent of the prior distribution of δ i and thus guarantees some robustness.
Following Foll and Gaggiotti , four data sets each comprising 50,000 (unlinked) SNPs with respectively 9300 (P = 0.186), 5000 (P = 0.1), 500 (P = 0.01) and 50 (P = 0.001) under selection (α i ≠ 0) were simulated under the inference model. To mimic characteristics of our samples, we considered 9 populations with parameters β j equal to the ones estimated on the real data (mean of the corresponding posterior distribution). For each pair of simulated locus/population, we subsequently simulated the η ij parameter by adding to the corresponding β j , a value of γ ij sampled from a Gaussian distribution (γ ij ~N(0,0.25)) and a value for α i equal to 0 if the locus was neutral or sampled from a Gaussian distribution (α i ~N(0,10)) if the locus was under selection. The frequency x i in the ancestral population (or migrant gene pool) was sampled from a Beta distribution with parameters equal to 0.7. This is equivalent to assume the ancestral population was at equilibrium with 4Nμ = 0.7 where N is the effective population size and μ the mutation rate . Based on these simulated values for x i and η ij we subsequently derived the p ij parameter of the binomial distribution from which the simulated observed a ij was sampled (taking n ij = 80 for all i and j). Only SNPs displaying MAF>0.01 in at least one simulated West-African taurine and one West-African zebu (identified according to the simulated β j ) were conserved which led to the exclusion of few SNPs (displaying α i too high in general). To a first approximation, this also allowed the mimicking of the ascertainment bias introduced by the SNP selection criterion. The posterior distributions of the different parameters of interest were then estimated as described above from 1,000 post burn-in samples (with a burn-in period of 3,000 iterations) and a thinning interval of k = 25. BF i for each simulated loci were then computed (see above) and used to compute the (true) FDR and FNR for a given threshold. To investigate relationships between the coefficient of selection and the magnitude of the locus effect α i we performed simulations under a simplified version (without mutation) of the Wright-Fisher island model with selection proposed by Beaumont and Balding . Briefly, 10 random mating populations with a constant size of 1,000 individuals and deriving from the same ancestral population were reproduced for 500 generations assuming a global F ST of 0.15. In total, 5,000 neutral SNPs and 500 SNPs for each investigated coefficient and mode of selection were simulated. Initial allele frequencies were drawn from a Beta distribution with both parameters equal to 0.7. At the end of the process, SNPs displaying a MAF <0.01 in at least 2 populations were discarded. Locus specific F ST were then computed following Weir and Cockerham . Following Foll and Gaggiotti , for each type of selection considered, the corresponding average F ST was compared to the one obtained under neutrality (α i = 0) to estimate the average SNP effect.
Because the annotation of the bovine genome is still sparse, the gene content information was derived from the TransMap cross-species alignments available in the UCSC Genome Browser http://genome.ucsc.edu/. For closer evolutionary distances, the alignments are created using syntenically filtered BLASTZ alignment chains, resulting in a prediction of the orthologous genes in cow. In total, 46,598 different RefSeq identifiers were anchored to the Btau_4.0 bovine genome assembly http://genome.ucsc.edu/. Considering that most consecutive SNPs on the map were separated by more than 20 kb and the relatively high level of LD at shorter distance (see Results and ) a SNP was considered as representative of a gene if it was localized within the boundary positions of the gene that extended by 15 kb upstream and downstream. According to this criterion, 15,360 out of the 36,320 SNPs were representative of 17,190 different TransMap RefSeq identifiers. Subsequent annotation and analyses were carried out with the Ingenuity Pathway Analysis (IPA) software v7.0 (Ingenuity Systems Inc., USA, http://www.ingenuity.com/). Among the 17,190 different TransMap RefSeq identifiers, 17,151 identifiers (99.7%) were represented in the Ingenuity Pathway Knowledge Base (IPKB) and corresponded to 7,177 different genes further considered as the reference set. Finally, 15,336 SNPs were located within a gene, leading to an average of 2.21 SNPs per annotated gene. Notice that 421 of these SNPs were representative of more than one gene.
Functional and Networks analyses were carried out using an approach similar to the one described in . Based on the information contained in the IPKB, IPA allows performing both functional and network analyses. The functional analysis aims at identifying the most significant biological functions represented in the candidate gene list from the reference set. Note that among the genes under selection, only those associated with at least one functional annotation in IPKB were eligible for functional analysis. The most significant functions are then obtained by comparing functions associated with eligible genes against functions associated with all the genes in the reference set using the right-tailed Fisher's exact test. The network analysis aims at searching for interactions (known from the literature) between candidate genes and all other IPKB molecules (genes, gene products or small molecules) and result in the definitions of networks which contains at most 35 molecules (including candidate genes). For each network, a score S is computed based on a right-tailed Fisher exact test for the overrepresentation of candidate genes (S = -log(p-value)). In our study, networks were considered as significant when S>3 and their associated top functions and diseases were further reported.
In order to identify regions under selection (with an unexpectedly high proportions of SNPs subjected to selection), we followed the locally adaptive procedure which allows to account for variations in distance between the different tested positions [81, 82]. Individual SNP BF values were first smoothed over each chromosome with a local variable bandwidth kernel estimator . We further performed 250,000 permutations to estimate local p-values which were corrected for multiple testing by computing q-values (see above).
Throughout the text, gene symbols are given according to HUGO standard nomenclature.
N'Dama (first population)
N'Dama (second population)
White Fulani Zebu
Zebu from Madagascar.
Allele Sharing Distance
Bos TAurus chromosome
Ingenuity Pathway Analysis
Ingenuity Pathway Knowledge database
False Discovery Rate
False Negative Rate
Minor Allele Frequency
Major Histocompatibility Complex
Principal Component Analysis
Single Nucleotide Polymorphism.
The authors wish to thank the farmers, the persons in charge of breeding programmes and the CIRDES (Burkina-Faso) team for their precious help in providing the samples from Burkina-Faso, Togo and Benin and Pr Lahoussine Ouragh from the Rabat Institute of Morocco for providing samples from Oulmès-Zaer breed. They also wish to thank anonymous referees for helpful corrections. The genotyping was supported by the French Ministry of Education and Research and the National Institute of Agronomic Research (INRA AIP "Bio-Ressources").
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.