Open Access

Association analysis for udder health based on SNP-panel and sequence data in Danish Holsteins

  • Xiaoping Wu1, 2,
  • Mogens S Lund1,
  • Goutam Sahana1,
  • Bernt Guldbrandtsen1,
  • Dongxiao Sun2,
  • Qin Zhang2 and
  • Guosheng Su1Email author
Genetics Selection Evolution201547:50

Received: 11 September 2014

Accepted: 21 May 2015

Published: 19 June 2015



The sensitivity of genome-wide association studies for the detection of quantitative trait loci (QTL) depends on the density of markers examined and the statistical models used. This study compares the performance of three marker densities to refine six previously detected QTL regions for mastitis traits: 54 k markers of a medium-density SNP (single nucleotide polymorphism) chip (MD), imputed 777 k markers of a high-density SNP chip (HD), and imputed whole-genome sequencing data (SEQ). Each dataset contained data for 4496 Danish Holstein cattle. Comparisons were performed using a linear mixed model (LM) and a Bayesian variable selection model (BVS).


After quality control, 587, 7825, and 78 856 SNPs in the six targeted regions remained for MD, HD, and SEQ data, respectively. In general, the association patterns between SNPs and traits were similar for the three marker densities when tested using the same statistical model. With the LM model, 120 (MD), 967 (HD), and 7209 (SEQ) SNPs were significantly associated with mastitis, whereas with the BVS model, 43 (MD), 131 (HD), and 1052 (SEQ) significant SNPs (Bayes factor > 3.2) were observed. A total of 26 (MD), 75 (HD), and 465 (SEQ) significant SNPs were identified by both models. In addition, one, 16, and 33 QTL peaks for MD, HD, and SEQ data were detected according to the QTL intensity profile of SNP bins by post-analysis of the BVS model.


The power to detect significant associations increased with increasing marker density. The BVS model resulted in clearer boundaries between linked QTL than the LM model. Using SEQ data, the six targeted regions were refined to 33 candidate QTL regions for udder health. The comparison between these candidate QTL regions and known genes suggested that NPFFR2, SLC4A4, DCK, LIFR, and EDN3 may be considered as candidate genes for mastitis susceptibility.


Mastitis, i.e. inflammation of the mammary gland, is a common and costly disease [1, 2] that is particularly problematic in the dairy industry. It adversely affects both animal and human health, since milk from affected cattle can enter the food supply and pose a health risk [3]. Mastitis also leads to economic losses owing to reduction in milk yield and quality, increased costs associated with treatment of affected cattle, discarded milk, and culling of severely affected cattle [4]. Therefore, prevention of mastitis is an extremely important breeding goal for dairy cattle. However, the heritability of clinical mastitis is low [59], and its genetic correlation with production traits is unfavorable [6, 10, 11]. Identification of specific genes associated with clinical mastitis would be extremely useful in breeding programs to reduce mastitis incidence. To that end, quantitative trait loci (QTL) mapping is a useful tool to identify genomic regions that are associated with clinical mastitis. Genome-wide association studies (GWAS) have been successfully used to identify QTL regions for a variety of traits in livestock [12]. Several studies have also attempted to detect QTL that affect clinical mastitis, with varying success [1316].

Many factors influence the efficiency of GWAS to detect QTL. One of the most important factors is marker density. Currently, a medium-density (MD) SNP chip with ~54 000 markers is widely used for GWAS in dairy cattle [1719]. In recent years, two high-density (HD) SNP chips with 777 962 SNPs from Illumina Inc. [20] and 648 874 SNPs from Affymetrix Inc. [21], and whole-genome sequencing (SEQ) data [22] have become available. Higher marker densities mean that markers are in stronger linkage disequilibrium (LD) with QTL that affect the trait of interest. In addition, genome sequencing data includes causative genetic variants. Thus, it is expected that using HD or SEQ markers will lead to the detection of more and more accurately localized QTL.

In addition to marker density, the statistical model selected for analysis has an effect on detection sensitivity. Several models have been proposed for GWAS, such as single-marker tests (SMT) [23], mixed model analysis (MMA) and linear models (LM) [24, 23], haplotype models (HM) [25], genealogy-based mixed-model (GENMIX) [25], and Bayesian variable selection models (BVS) [26]. Some studies have carried out model comparisons using simulated data. For example, Sahana et al. [23] compared the SMT, HM, LM, and BVS models, and concluded that the BVS model performed best. Dashab et al. [25] compared LM, random HM, GENMIX, and BVS models, and showed that LM and BVS were better than the other methods in terms of power and type-I error rate. However, there are very few reports on model comparisons based on real data from livestock [27, 28].

In general, LM models are performed in single-marker test analyses, for which each SNP is fitted separately in the model. In contrast, BVS models estimate the effects of all SNPs simultaneously. Usually, BVS models are implemented via Markov chain Monte Carlo (MCMC) algorithms [2931]. For both models, false positive associations due to population structure, such as family relationships can be controlled [23, 24] by including systematic factors and taking polygenic effects into account. For LM models, significant associations are commonly established using a t-test with Bonferroni correction, assessment of false discovery rate, or permutation testing [32]. For BVS models, the number of QTL can be determined by post-MCMC analysis using Bayes factors [33] or QTL intensity profiles [29]. Previous studies have claimed that BVS models yield higher power than LM and are a solution to the problem of establishing significance of multiple-testing in simulation studies [25, 23]. Sahana et al. [23] compared the marginal posterior probability for single markers (BAYSM) and the joint posterior probabilities for intervals of 11 markers (BAYINT) to infer the presence of QTL using a BVS model based on simulated data. They reported that BAYINT resulted in higher power to detect QTL, while BAYSM was more precise in estimating QTL position.

Complex traits are likely influenced by multiple QTL that have small individual effects but a large collective effect on the phenotype. Conversely, the effects of a given QTL can be distributed over several markers that are in LD with the QTL by using a Bayesian model which estimates effects of all SNPs simultaneously. The QTL intensity profile, which summarizes the marginal posterior probabilities of markers in a small region, was proposed to detect QTL regions [34]. In this situation, a small QTL can be detected via a marked peak, and multiple QTL can be separated [34]. QTL intensity profiles based on post-MCMC analysis provide a suitable approach to identify QTL regions using real data.

The first objective of this study was to refine six previously detected QTL regions for udder health in Danish Holstein cattle using three marker densities, i.e. MD, HD and SEQ datasets. The second objective was to compare the performance of the LM and BVS models to detect and separate QTL that are closely located.


Phenotype and genotype data

A total of 4496 Danish Holstein bulls with de-regressed estimated breeding values (DRP) for the udder health index were used. The index for udder health is a measure of the genetic value of the cow’s resistance to mastitis. The index was calculated based on estimated breeding values (EBV) of clinical mastitis (CM) from the 1st to 3rd lactations, i.e.:
$$ \mathrm{Index}\ \mathrm{of}\ \mathrm{udder}\ \mathrm{health} = 0.25*{\mathrm{CM}}_{11} + 0.25*{\mathrm{CM}}_{12} + 0.30*{\mathrm{CM}}_2 + 0.20*\mathrm{C}\mathrm{M}3, $$

