Background {#Sec1}
==========

The cow´s ability to reproduce is essential for milk production. Impaired reproduction will result in additional inseminations, higher replacement rate and increased culling rate. Recently genome sequencing technologies and bioinformatics analysis approaches have advanced tremendously. The availability of full genome sequence data can help to identify causal mutations underlying variation in female fertility. Application of these tools has led to remarkable increases in the numbers of trait markers available and thus enhanced precision of QTL mapping.

Genome-wide association studies (GWAS) examine common genetic variants in large numbers of individuals to determine whether an association with quantitative traits exists. GWAS have identified thousands of single nucleotide polymorphisms (SNPs) across the cattle genome associated with economically important traits in cattle (e.g. \[[@CR1]\]). However, most of the genetic variants detected by these studies are not causal for the traits themselves. Instead they are in linkage disequilibrium (LD) with the causal polymorphisms. Most GWAS studies in cattle have used data from one breed only. In most cases this has been the Holstein Friesian (HF) where the largest datasets are available. The low effective population size in HF has resulted in long-range LD \[[@CR2], [@CR3]\]. This limits the studies' ability to distinguish between causal factors and markers in strong LD with causal factors. The limitation can sometimes be overcome by using data from multiple breeds. A second limitation of previous mapping studies is that the marker panels used only represent a small fraction of the variants segregating in the population. Using a panel of individuals with whole genome sequence, in principle all SNPs can be imputed for all individuals in the mapping population \[[@CR4]\]. The data thereby comes to include the causal variants. Causal effects thereby become identifiable.

In addition to revealing the genetic architecture that underlies the physiological and biological process of female reproduction, this information could also be practically applied to genomic selection schemes. Genomic prediction helps to select breeding animals for the next generation more efficiently. Introducing high density SNP arrays (777 k) did not substantial increase accuracy of genomic predictions in cattle (0.5 to 1%) as compared to medium density 50 k SNP arrays \[[@CR5]\]. One of the reasons discussed was the increase in the number of unknown parameters to be estimated with high density data. With the availability of full genome sequence data this problem increases many folds. However, if the causal mutations underlying female fertility are identified and not only markers in LD with the causal mutation, the information could be included in genomic prediction models where additional weight can be put on certain genomic regions/variants which influence female fertility. This would in particular improve predictions over generations and across breeds.

Previously QTL areas influencing female fertility traits have been identified in the Nordic Holstein cattle population using the 50 k SNP array \[[@CR6]\]. In this study two of these QTL regions previously identified \[[@CR6]\] (BTA04 and BTA13) were first analyzed using the HD 777 k SNP array to further refine the regions in the NH. In addition, the Danish Jersey (DJ) and Nordic Red (NR) breeds were screened with the bovine HD array to investigate if the QTL on BTA04 and BTA13 were segregating in these breeds as well. To further refine the QTL regions in the NH, these regions were re-analyzed with the imputed full sequence data.

Results {#Sec2}
=======

HD SNP analysis {#Sec3}
---------------

BTA04 and BTA13 were screened with the bovine HD SNP array in the NH, NR and DJ breeds, based on the positions on the genome which were chosen due to high significance from Höglund et al. \[[@CR6]\]. BTA04 was screened for NRRH, IFLH, NRRC, ICF, IFLC, AISH, AISC in a 2,633,026 bp region spanning positions 38,257,758 to 40,890,784 bp (Figure  [1](#Fig1){ref-type="fig"}: left), whereas BTA13 was screened for ICF only, in the region ranging from 21,236,959 to 46,150,079 bp (Figure  [2](#Fig2){ref-type="fig"}: left). ICF was the only trait which showed significance in this position on BTA13 in the previous study \[[@CR6]\]. Correlations of the minor allele frequencies of the most significant markers in the NH, DJ and NR pointed in the direction that the QTL for ICF is breed specific on BTA13. However the results of correlations between the minor allele frequencies of the most significant marker in the three breeds on BTA04 were inconclusive whether the QTLs were breed specific. SNP markers in the NH breed reached significance levels on both BTA04 and BTA13 (Figure  [1](#Fig1){ref-type="fig"}: center and Figure  [2](#Fig2){ref-type="fig"} center; Additional files [1](#MOESM1){ref-type="media"} and [2](#MOESM2){ref-type="media"}). However, no markers in the DJ and NR breeds reached significance analyzing HD SNP data.Figure 1**Manhattan plots of GWAS on BTA04 for 5 female fertility traits with three densities of markers across successively narrower regions.** The y axes show --log~10~(p-value) of single-marker associations tests. The x axes show marker positions in base pairs (bp). Left panel: tests across three breeds with BovineSNP50 BeadChip \[[@CR6]\]. Vertical lines indicate the region shown in the middle panel. Center panel: tests in Holstein with the HD BeadChip. Vertical lines indicate the region shown in the right panel. Right panel: tests using sequence data for the Nordic Holstein breed.Figure 2**Manhattan plots of the GWAS results on BTA13 for ICF with three densities of markers across successively narrower regions.** The y axes show --log~10~(p-value) of single-marker associations tests. The x axes show marker positions in base pairs (bp). Left panel: GWAS result of the Nordic Holstein breed with the BovineSNP50 BeadChip \[[@CR6]\]. The red dots indicate the region shown in the center panel. Center panel: test result with the HD BeadChip. The blue dots indicate the region shown in the right panel. Right panel: test result using sequence data.

Sequence analysis {#Sec4}
-----------------

### BTA04 {#Sec5}

Based on the results using the HD array (Figure  [1](#Fig1){ref-type="fig"}: middle) an area from 38,721,162 -- 40,890,784 bp was selected for further analysis with full sequence data for this region. Based on the sequence data analysis (Figure  [1](#Fig1){ref-type="fig"}: right; Additional file [3](#MOESM3){ref-type="media"}) the QTL position was narrowed down to two regions. Region 1 spanned positions between 39,210,643 -- 39,498,957 bp (=288,314 bp) and Region 2 spanned 39,700,194 -- 40,890,784 bp (=1,190,590 bp). The two regions were searched for candidate genes/polymorphisms (Table  [1](#Tab1){ref-type="table"}). In total 9 genes were identified in these regions of which 5 were annotated, three were described as uncharacterized proteins and one was not annotated (Table  [1](#Tab1){ref-type="table"}). The annotated markers showing association with AISC, IFLC, NRRC and ICF are shown in Table  [2](#Tab2){ref-type="table"}. The annotated markers were not the markers with the highest --log~10~(P-value). The highest --log~10~(P-values) were 15.22 for AISC (Chr4:40298743), 16.14 for IFLC (Chr4:39213491), 7.18 for NRRC (Chr4:40728978) and 5.39 for ICF (Chr4:40599222).Table 1**Genes located on BTA04 in region 2 (39,700,194- 40,890,784 bp)**Ensembl gene IDDescriptionGene start (bp)Gene end (bp)Gene nameENSBTAG00000014800Known pseudogene4006766440068936ENSBTAG00000006138Semaphorin-3C4014049440345588SEMA3CENSBTAG00000047646Uncharacterized protein4043290640457597ENSBTAG00000014220Uncharacterized protein4045745040560780ENSBTAG00000046905Uncharacterized protein4054170140542459ENSBTAG00000017866Platelet glycoprotein 44058148440643369CD36ENSBTAG00000008641guanine nucleotide-binding protein G(t) subunit alpha-34077588540830619GNAT3ENSBTAG000000292925S ribosomal RNA40887316408874335S rRNAENSBTAG00000044967U6 spliceosomal RNA3975314939753244U6Table 2**Annotation of the markers on BTA04 associated with AISC, IFLC, NRRC and ICF in region 2 (39,700,194- 40,890,784 bp)**MarkerDescriptionPosition (bp)AISC ^1^ -log(p-value)IFLC ^1^ -log(p-value)NRRC ^1^ -log(p-value)ICF ^1^ -log(p-value)Gene namers43384664Non coding exon variant397532162.46052.81910.85070.7718rs43383647Non coding exon variant400677252.50252.58391.85740.2769rs43383646Non coding exon variant400677634.33963.41172.44950.7627rs43383645Non coding exon variant400677742.50292.58451.85770.2770rs43383644Non coding exon variant400678342.52732.60061.88790.2715rs43383643Non coding exon variant400679783.93035.17061.82901.7758rs109658404Non coding exon variant400680592.07452.07251.30270.2634rs43383637Non coding exon variant400686663.92275.15471.82721.7656rs43383636Non coding exon variant400686993.92275.15471.82721.7656rs109959240Non coding exon variant400687662.15162.14731.34160.2780rs136871729Non coding exon variant400689032.15592.15521.3450.280740322885C/ASynonymous variant403228852.74874.11290.27734.5324SEMA3C40438940A/GMissense variant404389401.42991.26070.40970.3905UPrs136410227Synonymous variant404550176.71357.36942.84452.655240585702G/ASynonymous variant4058570212.839914.01475.30215.3908CD3640599222G/ASynonymous variant4059922213.147614.23555.49045.3956CD3640614608C/TMissense variant4061460813.225914.27775.52815.3700CD3640614675G/A5 prime UTR variant4061467513.225514.28285.52495.3764CD3640807452A/GMissense variant408074529.27228.55494.14282.0438GNAT3rs110078696Synonymous variant408075769.3588.63514.20662.0567GNAT3rs109903966Synonymous variant408075949.35728.63464.20612.0566GNAT3rs108969608Synonymous variant408139779.42768.75994.21752.1207GNAT3^1^AISC: Number of inseminations per conception; IFLC: days from first to last insemination; NRRC: 56-day non-return rate; ICF: the length in days of the interval from calving to first insemination.

### BTA13 {#Sec6}

Based on the results of the HD array (Figure  [2](#Fig2){ref-type="fig"}: middle) an area from 33,201,457 bp to 34,393,006 bp was analyzed using sequence data (Additional file [4](#MOESM4){ref-type="media"}). The highest --log~10~(P-value) was 19.41 at 33,903,159 bp. Based on the analysis using sequence data two regions were identified (Figure  [2](#Fig2){ref-type="fig"}: right): Region 1: 33,900,143-33,908,994 bp (=8,851 bp) and Region 2: 34,051,815-34,056,728 bp (=4,913 bp). Linkage disequilibrium (LD) analysis showed that these two regions are in almost complete LD (Figure  [3](#Fig3){ref-type="fig"}). BTA13 was screened for candidate genes/polymorphisms in the range of 33 Mb to 34.4 Mb. In total 8 genes were identified in the area from 33.2 Mb to 34.4 Mb in of which 6 were annotated, one was described as uncharacterized protein and one was described as a pseudogene (Table  [3](#Tab3){ref-type="table"}). These genes were not located in Region 1 and Region 2 (Additional file [5](#MOESM5){ref-type="media"}: Figure S1). The SNPs in Region 1 and Region 2 as well as the region in-between were annotated by Variant Effect Predictor as "intergenic" i.e. no candidate genes were identified.Figure 3**Linkage disequilibrium plot of the pairwise D' for the identified regions for the length in days of the interval from calving to first insemination (ICF) on BTA13.**

###### 

**Genes located on BTA13 in the range of 33,000,000 bp to 34,400,000 bp**

  Ensembl gene ID      Description                                          Gene start (bp)   Gene end (bp)   Gene name
  -------------------- ---------------------------------------------------- ----------------- --------------- --------------
  ENSBTAG00000022000   calcium channel, voltage-dependent, beta 2 subunit   33005028          33285264        CACNB2
                       uncharacterized LOC100847770                         33281429          33283899        LOC100847770
  ENSBTAG00000022003   NOP2/Sun domain family, member 6                     33291604          33369610        NSUN6
  ENSBTAG00000001552   enhancer of polycomb homolog 1 (Drosophila)          33451679          33518736        EPC1
  ENSBTAG00000002280   kinesin family member 5B                             33619259          33665129        KIF5B
  ENSBTAG00000012677   Rho GTPase activating protein 12                     33713676          33830942        ARHGAP12
  ENSBTAG00000020053   zinc finger E-box binding homeobox 1                 34061282          34261303        ZEB1
                       deoxyuridine triphosphatase pseudogene               34103073          34103547        LOC100141025

Discussion {#Sec7}
==========

By analyzing the results of the significant QTL presented earlier \[[@CR6]\] using the HD data, it became evident that BTA04 and BTA13 showed the most significant results.

In this study we have re-analyzed two regions which were previously associated to female fertility traits and validated in different breeds of dairy cattle \[[@CR6]\]. A single marker analysis was performed using the HD array and subsequently DNA sequence for BTA04 and BTA13, respectively. The physical location of the SNPs has been used to search for candidate genes on the bovine genome.

BTA04 {#Sec8}
-----

In the QTL region 9 genes were annotated (Table  [1](#Tab1){ref-type="table"}). Based on the --log~10~(P-value) (Table  [2](#Tab2){ref-type="table"}), CD36 is the most obvious candidate gene. CD36 has a missense variant. A missense variant is a sequence variant that changes one or more bases, resulting in a different amino acid sequence but where the length is preserved (<http://www.sequenceontology.org/miso/current_release/term/SO:0001583>). CD36 has been suggested to be involved in numerous potential physiological functions; cell adhesion processes, binds long chain fatty acids and may function in the transport and/or as a regulator of fatty acid transport and act as receptor for thombospondins, THBS1 AND THBS2, mediating their antiangiogenic effects. This in turn involves many physiological processes \[[@CR7]\].

Another interesting candidate on BTA04 based on its suggested biological function is Semaphorin-3C (SEMA3C). SEMA3C binds to members of the plexin family and plays an important role in the regulation of developmental processes. Furthermore it is required for normal cardiovascular development during embryogenesis and plays an important role in axon growth and axon guidance by functioning as attractant for growing axons (<http://www.uniprot.org/uniprot/A7MB70>). However, the SNP markers located in these genes did not reach genome wide significance for the traits AISC, IFLC, and NRRC (Table  [2](#Tab2){ref-type="table"}). At this stage annotation of this region on BTA04 is not sufficient to reach a firm conclusion about the causal polymorphism. Even though some of the SNP markers were located within genes, the SNP markers with the most significant associations to AISC and IFLC were not annotated. Therefore, a better annotation might help to identify the nature of the causal mutation.

BTA13 {#Sec9}
-----

The HD analysis revealed a sharp peak on BTA13 in the range of 33.2 Mb to 34.4 Mb. Searching the chromosome for genes revealed a low number of annotated genes (918 transcripts on the entire chromosome). Eight of these genes were in the 33.2 Mb to 34.4 Mb range. The sequence analysis revealed two regions within the 33.2 Mb to 34.4 Mb range. Even though we have used sequence data it is still a challenge to disentangle whether the peaks harbor the same genetic variation influencing ICF, or are in reality a combination of two separate genetic variations. Therefor the two peaks were analyzed further and were determined to be in high LD. This suggests that it is the same genetic variation causing the two peaks. The SNPs annotated in these two regions were all annotated as intergenic variants. It is possible that BTA13 at this stage is poorly annotated such that the associated polymorphisms are located in as-yet undiscovered genes.

Fertility traits are complex traits as many different biological and physiological factors determine whether a cow is fertile, another challenge is the way the phenotypes are defined, they do not distinguish between different biological functions as they are defined more for a breeding purpose than from a biological point of view. Therefore it is not expected that there is a simple explanation with an obvious candidate gene, it is more likely a network of genes and intragenic variants that explain the variation of these traits. A number of genes and their potential functions in the vicinity of the two peaks are presented in Table  [3](#Tab3){ref-type="table"}.

Selection of markers {#Sec10}
--------------------

The two regions on BTA4 and BTA13 were chosen based on previous studies by the authors \[[@CR6], [@CR8]\]. The SNP markers for BTA04 and BTA13 from the bovine HD array have been tested in the three breeds. Even though the number of sires in the NR and DJ breeds are much lower compared to the NH breed the results suggest that the QTL for female fertility traits detected on both chromosomes are specific for the NH breed. What appeared to be the same hits on several chromosomes in a previous manuscript in several breeds when using the 50 k SNP array, clearly was not when using the 700 k SNP array.

We have been able to narrow down the regions on BTA04 and BTA13 by adding many more markers. However in our analysis we could not with certainty pinpoint the causal polymorphism for the fertility QTL. This might be due to the fact that half of the total genetic variants identified in the whole genome sequencing (WGS) were filtered out for various reasons. All the variants which were not bi-allelic were dropped due to limitations in the imputation software. Therefore, the actual causative polymorphism may be missing from the data analyzed here. Furthermore, there were hundreds of SNPs with very high --log~10~(P-values) due to high LD among themselves. Therefore, it is not possible to pick a few from them based on this analysis.

Effect of recessive lethals on fertility in Dairy Cattle {#Sec11}
--------------------------------------------------------

Recently it has been argued that an increase in recessive lethals play a part in the decline of fertility as elite sires obtain a very large number of offspring in cattle breeding, and the effective population size is around 50 in the Holstein breed \[[@CR9]\]. The development of increased information on the cattle genome has enabled us to detect putative recessive lethal mutations by the absence (or near absence) of homozygous recessive individuals in the population \[[@CR10]\]. Some recessive lethal alleles are known to affect the developing embryo. A recent study in the Nordic Holstein population \[[@CR10]\] has identified a number of recessive lethal haplotypes which appear to act in early pregnancy. The carrier frequency of the recessive lethal alleles is up to 20% in the Nordic Holstein population. Also, VanRaden et al. \[[@CR11]\] and Fritz et al. \[[@CR12]\] have mapped recessive lethals genes by homozygosity mapping in Holstein populations. However in the present study it does not seem to be the same effects that are picked up as the recessive lethals are not located on BTA04 or BTA13.

Annotation of the bovine genome {#Sec12}
-------------------------------

Tools and resources for annotation and gene discovery in the bovine genome are available (e.g. <http://www.ncbi.nlm.nih.gov/assembly/GCF_000003055.4/>; \[[@CR1]\]). Even though many genes have already been annotated on the bovine genome (<http://www.ncbi.nlm.nih.gov/assembly/GCF_000003055.4/>), the known set sites of gene transcription, initiation, termination as well as differential splicing remains incomplete. Therefore, information on genomic structure of organisms which are better annotated like mouse and human, are used as an information source to cover these gaps in the knowledge of the genome structure. However, identifying regulatory elements and non-protein coding regions annotation remains more challenging \[[@CR8], [@CR9]\].

Regulatory elements in the genome {#Sec13}
---------------------------------

It is a challenge to identify all regulatory elements in the genome also those which control gene expression. The question is how to relate genes and their products function. Even though we have access to the cattle sequence which facilitates the comprehensive identification of these transcriptional regulatory factors there is still a long way to go. The expressions of eukaryotic protein coding genes are regulated in different steps which include elements like; transcription initiation and elongation, mRNA processing, transport, translation and stability. Most regulation is believed to occur at the level of transcription \[[@CR13]\].

Directly measured phenotypes and genotypes of cows are necessary to provide a more direct link between phenotype and genotype, it can also enable mapping of non-additive effects which affect fertility traits in cattle \[[@CR14]\]. This was not investigated in the present study but could potentially shed more information regarding fertility in cattle.

Identification and analysis of phenotypic measures that reflect more directly the physiologic background of the reproduction traits could also be helpful in determining the precise physiological background represented by a specific QTL. Analyzing expression data from reproductive organs might help to identify the temporal-specific aspects of gene expression.

Conclusions {#Sec14}
===========

The results from screening BTA04 and BTA13 of NH, NR and DJ using the bovine HD SNP array for female fertility traits suggested that the QTL for female fertility located on these chromosomes were specific for the NH breed. The subsequent screening using the imputed sequence variants for NH narrowed down the QTL on BTA04 into two regions of 288,314 bp and 1.1 Mb respectively, while on BTA13 the QTL region for ICF was narrowed down into two regions of 8,851 bp and 4,913 bp, respectively. A candidate gene search for these QTL regions revealed that a missense mutation in the Platelet glycoprotein 4 (CD36) gene showed the strongest association with fertility traits and therefore is a strong candidate for the QTL on BTA04, whereas the annotated SNPs on BTA13 were all intergenic variants. Fertility traits are complex traits as many different biological and physiological factors determine whether a cow is fertile. Therefore it is not expected that there is a simple explanation with an obvious candidate gene but it is more likely a network of genes and intragenic variants that explain the variation of these traits.

Methods {#Sec15}
=======

Animal ethics statement {#Sec16}
-----------------------

Phenotypic data were obtained from routine records of cattle farms. Semen and blood samples for genotyping were collected in previous studies \[[@CR15], [@CR16]\] using standard procedure for breeding purposes by veterinary or authorized personnel, and all efforts were made to minimize suffering.

No data was collected for the purposes of this study. All DNA data were obtained through the analysis of materials collected as part of routine operation of cattle breeding programs. Farm management and breeding programs in the EU and EEA are subject to the "European Convention for the Protection of Animals kept for Farming Purposes" as implemented in national law.

Animal population {#Sec17}
-----------------

The animal population has been described before in Höglund et al. \[[@CR6]\]. In short, a total of 3,475 Nordic Holstein (NH) sires from Denmark, Sweden and Finland with official breeding values for female fertility traits were used to discover associations. We used 4,998 Nordic Red (NR) animals and 1,225 Danish Jersey (DJ) animals to check for segregation for QTLs in these populations.

Phenotypes {#Sec18}
----------

The traits evaluated included: number of inseminations per conception (AIS), 56-day non-return rate (NRR), days from first to last insemination (IFL), and the length in days of the interval from calving to first insemination (ICF). With exception of ICF, single trait breeding values (STBVs) from the national evaluation were available for both 1^st^ parity animals (heifers, suffixed H) and 2^nd^ and 3^rd^ parity animals (cows, suffixed C). For details regarding the phenotypes recorded and models used in routine breeding value prediction, see <http://www.nordicebv.info>.

Sequencing {#Sec19}
----------

The sequences for the reference population used for imputation of Nordic animals consisted of the whole genome sequence carried out at Aarhus University and in the 1,000 Bull Genome project \[[@CR4]\]. The sequencing of Nordic bulls at Aarhus University, Foulum was done using Illumina sequencers at Beijing Genomics Institute, Shenzhen, China. Sequencing was shotgun paired-end sequencing with a read length of 91 base pairs. Fastq data were converted from Illumina to Sanger quality encoding using a patched version of maq \[[@CR17]\]. They were aligned to the UMD3.1 assembly of the cattle genome \[[@CR18]\] using bwa version 0.6.2 \[[@CR19]\]. They were converted to raw BAM files using samtools \[[@CR20]\]. Quality scores were re-calibrated using the Genome Analysis Toolkit \[[@CR21]\] version 1.6's following the Human 1000 Genome guidelines incorporating information from dbSNP version 133 \[[@CR22]\]. Sequences were realigned around insertion/deletions using the Genome Analysis Toolkit version 1.6. Variants were called using the Genome Analysis Toolkit version 1.6's UnifiedGenotyper. Genomes for the 1,000 Bull Genomes project were sequenced in a number of laboratories and collected in the Department of Primary Industries, Victoria, Australia. Data processing was standardized. Sequences were aligned to the same reference genome as used at Aarhus University using versions of bwa \[[@CR19]\]. Variant calling was done using samtools's mpileup function. Variant Call Files from Aarhus University and the 1,000 Bull Genomes project were combined using the Genome Analysis Toolkit's CombineVariants with precedence given to calls from the Nordic dataset for animals appearing in both datasets.

Imputation HD and Sequence data {#Sec20}
-------------------------------

The sires used in this study were genotyped using the 50 k SNP array as described earlier \[[@CR6]\]. These 50 k SNP typings were the basis for the imputation to HD and sequence level. The imputation of 50 k SNP to the full sequence was done in two steps. First in another study (N.K. Kadri, *pers. comm*.), the 50 k genotypes (46,702 SNPs after quality control) for 12,322 Nordic bulls were imputed to HD genotypes (734,077 SNPs) using the software IMPUTE2 \[[@CR23]\]. The reference population with HD genotypes was available for 2,036 bulls (902 Holstein, 735 Nordic Red and 399 Danish Jersey).

In the second step of imputation, the 12,322 bulls imputed to HD genotypes were further imputed to the full sequence level, using a reference of 242 sequenced dairy bulls (132 Holstein, 42 Jersey, 52 Nordic Red and 16 Brown-Swiss). The 242 dairy cattle sequences originated from a combination of sequences processed at Aarhus University and sequences from the 1,000 Bull Genomes dataset. Only polymorphisms identified both in the Aarhus University dataset and the 1,000 Bull Genomes dataset were included. For positions containing both a SNP and an INDEL, the INDEL was deleted. SNPs at positions with disagreements between alleles from sequence and HD data were deleted. The reference data was pre-phased with BEAGLE v3.3.2 \[[@CR24]\] and all markers with an r^2^ value below 0.9 were removed. This left a total of 16,374,063 sequence markers and 629,028 HD markers for chromosomes 1--29. Chromosomes were divided into chunks of about 20,000 consecutive sequence markers with an overlap of 500 markers to minimize imputation error at ends of the chunks. Imputation was done using BEAGLE v3.3.2 \[[@CR24]\].

Statistical method for association analysis {#Sec21}
-------------------------------------------

The association between the SNP and the phenotype was assessed by a single-locus regression analysis for each SNP separately, using a linear mixed model \[[@CR25]\]. The model was as follows:

where **y** is the vector of phenotypes (de-regressed EBV), 1 is a vector of ones with length equal to the number of observations, *μ* is the general fixed mean, ***m*** is a vector of genotypic indicators (-1, 0, or 1) associating records to the marker effect, *g* is a scalar value of the additive fixed allele substitution effect of the SNP, ***Z*** is an incidence matrix relating phenotypes to the corresponding random polygenic effect, ***u*** is a vector of the random polygenic effect with a multivariate normal distribution N where ***A*** is the additive relationship matrix and is the polygenic variance, and *e* is a vector of random environmental deviates with a normal distribution ) where is the error variance and ***W*** is the diagonal matrix containing weights of the de-regressed estimated breeding values. The weight for the i^th^ animal was estimated as *w*~*i*~ = *r*^2^/(1 - *r*^2^), where r^2^ was the reliability of the de-regressed EBV of the i^th^ animal and *r*^2^ \> 0.98 was set to 0.98 to avoid very large weights for sires with very large number of progeny records. The model was fitted by REML using the software DMU \[[@CR26]\]. The standard error of the fixed effect estimates was obtained from DMU. Testing for the presence of an effect of a marker was done using a Wald test against a null hypothesis of H~0~: g = 0.

Significance levels {#Sec22}
-------------------

The significance level using genome wide bonferoni correction was 0.05/38,545 \~ 1.3 × 10^-6^ in the 50 k data (previous study). It was 0.05/734,077 \~ 6.8 × 10^-5^ SNPs in the 700 k data and 5.5 × 10^-8^ in the WGS data. The significance level using region wise bonferoni correction was 0.05/3716 \~ 1.3 × 10^-6^on BTA04 and on BTA13 0.05/3944 \~ 1.3 × 10^-4^.

Variant annotation {#Sec23}
------------------

Variants were annotated using Variant Effects Predictor version 2.8 \[[@CR27]\]. The underlying databases correspond to ENSEMBL databases version 70.

Linkage disequilibrium {#Sec24}
----------------------

To distinguish whether the two peaks in the sequence analysis on BTA13 were due to two separate mutations or one mutation the linkage disequilibrium of BTA13 was examined for these peaks. The top 100 markers in each peak were examined using HAPLOVIEW \[[@CR28]\].

Availability of supporting data {#Sec25}
-------------------------------

All DNA sequences used were taken from a publicly available assembly. The assembly is available for download (<ftp://ftp.ensembl.org/pub/release-73/fasta/bos_taurus/dna>). All variations used in the mapping study have been submitted by the 1000 Bull Genomes project for inclusion in dbSNP (<http://www.ncbi.nlm.nih.gov/SNP>). All annotations were obtained from a publicly available source (<http://www.ensembl.org>) downloadable including through Variant Effect Predictor (<http://www.ensembl.org/info/docs/tools/vep/script/index.html>).

Samples were collected in the context of previous studies \[[@CR15], [@CR16]\]. The pool of HD chip genotypes used as reference for imputations was created from a number of sources \[[@CR29]\]. Contributions were obtained through exchange agreements with other research institutions in Sweden, Finland, Germany, France, The Netherlands and Spain. Other genotypings were done commercially on behalf of Aarhus University. Access to HD chip typing data can be granted by each exchange partner individually.

Electronic supplementary material
=================================

 {#Sec26}

###### 

Additional file 1: **Results of BTA04 in the Danish Holstein breed for number of inseminations per conception (AIS), 56-day non-return rate (NRR), days from first to last insemination (IFL), and the length in days of the interval from calving to first insemination (ICF) using the bovine HD SNP array.** (TXT 148 KB)

###### 

Additional file 2: **Results of BTA13 in the Danish Holstein breed for the length in days of the interval from calving to first insemination (ICF) using the bovine HD SNP array.** (TXT 242 KB)

###### 

Additional file 3: **Results of the single-marker analysis of BTA04 in the Danish Holstein breed for number of inseminations per conception (AIS), 56-day non-return rate (NRR), days from first to last insemination (IFL), and the length in days of the interval from calving to first insemination (ICF) using sequence data.** Only significant P values are mentioned for each trait. NA: marker not significant. (TXT 70 KB)

###### 

Additional file 4: **Results of the single-marker analysis of BTA13 in the Danish Holstein breed for the length in days of the interval from calving to first insemination (ICF) using sequence data.** Only region-wide significant P-values are presented. (TXT 91 KB)

###### 

Additional file 5: **Manhattan plots of the GWAS results using imputed DNA sequence on BTA13 for ICF in the range of 33.2 Mb to 34.4 Mb.** The y axes show --log~10~(p-value) of single-marker associations tests. The x axes show marker positions in base pairs (bp). The blue dots indicate test results using sequence data. The red triangles on the x-axis indicate the position of the annotated genes in the region. (PDF 26 KB)

**Competing interests**

Authors declare no competing interests.

**Authors' contributions**

Conceived and designed the experiment: JKH, GS, BG, MSL. Analyzed the data: JKH, GS, RFB, BB, BG. Contributed reagents/materials/analysis tools: MSL, BG, GS. Wrote the paper: JKH. All authors contributed to the discussion of the results, read and approved the final manuscript.

We are grateful to the Danish Cattle Federation/NAV for providing the phenotypic data used in this study. This work was supported by a grant (No. 3405-10-0137) funded jointly by the Green Development and Demonstration Program of the Danish Ministry of Food, Agriculture and Fisheries, The Milk Levy Fund, Viking Genetics, and Nordic Cattle Genetic Evaluation. Semen samples were kindly provided by the Swedish Farmers Foundation for Agricultural Research in conjunction with Viking Genetics. The 1,000 Genomes Project is kindly acknowledged for sharing data to impute the genome sequence.
