Signatures of selection in five Italian cattle breeds detected by a 54K SNP panel

In this study we used a medium density panel of SNP markers to perform population genetic analysis in five Italian cattle breeds. The BovineSNP50 BeadChip was used to genotype a total of 2,935 bulls of Piedmontese, Marchigiana, Italian Holstein, Italian Brown and Italian Pezzata Rossa breeds. To determine a genome-wide pattern of positive selection we mapped the Fst values against genome location. The highest Fst peaks were obtained on BTA6 and BTA13 where some candidate genes are located. We identified selection signatures peculiar of each breed which suggest selection for genes involved in milk or meat traits. The genetic structure was investigated by using a multidimensional scaling of the genetic distance matrix and a Bayesian approach implemented in the STRUCTURE software. The genotyping data showed a clear partitioning of the cattle genetic diversity into distinct breeds if a number of clusters equal to the number of populations were given. Assuming a lower number of clusters beef breeds group together. Both methods showed all five breeds separated in well defined clusters and the Bayesian approach assigned individuals to the breed of origin. The work is of interest not only because it enriches the knowledge on the process of evolution but also because the results generated could have implications for selective breeding programs.


Introduction
Present day cattle breeds are the result of years of human selection, adaptation to different environments and demographic effects as domestication, migration and selection, all contributing to the actual patterns of genetic diversity [1,2]. During the domestication process, breeds were selected for productivity traits as, for example, milk yield [3]. Moreover, animal and semen exchange, carried out to improve production characteristics, have affected the genetic features of the breeds. This anthropic selection has influenced the genetic structure of cattle breeds, therefore a high percentage of loci purposely chosen for influencing potentially selected traits could result under selection [4,5].
Recently, the availability of high density SNP panels has given the possibility of performing population genetic studies in cattle populations using thousands of markers distributed across the entire genome. Medium density SNP panels have been used for example to analyze the genetic structure of cattle populations [6][7][8] to study past effective population size [9], to detect selection signatures [10], and to discover copy number variation (CNV) suitable for understanding genetic features and accelerating genetic improvement for complex traits [11].
We considered a total of 2,935 bulls belonging to two dairy (Italian Brown and Italian Holstein), two beef (Piedmontese and Marchigiana) and one double purpose (Italian Pezzata Rossa) breeds. The Italian Holstein derives from Dutch and North America Holstein breeds imported in Italy in the late XX century and it is currently the most common dairy breed. Piedmontese is mainly located in Northern Italy and it was in the past a dual purpose breed, while today it is selected for beef traits mainly exploiting a private myostatin mutation [12]. Marchigiana is a beef breed from central Italy derived from very ancient breeds like Chianina and Romagnola breeds. The Italian Brown was originally a multi purpose breed reared in the Alps; it was selected from 1950 as a dairy breed by importing Swiss Brown bulls from the U.S. Pezzata Rossa, Simmental, is a beef/dairy breed imported from Swiss/Austria and herded mostly in North East Italy.
The aim of this study was to identify genomic regions potentially under selection in the above five Italian cattle breeds using a 54K medium-density SNP panel. Our results could have implications for selective breeding programs by identifying signatures of artificial selection in gene involved in milk, meat or functional traits. We analyzed also the genetic structure of the breeds by classical multidimensional scaling and by Bayesian inference methods.

Samples and high throughput genotyping
The initial sample was formed by 2,935 bulls: 761 Italian Brown, 899 Italian Holstein, 323 Piedmontese, 464 Marchigiana and 488 Italian Pezzata Rossa bulls, respectively. The sample represents almost all the bulls available in Italy for all breeds but Holstein, where the bulls analyzed correspond to slightly less than a half of the available ones. Genomic DNA was extracted from semen using the NucleoSpin Tissue kit (Macherey-Nagel) according to manufacturer's instruction. DNA was checked for quality on agarose gel and quantified using a DTX microplate reader (Beckman Coulter) after staining with Picogreen (Invitrogen). Samples were genotyped using BovineSNP50 Genotyping BeadChips (Illumina, San Diego, CA, USA). Genotyping was outsourced to Geneseek (www.geneseek.com). The 50K SNP array contains 54,001 SNPs distributed across the entire genome, with an average SNP spacing of 51 Kb and a proportion of known chromosome positions of about 97 %; SNP positions within each chromosome were based on the Bos taurus genome assembly Btau_4.0 [13].
Data editing and genome-wide analysis Data were initially filtered using the GenABEL R package (http://www.r-project.org, http://mga.bionet.nsc.ru/*yurii/ ABEL/GenABEL/). Only autosomal markers were used and SNP with complete map information were used. Sires and markers with a call rate under 99 % were discarded as well as SNPs having a minor allele frequency (MAF)\5 % according to the currently employed thresholds [14]. Sires were checked for abnormally high autosomal heterozygosity and discarded when showing a false discovery rate (FDR)\1 % [15]. Then, sires of each breed were separated and Hardy-Weinberg equilibrium (HWE) was checked within each breed setting a threshold of P \ 0.01 in the filtered data set [16]. Finally, data were pooled again and filtering criteria explained above were applied once more. Kinship among sires was estimated directly from genomic data as proposed by Astle and Balding [17].
To determine a genome-wide pattern of positive selection, the F ST at each locus was calculated [18]. The loci under selection are expected to show an allele frequency that deviate from that of neutral loci, leading to an increased level of genetic differentiation. F ST values were then plotted against genome location. Signatures of selection can be recognized when adjacent SNPs in a region show high F ST [19] thus we used a sliding window approach, with a window of eight SNPs. A region with high F ST implies divergent selection between breeds, whereas low F ST imply balancing selection between breeds. Fixation index was calculated using the method proposed by Nei and Chesser [20] using in-house written R codes; this method was chosen because the sample includes (almost) all sires in the national herdbooks for three breeds (Marchigiana, Piedmontese and Pezzata Rossa) and thus fixed effects errors of sampling, i.e. effects unbounded by a prior distribution, seemed more important. Graphs were obtained using matplotlib (http://matplotlib.sourceforge. net/).
Genetic distance between breeds and sires was estimated calculating the matrix d ij = (0.5k ij ) where d is the distance and k is genetic kinship coefficient for sires i and j and then applying classical multidimensional scaling to the complete matrix. STRUCTURE software v. 2.3.2.1 [21] was used to analyse population structure. A total of 15,000 Markov chain Monte Carlo (MCMC) iterations (5,000 burn-in and 10,000 sampling) were performed for each tested K using the admixture model, considering allele frequencies correlated among populations and including no informative prior about individual membership; K values from 2 to 5 were used. Five independent runs for each tested K value were performed. The number of steps was chosen following [22], although for each K a single run of 50,000 iterations to test the effects of longer runs was performed. Evanno et al. [23] reported that in most cases, the estimated 'log probability of data' did not provide a correct estimation of cluster number (K value), and argued that an ad hoc statistic DK based on the rate of change in the log probability of data between successive K values could accurately detect true K. The suggested statistics was: where L(K) represents the Kth LnP(D), m is to the mean of 10 runs and s their standard deviation. We used the method of Evanno et al. [23] to estimate the number of populations. Graphical visualization of STRUCTURE results was performed by means of the DISTRUCT package [24].  (Figs. 1b, 2b). Three of the four SNPs identified in this region are located in the platelet-derived growth factor receptor, alpha polypeptide (PDGFRA) gene involved in the reproduction process and in the regulation of calcium level and near the KIT (gene v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog) gene expressed in the lactating bovine mammary gland and implicated in determining coat colour [5]. In dairy breeds a weaker peak, formed by four markers, is also observable around 38 Mb on BTA6 where ATP-binding cassette, sub-family G (ABCG2) and polycystic kidney disease 2 (PKD2) genes are located. The two genes play a role in the regulation of bovine lactation [25] and in calcium homeostasis, respectively [26].

Data editing
In beef breeds the peak on BTA6 consists of six SNPs spanning from 37.32 to 38.76 Mb (Figs. 1c, 2c). This interval contains 15 genes including LAP3 (leucine aminopeptidase 3) non-SMC condensin I complex subunit G (NCAPG) and ligand-dependent nuclear receptor corepressor-like protein (LCORL), genes involved in calving ease [27]. Figure 2a shows F ST values for markers on BTA6 only, without any averaging. Several markers that have F ST above the 99 % quantile can be observed at 18 Mb, between 36 and 39 Mb and at 95 Mb, but the strongest signal (considering either the number of SNPs with high F ST or the maximum value) can be observed at *38 Mb.
Moreover we calculated the F ST values for BTA6 comparing each breed against the remaining four (Fig. 3). In Italian Brown breed we can observe a signals at 72 and 37 Mb, in Italian Holstein and Pezzata Rossa breeds the strongest signal can be observed at 72 Mb, while in Marchigiana a peak is located around 37 Mb. The Piedmontese breed showed two signals one around 18 Mb and the other around 68 Mb.
The F ST calculated only for dairy breeds (Fig. 1b) shows a peak formed by four SNPs spanning from 47.23 to 48.30 Mb on BTA13. In this region we identified 14 genes among which CDS2 (CDP-diacylglycerol synthase) could affect milk fat composition. Four significant SNPs on chromosome 13 at positions 67197635, 67464116, 67490718, 67766784 were identified in beef breeds (Fig. 1c). Of the 8 genes in this region, SRC (v-src sarcoma viral oncogene homolog avian) and CTNNBL1 (catenin, beta like 1) may be related to muscle formation and to body weight, respectively.
Multidimensional scaling Figure 4 shows the first three components of the multidimensional scaling decomposition of the genetic distance  Fig. S1, Marchigiana and Piedmontese in Fig. S2) are located in the C1-C2 plane while the remaining population (Holstein in Fig. S1 and Pezzata Rossa in Fig. S2) is located in a different quadrant in the C1-C2 plane and is highly dispersed along C3.

Bayesian inference
To estimate the number of genetic clusters within the 2,797 cattle samples and 29,848 SNPs, a parametric genetic mixture analysis implemented in the STRUCTURE software was performed. Between 2 and 5 clusters (K values) were tested using the admixture model, considering allele frequencies correlated. Consistent results across runs were obtained and a clear clustering of breeds was observed for any K tested (Fig. 5). With K = 2 Brown and Holstein individuals are assigned to different clusters while sires of the other three breeds belong to each cluster with probabilities near 50 %. The separation becomes sharper when

Discussion
Recently many innovative tools, such as medium or high density SNP chips, have been developed for various domesticated species. In this study we presented two applications of population genetic analysis in five Italian cattle breeds using 50K bovine SNP chips. We first investigated the potential of SNP markers to identify selection signatures peculiar of each breed and then we analyzed the genetic structure in the same samples. By mapping the F ST values against genome location we identified genes showing signatures of positive selection involved in biological processes such as reproduction, metabolism of lipids, calving ease. A strong selection signal was observed on BTA6 when considering F ST across all cattle breeds. Interestingly, the F ST calculated only in dairy breeds revealed the evidence for selection in the region located at 72.45 Mb on chromosome 6, far from the caseins cluster, for which selection is carried out in some breeds. The peak is located near the PDGFRA gene which is associated with b-estradiol and implicated in the reproduction process and in the regulation of calcium level and near the KIT gene expressed in the lactating bovine mammary gland and implicated also in determining coat colour. These results are consistent with those obtained by Flori et al. [28] and Stella et al. [5] who found a positive selection signature in the same region in dairy cattle breeds. Flori et al. [28] used the F ST approach to detect the selection signatures in three French dairy cattle breeds and highlighted 13 significant signatures including the PDG-FRA gene which is proposed as candidate gene. Stella et al. [5] reported the largest composite log likelihood (CLL) in the same location on BTA6 within the KIT gene which is responsible for the piebald phenotype in four of the five dairy breeds analyzed. The peak around 38 Mb falls near the ABCG2 and PKD2 genes. Several studies identified a QTL affecting milk yield and milk composition on chromosome 6 in a region around 38 Mb containing ABCG2 gene [25,29,30]  Recently the results obtained by Wei et al. [25] suggested that ABCG2 plays a role in mammary epithelial cell proliferation and that the polymorphisms in this gene may influence milk production.
The other interesting gene in the region is PKD2 gene that could be related with the content of water in the milk since it is involved in calcium homeostasis [26].
The peak located at *38 Mb on BTA6 in beef breeds is near LAP3, NCAPG and LCORL. LAP3 encodes for a leucine aminopeptidase, which is responsible of the oxytocin hydrolysis [31]. Recent studies demonstrated the role of LAP3 in calving ease in Norwegian Red cattle [32]. Moreover, Bongiorni et al. [27] found a strong association between LAP3, NCAPG, LCORL and calving ease trait in Piedmontese.
These results are in agreement with the F ST values calculated for each breed against the other in BTA6.
The signal in dairy breeds is due mainly by Bruna, Holstein and Pezzata Rossa, while the signal around 37 Mb identified in beef breeds is due mainly to Marchigiana breed (Fig. 3). It is worth to notice that the Piedmontese breed shows a peak spanning between 17 and 18 Mb on BTA6. This region contains five genes including COL25A1 collagen gene.
A cluster of signals reflecting strong evidence of selection was observed also in BTA13. When we analyzed separately the F ST for beef breeds (Fig. 1c) a strong peak at position 67 Mb could be observed. Two interesting genes are located in this region: SRC and CTNNBL1. SRC is involved in the regulation of actin cytoskeleton and in the focal adhesion pathway [33]. Some studies reported the role of focal adhesion pathway for muscle formation in cattle [34] and muscle strength and integrity in racing horses [35]. The other gene, CTNNBL1, is associated with body weight and height [36]. In human it has been showed to be involved in the W nt /beta-catenin-signaling pathway and associated with obesity [36]. In dairy breeds we observed a peak spanning 47-48 Mb on BTA13 near the CDS2 gene which could influence milk fat composition; the gene is in fact involved in the phospolipid biosynthetic process. The Table 1 shows a list of candidate genes for genomic regions presenting the most extreme peaks in dairy and beef breeds.
No F ST peaks have been detected on or near the few genes today known to influence dairy or beef traits, like DGAT1, caseins, myostatin, leptin. This may be due to the large genetic network that influences the complex traits under selection, as well as the changes of the selection policies cross time: i.e. at least in Italy and many European countries, in the early days, the main selection objective was milk yield, afterwards it was protein and fat percentage, now sustainability traits are included in the selection index.
Regarding the assessment of the genetic structure we used Bayesian and multidimensional scaling approaches. Multidimensional scaling separated each breed in five well  defined clusters. Piedmontese formed the most compact cluster indicating that the breeding policy in this breed tends towards a narrower genetic basis. It has to be noted that in the past 50 years the selection has been oriented to strongly select for the double muscling trait and culling all non carrier subjects. It is worth to observe that the double purpose Pezzata Rossa is more distant from dairy breeds than beef breeds, suggesting in this breed a different management of selection and lack of admixture with other dairy breeds. Interesting, the Italian Brown showed a small group of outlying bulls in the MDS plot suggesting a potential substructure maybe due to the double type of exploitation of this breed both in high producing farms in the valleys and in harsher conditions in the mountains. The genetic isolation and lack of admixture among the two dairy breeds are confirmed by Bayesian analysis in which the breeds do not cluster on the basis of their purpose even for K = 2 (Fig. 5). This means that the differentiation pre-dates the selection for different purposes. Beef and dual purpose breeds tend to cluster together up to K = 3 and at K = 4 Piedmontese is assigned almost equally to the other two beef breeds. We hypothesize in this case a possible convergent artificial selection for beef breeds.