In theory, neutrally evolving alleles form the majority of the genome variation within and between species, and events such as population contraction, expansion, migration, isolation and admixture are responsible for random drift of such alleles. While the demographic history of a population determines its neutral allele frequency spectrum, alleles that impact fitness, survival, reproduction, or traits of human interest are subjected to natural or artificial selection, such that their frequencies deviate from the distribution of neutral alleles. Therefore, mapping loci under selection implicates searching for outlier alleles that substantially differ from the genome background.

The idea of outlier loci is dependent upon knowledge about the frequency distribution of neutrally evolving alleles. As different populations have distinct demographic histories, their allele frequency spectrum also differ. Consequently, statistical methods designed to map loci under selection must be calibrated to either a demographic model for the population under study or to the empirical distribution of scores across the genome, assuming most of the genome is evolving under neutrality. As most of the events that shaped cattle diversity at the species and breed level are still obscure, methods that are robust across different demographic scenarios that take advantage of genome-wide scores to detect outlier loci are more appropriate to analyze cattle data. Nevertheless, one can use the demographic model for cattle differentiation, domestication and expansion provided in this section and combine with specific models for specialized cattle breeds to simulate the distribution under neutrality and compare with empirical data.

## **METHODS FOR DETECTING POSITIVE SELECTION IN THE CATTLE GENOME**

