The neighbor-joining tree based on Allele Sharing Distances (ASD) [12] allowed us to unambiguously separate individuals according to their population of origin (Figure S2 in additional file 1). West African individuals branched in their expected intermediary position relatively to the taurine (OUL and AUB) and zebu (ZMA) outgroups. Moreover, West African zebus (ZFU and ZCH) were closer to ZMA while West African taurines (LAG, ND1, ND2, BAO and SOM) were closer to OUL and AUB. KUR and BOR individuals branched between West African taurines and zebus. To go further in the characterization of the population relationships we carried out a principal component analysis (PCA) [13] based on all available SNP information [14]. As shown in Figure ​Figure2A2A and in agreement with previous published studies [1,4,15,16], the first component which accounted for 7.88% of variation separated West African populations according to a zebu/taurine gradient while the second (accounting for 4.58% of the total variance) could be associated to an Africa/Europe gradient, the North African OUL being closer to the European AUB. A similar PCA performed after removing the three outgroup populations (ZMA, OUL and AUB) allowed to distinguish West African populations according to the same zebu/taurine gradient on the first axis (which accounted for 8.25% of the total variance). The second and third components (which accounted respectively for 2.59% and 1.41% of the total variance) separated respectively LAG and ND2 from the other West African taurines (BAO, ND1 and SOM) (Figure ​(Figure2B).2B). The position of ND1 along the third axis suggested a certain level of admixture with shorthorn populations, most probably of BAO origin, according to their sampling area (Figure ​(Figure1),1), although pairwise F_ST divergence was smaller when compared to SOM (Table S1 in additional file 2). Finally, the fourth component (which accounted for 1.10% of the total variance) separated both West-African hybrids (BOR and KUR) and zebus (ZFU and ZCH) (Figure ​(Figure2C).2C). Overall, these results demonstrated a clear partition of the West African populations considered. F_ST between pairs of populations (Table S1 in additional file 2) were all found significantly non null (P << 0.00001) and ranged between 0.013 (for ZFU compared to ZCH) and 0.28 (for LAG compared to ZFU).

We further studied the differentiation among SNPs and populations based on a Bayesian hierarchical model derived from [7] 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 [17]. 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) [18]. Recently, efficient algorithms were proposed to account for ascertainment biases [19], 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 [18]. 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).

Whole genome maps of adaptively differentiated loci in West African populations are given in Figure ​Figure44 for four different thresholds on BF value. In total, 2,054 (5.7%), 1,119 (3.1%), 619 (1.7%) and 375 (1.0%) of the 36,320 SNPs displayed a BF above 15, 20, 25 and 30 respectively (Table S3 in additional file 2). As expected (see above) most of the SNPs displaying high evidence of selection were highly differentiated suggesting they were subjected to positive selection (Figures ​(Figures33 and ​and4).4). More precisely, among the SNPs with a BF value above 15, 537 were under balancing selection with a F_ST < 0.011 while 1517 were clearly under positive selection with a F_ST > 0.28. At a more stringent threshold on BF of 30, 3 SNPs remained under balancing selection (F_ST < 0.0013) and 372 under positive selection (F_ST > 0.33). Interestingly, the whole genome distribution of SNPs declared under selection appeared less uniform as the threshold on BF increased (Figure ​(Figure4).4). This might be explained by the effect of selection which tends to increase LD between markers [21]. Yet, we did not include in the model any spatial structure, which may take into account the LD among SNPs, as recently proposed ().().(.)[22]. Within the different West African populations and at the distance of 70 kb corresponding to our average marker spacing (Figure S6 in additional file 1 and Methods), r² was close to the asymptotic value (<0.1 in all populations except ND2) (Figure S7 in additional file 1) which is related to the current effective population size [23]. As a result, the level of association between most SNP pairs was expected to be of similar magnitude to the one between unlinked markers as illustrated by the correlations of the estimated F_ST and BF between pairs of SNPs which quickly dropped towards 0 as they were more distantly related (Table S4 in additional file 2).

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 [24]. 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 ​(Table3)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 [27] participates to both the host anti-viral response through the modulation of IRF3 functions [28] 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 [31]. 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 [32]. 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 ​Table3)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 [33]. These two molecules were shown to act as co-receptors for HIV entry in human cells [34]. A strong signature of balancing selection was reported in the 5' cis-regulatory region of CCR5 in human [35], this region displaying as well a selective sweep in chimpanzee [36]. We also recently observed a significant signal of selection on a microsatellite located less than 1 Mb from CXCR4 in West African taurine [10] and defining the boundary of a QTL underlying the trypanotolerance trait in cattle [37]. Similarly, RFX3 (43.235-43.563 Mb on BTA08) is located within another QTL identified in this latter study [37]. 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 [41]. 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 [42].

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 [71]. We subsequently performed a PCA based on all available SNP information using the SMARTPCA software package [14]. As suggested by Patterson et al. [14], 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 [72]. The global F-statistics F_IT, F_ST and F_IS were estimated respectively in the form of F, θ and f [73] using the program GENEPOP 4.0 [74]. 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 [14] and GENEPOP [74] 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. [8], of the initial Bayesian hierarchical model developed by Beaumont and Balding [7]. 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² + F_IS^jp_ij; 2(1-F_IS^j)p_ij(1-p_ij) and (1-F_IS^j)(1-p_ij)² + F_IS^j(1-p_ij) [18]. 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(λ_ijx_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_iand Var(p_ij) = x_i(1-x_i)(1+λ_ij)^-1. Actually (1+λ_ij)^-1 can be identified as a F_ST^ij coefficient [17]. 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. [78]). Following Beaumont and Balding [7], 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, σ_γ²~_iidN(α_i + β_j, σ_γ²) with β_j | μ_β, σ_β²~_iidN(μ_β, σ_β²) and α_i | σ_α²~_iidN(0, σ_α²). Introduction of additional levels by defining priors on the variance components σ_α², σ_β² and σ_γ² 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): σ_α² = 1, σ_β² = 3.24 and σ_γ² = 0.25. The prior mean μ for the β_j was also taken as -2.0. Recently, Riebler et al. [8] 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 * | σ_α*²~_iidN(0, σ_α*²). 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) [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 σ_α*² can be derived from σ_α² since by construction σ_α² = E(P) σ_α*² (hence σ_α*² = 10).

In order to identify outlier loci, we decided to base the decision rule on a Bayes Factor BF_i defined as the ratio of the posterior to the prior odds that locus i is selected: