Using expression data to fine map QTL associated with fertility in dairy cattle

Abstract

Background: Female fertility is an important trait in dairy cattle. Identifying putative causal variants associated with fertility may help to improve the accuracy of genomic prediction of fertility. Combining expression data (eQTL) of genes, exons, gene splicing and allele specific expression is a promising approach to fine map QTL to get closer to the causal mutations. Another approach is to identify genomic differences between cows selected for high and low fertility and a selection experiment in New Zealand has created exactly this resource. Our objective was to combine multiple types of expression data, fertility traits and allele frequency in high- (POS) and low-fertility (NEG) cows with a genome-wide association study (GWAS) on calving interval in Australian cows to fine-map QTL associated with fertility in both Australia and New Zealand dairy cattle populations.

Results: Variants that were significantly associated with calving interval (CI) were strongly enriched for variants associated with gene, exon, gene splicing and allele-specific expression, indicating that there is substantial overlap between QTL associated with CI and eQTL. We identified 671 genes with significant differential expression between POS and NEG cows, with the largest fold change detected for the CCDC196 gene on chromosome 10. Our results provide numerous candidate genes associated with female fertility in dairy cattle, including GYS2 and TIGAR on chromosome 5 and SYT3 and HSD17B14 on chromosome 18. Multiple QTL regions were located in regions with large numbers of copy number variants (CNV). To identify the causal mutations for these variants, long read sequencing may be useful.

Conclusions: Variants that were significantly associated with CI were highly enriched for eQTL. We detected 671 genes that were differentially expressed between POS and NEG cows. Several QTL detected for CI overlapped with eQTL, providing candidate genes for fertility in dairy cattle.

Fig. 1

Overview of data used for analyses

Fig. 2

Volcano plot of differentially expressed genes. ‘Down’ indicates variants with a FDR ≤ 0.05 and a log2 fold change (logFC) smaller than 1, ‘up’ indicates variants with a FDR ≤ 0.05 and a log2 fold change (logFC) larger than 1, ‘no’ indicates all other variants

Fig. 3

Overlap between QTL for calving interval and eQTL. Y-axis = − log10(p) for the meta-analysis of calving interval (CI), red indicates variants that are significant (p ≤ 10^–6) for both CI and (a) gene expression (gene), (b) exon expression (exon), (c) gene splicing (splicing) and (d) allele specific expression (ASE)

Fig. 4

Overlap between a QTL for calving interval on chromosome 5 and GYS2 splicing. a y-axis = − log10(p) for the meta-analysis of calving interval (CI), and (b) y-axis = − log10(p) for the intron in GYS2 located between 88,636,100 and 88,647,589 bp, red indicates variants that are significant (p ≤ 10^–6) for both CI and gene splicing

Fig. 5

Overlap between a QTL for calving interval on chromosome 5 and TIGAR expression. a y-axis = − log10(p) for the meta-analysis of calving interval (CI), and (b) y-axis = − log10(p) for TIGAR expression, red indicates variants that are significant (p ≤ 10^–6) for both CI and gene expression

Fig. 6

Overlap between a QTL for calving interval on chromosome 6 and multiple types of eQTL. a y-axis = − log10(p) for the meta-analysis of calving interval (CI), (b) y-axis = − log10(p) for GC expression, (c) y-axis = − log10(p) for exon expression of exon of GC located between 86,968,870 and 86,968,921 bp, and (d) y-axis = − log10(p) for allele expression of tSNP located at 86,989,953 bp, red indicates variants that are significant (p ≤ 10^–6) for both CI and expression phenotype

Fig. 7

Overlap between a QTL for calving interval on chromosome 15 and multiple types of eQTL. a y-axis = − log10(p) for the meta-analysis of calving interval (CI), b y-axis = − log10(p) for OR4B1GP expression, c y-axis = − log10(p) for exon expression of exon of ARHGAP1 located between 76,627,377 and 76,627,508 bp, d y-axis = − log10(p) for the splice region between 77,959,141 and 77,997,756 bp, and (e) y-axis = − log10(p) for allele expression of tSNP located at 76,623,841 bp, red indicates variants that are significant (p ≤ 10^–6) for both CI and expression phenotype

Fig. 8

Overlap between a QTL for calving interval on chromosome 18 and multiple types of eQTL. a y-axis = − log10(p) for the meta-analysis of calving interval (CI), b y-axis = − log10(p) for SYT3 expression, c y-axis = − log10(p) for exon expression of exon of ENSBTAG00000038903 located between 58,342,620 and 58,343,553 bp, and (d) y-axis = − log10(p) for the splice region between 55,443,019 and 55,445,791 bp in in HSD17B14, red indicates variants that are significant (p ≤ 10^–6) for both CI and expression phenotype. Expression associations were only tested for variants within 1 Mb of the gene/exon/splice region; the gaps on the gene/exon/splice graphs fall outside these boundaries