where CM11, CM12, CM2, CM3 are EBV for clinical mastitis in lactation 1 from days in milk (DIM) -15 (15 days before calving) to 50, lactation 1 from 51 to 305 DIM, lactation 2 from −15 to 150 DIM, and lactation 3 from −15 to 150 DIM, respectively. Cows with and without clinical mastitis were recorded as 1 and 0, respectively. Breeding values of clinical mastitis for each lactation were estimated using a multi-trait random regression test-day animal model. Somatic cell counts for the 1st, 2nd and 3rd lactation, fore udder attachment, and udder depth from the 1st lactation were included as correlated traits to improve the accuracy of EBV for clinical mastitis [35], but EBV of these correlated traits were not included in the index of udder health. Udder health DRP [36, 37] were derived from the index of udder health evaluated in November 2010.

An association study for udder health was carried out using three marker datasets: (1) Illumina BovineSNP50 BeadChip (MD), (2) Illumina BovineHD BeadChip (HD), and (3) genome sequence data (SEQ). Six chromosome regions with significant effects on clinical mastitis that were reported in our previous study [38], were used for further association analysis. These six regions are located between 84 and 95 Mb on BTA5 (BTA for Bos taurus chromosome), 88 and 96 Mb on BTA6, 57 and 63 Mb on BTA13, 48 and 55 Mb on BTA16, 55 and 58 Mb on BTA19, and 32 and 40 Mb on BTA20. SNP array genotyping, sequence analysis, and imputation were performed as described by Sahana et al. [38]. Briefly, 5214 bulls were genotyped with the 50 k SNP chip. SNP chip data were edited using a number of criteria. First, individuals with a call rate higher than 85 % and SNPs with a call rate higher than 95 % were kept. SNPs with a minor allele frequency below 5 % and those that deviated from the Hardy-Weinberg equilibrium (P < 0.00001), or with average GenCall scores (Illumina Inc.) below 0.65 were excluded. After this quality control filtering, 43 415 SNPs remained. The MD data of 5214 Holstein bulls were imputed to HD data with BEAGLE [39], using 557 HD genotyped bulls (including 138 Nordic Holstein bulls, which represent a substantial contribution of the population) in the EuroGenomics project [40] as reference data. After quality control (using criteria similar to those described above for the 50 k chip), 648 219 HD SNPs remained. Finally, six regions with 790 SNPs from the MD chip and 8260 SNPs from the HD chip were targeted for analysis. Linkage disequilibrium (average r2 between adjacent marker pairs) estimates were equal to 0.21 for the MD and 0.56 for the HD data in the six targeted regions.

Sequence data analysis was performed at the Center for Quantitative Genetics and Genomics at Aarhus University. Details on sequencing, sequence analysis, and quality control are in [41, 42, 38]. Using the 89 sequenced bulls (including 32 Nordic Holstein bulls, which represent a substantial contribution of the present population) as reference data, 5214 Holstein bulls with imputed HD data for six targeted regions were imputed to SEQ data by BEAGLE [39]. A total of 196 882 SNPs were located within the six regions for the SEQ data. In the end, using all available phenotypic and genotypic data, 4496 Holstein bulls with MD, imputed HD, and imputed SEQ data were used in the analyses. These animals were the progeny of 373 bulls with 3169 dams.

The SNPs in each of the six regions were checked for quality as follows. First, SNPs that exclusively appeared in lower density data but disappeared in higher density data were discarded, removing 201 SNPs that were in MD but not in the HD and SEQ data as well as 308 SNPs that were in the HD but not in the SEQ data. Second, SNPs with a minor allele frequency less than 0.01 were discarded, removing one SNP from the HD dataset and 41 304 SNPs from the SEQ dataset. Finally, if two markers were in complete LD, then either the SNP that was absent from the lower marker density dataset was deleted or the more downstream SNP was deleted if both markers were absent from the lower density data. This removed two SNPs in the MD, 126 SNPs in the HD, and 76 722 SNPs in the SEQ data. After these edits, 587, 7825, and 78 856 SNPs in the six regions remained in the MD, HD, and SEQ data, respectively. The average distances between adjacent SNPs were approximately 70.0, 5.5, and 0.5 kb for the MD, HD, and SEQ data, respectively.

Statistical analyses

A linear mixed model and a Bayesian variable selection model were used to detect associations of SNPs with phenotypic data, as described below.

Linear mixed model (LM)

The LM model [24] used in this study is a single-SNP regression model. The model includes a fixed regression of phenotypes using SNP genotypes as a measure of the SNP effects. In addition, a random polygenic effect that accounts for shared genetic effects of related individuals was included. The statistical model is described by the formula:
$$ \mathbf{y}=\mathbf{1}\mu + \mathbf{x}\mathrm{g} + \mathbf{Z}\mathbf{u}+\mathbf{e}, $$

where y is the vector of phenotypes (DRP), 1 is a vector of ones, μ is the overall mean, g is the additive genetic effect of the analyzed SNP, x is a vector of SNP genotypes (coded as 0, 1 or 2), and u is a vector of random polygenic effects, which are normally distributed u ~ N(0, A σ u 2 ), where A is the pedigree-based additive relationship matrix, σ u 2 is the polygenic variance, Z is an incidence matrix relating phenotypes to the corresponding random polygenic effects, and e is a vector of residual effects, which are normally distributed e ~ N(0, Dσ e 2 ), where D is a diagonal matrix with elements d ii  = (1 − r DRP 2 )/r DRP 2 to account for heterogeneous residual variances due to different reliabilities of DRP (r DRP 2 ), and σ e 2 is the residual variance.

Significance testing of SNP effects was performed using a two-sided t-test. Our aim was to refine QTL regions from the selected large regions using the same population that was used to identify these regions. Although the segregation of QTL in the targeted regions was known, regions were large (up to 11 Mb). Therefore, we assumed that, except for a few, most of the SNPs were not associated with the trait. Accordingly, our null hypothesis was g = 0. A Bonferroni correction was applied to control for false positive associations. We declared a SNP significant if the P value was less than 0.05/N, where N is the number of SNPs. Therefore, the significant threshold value for − log10(P) was 4.07, 5.19, and 6.20 for the MD, HD, and SEQ data, respectively. Analyses were performed using the DMU package [43].

Bayesian variable selection model (BVS)

The BVS model [26, 44] used in this study describes SNP effects as a mixture distribution that estimates the effects of all SNPs simultaneously. The model is given by the following formula:
$$ \mathbf{y} = \mathbf{1}\upmu + {\displaystyle {\mathbf{\sum}}_{k=1}^m}{\mathbf{x}}_k{g}_k + \mathbf{Z}\mathbf{u} + \mathbf{e} $$
where y, 1, μ, u, Z, and e are defined as in the LM model. The term \( {\displaystyle {\mathbf{\sum}}_{k=1}^m}{\mathbf{x}}_k{g}_k \) fits additive SNP association effects, x k is a vector of SNP genotypes (0, 1 or 2), and g k is the effect of SNP k, which was assumed to follow a mixture distribution that assumes that most SNPs have a small phenotypic effect and some SNPs have a large effect:
$$ {\mathrm{g}}_k\sim \left\{\begin{array}{c}\hfill \mathrm{N}\left(0,\ {\sigma}_{g_0}^2\right)\ \mathrm{with}\ \mathrm{probability}\ {\pi}_0\hfill \\ {}\hfill \mathrm{N}\left(0,\ {\sigma}_{g_1}^2\right)\ \mathrm{with}\ \mathrm{probability}\ {\pi}_1=1-{\pi}_0\hfill \end{array}\right., $$

where π 0 is the probability of a SNP having a small effect (i.e., drawn from a distribution with a small variance) and π 1 is the probability of a SNP having a large effect (i.e., drawn from a distribution with a large variance). We assumed that the proportions in the two mixture distributions had a Bernoulli distribution and the prior distribution of π 0 and π 1 is a Beta (100, 1) distribution. We set π 0 = 0.98, 0.998, and 0.9998 as priors for the MD, HD, and SEQ data, respectively. We also assumed that the priors μ, \( {\sigma}_{g_0}^2 \), and \( {\sigma}_{g_1}^2 \) followed a uniform distribution. By assuming a small variance instead of zero for the first distribution (\( \mathrm{N}\left(0,\ {\sigma}_{g_0}^2\right) \)), the MCMC implementation is straightforward with recognizable conditional distributions for all model parameters, as described elsewhere [26, 44]. Analysis of the BVS model was performed using the Bayz software [45]. The Gibbs sampler was run as a single chain with 52 000 samples. The first 20 000 iterations were discarded as burn-in. After this, every 20th sample of the remaining 32 000 was saved for posterior analysis.

Within the Bayesian framework, the association of each SNP with the trait can be evaluated by a Bayes factor (BF). In this study, the BF for SNP i was calculated as the ratio between the estimated posterior probability and the average of posterior probabilities of SNPs with large effects (large variance, \( {\mathrm{g}}_k\sim \mathrm{N}\left(0,\ {\sigma}_{g_1}^2\right) \)) [46]:
$$ \mathrm{B}{\mathrm{F}}_i = \frac{{\hat{p}}_i/\left(1-{\hat{p}}_i\right)}{{\hat{\pi}}_1/\left(1-{\hat{\pi}}_1\right)}, $$

where \( {\hat{p}}_i \) is the posterior probability of the effect for SNP i being drawn from the distribution with large effects, and \( {\hat{\pi}}_1 \) is the average of posteriors probabilities of SNPs with large effects. According to commonly used criteria [33], a BF between 3.2 and 10 was considered moderate and a BF greater than 10 as strong evidence for the segregation of a QTL affecting a trait.

QTL regions

A QTL region was detected using the posterior QTL intensity function [29]. In a Bayesian analysis, each chromosome is divided into many small intervals of equal length (bins), i.e. for instance, 1 or 2 cM in simulation studies [29, 47]. In our study, we defined the QTL intensity for each SNP on interval Δ i as follows [29]:
$$ \mathrm{I}(i) = \frac{1}{N}\left[{\displaystyle {\sum}_{k=1}^N}\left({\displaystyle {\sum}_{q=1}^{N_{snp}^{(k)}}}{1}_{\left\{sn{p}_q^{(k)}\in \mathrm{T}\mathrm{S}\mathrm{D}\right\}}/{\Delta}_i\right)\right], $$

where i is the i th bin, N is the number of MCMC samples kept for analysis, and \( {1}_{\left\{sn{p}_q^{(k)}\in \mathrm{T}\mathrm{S}\mathrm{D}\right\}} \) is an indicator function with value 1 if the SNP in the k th sample falls in the distribution with large variance (TSD), and Δ i is the size of bin (constant number of markers in the current study) beginning from the i th position. To avoid having too many SNPs diluting the QTL effect and to balance different marker densities, the bin size was set to 14, 18 and 18 SNPs, which correspond to 1, 0.1, and 0.01 Mb for the MD, HD, and SEQ data, respectively. The QTL intensity for each bin was calculated as the average posterior probability of SNPs having a large effect within the bin.

The QTL intensity profile within a region is expected to show a peak if the region contains a QTL. Xu et al. [34] declared that the intensity profile is not able to distinguish two linked QTL regions and that the QTL effect profile is not able to accurately locate the causative mutation within an QTL region, while an intensity profile [48] that is weighted by estimated SNP effects can distinguish intervals with QTL and also show sharp peaks within intervals. Therefore, we also used the weighted QTL intensity profile that has been used in several previous studies [49, 48, 34]. The weighted QTL intensity is defined as:
$$ {\mathrm{I}}_{\mathrm{w}}(i) = \mathrm{I}(i)*\mathrm{w}(i) $$

where I(i) is the QTL intensity of i th bin, as given in Equation (5), and w(i) is the average SNP effect for the i th bin.

To detect QTL regions based on the weighted QTL intensity profile, bins that were three standard deviations above or below the mean of total bins were deleted as outliers (assuming that the bins are in a QTL region) and the mean and standard deviation of the remaining bins were re-calculated for the purpose of the t test. Then, a multiple t-test \( \Big(\mathrm{t} = \left(\mathrm{I}(i)\overline{-\mathrm{I}(i)}\right)/s \), where s is the standard deviation of the weighted QTL intensity of the bins) with a Bonferroni correction, was used to identify significant peaks among all the bins, including those outlier bins. Each significant peak represented a QTL region. A QTL region was specified as a region where the sum of the QTL intensities for the bins around the significant peak exceeded a predefined threshold (0.95). Finally, the average r2 (LD) between significant SNPs (BF ≥ 3.2) within a QTL region was calculated for each significant peak and between adjacent significant peaks. Adjacent peaks with an average r2 larger than 0.5 were merged into one region.

The position of each QTL region was defined according to the Bos taurus genome assembly UMD3.1 [50]. Genes that were located within or overlapped with the QTL regions were determined using information from the National Center for Biotechnology Information [51].

In the following, LM models with MD, HD, and SEQ data were designated LMMD, LMHD, and LMSEQ; BVS models with MD, HD, and SEQ data were designated BVSMD, BVSHD, and BVSSEQ; and BVS models using the QTL intensity profile with MD, HD, and SEQ data were designated BVSINTMD, BVSINTHD, and BVSINTSEQ, respectively.


Analysis using two models based on various marker densities

Figure 1 shows the associations of SNPs with mastitis using the LM and BVS models. For each model, the association patterns were similar for the three marker densities, although the number of significant SNPs decreased with decreasing marker densities, and both models showed peaks at similar locations. However, the BVS model presented clearer signals for QTL regions. In addition, with increasing marker densities, the peaks of putative QTL became sharper for the BVS model and the boundaries of adjacent QTL regions became more obvious for the LM model.
Fig. 1

Genome-wide plot for association of SNPs with udder health. (a) Plots generated by the LM model based on MD, HD, and SEQ data. (b) Plots generated by BVS model based on MD, HD, and SEQ data. Six chromosome regions between 84 and 95 Mb on BTA5, 88 and 96 Mb on BTA6, 57 and 63 Mb on BTA13, 48 and 55 Mb on BTA16, 55 and 58 Mb on BTA19, and 32 and 40 Mb on BTA20 are marked in alternating colors for clarity. The vertical axis is − log10(P) for the LM and BF for the BVS model, respectively; the horizontal dotted line indicates the genome-wide significance levels (for BVS, 3.2 is considered as putative and 10 is considered as strong evidence); the text on the subfigures presents the significant SNPs that were detected by the BVS model with a BF greater than 25

Details on the SNPs detected by LMMD, LMHD, LMSEQ, BVSMD, BVSHD, and BVSSEQ are in Table S1 [See Additional file 1: Table S1]. Table 1 presents the numbers of significant SNPs that were detected with the LM and BVS models. A total of 26 (MD), 75 (HD), and 465 (SEQ) SNPs were significant for both models. The total genetic variance estimated from the BVS model was equal to 110.20 for the MD data, and the variance explained by the SNPs of the six regions was equal to 13.10 (0.32 from SNPs with small effects and 12.77 from SNPs with large effects). Based on the HD data, the total genetic variance was equal to 113.6, and the variance explained by the SNPs of the six regions was equal to 20.54 (0.35 from SNPs with small effects and 20.20 from SNPs with large effects). Based on the SEQ data, the total genetic variance was equal to 123.70, and the variance explained by the SNPs of the six regions was equal to 30.50 (0.32 from SNPs with small effects and 30.18 from SNP with large effects). Using the LMHD, there were 5, 8, 3, 19, 17, and 27 significant SNPs in the six targeted regions with a higher test statistic (− log10(P)) than the highest test statistic obtained with the LMMD. Similarly, using LMSEQ, there were 0, 6, 0, 7, 2, and 0 SNPs with a − log10(P) greater than the highest − log10(P) achieved with the LMHD analysis for the six targeted regions.
Table 1

Number and percentage of significant SNPs






120 (20.4 %)

43 (7.3 %)


967 (12.4 %)

131 (1.7 %)


7209 (9.1 %)

1052 (1.3 %)

The analysis was implemented by the linear mixed model (LM) and the Bayesian variable selection model (BVS); figures in brackets are percentages of significant SNPs out of the total number of SNPs used in the analysis

The BVS model detected fewer significant SNPs in the six targeted regions than the LM model. Figure 2 graphically represents the position of significant SNPs on BTA6 using both models and the three marker densities. The LM model included one SNP for each run. Therefore, the − log10(P) values were consistent for the three marker densities with those of the LM model. However, when using the BVS model, the results among the three SNP datasets were not consistent. We observed three cases: (1) seven SNPs were detected by both BVSMD and BVSHD, while 70 SNPs were detected by both BVSHD and BVSSEQ; (2) some SNPs were only detected by BVSSEQ, e.g. the SNP at position 88 326 909 bp on BTA6, which is due to the fact that these SNPs were not present in MD and HD data; and (3) some SNPs were only detected by BVSMD, e.g the SNP at position 88 656 290 bp on BTA6, which may be explained by the QTL effect being distributed over several nearby SNPs, probably because the QTL was in strong LD with many markers in the HD and SEQ data.
Fig. 2

Positions of the significant SNPs on chromosome 6 between 88 and 89 Mb. a, b and c refer to the plots obtained with the LM model based on MD, HD, and SEQ data, respectively, while d, e and f refer to the plots obtained with the BVS model based on MD data, HD, and SEQ data, respectively. The horizontal axis indicates the position of the SNPs. The filled diamond indicates the genome-wide significant SNPs obtained with the LM and BVS (BF > 3.2) models. The SNP with a BF or − log10(P) greater than 25 was set at 25

Analysis of QTL intensity

QTL intensity profiles from the analysis based on the three marker datasets are in Fig. 3. The corresponding weighted QTL intensities are in Fig. 4. The association patterns of the weighted QTL intensities were consistent across the three marker densities. BVSINTMD, BVSINTHD, and BVSINTSEQ detected one, 16, and 33 QTL intensity peaks, respectively (Fig. 5). Among the SNPs that were significant with the BVS models (BF > 10), 36.4, 87.5, and 86.7 % were within the QTL regions identified by BVSINTMD, BVSINTHD, and BVSINTSEQ, respectively. The positions and intervals of the detected QTL intensity peaks by BVSINTMD and BVSINTHD are in Table 2. BVSINTMD detected only one QTL intensity peak on BTA6. Table 3 shows the 51 genes that were located within or adjacent to the QTL intensity peaks detected by BVSINTSEQ. Among these QTL intensity peaks, 27 were located within or overlapped with known genes, while the others were 5 to 165 kb away from the nearest known gene. The average LD (r2) of the 33 QTL intensity regions was equal to 0.68. The average interval length of the QTL intensity peaks became smaller as marker densities increased and was approximately 1.39, 0.22, and 0.10 Mb for BVSINTMD, BVSINTHD, and BVSINTSEQ, respectively.
Fig. 3

QTL intensity profiles for different marker densities. From top to bottom are shown the plots based on MD, HD, and SEQ data. The six chromosome regions are between 84 and 95 Mb on BTA5, 88 and 96 Mb on BTA6, 57 and 63 Mb on BTA13, 48 and 55 Mb on BTA16, 55 and 58 Mb on BTA19, and 32 and 40 Mb on BTA20; the regions are separated by vertical dotted lines

Fig. 4

Weighted QTL intensity profiles for different marker densities. From top to bottom are shown the plots based on MD, HD, and SEQ data. The six chromosome regions are between 84 and 95 Mb on BTA5, 88 and 96 Mb on BTA6, 57 and 63 Mb on BTA13, 48 and 55 Mb on BTA16, 55 and 58 Mb on BTA19, and 32 and 40 Mb on BTA20; the regions are separated by vertical dotted lines

Fig. 5

Weighted QTL intensity peaks detected by multiple t-tests for different marker densities. Blue solid circles are the weighted QTL intensities that are significant with the t test. From top to bottom are shown the plots based on MD, HD, and SEQ data, respectively. The six chromosome regions are between 84 and 95 Mb on BTA5, 88 and 96 Mb on BTA6, 57 and 63 Mb on BTA13, 48 and 55 Mb on BTA16, 55 and 58 Mb on BTA19, and 32 and 40 Mb on BTA20; the regions are separated by vertical dotted lines

Table 2

QTL regions detected by QTL intensity analysis using the MD and HD data



Start_pos (bp)

End_pos (bp)

Reg_size (bp)

















































































Chr = chromosome; Start_pos = start position of the QTL region; End_pos = end position of the QTL region; Reg_size = size of the QTL region

Table 3

Genes located within or closest to the QTL regions





Gene symbol























































































































SLMO2, LOC100140850













































































SAP30BP, RECQL5, SMIM5, LOC100302389,MYO15B, LOC101907888






































QTL regions were detected by analysis of QTL intensity using SEQ data; Chr = chromosome; Start_pos = start position of the QTL region; End_pos = end position of the QTL region; Reg_size = size of the QTL region; Distance = distance between gene and QTL region; r2 = linkage disequilibrium between significant SNPs (BF ≥ 3.2) within each region


To our knowledge, only two studies have performed sequence-based association studies for clinical mastitis [52, 38]. In this study, based on SEQ data and a BVS model for QTL intensity analysis, 33 QTL regions with an average length of 0.069 Mb were detected in the six pre-selected chromosome regions.

Comparison of the two models

The association patterns between SNPs and mastitis were similar for the LM and BVS models. However, signals were much more distinct when using the BVS model (Fig. 1). This indicates that the BVS model is able to identify and narrow down regions of putative QTL and to distinguish linked QTL.

According to the significance level used to detect SNPs in this study, the number of significant SNPs identified with the BVS model was smaller than with the LM model (see Table 1). The difference in the percentage of significant SNPs between the two models increased with increasing marker densities. This is because the LM model is a single-SNP model and, thus, a QTL effect is determined by a single SNP in each run. Consequently, many SNPs in LD with the QTL can have a significant effect on the model’s output. In contrast, the BVS model estimates the effects of all SNPs simultaneously, and the effect of a QTL may be represented by a single SNP or distributed over several SNPs that are in strong LD with the QTL [53] (i.e., several SNPs could together represent the effect of a single QTL).

Comparison of the three marker densities

In a previous study, using a simulated dataset, Cleveland et al. [53] reported that Bayesian analysis showed clear QTL signals when SNPs were clustered in a 0.5 to 5.0 cM region. In this study, three marker densities were analyzed, with average distances between markers of 70 kb (for MD), 5.5 kb (for HD), and 0.5 kb (for SEQ).

Our results (Fig. 1) showed that peak locations were largely consistent for the three marker densities. Spencer et al. [54] reported that marked differences in genome coverage may not translate into appreciable differences in power to detect causative variants when using an additive model. However, it is obvious that there were fewer markers near a QTL in the MD data compared to the HD data, and therefore, there were few markers in LD with a QTL [53].

Since, in this study, a single-SNP regression was performed using the LM model and all lower-density SNPs were included in the higher-density dataset, all significant SNPs that were identified with LMMD were also detected with LMHD, and all significant SNPs that were identified with LMHD were also detected with LMSEQ if the same significance threshold was used. However, some SNPs that were detected with LMMD did not reach significance with LMHD, and some SNPs detected with LMHD did not reach significance with LMSEQ, because of the different significance thresholds for the three marker densities as a result of Bonferroni multiple-testing correction. In addition, as density increased, the LD within a given region became stronger and a larger number of markers around a QTL could show an effect on the trait. Therefore, increasing marker density does not necessarily clarify the QTL region boundaries when using LM analysis. In addition, there were 79 significant SNPs for LMHD with a − log10(P) that was larger than the highest − log10(P) obtained with LMMD for the six regions, but only 15 significant SNPs for LMSEQ with a − log10(P) that was larger than the highest − log10(P) obtained with LMHD for the six regions. Thus, the relative increase in associated SNPs from the HD to the SEQ data was less than the increase from the MD to the HD data. This may be due to imperfect imputation of sequence data due to the small reference size or may reflect that the density of the HD data was sufficient for QTL detection and further density increases only give small improvements. However, it should also be noted that both the MD and HD data have ascertainment biases towards common variants, which are more suitable for GWAS.

In contrast, some significant SNPs in the BVSMD set were not detected when using BVSHD and BVSSEQ data, and some significant SNPs in the BVSHD set were not detected in the BVSSEQ set, although the same BF criterion was used for the three maker densities. Moreover, increasing marker density did not result in a higher BF. Since the BVS model fits all markers simultaneously, it is possible that when markers are very close to a QTL but are not part of low-density data, the QTL effect may be shifted to nearby markers [53], or a QTL effect may be represented by a single SNP when using low-density data but shared by many SNPs when using high-density data. This suggests that instead of a single SNP, a summary statistic across a small region is necessary to detect a QTL when using the BVS model.

In spite of the interesting properties of high-throughput sequencing, it is necessary to take some of the limitations of imputed sequence data into consideration. Only 89 sequenced animals from three breeds were used as imputation reference population and, therefore, imputed sequences were expected to have relatively low imputation accuracy, especially for SNPs with low minor allelic frequencies. Brøndum et al. [42] reported an imputation accuracy of around 0.90 when the reference population included 242 individuals, which is much lower than the accuracy obtained (0.97) in imputation of HD markers from 50 k data for Nordic Holstein cattle [55]. In addition, the errors in calling variants may be higher with high-throughput sequencing than with SNP-array genotyping data [42]. Daetwyler et al. [22] also pointed out that results of association studies that are based on imputed sequences should be interpreted with caution, since SNPs with slightly higher P values than the most significant SNPs from multiple-testing can also be considered as potential causative mutations, particularly if there is strong supporting functional evidence. In addition, with the very high density of SNPs in sequence data, it is expected that a large number of SNPs will be in high or near perfect LD with the causal variants. This makes it difficult to distinguish causal variants from SNPs in high or near perfect LD with the causal variants using a Bayesian approach since the QTL effect may be diluted across many neighboring SNPs. Therefore, in this study, we removed one of the SNPs from a pair in which two SNPs were in complete LD although this may cause the loss of causative variants. However, it is most likely that adjacent SNPs that are in complete LD are close to each other. In other words, if a causative variant is removed by LD-based pruning, it will be replaced by another SNP that is closely associated with the causal variant. Thus in these conditions, the putative loss of causative variants should have little influence on the identification of QTL peaks. In fact, if a neutral SNP is in complete LD with a causal SNP, the latter cannot be differentiated using either the single-SNP linear mixed model or the Bayesian model. However, the Bayesian approach is more appropriate for narrowing QTL regions and distinguishing multiple QTL regions that segregate closely. Post-analysis based on bins using the Bayesian approach without pruning the SNPs may avoid the possible loss of causal variants and dilution of QTL effects.

QTL regions

In the BVS model, the QTL effect may be distributed over several markers that are in LD with the QTL. Therefore, combining the posterior probabilities of closely located markers can result in higher power when inferring the presence of a QTL [23]. Detection of QTL regions based on QTL density profiles was initially proposed by Sillanpaa and Arjas [29]. Using simulated data, QTL regions are easy to determine because of the simplicity of the simulation [56]. However, when using real data, the significance threshold for QTL intensity peaks is not as well defined because of the low signal to noise ratio (Fig. 4). Previous studies have detected QTL peaks by using QTL intensity, but only a few defined the interval and specific position of QTL regions [34, 56, 48]. In this study, we used multiple t-tests based on weighted QTL intensities to determine QTL regions.

Based on the present data, the association patterns of QTL intensities (Fig. 4) and BF (Fig. 1) were similar. However, the QTL intensities presented clear peaks. Analysis of the QTL intensity profiles detected only one QTL peak on BTA6. This indicates that studies with higher marker densities have greater power to detect QTL regions by QTL intensity profile analysis. In addition, the average length of the detected QTL regions became shorter as marker densities increased (Tables 2 and 3), which suggests that QTL regions can be refined by using high-density markers and applying the BVS model together with a QTL intensity profile.

Based on SEQ data, the average r2 (LD) of significant SNPs between a pair of QTL regions was lower than 0.55, while the average r2 of significant SNPs within each QTL region (average 0.68) was higher than the r2 among different QTL regions. In addition, comparisons between QTL intensity profiles and BF (>10) showed that 87 % of the significant SNPs detected by BF were located within the QTL regions determined by QTL intensity profiles based on the HD and SEQ data. Among the 33 QTL regions detected by QTL intensity profiles using the SEQ data, 27 regions contained a known gene. For the remaining six regions, the closest gene was at most 165 kb away.

Five of the detected regions have been reported by previous studies or contain genes that are known to be functionally associated with mastitis traits. One of these QTL regions (between 88 056 115 and 88 164 130 bp on BTA6) overlaps with the deoxycytidine kinase (DCK) gene. Abdel-Shafy et al. [57] reported that this gene is associated with somatic cell score, which is a mastitis indicator [58]. Another QTL region (between 88 355 595 and 88 476 214 bp on BTA6) overlapped with the solute carrier family 4, sodium bicarbonate co-transporter, and member 4 (SLC4A4) gene belonging to the SLC4 family. Sodeland et al. [59] reported a high correlation of SNPs surrounding the SLC4A4 gene with clinical mastitis. A third QTL region between 88 984 167 and 89 072 556 bp on BTA6, overlapped with a previously reported QTL region between 89 and 91 Mb for mastitis [60]. This region contains a part of the neuropeptide FF receptor 2 (NPFFR2) gene that encodes a member of a G-protein-coupled neuropeptide receptors subfamily that is activated by the neuropeptides A-18-amide (NPAF) and F-8-amide (NPFF) [61]. Sun et al. [61] showed that NPFF is involved in anti-inflammatory effects, both in vitro and in vivo, and Minault et al. [62] showed that it modulates the proliferation of human T lymphocytes. In addition, Sodeland et al. [59] detected a QTL region around 35.5 Mb on BTA20 which is associated with clinical mastitis. A fourth QTL region, between 35 900 432 and 36 010 180 bp on BTA20, was detected and overlapped with the leukemia inhibitory factor receptor alpha (LIFR) gene. LIFR is a breast cancer metastasis suppressor upstream of the Hippo-YAP pathway and a prognostic marker [63]. Finally, a QTL region between 57 509 023 and 57 581 780 bp on BTA13 overlapped with the QTL at 57.54 Mb reported by Sahana et al. [64]. This region overlaps with the endothelin 3 (EDN3) gene, which influences neutrophil activation [65]. In blood, neutrophils are the major leukocytes that respond to inflammatory stimuli. Therefore, NPFFR2, SLC4A4, DCK, LIFR, and EDN3 are candidate genes for susceptibility to mastitis.


The power of QTL detection can be increased by increasing marker densities and the BVS model outperforms the LM model in refining QTL locations with clear boundaries between linked QTL. Based on the results obtained with the SEQ data, six preselected regions were refined into 33 candidate QTL regions for udder health. Furthermore, the comparison between these candidate QTL regions and known genes suggests that NPFFR2, SLC4A4, DCK, LIFR, and EDN3 may be considered as candidate genes for mastitis susceptibility. Further studies are required to validate the causative loci that underlie these QTL and to investigate the function of the candidate genes that affect udder health.



This work was supported by the research project “Genomic Selection – from function to efficient utilization in cattle breeding (grant no. 3412-08-02253, Denmark)” and the Chinese Scholarship Council (CSC) that provided a scholarship for the first author.

Authors’ Affiliations

Center for Quantitative Genetics and Genomics, Department of Molecular Biology and Genetics, Aarhus University
Key Laboratory of Animal Genetics and Breeding of Ministry of Agriculture, National Engineering Laboratory for Animal Breeding, College of Animal Science and Technology, China Agricultural University


  1. Hogeveen H, Pyorala S, Waller KP, Hogan JS, Lam TJGM, Oliver SP, et al. Current status and future challenges in mastitis research. In Proceedings of the National Mastitis Council 50th Annual Meeting: 23–26 January 2011; Arlington: the National Mastitis Council 50th Annual Meeting; 2011:36–48.Google Scholar
  2. Awale M, Dudhatra G, Avinash K, Chauhan B, Kamani D, Modi C, et al. Bovine mastitis: a threat to economy. Sci Rep. 2012;1:1–10.Google Scholar
  3. Hameed KGA, Sender G, Korwin-Kossakowska A. Public health hazard due to mastitis in dairy cows. Anim Sci Pap Rep. 2007;25:73–85.Google Scholar
  4. Bennedsgaard TW, Enevoldsen C, Thamsborg SM, Vaarst M. Effect of mastitis treatment and somatic cell counts on milk yield in Danish organic dairy cows. J Dairy Sci. 2003;86:3174–83.PubMedView ArticleGoogle Scholar
  5. Nash DL, Rogers GW, Cooper JB, Hargrove GL, Keown JF, Hansen LB. Heritability of clinical mastitis incidence and relationships with sire transmitting abilities for somatic cell score, udder type traits, productive life, and protein yield. J Dairy Sci. 2000;83:2350–60.PubMedView ArticleGoogle Scholar
  6. Lund T, Miglior F, Dekkers JCM, Burnside EB. Genetic-relationships between clinical mastitis, somatic-cell count, and udder conformation in Danish Holsteins. Livest Prod Sci. 1994;39:243–51.View ArticleGoogle Scholar
  7. Hansen M, Lund MS, Sorensen MK, Christensen LG. Genetic parameters of dairy character, protein yield, clinical mastitis, and other diseases in the Danish Holstein cattle. J Dairy Sci. 2002;85:445–52.PubMedView ArticleGoogle Scholar
  8. Heringstad B, Klemetsdal G, Ruane J. Selection for mastitis resistance in dairy cattle: a review with focus on the situation in the Nordic countries. Livest Prod Sci. 2000;64:95–106.View ArticleGoogle Scholar
  9. Rupp R, Boichard D. Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins. J Dairy Sci. 1999;82:2198–204.PubMedView ArticleGoogle Scholar
  10. Lund MS, Jensen J, Petersen PH. Estimation of genetic and phenotypic parameters for clinical mastitis somatic cell production deviance, and protein yield in dairy cattle using Gibbs sampling. J Dairy Sci. 1999;82:1045–51.PubMedView ArticleGoogle Scholar
  11. Lescourret F, Coulon JB. Modeling the impact of mastitis on milk-production by dairy cows. J Dairy Sci. 1994;77:2289–301.PubMedView ArticleGoogle Scholar
  12. Georges M. Mapping, fine mapping, and molecular dissection of quantitative trait Loci in domestic animals. Annu Rev Genomics Hum Genet. 2007;8:131–62.PubMedView ArticleGoogle Scholar
  13. Sahana G, Lund MS, Andersson-Eklund L, Hastings N, Fernandez A, Iso-Touru T, et al. Fine-mapping QTL for mastitis resistance on BTA9 in three Nordic red cattle breeds. Anim Genet. 2008;39:354–62.PubMedView ArticlePubMed CentralGoogle Scholar
  14. Holmberg M, Andersson-Eklund L. Quantitative trait loci affecting health traits in Swedish dairy cattle. J Dairy Sci. 2004;87:2653–9.PubMedView ArticleGoogle Scholar
  15. Klungland H, Sabry A, Heringstad B, Olsen HG, Gomez-Raya L, Vage DI, et al. Quantitative trait loci affecting clinical mastitis and somatic cell count in dairy cattle. Mamm Genome. 2001;12:837–42.PubMedView ArticleGoogle Scholar
  16. Kühn C, Bennewitz J, Reinsch N, Xu N, Thomsen H, Looft C, et al. Quantitative trait loci mapping of functional traits in the German Holstein cattle population. J Dairy Sci. 2003;86:360–8.PubMedView ArticleGoogle Scholar
  17. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high density SNP genotyping assay for cattle. PLoS One. 2009;4, e5350.PubMedView ArticlePubMed CentralGoogle Scholar
  18. Cole JB, Wiggans GR, Ma L, Sonstegard TS, Lawlor Jr TJ, Crooker BA, et al. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics. 2011;12:408.PubMedView ArticlePubMed CentralGoogle Scholar
  19. Mai MD, Sahana G, Christiansen FB, Guldbrandtsen B. A genome-wide association study for milk production traits in Danish Jersey cattle using a 50 K single nucleotide polymorphism chip. J Anim Sci. 2010;88:3522–8.PubMedView ArticleGoogle Scholar
  20. Matukumalli LK, Schroeder S, DeNise SK, Sonstegard T, Lawley CT, Georges M, et al. Analyzing LD blocks and CNV segments in cattle: novel genomic features identified using the BovineHD BeadChip. San Diego: Illumina Inc; 2011. Pub. No. 370-2011-002.Google Scholar
  21. Rincon G, Weber KL, Eenennaam AL, Golden BL, Medrano JF. Hot topic: performance of bovine high-density genotyping platforms in Holsteins and Jerseys. J Dairy Sci. 2011;94:6116–21.PubMedView ArticleGoogle Scholar
  22. Daetwyler HD, Capitan A, Pausch H, Stothard P, van Binsbergen R, Brondum RF, et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. Nat Genet. 2014;46:858–65.PubMedView ArticleGoogle Scholar
  23. Sahana G, Guldbrandtsen B, Janss L, Lund MS. Comparison of association mapping methods in a complex pedigreed population. Genet Epidemiol. 2010;34:455–62.PubMedView ArticleGoogle Scholar
  24. Yu JM, Pressoir G, Briggs WH, Vroh Bi I, Yamasaki M, Doebley JF, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat Genet. 2006;38:203–8.PubMedView ArticleGoogle Scholar
  25. Dashab GR, Kadri NK, Shariati MM, Sahana G. Comparison of linear mixed model analysis and genealogy-based haplotype clustering with a Bayesian approach for association mapping in a pedigreed population. BMC Proc. 2012;6:S4.PubMedView ArticlePubMed CentralGoogle Scholar
  26. George EI, McCulloch RE. Variable selection via Gibbs sampling. J Am Stat Assoc. 1993;88:881–9.View ArticleGoogle Scholar
  27. Sahana G, Kadlecova V, Hornshoj H, Nielsen B, Christensen OF. A genome-wide association scan in pig identifies novel regions associated with feed efficiency trait. J Anim Sci. 2013;91:1041–50.PubMedView ArticleGoogle Scholar
  28. Pausch H, Jung S, Edel C, Emmerling R, Krogmeier D, Gotz KU, et al. Genome-wide association study uncovers four QTL predisposing to supernumerary teats in cattle. Anim Genet. 2012;43:689–95.PubMedView ArticleGoogle Scholar
  29. Sillanpää MJ, Arjas E. Bayesian mapping of multiple quantitative trait loci from incomplete inbred line cross data. Genetics. 1998;148:1373–88.PubMedPubMed CentralGoogle Scholar
  30. Satagopan JM, Yandell YS, Newton MA, Osborn TC. A Bayesian approach to detect quantitative trait loci using Markov chain Monte Carlo. Genetics. 1996;144:805–16.PubMedPubMed CentralGoogle Scholar
  31. Heath SC. Markov chain Monte Carlo segregation and linkage analysis for oligogenic models. Am J Hum Genet. 1997;61:748–60.PubMedView ArticlePubMed CentralGoogle Scholar
  32. Bush WS, Moore JH. Chapter 11: Genome-wide association studies. PLoS Comput Biol. 2012;8, e1002822.PubMedView ArticlePubMed CentralGoogle Scholar
  33. Kass RE, Raftery AE. Bayes Factors. J Am Stat Assoc. 1995;90:773–95.View ArticleGoogle Scholar
  34. Xu CW, Wang XF, Li ZK, Xu SZ. Mapping QTL for multiple traits using Bayesian statistics. Genet Res (Camb). 2009;91:23–37.View ArticleGoogle Scholar
  35. Nordic cattle genetics evaluation. Accessed 20 March 2013.
  36. Goddard M. A method of comparing sires evaluated in different countries. Livest Prod Sci. 1985;13:321–31.View ArticleGoogle Scholar
  37. Schaeffer LR. Model for international evaluation of dairy sires. Livest Prod Sci. 1985;12:105–15.View ArticleGoogle Scholar
  38. Sahana G, Guldbrandtsen B, Thomsen B, Holm LE, Panitz F, Brondum RF, et al. Genome-wide association study using high-density single nucleotide polymorphism arrays and whole-genome sequences for clinical mastitis traits in dairy cattle. J Dairy Sci. 2014;97:7258–75.PubMedView ArticleGoogle Scholar
  39. Browning BL, Browning SR. A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet. 2009;84:210–23.PubMedView ArticlePubMed CentralGoogle Scholar
  40. Lund MS, de Roos APW, de Vries AG, Druet T, Ducrocq V, Fritz S, et al. A common reference population from four European Holstein populations increases reliability of genomic predictions. Genet Sel Evol. 2011;43:43.PubMedView ArticlePubMed CentralGoogle Scholar
  41. Hoglund J, Sahana G, Brondum R, Guldbrandtsen B, Buitenhuis B, Lund M. Fine mapping QTL for female fertility on BTA04 and BTA13 in dairy cattle using HD SNP and sequence data. BMC Genomics. 2014;15:790.PubMedView ArticlePubMed CentralGoogle Scholar
  42. Brondum RF, Guldbrandtsen B, Sahana G, Lund MS, Su G. Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle. BMC Genomics. 2014;15:728.PubMedView ArticlePubMed CentralGoogle Scholar
  43. Madsen P., Jensen J. A user's guide to DMU. University of Aarhus, Faculty Agricultural Sciences. 2010. Accessed 20 June 2013.
  44. Meuwissen THE. Accuracy of breeding values of 'unrelated' individuals predicted by dense SNP genotyping. Genet Sel Evol. 2009;41:35.PubMedView ArticlePubMed CentralGoogle Scholar
  45. Heuven HC, Janss LL. Bayesian multi-qtl mapping for growth curve parameters. BMC Proc. 2010;4 Suppl 1:S12.PubMedView ArticlePubMed CentralGoogle Scholar
  46. Janss L. Genomic models in bayz. Accessed 20 June 2013.
  47. Wang H, Zhang YM, Li XM, Masinde GL, Mohan S, Baylink DJ, et al. Bayesian shrinkage estimation of quantitative trait loci parameters. Genetics. 2005;170:465–80.PubMedView ArticlePubMed CentralGoogle Scholar
  48. Xu S, Yi N. Mixed model analysis of quantitative trait loci. Proc Natl Acad Sci U S A. 2000;97:14542–7.PubMedView ArticlePubMed CentralGoogle Scholar
  49. Fang M. Bayesian shrinkage mapping of quantitative trait loci in variance component models. BMC Genet. 2010;11:30.PubMedView ArticlePubMed CentralGoogle Scholar
  50. Zimin AV, Delcher AL, Florea L, Kelley DR, Schatz MC, Puiu D, et al. A whole-genome assembly of the domestic cow, Bos taurus. Genome Biol. 2009;10:R42.PubMedView ArticlePubMed CentralGoogle Scholar
  51. National Center for Biotechnology Information. Accessed 20 July 2013.
  52. Qanbari S, Pausch H, Jansen S, Somel M, Strom TM, Fries R, et al. Classic selective sweeps revealed by massive sequencing in cattle. PLoS Genet. 2014;10, e1004148.PubMedView ArticlePubMed CentralGoogle Scholar
  53. Cleveland MA, Deeb N. Evaluation of a genome-wide approach to multiple marker association considering different marker densities. BMC Proc. 2009;3:S5.PubMedView ArticlePubMed CentralGoogle Scholar
  54. Spencer CCA, Su Z, Donnelly P, Marchini J. Designing genome-wide association studies: sample size, power, imputation, and the choice of genotyping chip. PLoS Genet. 2009;5, e1000477.PubMedView ArticlePubMed CentralGoogle Scholar
  55. Brondum RF, Ma P, Lund MS, Su G. Short communication: genotype imputation within and across Nordic cattle breeds. J Dairy Sci. 2012;95:6795–800.PubMedView ArticleGoogle Scholar
  56. Liu JF, Liu YJ, Liu XG, Deng HW. Bayesian mapping of quantitative trait loci for multiple complex traits with the use of variance components. Am J Hum Genet. 2007;81:304–20.PubMedView ArticlePubMed CentralGoogle Scholar
  57. Abdel-Shafy H, Bortfeldt RH, Tetens J, Brockmann GA. Single nucleotide polymorphism and haplotype effects associated with somatic cell score in German Holstein cattle. Genet Sel Evol. 2014;46:35.PubMedView ArticlePubMed CentralGoogle Scholar
  58. Heringstad B, Gianola D, Chang YM, Odegard J, Klemetsdal G. Genetic associations between clinical mastitis and somatic cell score in early first-lactation cows. J Dairy Sci. 2006;89:2236–44.PubMedView ArticleGoogle Scholar
  59. Sodeland M, Kent MP, Olsen HG, Opsal MA, Svendsen M, Sehested E, et al. Quantitative trait loci for clinical mastitis on chromosomes 2, 6, 14 and 20 in Norwegian Red cattle. Anim Genet. 2011;42:457–65.PubMedView ArticleGoogle Scholar
  60. Sodeland M, Grove H, Kent M, Taylor S, Svendsen M, Hayes BJ, et al. Molecular characterization of a long range haplotype affecting protein yield and mastitis susceptibility in Norwegian Red cattle. BMC Genet. 2011;12:70.PubMedView ArticlePubMed CentralGoogle Scholar
  61. Sun YL, Zhang XY, Sun T, He N, Li JY, Zhuang Y, et al. The anti-inflammatory potential of neuropeptide FF in vitro and in vivo. Peptides. 2013;47:124–32.PubMedView ArticleGoogle Scholar
  62. Minault M, Lecron JC, Labrouche S, Simonnet G, Gombert J. Characterization of binding-sites for neuropeptide Ff on T-lymphocytes of the Jurkat cell-line. Peptides. 1995;16:105–11.PubMedView ArticleGoogle Scholar
  63. Chen DH, Sun YT, Wei YK, Zhang PJ, Rezaeian AH, Teruya-Feldstein J, et al. LIFR is a breast cancer metastasis suppressor upstream of the Hippo-YAP pathway and a prognostic marker. Nat Med. 2012;18:1511–7.PubMedView ArticlePubMed CentralGoogle Scholar
  64. Sahana G, Guldbrandtsen B, Thomsen B, Lund MS. Confirmation and fine-mapping of clinical mastitis and somatic cell score QTL in Nordic Holstein cattle. Anim Genet. 2013;44:620–6.PubMedView ArticleGoogle Scholar
  65. Elferink JGR, De Koster BM. The involvement of protein kinase G in stimulation of neutrophil migration by endothelins. Eur J Pharmacol. 1998;350:285–91.PubMedView ArticleGoogle Scholar


© Wu et al. 2015

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.