As proposed by Darwin and Wallace (1858), positive selection is the phenomenon where phenotypes that increase the likelihood of survival and reproduction (i.e., that increase an individual's fitness) become more prevalent in populations over time. In the context of the genome sequence, if a specific allele confers advantage, its carrier is more likely to thrive and leave more offspring than non-carriers, causing the haplotype containing that beneficial allele to spread quickly and increase in frequency in the population (Sabeti et al., 2002, 2006).

The majority of the genetic variation found within and between populations is deemed to have little or no effect on fitness, so that haplotype frequencies follow Hardy-Weinberg expectations (Kimura, 1968; Hellmann et al., 2003). Therefore, positive selection leaves distinctive tractable patterns of genetic variation that differ from the neutrally evolving background DNA sequence. These patterns are broadly referred as "signatures" or "footprints" of selection (Oleksyk et al., 2010).

In this section, we describe classes of signatures in terms of different signals. In each class, whenever convenient, we describe the basis of some popularly adopted methods to highlight their strengths and weaknesses for cattle studies. For a broader overview, we suggest the reader to consult other previously published reviews on the topic (Sabeti et al., 2006; Oleksyk et al., 2010; Vitti et al., 2013; Qanbari and Simianer, 2014). Additionally, we provide examples of studies reporting putative loci under selection in different cattle breeds.

#### **LOCAL GENETIC DIVERSITY DEPRESSION**

At each generation, recombination shuffles and breaks haplotypes down, producing linkage equilibrium. When selection increases the frequency of an advantageous allele, neighbor variants "hitchhike" and rise in frequency together so quickly that recombination does not prevent linkage disequilibrium, causing an entire chromosome segment to lose diversity. Therefore, positive selection can be probed by searching for chromosome regions where heterozygosity is much lower than expected under neutrality (Oleksyk et al., 2008).

Ramey et al. (2013) scanned Illumina® BovineSNP50 (50 k) genotypes of over 6000 animals from 13 taurine and one zebu breed using an *ad-hoc* method to identify selective sweeps reaching fixation. Briefly, they considered a candidate locus under selection if at least five contiguous SNPs presenting a minor allele frequency (MAF) > 0.01 spanned a chromosome segment of 200 kb or more. Although the strategy was successful to find validated selective sweeps, such as the *POLL* locus that controls horn development (Georges et al., 1993), they showed that the SNP ascertainment bias of the 50 k assay incurred in false positives in breeds that are genetically distant from the SNP discovery breeds. This is not unexpected, as the influence of the 50 k bias on population genetics parameters, such as heterozygosity, genetic structure and differentiation, is well documented (Decker et al., 2009). Therefore, methods relying on heterozygosity are largely prone to type I errors generated from the SNP discovery process. Ideally, genotype data should be generated from less biased commercial arrays, such as the Illumina® BovineHD (HD), or from re-sequencing efforts.

Another approach to search for deficit of heterozygosity is the identification of islands of runs of homozygosity (ROH). As long stretches of consecutive homozygous genotypes indicate identical-by-descent haplotypes (Gibson et al., 2006; Lencz et al., 2007), ROH have been recently used to characterize genomewide inbreeding in cattle (Purfield et al., 2012; Ferencakovi ˇ c et al., ´ 2013; Kim et al., 2013). However, as recombination is random, the distribution of ROH across samples is expected to be highly heterogeneous under neutrality, so genomic hotspots of ROH can be a signal of selective sweep (Curik et al., 2014). Interestingly, the length of the run is negatively correlated with the number of generations back to the selective pressure or inbreeding event (Howrigan et al., 2011), so ROH size can be used to date the age of the selective sweep.

Focusing on the identification of ROH islands produced by recent artificial selection in U.S. Holstein cattle, Kim et al. (2013) proposed a new statistic, namely *FRL*, that compares local homozygosity between a population under artificial selection and a control population. Briefly, for each SNP, animals are scored as 1 if the locus is encompassed by a ROH or 0 otherwise. Then, *FL* is computed as the proportion of animals with scores equal to 1. Finally, *FRL* is obtained as:

$$F\_{RL} = \ln\left(\frac{F\_{L(selected)}}{F\_{L(control)}}\right).$$

All scores are standardized (i.e., subtracting each value by the average score and dividing by the standard deviation) to produce a distribution with mean zero and variance one. Extreme positive values of *FRL* represent changes in allele frequency in the artificially selected population in comparison to the control group, and therefore reflect homozygosity attributable to recent inbreeding, drift or selection at the analyzed locus. As for the approach adopted by Ramey et al. (2013), SNP density and ascertainment bias are important confounders that can produce false positive ROH (Ferencakovi ˇ c et al., 2013 ´ ), so commercial SNP assays must be used with caution here.

Methods that search for local absence of heterozygosity are most powerful to detect haplotypes that hitchhiked to fixation. However, allowing for some heterozygosity can be of help to detect partial sweeps, i.e., haplotypes under ongoing selection that have not reached fixation yet. For a given bi-allelic site, let *nREF* be the number of observed reference alleles, *nALT* the number of observed alternative alleles, and p and q their respective frequencies. The expected frequency of heterozygote genotypes under Hardy-Weinberg Equilibrium is given by:

$$H = 2pq = \frac{2n\_{REF}n\_{ALT}}{\left(n\_{REF} + n\_{ALT}\right)^2}$$

This expectation can be easily extended for a chromosomal segment containing multiple bi-allelic sites as:

$$H = \frac{2\Sigma n\_{\rm REF} \Sigma n\_{\rm ALT}}{\left(\Sigma n\_{\rm REF} + \Sigma n\_{\rm ALT}\right)^2}$$

Rubin et al. (2010) proposed running sliding windows across the genome, calculating *H* for these windows and then standardizing the obtained values. The method named as *ZHp* produces standard deviation scores, and extremely negative values represent chromosome windows where the regional diversity is substantially lower than the average genome diversity. Although the method has been successfully used to map candidate variants under selection in chickens (Rubin et al., 2010), dogs (Axelsson et al., 2013) and pigs (Rubin et al., 2012), no studies applying *ZHp* to cattle data have been published to date.

#### **CHANGE IN THE ALLELE FREQUENCY SPECTRUM**

After a complete selective sweep (i.e., the haplotype containing the selected variant reaches fixation), new mutations slowly restore local diversity in the course of many generations. As mutations are generally rare and may take a large number of generations to drift to high frequency under neutral evolution, the local heterozygosity depression signal is deemed to remain robust for several generations after the selective pressure has occurred.

Newly acquired mutations or derived alleles (i.e., variants that differ from the original or ancestral allele) occur in lower frequency in comparison to ancestral alleles under neutrality, but when they arise within a selective sweep they will hitchhike to high frequency quickly in the selected population. Therefore, another class of signals that emerge after the depression in local diversity is the enrichment for moderate to high frequency derived alleles. Several methods have been proposed in this category, including Tajima's (1983), Fay and Wu's (2000), and -*DAF* (Grossman et al., 2010).

The limitation of the use of these methods in cattle data is that inference of ancestral alleles should be preferably performed by comparison of domesticated cattle genomes with wild type genomes. As these are not available, ancestral allele information for cattle SNP assays were derived from genotypes of outgroup species, such as Gaur (*Bos gaurus*), Buffalo (*Bubalus bubalis*) and Yak (*Bos grunniens*), which are assumed to descend from a common founder *Bovinae* species (Matukumalli et al., 2009; Utsunomiya et al., 2013). However, as not all SNP probes hybridize against the genomes of these outgroup species, ancestral allele information for cattle is incomplete. Future availability of genome assemblies for these outgroup species may help to better infer ancestral status for common SNPs.

One way to bypass the limitation of ancestral allele information is by using samples from several genetically distinct populations to estimate average allele frequencies that could represent the spectrum in the common ancestral population. In this approach, instead of searching for enrichment of high frequency derived or rare alleles, one may look for a shift in the allele frequency spectrum in one population in comparison to the average across populations. Stella et al. (2010) successfully applied the negative composite log-likelihood (*CLL*) approach, an extension of the composite likelihood ratio (*CLR*) test (Kim and Stephan, 2002; Nielsen et al., 2005) to 13 taurine, 4 zebu and 2 synthetic breeds, and reported several candidate loci under selection. These included *KIT* (mast/stem cell growth factor receptor gene), responsible for the piebald and color sidedness phenotype (Durkin et al., 2012), and *MC1R* (melanocortin 1 receptor gene), incriminated in the black/red coat color in Holstein and Angus (Klungland et al., 1995). These loci were further confirmed using whole genome sequence data of German Fleckvieh (Qanbari et al., 2014).

Briefly, *CLL* is computed for each SNP window as:

$$CLL = -\sum\_{i=1}^{k} \ln\left[1 - Pr\left(d < |p\_i - \overline{p}\_i| \mid \mu\_i\right)\right],$$

Where, relative to SNP *i*, *d* is any random value from the theoretical distribution of allele frequencies with mean μ*i* = *pi* , *pi* is the reference allele frequency averaged across populations, and *pi* is the allele frequency in the population being investigated. The theoretical distribution of allele frequencies can be modeled as a binomial or a normal approximation to the binomial distribution. As calculations are very straightforward and only require a dataset with multiple breeds instead of ancestral allele information, it is general enough for cattle SNP data.

Another extension of the *CLR* statistic was proposed by Chen et al. (2010), namely *XPCLR* (cross-population composite likelihood ratio test), which attempts to detect a shift in allele frequency in a target population in respect to a reference population. Lee et al. (2014) applied *XPCLR* to Holstein and Hanwoo whole genome sequence data and reported that the chromosome segment encompassing the kappa-casein gene (*CSN3*) exhibited high *XPCLR* scores in Holstein cattle.

#### **LONG-RANGE HAPLOTYPES**

The concept of Extended Haplotype Homozygosity (*EHH*), first introduced by Sabeti et al. (2002), attempts to identify haplotypes that increased so rapidly in frequency that recombination could not substantially break them down, so linkage disequilibrium presents a long-range persistency. Briefly, consider *N* chromosomes in a sample, and *G* unique haplotypes extending from a core SNP site to an upstream or downstream position *x*, with each group *g* having *ng* observations. The *EHH* score for the entire sample is calculated as:

$$EHH = \frac{\sum\_{\mathfrak{g}=1}^G \left(\frac{n\_{\mathfrak{g}}}{2}\right)}{\left(\frac{N}{2}\right)}$$

This score serves as a proxy for the probability of identity-bydescent of haplotypes within the chromosome segment being investigated. Generally, *EHH* is calculated at varying distances from the core SNP position, so that the decay of *EHH* as a function of physical distance can be assessed to determine the extension of the haplotype homozygosity. From this seminal concept, a family of statistical methods was developed in order to scan entire genomes in the search for evidence of selection.

Voight et al. (2006) proposed to measure how rapidly *EHH* decays from a core SNP site by calculating the area under the *EHH* curve,

$$iHH = \int\_{a}^{b} EHH(\mathbf{x})d\mathbf{x}$$

where *iHH* represents the definite integral of *EHH* evaluated over the domain of the chromosome segment delimited by upstream position *a* and downstream position *b* where *EHH* decays to some arbitrary small value (originally 0.05). As the area under the curve is not tractable analytically, a trapezoid quadrature with nonuniform grid can be adopted as a deterministic approximation:

$$iHHI \sim \sum\_{k=1}^{K} \left( \chi\_{k+1} - \chi\_k \right) \frac{(EHH\_{k+1} + EHH\_k)}{2}.$$

A within population score, namely Integrated Haplotype Score (*iHS*), for a given site *i*, was introduced by Voight et al. (2006) as the log-ratio between the integrated *EHH* for the haplotypes containing the ancestral allele (*iHHA*) and the derived allele (*iHHD*):

$$iHS\_i = \ln\left(\frac{iHH\_{A,i}}{iHH\_{D,i}}\right)$$

These scores are then standardized to have mean zero and variance one.

Extremely negative standardized *iHS* values have been of particular interest in human genetics, as they represent a recently acquired mutation that increased very rapidly in frequency (i.e., there is a partial sweep due to ongoing selection) or a haplotype that hitchhiked to fixation and then became enriched for derived alleles (Voight et al., 2006). However, a sweep can also produce large positive *iHS* values at nearby SNPs if ancestral alleles hitchhike with the selected site, so the chromosome region surrounding the selected variant typically exhibits a cluster of extreme positive and negative *iHS* values. Furthermore, in the context of cattle data, artificial selection and domestication probably favored "beneficial" alleles in the sense of human interest, regardless if it is ancestral or derived. Therefore, both positive and negative values should be investigated in cattle data. This implicates that the absolute value of standardized *iHS* scores should be preferred over the signed values, or, equivalently, that a two tailed hypothesis test should be assumed. As only partial ancestral allele information is available for cattle SNP assays, and the search for footprints of selection by *iHS* in cattle should disregard the direction of the sweep, a more appropriate generalized version of *iHS* can be postulated as the log-ratio between the integrated *EHH* for an arbitrary reference allele (*iHHREF*) and for the alternative allele (*iHHALT*).

One of the limitations of this method is that if a given marker presents a nearly or completely fixed allele in the population being analyzed, this allele will have no integral to be calculated or an integral close to zero, so the log-ratio will result in a positive or negative infinite value. In this scenario, the calculation of *iHS* must be conditioned by *iHHREF* > 0 and *iHHALT* > 0, which indirectly leads to a minor allele frequency (MAF) constraint. This limitation renders *iHS* underpowered to detect very recent nearly fixed selective sweeps, which are of primary interest in the cattle community. However, as discussed earlier, a crucial point to be considered is that contiguous chromosomal segments containing SNPs with MAF = 0 can also result from SNP chip ascertainment bias, which may produce false positive signals.

Tang et al. (2007) and Sabeti et al. (2007) have independently developed equivalent methods, *Rsb* and *XPEHH*, respectively, which attempt to compare long-range haplotypes between populations in order to increase the power of selective sweep detection. The most crucial improvement is that, for each population being analyzed *iHH* is calculated for the entire sample, instead of being partitioned between derived and ancestral alleles. This eliminates the MAF constraint and recovers the power to detect sweeps reaching fixation. The comparison with a population where the selective sweep may not have occurred adds extra power to the method. Calculations are performed as follows:

$$XPEHH\_i = Rsb\_i = \ln\left(\frac{iHH\_{pop1,i}}{iHH\_{pop2,i}}\right)$$

Where, relative to SNP *i*, *iHHpop*1,*i* is the integrated *EHH* in the first population and *iHHpop*2,*i* is the integrated *EHH* in the second population. Scores are also standardized to produce a distribution of standard deviates. Positive values indicate selective sweeps in the population used in the numerator, while negative values indicate selection in the population used in the denominator. Here, it is easy to keep track of the signals by using one-tailed hypothesis tests.

Studies applying *EHH*-based methods to cattle data are numerous (for instance, Hayes et al., 2008; Gautier and Naves, 2011; Qanbari et al., 2011, 2014; Flori et al., 2012; Utsunomiya et al., 2013; Huson et al., 2014). The reported loci are deemed to be genome responses to a variety of different selective pressures, such as milk and meat production, coat color, heat stress, and reproductive performance. Among these, one particularly interesting selective sweep, most likely related to adaptation to heat stress, has been reported in Creole cattle, including Senepol, Carora, Romosinuano, and cross-bred lineages (Flori et al., 2012; Huson et al., 2014). These cattle breeds present the slick hair coat phenotype, a dominant trait associated to heat tolerance in tropically adapted cattle that descend from Spanish cattle introduced to the New World. The chromosome segment containing the selective sweep ranges from 37.5 to 39.6 Mb on chromosome 20, with a variable peaking position (39.5 or 37.7 Mb) depending on the SNP panel (BovineSNP50 or BovineHD) and dataset analyzed (Flori et al., 2012; Huson et al., 2014). The disputed positional candidate genes are the retinoic acid induced 14 (*RAI14* or *NORPEG*), prolactin receptor (*PRLR*), and S-phase kinaseassociated protein 2 (*SKP2*). A strong candidate mutation has been recently proposed for *PRLR*, a single base deletion in exon 10 (ss1067289408) predicted to cause a frameshift that introduces a premature stop codon (p.Leu462∗) and consequent loss of 120 C-terminal amino acids from the long isoform of the prolactin receptor (Littlejohn et al., 2014).

#### **POPULATION DIFFERENTIATION**

Following the same principle as in *Rsb*/*XPEHH*, *XPCLR* and *FRL*, although positive selection may act across populations sharing geographical proximity, environmental factors or a common phenotype, outgroup populations may not share the same selective pressure. Therefore, changes in allele frequency promoted by selection in one group will not be detectable in the other, and large differences in allele frequencies between populations will be observed.

The fixation index *FST* (Wright, 1950; Weir and Cockerham, 1984) and its abundant estimators is the gold standard method for detecting highly differentiated loci between populations. In essence, given the average allele frequency *p* across subpopulations, *FST* is simply the ratio between the variance in the allele frequency in different subpopulations σ2 *S* = *k j* = 1 *pj* − *p* 2 and the variance in the total population σ2 *T* = *p*(1 − *p*). In pair-wise comparisons, calculations simplify to:

$$F\_{ST} = \frac{(p\_1 - p\_2)^2}{\left(p\_1 + p\_2\right)\left(2 - p\_1 - p\_2\right)}$$

Scores can be averaged across SNP windows or smoothed against genomic positions using a local variable bandwidth kernel estimator (Porto-Neto et al., 2013).

Flori et al. (2009) applied *FST* to three French dairy cattle breeds (Holstein, Normande, and Montbéliarde), and found that some of the putative loci under selection in that study overlapped genes that strongly affect milk production traits, such as the growth hormone receptor (*GHR*), and coat color, for instance *MC1R*. Porto-Neto et al. (2013) generated a comprehensive map of divergent loci between taurine and zebu cattle using over 777,000 SNPs and 13 cattle breeds, and reported that the highest scoring locus in the *FST* analysis maps to chromosome 7:47.2– 53.7 Mb, which shelters a cluster of immune-related genes and *SPOCK1*, a gene previously implicated in puberty (Fortes et al., 2010).

Extensions of *FST* led to the development of the *FLK* test (Lewontin and Krakauer, 1973; Bonhomme et al., 2010; Fariello et al., 2013), which not only account for effective population size and hierarchical structure among populations, but also have known distributional properties under neutrality. Briefly, let *qi* be the vector of reference allele frequencies at marker *i* for the populations being compared, and *q*0 the ancestral allele frequency for the same allele. The *FLK* method relies on the linear model:

$$q\_i = 1\_n q\_0 + e$$

where the residual term *e* is assumed *N* (0, *Vi*), and *Vi* is the expected variance-covariance matrix for vector *qi*. This matrix is modeled as:

$$V\_i = \mathcal{F}q\_0(1 - q\_0)$$

where *F* is a kinship matrix. The diagonal elements of *F* represent the expected inbreeding coefficients in each respective population, and off-diagonal elements represent the amount of drift accumulated on the different branches of the population tree. A weighted least squares estimate of *q*0 is then obtained as:

$$\hat{q}\_0 = \left(\mathbf{1}\_n^\prime \mathcal{F}^{-1} \mathbf{1}\_n\right)^{-1} \mathbf{1}\_n^\prime \mathcal{F}^{-1} q\_i$$

The *FLK* score is a measure of goodness-of-fit of this model, and is calculated as the deviance (i.e., residual sum of squares):

$$FLK\_i = \left(q\_i - \hat{q}\_0 \mathbf{1}\_n\right)' V\_i^{-1} (q\_i - \hat{q}\_0 \mathbf{1}\_n)'$$

Under the assumption of a star-like tree-pure drift model (i.e., the populations being compared evolved in parallel from a single ancestral population, with no mutations or admixture), matrix *F* can be simplified to *F* = *InFST*, where *In* is an identity matrix and *FST* is the average *FST* over all SNP loci. In this case, the average allele frequency across populations is an unbiased estimator of *q*0, and the deviance is simplified to *FST*(*n* − 1)/*FST*, which gives a test statistic that is linearly correlated with *FST*. As discussed later, while *FST* has no known theoretical distribution under the neutral model, *FLK* can be modeled as a chi-squared distribution. Moreover, *FLK* outperformed *FST* in simulations, especially when scores were computed based on haplotypes (*hapFLK*) instead of single markers (Fariello et al., 2013). This method is a suitable alternative to *FST* for cattle data.

#### **DISTRIBUTIONS UNDER THE NULL HYPOSTHESIS AND** *p***-VALUES**

Although different in formulation and assumptions, all methods presented so far attempt to address the same null and alternative hypothesis: *H*0 = locus is neutral; *H*1 = locus is not neutral. Under the hypothesis of neutrality, *EHH*-based methods are standard normal deviates (Voight et al., 2006; Sabeti et al., 2007; Tang et al., 2007). Hence, the probability that SNP *i* with *Rsb* or *XPEHH* score *xi* is neutral in the population used as numerator can be approximated by an upper tail *p*-value derived from the normal cumulative density function (CDF) :

$$\Pr\left(Rsh, XPEHH > \chi\_i \mid neutral\right) = 1 - \Phi\left(\chi\_i\right)$$

Likewise, the probability that SNP *i* with *iHS* score *xi* is neutral, given that both reference and alternative alleles are of interest, can be approximated by a two-tailed *p*-value:

$$\begin{aligned} \Pr\left( |iHS| > |\mathbf{x}\_i| \; | \; neutral \right) &= 1 - |\Phi\left(\mathbf{x}\_i\right) - \Phi\left(-\mathbf{x}\_i\right)| \\ &= 1 - 2 \left| \Phi\left(\mathbf{x}\_i\right) - \mathbf{0}.5 \right| \end{aligned}$$

Recently, concerns have emerged on the interpretation of *p*-values in signatures of selection analyses (Simianer et al., 2014). It is argued that, at least for the cases of *EHH*-based methods, scores standardized using genome-wide data are not test statistics in the classical sense but only deviates from an average (Voight et al., 2006), so *p*-values would only represent quantiles from an empirical distribution, rather than formal significance values. This implies that the probability of obtaining an arbitrary score given the locus is neutral cannot be exactly computed once the underlying true probability distribution may vary according to different demographic scenarios. The major caveat here is, for any given test discussed, the proposed asymptotic distribution under the hypothesis of neutrality is still largely based on coalescent simulations with demographic models calibrated for human data. However, if the majority of the genome is under neutrality and the null distribution is robust to a wide range of demographic models, one may expect that genome-wide distribution of scores in cattle data should mimic the simulations performed for the human model, serving as a control (Gianola et al., 2010). This has been observed in SS studies in cattle (Gautier and Naves, 2011; Flori et al., 2012; Utsunomiya et al., 2013), and therefore these theoretical distributions under neutrality are suitable for cattle data.

In the case of *FST*, although approximate distributions are available (e.g., exponential or beta), such approximations are sub-optimal. One advantage of *FLK* over *FST* is that scores are asymptotically distributed as χ2 ν under the null hypothesis, where degrees of freedom ν can be equivalently calculated as *number of subpopulations* − 1 or as the average *FLK* across all loci. Upper tail *p*-values can be computed as:

$$\Pr\left(FLK > \chi\_i \mid \mathit{neutral}\right) = 1 - F\_{\mathbb{V}}(\chi\_i)$$

where *F*ν is the χ2 CDF with ν degrees of freedom.

Although *FST* and *CLL* scores have unknown distributions under neutrality, both rely on population comparisons. The use of datasets consisting of multiple populations allows for permutation tests via random sampling of individuals or random sorting of population labels in order to compute an empirical null distribution. Nevertheless, permutation tests are computationally intensive and may be impracticable in re-sequencing data.

The *ZHp* method suffers from the same problem, with the additional limitation of not allowing data permutations when a single population is surveyed. Although a truncated normal distribution could be suggested to approximate its null distribution, this has not been properly assessed in practice. Another challenge is deriving significance values when scores are averaged across SNP windows. Following the assumption that the majority of the loci (single markers or SNP windows for that matter) are under neutrality, an empirical CDF derived from genome-wide scores should converge to the underlying true CDF, so probability values could be empirically obtained from a step function:

$$\begin{aligned} \Pr\left(\text{Test} > \text{x}\_{i} \mid neutral\right) &= 1 - ECDF\left(\text{x}\_{i}\right) \\ &= 1 - \frac{\left(number\ of\ scores\leq\text{x}\_{i}\right)}{number\ of\ observed\ scores} \end{aligned}$$

The probability that SNP *i* with score *xi* in a given test is selected is not as trivial to approximate. While the distributions under the hypothesis of neutrality seems to be robust across a wide range of demographic models, the distribution under the hypothesis of selection may vary widely depending upon the demographic history that shaped the allele frequency spectrum. Therefore, there is no unique theoretical distribution to represent selected variants, and coalescent simulations using well calibrated demographic models are required in order to generate empirical distributions. This has been successfully done for human data (Grossman et al., 2010), but is highly challenging in cattle as most of the demographic history of the bovine species is yet to be uncovered. Nevertheless, promising methods to infer demographic scenarios from the data, including estimates of population sizes and population separations over time, are now emerging (Schiffels and Durbin, 2014), which could be useful to elucidate cattle history and customize coalescent simulations based on empirical data.

#### **AVAILABLE SOFTWARE AND ANALYSIS BEST PRACTICES**

**Table 1** summarizes all the software mentioned below. The first step to be considered before starting a SS analysis is filtering out poor quality data and markers and samples that are not informative or that may confound the analysis. Discussing the particularities of measures of quality control is beyond the scope of the present article, but some of the best practices in SS studies can emerge from common sense. First, metrics such as Hardy-Weinberg Equilibrium (HWE) and MAF, highly used in GWA studies, should be applied with caution. Elimination of markers with extreme deviations from HWE expectations may counteract the whole SS enterprise, as in this type of study outlier loci are being sought. Likewise, MAF controls may cause signals that are reaching fixation to be completely lost. A situation where HWE and MAF thresholds can be benign is when only markers with extreme excess of heterozygotes are eliminated, and MAF constraints are applied to the pooled allele frequencies across all populations in the study. We suggest *PLINK* v1.07 (Purcell et al., 2007) or *PLINK* v1.90 for this first data screening. Other usual filters such as GenCall and GenTrain scores and call rates should follow the same guidelines as in GWA studies.

Another important issue is cryptic relatedness. Eliminating samples of closely related animals is of paramount importance to reduce false positive signals. We have previously reported an algorithm to find the maximum independent set based on identity-by-descent, i.e., maximize the number of samples while eliminating first and second degree relationships using SNP data (Utsunomiya et al., 2013). Other software such as *PLINK* v1.90 and *GCTA* v1.24.2 (Yang et al., 2011) also provide means to find the optimal set of independent individuals.

The next key step is producing high quality phased data. There are several methods and implementations for phasing, such as *fastPHASE* v1.2 (Scheet and Stephens, 2006), *Beagle* v3.3 or later (Browning and Browning, 2008), and *SHAPEIT2* (O'Connell et al., 2014). Although *fastPHASE* exhibits high phasing accuracy, it is orders of magnitude more computationally intensive than *Beagle* or *SHAPEIT2*. It is important to notice that phasing is not a straightforward procedure, and is highly prone to errors. Consequently, results from haplotype-based methods should be assessed with caution. The effect of haplotype errors on SS results remains underexplored and at some extent neglected.

Regarding the SS analysis *per se*, *EHH*-based methods can be computed using *Sweep*, *selscan* (Szpiech and Hernandez, 2014) or the *R* package *rehh* (Gautier and Vitalis, 2012). *FRL* has a dedicated software, as well as *CLR* (available at: and *XPCLR*. For ROH-based tests, runs can be computed using either *PLINK* or *SNP & Variation Suite v7.6.8* or later, and *FRL*can be easily computed with home grown scripts. In the cases of *FST*, *CLL*, *ZHp* and

#### **Table 1 | Available software for signatures of selection analysis.**

the method reported by Ramey et al. (2013), allele frequencies and genotype counts can be obtained from either *PLINK* or *SNP & Variation Suite*, and calculation of scores can be easily implemented in *R* (available at: http://www.r-project.org/) or other languages. Simulations of population genetics datasets under specific demographic models, including neutral and selected loci, can be performed using coalescent simulators such *cosi* (Schaffner et al., 2005), *cosi2* (Shlyakhter et al., 2014), or *MSMS* (Ewing and Hermisson, 2010).

## **COMBINING SELECTION SIGNALS**

The available methodologies to detect positive selection differ substantially from each other in terms of the pattern of genetic variation encrypting a "signal." However, all of them have a shared objective: to identify loci that have undergone positive selection. Indeed, at least for recent selective pressures (up to a few thousand generations back), a selected variant is expected to be in a chromosome segment where there has been loss of diversity, enrichment for derived or rare alleles, population differentiation, and highly frequent long-range haplotypes. Therefore, collecting evidence across different methodologies targeting distinct classes of signals may help in identifying loci under positive selection. This section explores the statistical properties and limitations of some of the available methods designed to combine different methods for signatures of selection.

#### **COMPOSITE OF MULTIPLE SIGNALS (CMS)**

In the original implementation of the method (Grossman et al., 2010), CMS is a local test designed to narrow down signals detected from *s* distinct genome-wide scans, and is defined as the approximate joint posterior probability that a given variant is selected:

$$\begin{split} \text{CMS}\_{\text{local}} &= \Pr\left(\text{selected} \mid \mathbf{x}\_{\vec{\text{ij}}}\right) \\ &= \prod\_{\vec{\text{ij}}=1}^{s} \frac{\Pr\left(\boldsymbol{\omega}\_{\vec{\text{ij}}} | \text{selected}\right) \Pr\left(\text{selected}\right)}{\Pr\left(\boldsymbol{\omega}\_{\vec{\text{ij}}} \mid \text{selected}\right) \Pr\left(\text{selected}\right)} \\ &\quad + \Pr(\boldsymbol{\omega}\_{\vec{\text{ij}}} | \text{neutral}) \Pr(\text{neutral}) \end{split}$$

The genome-wide extension (Grossman et al., 2013) focuses on the product of the Bayes Factors for each one of the tests to be combined. For each test, BF is computed as the ratio between the posterior and prior odds:

$$\text{CMS}\_{\text{GW}} = \prod\_{j=1}^{s} \text{BF}\_{\text{ij}} = \prod\_{j=1}^{s} \frac{\Pr(\text{x}\_{\vec{\text{ij}}} | \text{selected}) \Pr(\text{selected})}{\Pr(\text{x}\_{\vec{\text{ij}}} | \text{neutral}) \Pr(\text{neutral})}$$

In the absence of prior information on the number of loci under selection across the genome, CMS scores simplify to:

$$\text{CMS}\_{\text{GW}} = \prod\_{j=1}^{s} \frac{\Pr(\boldsymbol{\kappa}\_{ij}|s \text{selected})}{\Pr(\boldsymbol{\kappa}\_{ij}|neutral)}$$

A challenging aspect of the implementation of these methods is computing Pr(*xij*|*selected*), which requires simulations under clear demographic assumptions. In contrast, as discussed earlier, many methods designed to detect markers under positive selection allow for approximating Pr(*xij*|*neutral*) from asymptotic theoretical or empirical distributions. Therefore, composite tests considering only the distributions under neutrality are more appropriate for cattle data.

The original CMS score can be modified to take advantage of the assumed distributions under neutrality in order to relax demographic assumptions and avoid expensive simulations. The most essential modification involves reformulating the problem of detecting markers departing from neutrality. Instead of considering the assessment of whether a marker has been selected or not, one can look for support from the data against the null model, i.e., that the marker does not fit well to the neutral model. First, let the new CMS score be the approximate joint posterior probability of a given variant not being neutral:

$$\begin{aligned} \text{CMS}\_{\text{new}} &= \Pr\left(\text{not } neutral[\mathbf{x}\_{ij}] \\ &= \prod\_{j=1}^{s} \frac{\Pr(\mathbf{x}\_{ij}|not \text{ } neutral)\Pr(\text{not } neutral)}{\Pr\left(\boldsymbol{\chi}\_{ij}|neutral\right)\Pr(\text{neutral})} \\ &\quad + \Pr\left(\boldsymbol{\chi}\_{ij}|not \text{ } neutral\right)\Pr(\text{not } neutral) \end{aligned}$$

Here, Pr(*xijneutral* | ) is computed directly from its theoretical distribution, and Pr *xij not neutral* is computed as 1 − Pr(*xij*|*neutral*). Also, it can be assumed that the prior Pr(*not neutral*) = 1 − Pr(*neutral*). Therefore, the new CMS score can be re-written as:

$$\text{CMS}\_{\text{new}} = \prod\_{j=1}^{s} \frac{\left[1 - \Pr\left(\boldsymbol{\chi\_{ij}}[neutral]\right)\right] \left[1 - \Pr(neutral)\right]}{\Pr\left(\boldsymbol{\chi\_{ij}}[neutral]\right) \Pr(neutral)}$$

Likewise, the genome-wide modified CMS scores can be reformulated as:

$$\text{CMSG}\_{\text{GW}-\text{new}} = \prod\_{j=1}^{s} \frac{\Pr(\text{x}\_{\vec{\text{ij}}} | \text{not } neutral)}{\Pr(\text{x}\_{\vec{\text{ij}}} | neutral)} = \prod\_{j=1}^{s} \frac{1 - \Pr(\text{x}\_{\vec{\text{ij}}} | neutral)}{\Pr(\text{x}\_{\vec{\text{ij}}} | neutral)}$$

It is important to note that this modification does not allow for the same interpretation as the original CMS method: the composite likelihood does not indicate selection, but rather, that a marker does not fit well the neutral model.

#### **META ANALYSIS OF SELECTION SIGNALS (META-SS)**

Following the ideas expanded from the landmark publication of Grossman et al. (2010), given the probabilities under the null hypothesis for each test, our interest is to identify loci presenting consistent rejection of the neutral model across the different tests. For any given statistic, *p*-values are uniformly distributed in the interval between 0 and 1 under the null hypothesis. This property makes possible to use an inverse CDF, such as the Gaussian density, to produce scores for each test derived from a single theoretical distribution. Therefore, for each SNP *i* and test *j*, a new score can be computed as *Zij* = −1(1 − *Pij*), where *Pij* is the *p*-value. These Z-transformed *p*-values can be then averaged and standardized to produce a composite score.

We have previously described *meta* − *SS* (Utsunomiya et al., 2013), an abstraction of the Stouffer Z-transformation for combining different selection signals using the aforementioned framework. As the Stouffer method assumes the tests are uncorrelated under the shared null hypothesis and the use of pair-wise comparisons produce correlated scores, a weighted average was originally proposed to penalize dependent tests:

$$meta - SS = \frac{\sum\_{j=1}^{s} \alpha\_{j} Z\_{j}}{\sqrt{\sum\_{j=1}^{k} \alpha\_{j}^{2}}}$$

where ω*j* is the weight for test *j*. In this setting, a uniform penalization can be applied to control for the inflation of correlated tests. As this penalization does not incorporate the strength of correlations among tests, the *meta* − *SS* test can be modified to explicitly account for the magnitude of correlations between scores. Considering all scores are equally weighted, the corrected composite score can be computed as:

$$meta - \text{SS} = \frac{\sum\_{j=1}^{s} Z\_j}{\sqrt{k + 2R}}$$

where *R* is the sum of all pairwise Pearson's product-moment correlations. Under the hypothesis of neutrality, these composite scores are distributed as *N*(0, 1), so the higher is the Z-transformed value, the worse the marker fits the neutral model. Upper tail *p*-values can then be obtained from the standard normal CDF.

Obvious limitations from *meta* − *SS* and CMS is the inability to incorporate statistics for which *p*-values cannot be derived. Randhawa et al. (2014) proposed Composite Selection Signals (CSS), a nonparametric interpretation of *meta* − *SS*, where fractional ranking is used instead of *p*-values to combine different tests. Briefly, the vector of test statistics for method *j* is first sorted and then ranked, taking values 1, ..., *k*. Next, the vector of ranks is re-scaled by dividing all elements by *k* + 1, thus producing a variable ranging from 0 to 1. These re-scaled ranks are treated as they were *p*-values for the test statistics, and then combined as in the *meta* − *SS* approach. This strategy is equivalent to computing probabilities from an empirical CDF using a step function, as discussed earlier, which has an appealing feature: as fractional ranking can be generated for any particular test, signature blending is made feasible even if the theoretical distributions are unknown or if scores have been averaged in chromosome windows. However, a caveat is that the magnitude of the actual test statistics may be lost, so one may expect loss of power compared to the use of theoretical or simulated distributions.

#### **PRINCIPAL COMPONENTS ANALYSIS**

Simianer et al. (2014) proposed combining different tests by applying an eigendecomposition of the correlation matrix of the scores. The attractive feature of this method is the possibility of using standardized scores instead of approximate probabilities. However, as each principal component has heterogeneous loadings from each test, deriving a single synthetic score that summarizes all different tests remains a challenge in this framework.

## **CHALLENGES AND FUTURE DIRECTIONS**

In theory, genome-wide genotypes are a vast source of information that can be explored in the search for large effect mutations that underwent selection. However, the existing data and methods still suffer from power issues and confounding effects that can give rise to false positive and false negative signals.

Although simulations suggest that only marginal gains in power are obtained when the sample size is increased from tens to hundreds of unrelated samples, marker density and allele frequency spectrum seems to impact power dramatically (Lappalainen et al., 2010; Simianer et al., 2014). Genotypes derived from commercial SNP arrays have two important limitations in this context: (1) incomplete genome coverage by markers; and (2) ascertainment bias. The search for SS must be preferentially performed using high density SNP panels, although optimal average intermarker distances to detect a sweep may vary depending on the effective population size, extent of linkage disequilibrium and the nature of the signal. Regarding ascertainment bias, commercial SNP arrays are suitable for cattle populations that are closely related to the breeds used in the SNP discovery process, but there is no guarantee they will be informative in genetically distant populations. Indeed, with a few exceptions, little congruence has been reported between candidate selected loci identified using whole genome sequence and different commercial genotyping platforms in African humans not included in the HapMap data (Lachance and Tishkoff, 2013). Altogether, these arguments suggest that re-sequencing data is the optimal choice in SS studies in cattle. At some extent, the HD assay is appropriate, as it has a high-density coverage of the genome with SNPs that are less biased than competing panels.

Another important source of confounding comes from the methods available to detect SS. First, all methods assume that individuals have no recent relationships in their pedigrees, a condition that is hardly true and generally ignored. It is essential to filter the data for cryptic relationships and assure to include only samples that are unrelated for at least two generations. Second, most of the methods rely on haplotypes and SNP coordinates, so further improvement of phasing strategies and of the bovine reference genome assembly is crucial to assure high quality results. Third, variants can depart from neutrality not only due to positive selection, but also as consequence of demographic events such as bottlenecks, genetic drift and admixture. Distinguishing loci under selection from neutrally evolving loci remains a major challenge in the field, and will require refinement of existing methods and development of new tests. Nevertheless, combining signals across different methods seems to be a promising approach to mitigate the individual methodological limitations. Also, when available, the concomitant analysis of environmental data (e.g., temperature, humidity, precipitation, disease prevalence, etc.) may be of great help in distinguishing true positives and accelerating the link between signal and phenotype (Lv et al., 2014).

Well-planned study designs will be crucial to exploit the full potential of SS in the detection of large effect mutations favored by selection. The identification of common adaptive phenotypes, together with geographical information data, should be an important player in sampling and decisions of population comparisons. Cattle breeds that are not highly productive but that exhibit genetic local adaptation should be considered as priority targets, as their environmental fitness was probably forged by hundreds of years of natural and artificial selection. In the context of artificial selection for complex traits, as large cattle pedigree cohorts for genomic selection become available, it will be soon possible to actually assess rapid changes in allele frequency using historical data, rather than present date data only. First demonstrations of such ideas were presented by Decker et al. (2012), and are likely to be incorporated as routine monitoring tools of genomic resources in breeding programs in the future.

Recently, comparing candidate loci across GWA studies has been facilitated in cattle with the advent of the Animal QTLdb (Hu and Reecy, 2007). Similarly, results from a SS scan on the human 1000 Genomes data reported by Grossman et al. (2013) were made available through the CMS viewer tool (available at: http://www.broadinstitute.org/mpg/cmsviewer/). Pybus et al. (2014) have also launched a comprehensive database of SS in the 1000 Genomes data (available at: http://hsb.upf.edu/). The research community would highly benefit from the development of a SS database for livestock species, which would not only facilitate cross-referencing, but would also help researchers willing to dig deep into the functional meaning of the signals to select promising candidates emerging from multiple preexisting studies.

Finally, similarly to the recent developments in human SS (Kamberov et al., 2013), unraveling the functional relevance of the putative selected variants will demand interdisciplinary reasoning, compilation of a wide range of data types (e.g., transcriptomic, proteomic), and assemblage of an arsenal of *post-hoc* assays, such as genomic editing, culture, phenotyping and challenge of specific cell lines, production of knock-out models, and generation of cross-bred lines for confirmatory segregation analyses.

## **ACKNOWLEDGMENTS**

The first author is supported by São Paulo Research Foundation (FAPESP) - process 2014/01095-8. Mention of trade name proprietary product or specified equipment in this article is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the authors or their respective institutions.

#### **REFERENCES**

Ajmone-Marsan, P. (2010). On the origin of cattle: how aurochs became cattle and colonized the world. *Evol. Anthropol. Issues News Rev.* 19, 148–157. doi: 10.1002/evan.20267

highlights the importance of the slick locus in tropical adaptation. *PLoS ONE* 7:e36133. doi: 10.1371/journal.pone.0036133

in Northern European populations.*Eur. J. Hum. Genet*. 18, 471–478. doi: 10.1038/ejhg.2009.184

Wright, S. (1950). Genetical structure of populations. *Nature* 166, 247–249. doi: 10.1038/166247a0

Yang, J., Lee, S. H., Goddard, M. E., and Visscher, P. M. (2011). GCTA: a tool for genome-wide complex trait analysis. *Am. J. Hum. Genet.* 88, 76–82. doi: 10.1016/j.ajhg.2010.11.011

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 31 October 2014; accepted: 26 January 2015; published online: 10 February 2015.*

*Citation: Utsunomiya YT, Pérez O'Brien AM, Sonstegard TS, Sölkner J and Garcia JF (2015) Genomic data as the "hitchhiker's guide" to cattle adaptation: tracking the milestones of past selection in the bovine genome. Front. Genet. 6:36. doi: 10.3389/ fgene.2015.00036*

*This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics.*

*Copyright © 2015 Utsunomiya, Pérez O'Brien, Sonstegard, Sölkner and Garcia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Revisiting demographic processes in cattle with genome-wide population genetic analysis

Pablo Orozco-terWengel 1 \*, Mario Barbato1 , Ezequiel Nicolazzi 2 , Filippo Biscarini 2 , Marco Milanesi 3 , Wyn Davies 4 , Don Williams 5 , Alessandra Stella2 , Paolo Ajmone-Marsan3 and Michael W. Bruford1

1 School of Biosciences, Cardiff University, Cardiff, UK, 2 Parco Tecnologico Padano, Lodi, Italy, 3 Faculty of Agriculture, Università Cattolica del Sacro Cuore, Piacenza, Italy, 4 Dinefwr, National Trust, Llandeilo, UK, 5 Independent Researcher, Swansea, UK

#### Edited by:

Peter Dovc, University of Ljubljana, Slovenia

#### Reviewed by:

Ino Curik, University of Zagreb, Croatia Lingyang Xu, University of Maryland, USA

#### \*Correspondence:

Pablo Orozco-terWengel, School of Biosciences, Cardiff University, Sir Martin Evans Building, Museum Avenue, Cardiff CF10 3AX, UK<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics

Received: 12 February 2015 Paper pending published: 26 March 2015 Accepted: 12 May 2015 Published: 02 June 2015

#### Citation:

Orozco-terWengel P, Barbato M, Nicolazzi E, Biscarini F, Milanesi M, Davies W, Williams D, Stella A, Ajmone-Marsan P and Bruford MW (2015) Revisiting demographic processes in cattle with genome-wide population genetic analysis. Front. Genet. 6:191. doi: 10.3389/fgene.2015.00191 The domestication of the aurochs took place approximately 10,000 years ago giving rise to the two main types of domestic cattle known today, taurine (Bos taurus) domesticated somewhere on or near the Fertile Crescent, and indicine (Bos indicus) domesticated in the Indus Valley. However, although cattle have historically played a prominent role in human society the exact origin of many extant breeds is not well known. Here we used a combination of medium and high-density Illumina Bovine SNP arrays (i.e., ∼54,000 and ∼770,000 SNPs, respectively), genotyped for over 1300 animals representing 56 cattle breeds, to describe the relationships among major European cattle breeds and detect patterns of admixture among them. Our results suggest modern cross-breeding and ancient hybridisation events have both played an important role, including with animals of indicine origin. We use these data to identify signatures of selection reflecting both domestication (hypothesized to produce a common signature across breeds) and local adaptation (predicted to exhibit a signature of selection unique to a single breed or group of related breeds with a common history) to uncover additional demographic complexity of modern European cattle.

#### Keywords: cattle, Bos taurus, Bos indicus, phylogeography, adaptation, SNP array

## Introduction

Modern cattle comprise two species, Bos taurus and Bos indicus, both derived from the extinct aurochs (B. primigenius). The divergence between the two species has been dated ∼250,000 years ago on the basis of mitochondrial DNA haplotypes (Bradley et al., 1996). The two species are interfertile and a number of hybrid breeds have been established (e.g., the American Beefmaster that was established from a Shorthorn × Hereford × Brahman cross). The details of cattle domestication have been a contentious, and arguments in favor of two or three domestication events have been put forward. The hypothesis arguing for two domestication sites suggests that taurine cattle were domesticated in the Fertile Crescent, likely from B. primigenius, while indicine cattle were domesticated in the Indus valley, likely from B. primigenius namadicus (Loftus et al., 1994; Troy et al., 2001; Lewis et al., 2011). The hypothesis arguing for three domestication sites adds a Northeast African center, that gave rise to the African taurine breeds deriving form B. primigenius opisthonomous (Wendorf and Schild, 1994; Payne and Hodges, 1997). Mitochondrial DNA evidence suggests that two main lineages were domesticated, one now typical of taurine breeds, and the other now typical of indicines (Loftus et al., 1994; Bradley et al., 1996). The domestication process is thought to have originated sometime between 8000 and 10,000 years ago (Ajmone-Marsan et al., 2010; Gautier et al., 2010; Lewis et al., 2011). Several potential routes have been hypothesized to explain the dispersal of taurine cattle from the Fertile crescent, mediated by Neolithic farmers, including alternative routes into Europe, occurring on at least two occasions (a Danubian route and a Mediterranean route Payne and Hodges, 1997; Pellecchia et al., 2007; Ajmone-Marsan et al., 2010; Taberlet et al., 2011). As aurochs only became extinct recently (Machugh et al., 1998; Beja-Pereira et al., 2006), it is possible that both taurine and indicine cattle crossed with their wild ancestor during their expansion out of the domestication centers, although the evidence supporting it is much debated (Troy et al., 2001; Gotherstrom et al., 2005; Beja-Pereira et al., 2006; Bollongino et al., 2008; Perez-Pardal et al., 2010; Edwards et al., 2011).

Several studies have addressed the population structure and phylogenetic relationships between taurine and indicine cattle breeds. Initially, studies were based on allozymes, mitochondrial DNA sequences, and microsatellites (Loftus et al., 1994; Medjugorac et al., 1994; Bradley et al., 1996; Machugh et al., 1998; Kantanen et al., 2000; Hanotte et al., 2002; Gautier et al., 2007) describing patterns of genetic variation in cattle, e.g., the two main lineages and differentiation between taurine and indicine breeds, albeit largely focusing on taurines. Nevertheless, such studies identified a clear differentiation between breeds from the two species, and between African and European taurines (Gautier et al., 2007). Interestingly, this handful of markers separated the African taurine breeds (e.g., Lagune, N'Dama and Somba) living in tsetse fly endemic regions from those that did not co-occur with these flies (e.g., Kuri and Borgou), suggesting that demographic signals among these breeds are strong and possibly related to properties of the environment they occupy. Recently, surveys of genome-wide genetic variation in many breeds have been carried out using SNP arrays (Bovine Hapmap et al., 2009; Mctavish et al., 2013; Decker et al., 2014). These more comprehensive analysesincluded tens of thousands of SNPs distributed across all chromosomes genotyped in several breeds from around the world. With these data, the differentiation between taurine and indicine breeds has been measured at approximately 10%, with approximately 5% for the divergence between African and non-African taurines (Gautier et al., 2010; Mctavish et al., 2013; Decker et al., 2014). Analyses using this type of large-scale data further identified a split between African longhorn taurines (N'Dama) and shorthorn taurines (e.g., Somba and Lagune—Gautier et al., 2010). Such analyses have also proved powerful enough to identify the mixed indicine/taurine genetic component of recently established hybrid breeds such as Santa Gertrudis and Beefmaster (Gautier et al., 2010; Mctavish et al., 2013).

Cattle are the most common large livestock species in the world; the global population size of which is approximately 1400 million animals (Felius et al., 2011; Taberlet et al., 2011). Among these, 1311 breeds have been recognized, 209 of which are now extinct (Rischkowsky and Pilling, 2007). Since domestication, cattle have spread around the world, with farmers selecting animals with desirable characteristics and establishing many breeds. This process, while not passive, was until recently slow and resulted in animals becoming adapted to the local conditions, e.g., feed types, local weather, and diseases (Wood, 1973; Cobb, 2006; Russell, 2007). The last two centuries saw a dramatic acceleration in the process of artificial/directional selection, largely as a consequence of the development of the breed concept, artificial insemination, and the development of statistical approaches used to estimate breeding values used to choose sires and dams for breeding purpose. This has limited the intercrossing between animals of divergent genetic background in order to increase or maintain breed purity, especially in dairy cattle [in beef cattle, conversely, cross-breeding is often used to exploit also non-additive genetic variation i.e., heterosis (Simm, 1998); still, pure lines for cross-breeding need to be selected and maintained]. Additionally, as a consequence of artificial insemination and marketing policies, few "champion" males with characteristics of interest have been used extensively to fertilize large proportions of females. The above developments have led to: (i) intense within-breed (or line) selection; (ii) importation/exportation of elite breeding stock to different environments/countries (migration); (iii) strong founder effect of "champion/top ranked" males both in breed formation or in already established breeds whose effective population size has thus been shrunken; (iv) larger genetic drift as a consequence of breed isolation. Generally, by trying to keep breeds pure, gene-flow among them has been reduced or even halted (unless done to improve breeds by upgrading Felius et al., 2011), and by using relatively few males to fertilize many females, inbreeding is increased (Taberlet et al., 2011). These approaches are used to manipulate the properties of a breed at the expense of the decrease in the genetic variability (Bulmer, 1971) which may be accompanied by a decrease in genetic health of the breed (e.g., reduction in fertility in Holstein cattle via the increase in frequency of deleterious variants through the process of inbreeding) (Pryce et al., 1997, 2004; Biffani et al., 2005). Overall, this change in farming practices resulted in a change from limited selection throughout much of the domestic history of cattle (∼9800 years) to relatively stronger selection in the last ∼200 years.

The relatively rapid distribution of cattle around the world over the last ∼10,000 years suggests that populations have been potentially exposed to selective pressures deriving from the environmental variables that they had not been exposed previously, e.g., diseases occurring outside of the domestication centers (Hanotte et al., 2002; Felius et al., 2014; Xu et al., 2014). Consequently, while some cattle breeds represent recent introductions to new geographical areas such as the New World or Oceania, other breeds represent long established populations outside of the domestication center. The genetic makeup of such populations is a valuable source to understand the processes driving local adaptation, as the former represent genetic pools currently being shaped by local selective pressures (i.e., may serve as examples of selection in action) (Flori et al., 2012), while the latter represent established populations where the extant genetic variation reflects the outcome of local adaptation (Felius et al., 2014). In particular, long established local breeds are likely to represent the reservoirs of important genetic variation with adaptive potential that needs to be characterized (Taberlet et al., 2011; Felius et al., 2014) before it possibly disappears as a consequence of the widespread use of a handful of industrial breeds among others (Fao, 2007; Thornton, 2010; Herrero and Thornton, 2013; Perry et al., 2013).

Here, we collated a large-scale genome-wide dataset of SNP polymorphism in cattle breeds from around the world representing published taurine and indicine breeds, and new data for four European taurine breeds, all genotyped for a common set of ∼35,000 SNPs in over 1300 animals. We used these data to characterize the distribution of genetic variation among cattle breeds around the world, and the genealogical relationships between these breeds. The dataset described here represents multiple breeds from various continents (e.g., multiple European and African taurines, as well as Asian and African indicines) providing a unique experimental set up that allows searching for signatures of selection that are consistent across multiple breeds from different continents. This approach enables to separate signatures of selection that may be breed specific from selection signals common for several breeds occurring in the same geographic region. We extend the published dataset with new SNP data to a case study where we compare white cattle breeds of conservation concern in the UK (White Park and Chillingham) with a candidate Roman ancestor, the Chianina, to help resolve the commonly held belief of a Roman origin of the Welsh White Park cattle (Felius et al., 2011).

## Materials and Methods

## SNP Array Data

SNP array data was obtained for a total of 1346 animals representing 46 cattle breeds, and four species related to cattle (B. gaurus—Gaur, B. javanicus—Banteng, B. gruniensis— Yak, and Bison bison—Bison) used as outgroups (**Table 1**). The data comprises a combination of previously published datasets and new data collected for this study. The published data represents breeds published by Gautier et al. (2009); Gautier et al. (2010) and the Bovine HapMap project (Bovine Hapmap et al., 2009) genotyped with Illumina's BovineSNP50 v.1 and v.2 chip assay, respectively (Bovine Hapmap et al., 2009) (**Table 1**). Four new breeds were genotyped with the BovineHD Beach Chip for this project, namely, Welsh While Park (Dinefwr), Chianina, Romagnola and Chillingham. All SNP coordinates were converted to UMD3.1 bovine assembly (RefSeq: GCF\_000003055.5). The datasets were merged using PLINK v1.7 (Purcell et al., 2007). The merged dataset was filtered to keep only individuals with 95% of their SNPs called and SNPs with 95% call rate across all samples. The SNPs and individuals that did not pass these filters were discarded. The total SNPs left in the dataset was 36,503 (data is available upon request from the authors). The data for each breed was phased with fastPHASE v1.4 (Scheet and Stephens, 2006).

## Population Divergence

An analysis of population structure was carried out using the software Admixture v1.22 (Alexander et al., 2009), which uses a model-based estimation of indivdual ancestry for a range of prior values of K defined by the user. A cross-validation approach is used to determine the most likely number of populations (K) in the data. For each tested value of K, Admixture estimates the proportion of each individual's genotype deriving from each cluster. The values of K tested were in the range between 1 and 60 to accommodate potential population structure within breeds. For this analysis we excluded outgroup samples, breeds represented by fewer than 10 individuals, and SNPs with a minimum allele frequency lower than 0.01, linkage disequilibrium higher than 0.1 using a sliding window approach of 50 SNPs and a step size of 10 SNPs (Alexander et al., 2009). A principal components analysis was carried out using flashpca v1.2 (Abraham and Inouye, 2014) with default settings in order to investigate the ordinal relationships among breeds and individuals. Lastly, a NeighbourNet network was constructed using Reynold's distance among breeds using Splitstree v4.13.1 (Huson and Bryant, 2006), and a dendrogram using the same genetic distance and 100 bootstrap replicates to assess the statistical support of breed; clustering was performed using

Phylip v3.69 (Felsenstein, 1993). Within this framework, the Welsh White Park cattle breed was analyzed along the other cattle breeds to identify their relative position within the dendrogram and PCA explaining the similarities between cattle breeds. In particular, the common belief of a Roman origin of the Welsh White Park was addressed by comparing the position of this breed in relation to the Italian breeds, in particular the large Chianina breed, one of the oldest cattle breeds originating from the region of Valdichiana in central Italy, and Romagnola. Additionally, we estimated the relationships between Welsh White Park and other 17 breeds in the dataset using Treemix v1.12 (Pickrell and Pritchard, 2012) in order to determine the historical relationships between these breeds in terms of splits and migration (mixtures) between breeds. The breeds used for this analysis were Welsh White Park, the Italian Chianina, Romagnola (ROM), and Piedmontese, the British breeds Chillingham, Hereford, Angus, Jersey (JER), as well as five other non-African taurines (Holstein (HOL), Brown Swiss, Charolais, Normande, Vosgienne), the four African taurines Baoule, Lagune, Somba and N'Dama (NDA), and the indicine breed Brahman as outgroup. Treemix was run iteratively for values of the migration parameter between 0 and 12. The f index representing the fraction of the variance in the sample covariance matrix (W) accounted for by the model covariance b matrix (W) was used to identify the number of modeled migration events that best fitted the data (Pickrell and Pritchard, 2012).

## Demographic Estimation

Recent demographic history was measured by the trend in effective population size (Ne) change over the last 8000 years using the software SNeP v1.0 and default settings (Barbato et al., 2015). Loci with missing data were excluded from this analysis, as well as SNPs with a minimum allele frequency (MAF) smaller than 5% (Sved et al., 2008; Corbin et al., 2012). Linkage disequilibrium (LD) was calculated between each pair of SNPs separated by a minimum distance of 5 kbp and maximum distance of 1 Mbp using Hill and Robertson's r 2

#### TABLE 1 | Breed abbreviations and species identification.

(Continued)

#### TABLE 1 | Continued

The three letter code used to identify the breed throughout the manuscript is shown (Abb), along with the breed name (BN) and the taxonomic status indicating the species name or whether the breed is a B. taurus x B. indicus hybrid (Tax), the data set's origin (D), namely Bovine Hapmap Consortium (1), Gautier et al., 2010 (2), this study (3), the breed's sample size (N), the observed heterozygosity (Ho), and the inbreeding coefficient (F).

(Hill and Robertson, 1968). LD values were grouped in 30 distances bins, such that each bin represented approximately the same number of pairwise LD estimates used to estimate Ne (Barbato et al., 2015). Within each bin the mean r 2 was calculated and used to estimate Nt = 1/(4f(c)) (1/r 2 − 1), where f(c) = c [(1−c/2)/(1−2)2 ] (Sved, 1971); c is the linkage distance inferred from the physical distance between SNPs assuming 1 Mb ∼ 1 cM and Nt represents the effective population size estimate at t = 1/2c generations ago.

#### Signatures of Selection

As the dataset analyzed here comprised B. taurus breeds from Europe and Africa, as well as B. indicus breeds from Asia and Africa, we sought to identify signatures of selection in a replicated manner (**Figure 1**). Initially we compared four taurine breeds (two European—Gascon and Chianina; and two African—Baoule and N'Dama) with four indicine breeds (two African – Zebu Fulani and Zebu Bororo; and the Asian Gir and Brasilian Nelore; Nelore was recently derived from Indian Ongole cattle, therefore, we broadly refer to these breeds as Asian indicine) to identify signatures of selection between B. taurus and B. indicus. Within B. taurus, we randomly chose three European breeds to compare against three African breeds in order to identify selection signatures specific to breeds in these continents. In the same way, for B. indicus we compared three Asian breeds against three African breeds (**Figure 1**). In contrast to previous studies, this experimental design enabled us identifying potential signatures of selection reflecting adaption to the local environment, instead of breed specific signatures that potentially reflect the breed's particular history.

The test of selection used was the cross population extended haplotype homozygosity test (XP-EHH) calculated for each SNP in the dataset as implemented in Selscan v1.0 (Szpiech and Hernandez, 2014). Selscan was used with default settings except for the maximum gap in bp allowed between SNPs, which was set to 400,000 to address the BovineSNP50 SNP density. For the comparisons between taurine and indicine breeds the following pairwise combinations were used: N'Dama (African taurine)/Gir (Asiatic indice), Gascone (European taurine)/Zebu Bororo (African indicine), Baoule (African taurine)/Zebu Fulani

(African indicine) and Chianina (European taurine)/Nelore (Asian indicine). For the within-taurine comparison three pairwise tests were done between a European breed and an African breed, and for the within-indicine comparison three pairwise comparisons were carried out between an Asiatic breed and an African breed (**Figure 1**). For each pairwise comparison, the XP-EHH results where standardized for each chromosome separately, and the 5 and 95% quantiles of the standardized XP-EHH distribution for all markers where used as threshold to identify outlier SNPs, i.e., SNPs with XP-EHH values equal or smaller (more negative) or equal or larger (more positive) than the chosen thresholds. As multiple XP-EHH were carried out for each comparison (e.g., four comparisons between taurine and indicine breeds) this analysis was the equivalent of a replicated test. We took advantage of this design to identify SNPs under selection as those that showed for several pairwise comparisons an extreme XP-EHH with the same sign (e.g., all pairwise comparisons should an extreme positive XP-EHH).

For the between species comparison we defined the consistent SNPs under selection as those showing the same XP-EHH trend in at least three of the four pairwise comparisons. We did not restrict ourselves to only those SNPs where the four tests showed the same result to avoid excluding SNPs that failed to pass the 5% MAF threshold within a particular breed but which showed a consistent signature of selection across the remaining three comparisons. For the within species comparison of breeds occurring on different continents we identified consistent SNPs under selection as those showing an extreme XP-EHH with the same sign in each of the three comparisons made. The SNPs showing consistent signatures of selection across multiple pairwise comparisons where linked to neighboring genes using a window approach of 50 K base pairs, and a Gene Ontology analysis was carried out on these SNPs using Gorilla (www. cblgorilla.cs.technion.ac.il; Eden et al., 2009).

## Results

Of the 36,503 SNPs analyzed in 1345 animals representing taurine and indicine breeds, the observed heterozygosity varied between a minimum of 0.026 (Chillingham) and a maximum of 0.33 (the hybrid breed Beefmaster; **Table 1**). The comparison of expected heterozygosity (unbiased to sample size) between the taurine and indicine breeds was not significant, as well as the comparison between breeds in and out of Africa (all Welch t-test p-values > 0. 2). However, the average observed heterozygosity in taurine breeds was slightly larger than that in indicine breeds [0.28 standard deviation (sd) 0.05 and 0.20 sd 0.04, respectively]. The average expected heterozygosity across breeds in both datasets was similar, i.e., 0.34 sd 0.0001. The inbreeding coefficient (FIS) was significantly lower for taurine breeds (FIS: 0.18 sd 0.16) than for indicine breeds (FIS: 0.42 sd 0.11; Welch t-test p-value: 0.0011). Within B. taurus, the African breeds presented a significantly higher inbreeding coefficient than the non-African taurines (Welch t-test p-value: 1.54×10−7 ). Similarly, the non-African taurines also presented significantly lower FIS values than either of the groups of indicine breeds (Welch t-test p-value: 1.84 × 10−5 and 0.026, respectively). Among the non-African taurines, the Chillingham presented by far the highest inbreeding coefficient of the entire dataset (FIS: 0.92), on average five-fold larger than other taurines, and twofold larger than the indicine breeds. The Welsh White Park is thought to be related to the Chillingham breed, but presented an observed heterozygosity (0.245) only slightly smaller than the average value observed in the other non-African taurines. The FIS value for Welsh White Park (0.27) was less than a third of that of Chillingham but still almost twice the average across the other non-African taurine breeds (0.14 sd 0.14). With the exception of Chillingham, the four breeds genotyped here presented typical observed heterozygosity and FIS values for European taurine breeds (**Table 1**). Removing from the analysis the breeds with samples sizes smaller than 20 did not change the results (results not shown).

## Admixture and Genetic Relationships

Admixture analysis between taurine and indicine breeds was initially run for values of the number of clusters (K) between 1 and 60. The cross validation (CV) statistic used to choose the most suitable number of clusters had its lowest value at K = 60. However, the shape of the CV curve suggested that higher K values may present lower CV. Consequently, we carried out runs of Admixture for larger values of K (e.g., 80). However, the CV values continued to reduce the larger the K numbers tested (SM Figure 1). While it may be possible that some family structure within breeds explains this trend, the substantially longer run time needed by the algorithm for such K values combined with the small decrease in CV suggests that clustering solutions larger than K ∼60 potentially represent spurious results (a similar observation was made by Mctavish et al. (2013) who only achieved marginal increases on the likelihood scores using structure in their cattle dataset for values of K beyond 3).

The hierarchical clustering analysis with Admixture showed that the largest difference for K = 2 was between taurine and indicine breeds. Conditioning the data to K = 3 separated the taurine breeds with an African origin from the non-African, suggesting a large differentiation between taurine breeds in these two groups of samples (**Figure 2**). Contrastingly, for values up to K = 10, the African indicine breeds presented a remarkably similar distribution of genetic variation to that of the Asian indicine breeds, which did not form a separate cluster. Nonetheless, the African indicine breeds consistently presented a distinctive proportion of their genotype representing taurine ancestry throughout the analysis (a pattern almost entirely absent in the breeds of Asian origin). The B. taurus/B. indicus hybrids show a genetic make-up consisting of non-African taurine and B. indicus genetic components (Beefmaster and Santa Gertrudis), or African taurine and B. indicus (Borgou, Kuri, and Sheko).

The Admixture analysis was complemented with a principal component analysis (**Figure 3**). The first principal component (PC1) explained ∼10% of the total variation in the data and separated taurine form indicine breeds. Besides of PC1, the only other component explaining a substantial part of the variance was PC2 (∼5%), which separates the African taurine from the remaining taurine, as well as the African indicine from the Asian indicine breeds. PC3 to PC7 explained between less than 2 and 1.4% of the variance, while the remaining components explained less than 1% of the total variance, thus were not taken into consideration. PC3 separated Holstein from the remaining taurine breeds. PC4 identified a group of British breeds formed by Angus, Red Angus and Hereford, and PC5 identified a component of variation specific to Hereford separating it from the other taurines. PC6 identified the Chillingham cattle as well as the Jersey breed, and PC7 identified a specific component of variation to Chillingham separating it from the remaining taurine breeds. The combination of PC1 and PC2 identify the major groups of individuals described here in terms of species and geographic component (**Figure 3A**) with the hybrid breeds occurring between these. As expected from the Admixture analysis, the breeds Kuri, Borgou, and

FIGURE 2 | Analysis of Population Structure. Results of the analysis of population structure conditioning the dataset to 3 clusters (top row) and to 51 clusters (bottom row). Each animal is represented by a straight bar that is colored. The amount of a color reflects the individual's proportion of genetic variation originating in the cluster of

that color. Each breed is labeled in the center of its box on the bottom of each plot. In the top row the non-African taurine breeds are labeled in blue, the African taurine in red and the indicine (both African and not African in green). The new datasets produced for this project are highlighted in green.

Sheko are placed with these principal components between the group of African taurine and indicine. Similarly, Beef Master and Santa Gertrudis are placed between the group of non-African taurine and indicine breeds, although from this analysis it is not clear whether their indicine genetic component is more similar to African indicine or Asian indicine.

The methods described above identified differences between breeds but did not allow us to determine how these relate to each other. We therefore estimated a NeighbourNet network with these breeds using Reynold's distance. The torso of the network shows multiple pathways of connection between breeds reflecting the relatively recent divergence between many populations (**Figure 4**). However, the topology of the network resulted in the pattern expected from the previous analyses (Admixture and PCA), where a clear separation between taurine and indicine occurs, as well as the distinct separation between the African taurine breeds from the non-African taurines. Similarly to the PCA and in contrast to the Bayesian approach of Admixture, the network also resolved the African indicine breeds separately from their Asian counterparts. The dendrogram depicting the clustering of breeds using Reynold's distance recovered the grouping of breeds identified in the PCA with most branches (83%) supported with high statistical support (bootstrap values higher than 70%; SM Figure 2). The clustering of the outgroups was not well resolved, except for the grouping between bison and yak. The indicine breeds clustered together with the Asian indicine showing a greater similarity to each other, and then to the Madagascan zebu. The two African indicine clustered together and then clustered with the remaining indicines. The taurine breeds separated into two clusters, one representing the African taurines and the other the non-African taurines. The hybrids Beefmaster and Santa Gertrudis clustered together within the non-African taurine group reflecting their larger taurine genetic component. Contrastingly, the African hybrids (Sheko, Borgou, and Kuri) clustered with the indicine breeds suggesting a higher indicine genetic component in their genotype.

A comparison between the Welsh White Park cattle and the other British breeds (Chillingham, Hereford, and Jersey) and Italian breeds (Chianina, Romagnola, Piedmontese) was carried out to establish whether the UK white cattle have a genetic link with Roman white cattle, as is popularly believed. The Principal Components Analysis of this datasets identified the White Park breed clustering near the other British breeds (i.e., HFD, ANG, RGU), instead of closer to the Italian breeds (**Figure 3B**), as would be expected if White Park cattle had a strong Italian genetic

component due to its speculated Roman origin. Similarly, in the dendrogram, the Italian breeds clustered together and separately from the Welsh White Park and Chillingham breeds, the latter two not clustering with any of the other taurine breeds suggesting their genetic distinctiveness. The Treemix analysis revealed a similar topology with the Welsh White Park and Chillingham clustering together and then clustering to Hereford and Angus, rather than clustering with the Italian breeds (**Figure 5**, SM Figures 3, 4). Adding up to seven or eight migration edges to the tree topology improved the amount of the variance explained by the phylogenetic model, albeit marginally (f statistic ∼ 0.99965; **Figure 5**; Pickrell and Pritchard, 2012). The model without migration edges already presents an f statistic of 0.9981 (**Figure 5** inset), a value above the threshold described previously, and above which the phylogram was not better explained by adding additional migration events (Pickrell and Pritchard, 2012; Decker et al., 2014). However, more migration edges did not seem to further increase this variance as the f statistic reached an asymptote between 7 and 8 migration edges (**Figure 5** inset). The largest increase in variance occurred by adding the first migration edge to the graph from Brahman to the ancestor of Chianina and Romagnola, consistently with the Admixture results that suggest these Italian breeds have a minor (∼10%) indicine genetic component. With the exception of two edges linking Brahman, Angus and Somba, the remaining migration edges were between taurine breeds. All migration edges had a weight less than ∼0.2, with the exception of those from Jersey into Hereford, and from the ancestor of Hereford, White Park and Chillinghan, into Charolais, both of which have a weight above 0.4, suggesting the source population made a substantial genetic contribution to the recipient breeds.

## Demographic History

The trend in historic effective population sizes for each breed was estimated using the distribution of linkage disequilibrium across the genome. All breeds showed a declining effective population size over the last ∼2000 generations (i.e., approximately 8000 years assuming an average generation length of 4 years throughout these species history; **Figure 6**). Interestingly the African breeds showed an average higher effective population size for each species than the non-African breeds. Consistent with their hybrid nature, the five hybrid breeds presented a larger effective population size, probably reflecting the artificial increase in heterozygosity deriving from the admixture event (**Figure 6**).

## Signals of Selection

A replicated approach using XP-EHH was used to identify consistent SNPs showing signatures of selection in the dataset (**Figure 1**). Four pairwise comparisons between taurine and indicine breeds were carried out resulting in 10,150 SNPs that passed the MAF threshold within population and which were tested in at least three out of the four pairwise comparisons. Among these SNPs, 3029 (∼30%) showed consistent signatures of selection in the taurine breeds and 2385 (∼24%) in the indicine breeds. The comparison between African and non-African indicine breeds resulted in 9388 SNPs tested, out of which 1768 (∼19%) showed signatures of selection specific to the Asian breeds and 2049 (∼22%) showed signatures of selection specific to the African breeds. For the comparison between taurine breeds, a similar number of SNPs was tested, 9938 SNPs. Of these 1872 (∼19%) showed signatures of selection in European cattle breeds and 1739 (∼17%) in African breeds. The SNPs showing consistent signatures of selection were linked to neighboring

genes that occurred up to 50,000 bp downstream or upstream of the SNP position. No significantly enriched GO categories were found for any of the sets of breeds analyzed here. However, some patterns were observed, e.g., for the candidate SNPs under selection in Asian indicine several genes related to immunity were identified, e.g., LIPH (involved in platelet aggregation), Tespa1 (involve in thymocyte development), CCL14 (activates monocytes) (SM Table 1). However, interestingly, several RNA types were identified for several of the groups of breeds tested (**Table 2**). For the comparison between taurine and indicine breeds 6 snoRNA and 10 snRNA were identified among the taurine breeds, and two miRNA and 1 snoRNA in the indicine breeds. For the comparison between Asian and African indicine, a total of three miRNA, 1 snoRNA and 2 snRNA were identified for the Asian breeds and none for the African, while for the taurine comparison, the opposite was true, i.e., 1 snoRNA and 1 snRNA were identified for the African breeds and none for the European breeds. Overall this suggests that a substantial amount of selection signature occurred near non-coding RNAs that modify transcription of genes or of other RNAs.

This dataset represents one of the most comprehensive analyses of B. taurus and B. indicus cattle breed genetic data. Overall a higher observed genetic variation was observed in taurine breeds that originated or derived from European cattle breeds (∼0.29), than in the other cattle breeds. Both African taurine and indicine breeds presented an intermediate level of observed heterozygosity (∼0.22 in both sets) and the Asian indicine breeds presented the lowest observed heterozygosity (∼0.17). While the observed difference could have a generic demographic explanation, e.g., a larger effective population size in the taurine breeds of European descent, such pattern is not consistent with what has been described with other types of molecular markers such as microsatellites (Rothammer et al., 2013). The most likely explanation for this pattern is ascertainment bias (Gautier et al., 2009; Matukumalli et al., 2009). In particular, for those breeds not included in the panel of breeds used to identify SNPs in the BovineSNP50 chip assay, minimum allele frequencies are expected to be substantially lower and the genetic variation captured with the chosen SNPs may not be reflective of the true underlying genetic variation in those breeds. In this case, SNP discovery was based on the genome sequence reads of five taurine breeds (Holstein, Angus, Limousin, Jersey, Norwegian Red) and one indicine (Brahman) breed compared to the reference genome of a Hereford cow (Bovine Hapmap et al., 2009). Ascertainment bias deriving from the design of this SNP array was previously described for a subset of the dataset analyzed here, where up to 30% of the SNPs had a MAF > 0.3 in taurine breeds, while only 19% of the SNPs had an equivalent MAF in indicine breeds (Nelore, Brahman, and Gir) (Bovine Hapmap et al., 2009). Similarly, the analysis of Holstein, Angus and Brahman samples carried out with a subset of the SNPs in the BovineSNP50 chip assay ascertained in each of the breeds separately, showed how Brahman animals presented a lower heterozygosity using the taurine SNPs (pvalue < 0.05 and p-value = 0.07, respectively) than using the indicine SNPs (Neto and Barendse, 2010). Contrastingly, in the same analysis, for the Holstein and Angus derived SNPs, the taurine samples presented an observed heterozygosity typical to that of any other taurine breed. Furthermore, as no African taurine breeds and most indicine breeds were not included in the SNP ascertainment, the observed levels of genetic variation are likely to be an underestimate of the true genetic variation in these breeds (Matukumalli et al., 2009), and the associated inbreeding coefficient may consequently be an overestimate of the real ones. Thus, comparisons of genetic diversity discussed here are only likely to be meaningful within major genetic groups.

Among all cattle, the British Chillingham breed possessed unusually low observed heterozygosity (0.026) and an extremely high inbreeding coefficient (0.92). The Chillingham occurs in Northumberland, UK, and for the last ∼300 years has lived in a feral state with almost no human intervention and no immigration (Visscher et al., 2001). Additionally, this breed passed through a bottleneck in 1947 where the herd reduced to five males and eight females. Previous analysis of genetic variation using 25 microsatellites identified only one marker with

polymorphism (Visscher et al., 2001). Additionally, sequences of the full mitochondria of eight animals (sequencing depth ∼2935x using next generation sequencing) revealed a single haplotype in the herd that differentiated them from other taurine breeds by three mutations, and one mutation (likely to be a recurrent mutation) that related them to indicine breeds and yak (Hudson et al., 2012). These studies suggest that Chillingham probably survives because the population has purged most deleterious mutations. The results shown here represent the largest study of genetic variation in this breed so far and did find similarly low levels of genetic variation in the breed (∼900 polymorphic SNPs out of ∼35,000). An extended analysis including more individuals and genetic markers in this breed may identify genetic variation of evolutionary importance not captured in this study, as well as shedding light on the distribution of homozygous and heterozygous tracts in the breed that could be used to monitor its genetic health.

## Population Structure

The analyses of population structure easily identified the taurine and indicine genetic components in the dataset (**Figure 2**). However, when conditioning the number of clusters in the data to K = 2, the African taurine cattle presented a genetic composition similar to that of the American hybrid breeds Beef Master and Santa Gertrudis, the latter two are breeds of recent origin (less than 100 years old). Increasing the number of clusters to K = 3 clearly separates the African taurines in their own cluster, suggesting that despite the strong non-African taurine and indicine genetic component (as seen for K = 2), African taurines (both longhorn and shorthorn) have a substantial third genetic component specific to Africa which may derive from cross breeding with local aurochs (Mctavish et al., 2013). Consistent with previous observations, the African taurine N'Dama (ND1 and ND2) and Somba (SOM) breeds presented a minor indicine genetic component (Gautier et al., 2010). Among the hybrids, Beefmaster (BMA) and Santa Gertrudis (STG) presented a major genetic component of non-African taurine origin and the remaining of indicine origin, in accordance with their breed formation. Beefmaster and Santa Gertrudis are both hybrid breeds of North American origin resulting from the cross between European taurine cattle of British origin (Hereford and Shorthorn) and Brahman cattle (Mctavish et al.,

DNA sequences coding for molecules in close proximity to SNPs under selection were grouped as Protein Coding, RNA types (mi, sno, and sn), rRNA and pseudogenes. The number of SNPs under selection close to each of these classes is shown for each of the three analyses presented.

2013). On the other hand, the remaining three hybrid breeds, the Borgou, Kuri, and Sheko present a genetic make-up of African taurine and indicine origin, with the first two having a somewhat larger taurine proportion than the Sheko breed. Among the taurine breeds, the Moroccan Oulmès Zaer (OUL) breed of unknown origin, presented approximately ∼3/5 of its genotype of non-African taurine type and the remaining of African taurine type (**Figure 2**), consistent with previous results suggesting an influence of African longhorn taurine into the otherwise European taurine genetic background of this breed (Gautier et al., 2010; Decker et al., 2014). Interestingly, the non-African taurine cattle presented around 10% or less of their genetic make-up showing African taurine ancestry (**Figure 2**), possibly deriving form cross Mediterranean movement of cattle during the last millennia (Dürst, 1899; Do Vale, 1907; Felius et al., 2011). A further indicine genetic component was also detected in the French Charolais and Gascon breeds (<5%), and the Italian Romagnola and Chianina breeds (∼10%). These results support observations by Decker et al. (2014) that suggest indicine admixture in European taurine breeds (the bulk of non-African taurines in this study) is rare instead of widespread (Mctavish et al., 2013). The genetic variation of African taurine is clearly distinct from that of other cattle breeds, and the PCA analysis suggests that the three main cattle groups (African and non-African taurine, and indicine) are almost as different from each other. While these results do not provide evidence for or against the hypothesis of a third domestication event in Africa (Felius et al., 2011), it suggests that African taurine cattle present a unique characteristics in their gene-pool.

The analysis of population structure using Admixture resulted in asymptoting CV values reflecting problems in identifying the true number of clusters in the data, possibly due to the presence of family structure within the breeds, or relatively low divergence between the estimated clusters. Thus, we also carried out principal components analysis, and as expected from the problems with Admixture, there is relatively low variation among cattle groups, i.e., the largest principal component (PC1) only explains ∼10% of the variance in the dataset (separating taurine form indicine breeds), while PC2 explains ∼5% of the variance (separating African from non-African taurines). Consistent with the Admixture result, the combination of PC1 and PC2 places the hybrid breeds Kuri, Borgou and Sheko between the African taurine and African indicine clusters, with the Sheko much closer to the indicine cluster reflecting its larger indicine background than Kuri and Borgou. The American Beef Master and Santa Gertrudis are placed between the non-African taurine and Asian indicine, but somewhat closer to the taurine cluster reflecting the larger taurine component of their genetic variation.

Lastly, the NeighbourNet analysis and the dendrogram identify the Madagascan zebu as more similar to the Asian indicine than to the other two African indicine. While this placing may be biased due to the ascertainment bias in the SNP array, it is tempting to suggest that the Madagascan population share a more recent history with Asian indicine breeds than the continental African breeds (e.g., a more recent establishment of the Madagascan population).

## Roman Origin of Welsh White Park

Within the United Kingdom, traditional knowledge suggests that the Welsh White Park were brought to Britain by the Romans (Felius et al., 2011). While there is no evidence of the Romans having introduced cattle to Britain (Visscher et al., 2001), PCA places the White Park almost in the middle of the group of non-African taurine breeds (**Figure 3A**), consistently with the results by Decker et al. (2014) based on genotypes from five individuals. In particular, a comparison between Welsh White Park and the three Italian breeds in this dataset, Chianina, Romagnola and Piedmontese, revealed that the Italian breeds are more closely related to each other than to White Park, and the Welsh breed is more closely related to the other British breeds (**Figures 3B**, **5**). Nevertheless, while this suggests that Welsh White Park are not similar to the Italian breeds it still does not answer the question regarding their origin. Both the NeighbourNet (**Figure 4**) and the phylogenetic network (**Figure 5**) show that the Welsh White Park and Chillingham cluster near each other (**Figure 4**). Such observation could be the outcome of genetic drift or inbreeding (or an interplay of these two parameters) in these breeds driving the divergence of their allelic frequencies (particularly in the case of Chillingham). Consistently with this hypothesis, the phylogenetic network analysis modeling genetic drift suggests that Welsh White Park and Chillingham have been through substantial drift in comparison to other taurines (**Figure 5**). Additionally, the Treemix analysis suggests that no migration events connect either the Welsh White Park or Chillingham with any of the Italian breeds, unless an additional edge that marginally increases the f-statistic by 0.01% (i.e., from 99.96 to 99.97% of the variance explained) is added; however, that additional migration edge is from Welsh White Park into Romagnola (contrary to common believe of an Italian ancestry of White Park) and has a weight near 0 (SM Figure 4) indicating that a very small fraction (if at all) of the Romagnola allelic variation may be of Welsh origin. While such idea could be appealing to some, it does not necessarily reflect a real historic event as the inferred edge provides an almost negligible contribution to the variance in the phylogram and corresponds to a minimal contribution of Welsh allelic variants into Romagnola. Additionally, the question of whether the inferred patterns of migration by Treemix are somehow dependent on the data analyzed could be raised. Further studies on the genetic history of these iconic breeds should include additional breeds such as the Podolian Hungarian Gray, or breeds close to the domestication center in order to better characterize the origin of this animals.

### Demographic History

The demographic history of a population can be traced using the distribution of LD across the genome (e.g., Hill and Robertson, 1968), and if the cattle domestication event was marked by a strong bottleneck, it is expected that such a pattern should be visible with this dataset. However, it is important to note that the analyses presented here might provide a biased estimate of the true extent of demographic history, as LD calculated over SNP array markers is expected to decay at a lower rate than when the same estimates are obtained from whole-genome sequence data (Qanbari et al., 2014).

All the breeds analyzed present trends marked by a decline in their effective population size during the last ∼2000 generations (**Figure 6**). Assuming an average generation length of 4 years, that suggests a constant decline across all breeds since the onset of the domestication process approximately 8000 years ago (Gautier et al., 2010). These results are consistent with observations describing a general decrease in breeds such as Angus, Holstein, Hereford, N'Dama and Brahman over the same period of time (Gautier et al., 2007; De Roos et al., 2008; Flury et al., 2010). Interestingly, overall, the harmonic-mean effective population size of the non-African taurine breeds was smaller than that of the non-African indicine, but both of those means were much smaller than that of the African taurines and indicines (**Figure 6**). These results suggest that either crossing between African taurine and indicine breeds has increased the effective population size in these breeds by moving genetic variation between breeds, and/or that more benign management practices in Africa have resulted in a milder reduction in effective population size. Although the second hypothesis cannot be ruled out, consistent with the first hypothesis, the hybrid breeds present overall a higher average effective population size.

#### Signatures of Selection

The unique nature of the data presented here provided us with a potentially powerful statistical framework to detect signatures of selection. In contrast to other studies, we examined multiple breeds for each species, accounting for their distribution in and out of Africa to search for signatures of selection in the equivalent of a replicated experimental design, where instead of identifying selection in a particular breed, we used several breeds of each species in each continent as biological replicas. The test used is an extension of Sabeti's extended haplotype homozygosity (Sabeti et al., 2002), which is based on the principle that a site under positive selection will rapidly increase in frequency, bringing neighboring variants to a high frequency in a process called hitchhicking (Smith and Haigh, 1974), thereby removing genetic variation around the site under selection (a selective sweep) and leaving a track of extended homozygosity. The magnitude of the effect of the selected site on the neighboring variation largely depends on how strong and recent the selective pressure is. If the process is slow, because selection is weak, recombination events separate the selected site from neighboring variants breaking down linkage and consequently resulting in a short area of extended homozygosity (Sabeti et al., 2002). On the contrary, if the process is fast because selection is strong, the area of extended homozygosity will be large, as there is not enough time for recombination to break down linkage. Consequently, identifying old signatures of positive selection using the EHH or any of its derived versions (e.g., XP-EHH) is rather difficult as the expected pattern of linkage disequilibrium is likely to have disappeared (Sabeti et al., 2007; Chen et al., 2010). However, despite the ancient divergence between B. taurus and B. indicus, i.e., approximately 250,000 years old (Bradley et al., 1996; Bovine Hapmap et al., 2009), a consistent signature of selection was found for several thousands of SNPs (5385 SNPs) probably reflecting relatively recent selective events rather than ancient ones.

Consistent with the expectation for the comparison between breeds in both species, a similar number of SNPs in each species showed signatures of selection (i.e., 3029 in taurine and 2385 in indicine). For the African vs. non-African comparison of taurine breeds 1872 of the SNPs under selection were on European breeds and 1739 were in African breeds. In contrast, for the comparison of African vs. non-African indicine breeds, 1768 of the SNPs under selection were on the Asian indicine breeds and 2049 were in the African breeds. The difference between the comparison of taurine breeds and indicine breeds in and out of Africa may reflect the shorter period of time that the later have spent in Africa since their introduction to that continent. The first domestic cattle to have entered Africa were taurines approximately 7000 years go (Payne and Hodges, 1997; Hassan, 2000; Felius et al., 2014). These domestic animals probably entered via a route crossing through Egypt and moving into northern Africa. In contrast, indicine breeds entered Africa much later via routes leading to East Africa via Egypt some time between ∼3500 and 2500 years ago (Payne and Hodges, 1997; Hanotte et al., 2002). The timing of these migration events into Africa implies that taurine breeds have been in Africa between two and three times longer than indicine breeds.

At least three non-exclusive hypotheses can be derived from the taurine colonization history of Africa. Firstly, African taurine have become clearly differentiated from other non-African taurines as genetic drift exacerbated differences between the taurine propagules that spread from the domestication center into Europe and Africa. Secondly, African taurines had a longer time to hybridize with local aurochs, as possibly was the case for the European taurines, which hybridized with European aurochs (Decker et al., 2014). Consequently, African taurines are expected to carry unique genetic variation (deriving from local aurochs) that differentiates them from non-African taurines, as African aurochs were likely to genetically differ from European aurochs due to isolation by distance (Perez-Pardal et al., 2010; Mctavish et al., 2013). Thirdly, by having spent a longer time in Africa, African taurines are expected to have had a longer time to adapt to the local environment and although some of their signatures of local adaptation may have dissipated, they may have left a trace of genetic identity that helps differentiating African taurine form other cattle breeds. In contrast, B. indicus was introduced to Africa approximately ∼2000 years ago, and during this short period of time African indicine breeds may have not had enough time to differentiate from other indicine breeds, as much as the African taurines did. However, as the presence of indicine breeds in Africa is recent it is reasonable to expect that these carry more signatures of selection than African taurines (∼22 vs. ∼17%, respectively).

The analysis of selection carried out here sought to identify SNPs under selection using the equivalent of a replicated experiment. For that purpose it was necessary that several pairwise comparisons identified the same SNP under selection showing the same sign in the XP-EHH analysis (e.g., a SNP showing selection for the European taurine variant in the three pairwise comparisons). The rationale for this is that genetic drift can cause extreme differences in allele frequencies between populations, however, the chance that the same allele changes allele frequency in the same direction in multiple populations is low (e.g., A gets fixed by chance in the three replicates). While our approach is expected to be robust as it combines the power of the track of linkage disequilibrium left after positive selection and its replicability across populations occupying similar environments (e.g., three European breeds tested against three African breeds), it is possible that weaker signatures of selection were not identified. Nevertheless, it is interesting to note that for the pairwise comparisons where several genes were identified to be in close proximity to the SNPs under selection some chromosomes seems to be more frequent under selection than others. For example, among the Asian indicine, eight genes occurred on the first chromosome, and 11 on the 19th chromosome (SM Table 1), while three out of the four genes identified among the European taurine breeds occurred on the 10th chromosome. Lastly, it is expected that short term adaptive processes may be mediated through modifications of regulatory pathways (Karim et al., 2011; Visscher and Goddard, 2011; David et al., 2013; Furusawa and Kaneko, 2013). Consistently with such hypothesis we find that several SNPs under selection occur in close proximity to RNA coding sequences involved in transcription regulation (e.g., miRNAs) and modification of other RNAs (e.g., snoRNAs).

## References

The domestication of cattle in Eurasia had profound implications for society (Ajmone-Marsan et al., 2010; Groeneveld et al., 2010; Felius et al., 2014). Cattle are one of the major livestock species in the world and represent a large assemblage of widespread breeds as well as locally adapted breeds (Ajmone-Marsan et al., 2010; Felius et al., 2011, 2014). Studying the genetic resource of cattle breeds around the world is an on going work (e.g., Bovine Hapmap et al., 2009; Gautier et al., 2010; Mctavish et al., 2013; Decker et al., 2014) and approaches like the ones used here show the power that multiple comparisons have to derive particular trends and patterns that have characterized these breeds' histories. Additionally, identifying gaps representing regions where cattle has not yet been studied, as well as breeds that have not yet been genetically characterized is of utter importance to identify the genetic variation with adaptive value (e.g., local adaptation) and that should be the target of conservation efforts. Finally, whether Welsh White Park cattle were introduced to Britain by the Romans seems unlikely with the data described here, but cannot be entirely ruled out. The current results suggest that the breed presents a distinct genetic variation possible reflecting its ancient origin, but without some sort of calibration for the molecular phylogeny, an Approximate Bayesian computation approach that estimates times of divergence between breeds, or its comparison to local breeds close to the domestication center, this question is left partially unanswered (Csillery et al., 2010a,b; Sunnaker et al., 2013).

## Author Contributions

PO and MWB conceived the project; PO, MB, EN, and FB carried out analyses; WD, DW, MM, AS, PA contributed samples and reagents; PO wrote the manuscript; MB, EN, FB, and MWB revised the manuscript. All authors approved the final manuscript.

## Acknowledgments

The authors acknowledge the Chillingham Wild Cattle Association for the collection of samples for DNA analysis, as well as the NextGen consortium (http://nextgen.epfl.ch/) and Genomic-Resources (European Science Foundation; http:// www.esf.org/index.php?id=7009) for funding, and two reviewers for their comments.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fgene. 2015.00191/abstract

modern and ancient DNA. Proc. Natl. Acad. Sci. U.S.A. 103, 8113–8118. doi: 10.1073/pnas.0509210103

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Orozco-terWengel, Barbato, Nicolazzi, Biscarini, Milanesi, Davies, Williams, Stella, Ajmone-Marsan and Bruford. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Microsatellite genotyping of medieval cattle from central Italy suggests an old origin of Chianina and Romagnola cattle

#### *Maria Gargani 1 \*, Lorraine Pariset 1, Johannes A. Lenstra2, Elisabetta De Minicis 3, European Cattle Genetic Diversity Consortium† and Alessio Valentini 1*

*1 Department for Innovation in Biological, Agro-food and Forest systems, University of Tuscia, Viterbo, Italy*

*2 Faculty of Veterinary Medicine, Utrecht University, Utrecht, Netherlands*

*3 Department of Sciences of Cultural Heritage (DISBEC), University of Tuscia, Viterbo, Italy*

#### *Edited by:*

*Stéphane Joost, École Polytechnique Fédérale de Lausanne, Switzerland*

#### *Reviewed by:*

*Pablo Orozco ter Wengel, Cardiff University, UK Letizia Nicoloso, Università degli Studi di MIlano, Italy Florian Alberto, CNRS, France*

#### *\*Correspondence:*

*Maria Gargani, Department for Innovation in Biological, Agro-food and Forest systems, University of Tuscia, Via S. Camillo de Lellis, 01100 Viterbo, Italy e-mail<EMAIL_ADDRESS>
Analysis of DNA from archeological remains is a valuable tool to interpret the history of ancient animal populations. So far most studies of ancient DNA target mitochondrial DNA (mtDNA), which reveals maternal lineages, but only partially the relationships of current breeds and ancient populations. In this study we explore the feasibility of nuclear DNA analysis. DNA was extracted from 1000-years old cattle bone collected from Ferento, an archeological site in central Italy. Amplification of 15 microsatellite FAO-recommended markers with PCR products yielded genotypes for four markers. Expected heterozygosity was comparable with values of modern breeds, but observed heterozygosity was underestimated due to allelic loss. Genetic distances suggested a position intermediate between (1) Anatolian, Balkan, Sicilian and South-Italian cattle and (2) the Iberian, North-European and Central-European cattle, but also a clear relationship with two central-Italian breeds, Chianina and Romagnola. This suggests that these breeds are derived from medieval cattle living in the same area. Our results illustrate the potential of ancient DNA for reconstructing the history of local cattle husbandry.

**Keywords: microsatellite, ancient DNA, cattle, Ferento, NeighborNet**

## **INTRODUCTION**

The study of ancient DNA (aDNA) has developed since about 28 years from the analysis of short segments of mitochondrial DNA (mtDNA) to the spectacular whole-genome analysis of close relatives of the *Homo sapiens* (Higuchi et al., 1984; Pääbo, 1985; Reich et al., 2010). It complements archeological and historical investigations by providing ancestral links between ancient samples and present organisms (Vernesi et al., 2004; Caramelli et al., 2007). The majority of aDNA studies still target mtDNA, which is facilitated by its high copy number per cell. This allowed to determine the sequences of mtDNA control regions extracted from cattle remains excavated at several sites (Lenstra et al., 2014), but also the entire mitochondrial genome of two *Bos primigenius* specimens (Edwards et al., 2010; Lari et al., 2011). However, the genealogy inferred from the mtDNA not necessarily reflect that inferred with nuclear markers. Analysis of nuclear DNA from ancient remains is more challenging, but has been reported in several studies (Greenwood et al., 1999, 2001; Noonan et al., 2005; Pariset et al., 2007; Green et al., 2010; Reich et al., 2010). Ancient genotypes of nuclear microsatellites markers are informative for genetic diversity, subdivision and geographical origin (Edwards et al., 2003; Allentoft et al., 2011; Ishida et al., 2012; Nyström et al., 2012). Moreover, there are several studies carried out with the so called FAO list of microsatellites in almost every continent on a huge amount of individuals of several breeds (Ajmone-Marsan and The GLOBALDIV Consortium, 2010). This valuable database could be related to ancient remains in order i. to discover the closeness of ancient individuals to extant populations to understand the type of husbandry and production of the time and ii. to infer the origin, migration and admixture of present breeds.

Authentic cattle breeds in central and southern Italy are of the Podolian type (Felius, 1995). However, it is not yet clear if these are of ancient origin or have been imported (Felius et al., 2014). The aim of this study was to gain insight into history of Italian cattle by comparing 1000-year old cattle DNA

†The following members of the European Cattle Genetic Diversity Consortium contributed to this study: K. Moazami-Goudarzi, D. Laloë, INRA, Jouy-en-Josas, France; J. L. Williams (presently at the Parco Tecnologica Padano, Lodi, Italy), P. Wiener, Roslin Institute, UK; D. G. Bradley, Trinity College, Dublin, Ireland; G. Erhardt, Justus-Liebig Universität, Giessen, Germany; B. Harlizius, School of Veterinary Medicine, Hannover, Germany; I. Medugorac, Ludwig-Maximilians-Universität, Munich C. Looft, E. Kalm, Christian-Albrechts-Universität, Kiel, Germany; J. Cañón, S. Dunner, Universidad Complutense de Madrid, Spain; C. Rodellar, I. Martín-Burriel, Veterinary Faculty, Zaragoza, Spain; A. Sanchez, Universitat Autonoma de Barcelona; A. M. Martínez, Universidad de Córdoba, Córdoba, Spain; C. Ginja, Instituto Superior de Agronomia, Lisboa, Portugal; P. Ajmone-Marsan, Universitá Cattolica del S. Cuore, Piacenza, Italy; F. Pilla, A. Bruzzone, Universitá del Molise, Campobasso, Italy; D. Marletta, Università degli Studi di Catania, Catania, Italy; L. E. Holm, Danish Institute of Agricultural Sciences, Tjele, Denmark; G. Dolf, University of Berne, Switzerland.

of Italian origin with DNA from modern cattle. DNA was extracted from bones collected in Ferento, an archeological site near Viterbo in central Italy inhabited since Bronze Age, but developed mainly during Roman and Medieval ages. It rises on a triangular plateau made of tuff stone [Pianicara, IGM F.137 II NE] and delimited by two ditches: Guzzarella (or Vezzarella) and Acqua Rossa, subaffluents of the Tiber. The archeological investigations in Ferento revealed five areas, named assay I, assay II, assay III, assay IV, and assay V. The assay I comprises the area northward the Decumano in the medieval quarter; the assay II is localized on west area at the limit of the plateau, where an important necropolis was found; the assay III was occupied with a residential building on the southern side of the Decumano and with an important system of tanks, leant against the Theater; the assay IV, placed in the northwestern area of the plateau, and the assay V, near the assay I close to the fortification.

We compared the microsatellite genotypes obtained for four microsatellites with those of modern cattle breeds in order to uncover the origins of the Central-Italian medieval cattle.

## **MATERIALS AND METHODS**

#### **ANCIENT AND MODERN SAMPLES**

Thirty bone samples of medieval cattle (∼1000 BP) were collected at Ferento, an archeological site near Viterbo in central Italy (**Table 1**). The analyzed remains were recovered from five different excavation areas. The dating was made on the basis of pottery found in the same layer. Moreover, one samples were sent to CEDAD (CEntre for DAting and Diagnostics) belonging to the University of Salento (Department of Innovation Engineering) and were carbon-dated using high resolution mass spectrometry. The results confirmed the age estimated by pottery layers. Species identification, for the selection of the samples to be analyzed, was based on morphology and dimensions of the specimens.

In order to avoid contamination and degradation by sudden change in environmental conditions, samples were picked up using latex gloves, put immediately in vacuum-pack, stored at −20◦C (Pariset et al., 2007) and transferred to a lab dedicated to ancient DNA analysis. In a dedicated laboratory room, the outer surface was removed with sand paper and bone powder was collected by drilling. DNA was extracted according to Yang et al. (1998) and quantified using a DTX Multimode Detector

#### **Table 1 | Medieval samples.**

880 (Beckman) with Picogreen method (Quant-iT PicoGreen, Invitrogen). At least two independent DNA extractions were performed. One single sample per day was extracted with a negative control. Comprehensive datasets of modern samples of genotypes of cattle breeds from Europe, Africa and Asia (Table S1) have been described previously (European Cattle Genetic Diversity Consortium, 2006; Medugorac et al., 2009; Laloë et al., 2010; Felius et al., 2011; Delgado et al., 2012).

#### **MICROSATELLITE GENOTYPING**

Fifteen microsatellites were selected from the FAO list of recommended microsatellites (Table S2) (FAO, 2011). Only for three of these markers the published primers generated PCR products shorter than 160 bp. In order to optimize the amplification with aDNA templates, we redesigned primers for the other 12 microsatellites (Table S2) using Primer3 software. In a separate laboratory room, a first-round of PCR was carried out using a 20µl reaction volume containing PCR buffer with 2.5 mM MgCl2, 200µM of each dNTPs, 0.05µM of both the forward (fluorescent labeled dyes) and the reverse primer (Proligo), 0.2 units of Bio-x-act short Taq Polymerase (Bioline) and 10 ng of DNA. A 5 min denaturation step at 94◦C was followed by 30 cycles of 30 s denaturation at 94◦C, 1 min at the respective annealing temperatures (Table S2) and 1 min of extension at 68◦C, followed by a final 5 min extension step at 68◦C. Products were purified with exosap (Promega). A second round of PCR was performed with as template the product (1µl) of the first PCR reaction. During each PCR experiment both extraction and the PCR negative controls were checked for the absence of amplification products. At least four amplifications per sample were performed in order to monitor for allelic drop out. PCR products were separated by electrophoresis on a CEQ 8800 sequencer (Beckman Coulter) and sized with the standard CEQ™ DNA Size Standard Kit—400 (Beckman Coulter). The alleles were scored using the proprietary CEQ fragment analysis software.

Microsatellite dataset from modern cattle breeds were obtained from the EU Resgen project. In order to compare the allele lengths of the ancient cattle with those of the modern cattle database, we have included in all plates subject to PCR three modern samples (ITROM10, ITROM11, ITROM20) as a reference. The ancient cattle alleles were recalculated for comparison with the reference individuals. Genotypes were scored as heterozygous when the electropherogram showed a clear biallelic profile or when two different mono-allelic profiles were observed in different PCRs (Allentoft et al., 2011). The allelic profiles have been deposited in the Dryad database (http://dx*.*doi*.*org/10*.*5061/ dryad*.*d4500). Allele frequencies and observed and expected heterozygosity were calculated for each population using MSToolkit. Genetic distances were calculated by using the program MSAT2 (http://genetics*.*stanford*.*edu/hpgl/projects/microsat/) and visualized in NeighborNet graphs with SplitsTree4 program (http:// www*.*splitstree*.*org/). The NeighborNet graphs were simplified by combining several breeds in regional pools as approximations of ancestral gene populations. Thus, we made metapopulations of the dairy breeds from the Northwest-European Lowland, Central Brown cattle (including Swiss-Brown, Rendena, Bruna, and Cabannina), a West-Central cluster (Simmental and its derivatives); Iberian cattle, Sicilian Modicana and Cinesara, Balkan Busha, Anatolian cattle and Indopakistan zebu breeds, respectively. Principal component analysis of a matrix of F*ST* genetic distances was performed using the program Genalex (Peakall and Smouse, 2006). Unsupervised model-based clustering was done with the program Structure (Pritchard et al., 2000).

## **RESULTS**

Extraction of genomic DNA yielded for all 30 ancient samples quantifiable DNA and for 16 of these also amplification products. From the 15 microsatellites, four (INRA5, HEL1, MM12, CSRM60) could be amplified successfully in at least one sample. A total of 760 amplifications were performed on ancient samples. The success rate PCR per sample and per locus is shown in Table S3. Markers INRA5, HEL1, MM12, and CSRM60 yielded genotypes in 6, 8, 10, and 10 samples, respectively. All extraction and PCR controls for different DNA isolation, purification and amplification steps did not yield PCR products. In 35 PCR reactions we observed allelic dropout of 50% of heterozygous genotypes for HEL1, 60% for CSRM60, 75% for MM12, 73% for INRA5 and 64% averaged of the four loci, but never more than two alleles per marker in a given sample. The expected (He) and observed (Ho) heterozygosity values for a panel of European breeds calculated on the basis of 4 microsatellites ranged from 0.48 (Lagunaire breed) to 0.771 (Istrian breed) and from 0.43 (Lagunaire breed) to 0.76 (Istrian breed) respectively. The He and Ho values for pools of breeds ranged from 0.68 (Zebu breed pool) to 0.786 (Sicilian breed pool) and from 0.62 (Zebu breed pool) to 0.79 (Anatolian breed pool) respectively. These values were in the same range of those based on 30 microsatellites: He and Ho for panel of European breeds ranged from 0.544 (Lagunaire breed) to 0.741(Podolica breed) and from 0.51(Lagunaire breed) to 0.73 (Bohemian Red breed); He and Ho for pools of breed ranged from 0.67 (Zebu breed pool) and 0.780 (Anatolian breed pool) and from 0.59 (Zebu breed pool) and 0.72 (Busha breed pool) respectively (Table S1). The expected heterozygosity (0.698) in Ferento samples is within the range of the values found for modern breeds. We observed a clear heterozygote deficit (observed heterozygozity 0.588) related to allelic dropout for ancient samples.

Principal component analysis (Figure S1) as well as modelbased clustering of Ferento cattle (not shown) combined with a panel of European breeds typed with the same markers differentiated three clusters: (1) African cattle, (2) Northwestern-European Lowland dairy breeds, and (3) Central-European and Iberian cattle as well as the ancient samples from Ferento. However, PCA did not establish relationships between Ferento cattle and individual breeds, since relative positions of breeds could not be reproduced with other panels of four microsatellites (not shown). We are aware that four microsatellites may have a low power in resolving diversity in populations that diverged since very few generations, however, we could not increase the number of markers that successfully and repeatedly amplified ancient DNA within the FAO list that permits to relate ancient to extant breeds.

In order to get a more informative comparison of ancient Ferento genotypes with those of other cattle, we constructed NeighborNet graphs. **Figure 1** shows NeighborNet graphs of the clusters with several modern Podolian breeds and ancient

Ferento cattle. Although distances on the basis of only four microsatellites are not expected to be accurate, these distances reproduce the clustering that are generated by the 30 microsatellites. Interestingly, Ferento cattle is clearly linked with the Tuscanian Chianina and the related Romagnola cattle. The ancient and modern central Italian cattle is intermediate between (1) the Anatolian, Sicilian and Balkan Podolian breeds, and (2) the Iberian, Central European and Lowland dairy breeds. We observed that none of the individuals from modern breeds were as close to Ferento as the Chianina and Romagnola breeds.

## **DISCUSSION**