Prioritizing candidate genes for fertility in dairy cows using gene-based analysis, functional annotation and differential gene expression

Association analysis for fertility

Using a previously published approach [12], we performed GWAS and found seven independent QTL (−log₁₀(P) > 8.5) across Bos taurus autosomes (BTA) 1, 6, 10, 13, 17 and 24 (Table 1 and Fig. 1). In total, 5806 SNPs exceeded the genome-wide significance level. The most significant signal was on BTA13 with the lead SNP located at 32,852,133 bp (rs210238678). This QTL was close to previous reported QTL in Nordic Holstein cattle (at BTA13: 33,903,159 bp) [13]. The second strongest association signal was on BTA6 (95,867,927 bp, rs41567777) with GK2 as the closest gene. GK2 belongs to GO terms “carbohydrate metabolic process” and “galactose metabolic process”. On chromosome 1, we identified two QTL. The first had its lead SNP at 69,742,415 bp (rs208311936) which is an intronic variant in the UMPS gene. The second had its lead SNP at 140,785,028 bp (rs385628476), an intronic variant in the BRWD1 gene. UMPS belong to GO terms “female pregnancy” (The set of physiological processes that allow an embryo or foetus to develop within the body of a female animal that covers the time from fertilization of a female ovum by a male spermatozoon until birth) which make it a good candidate for fertility. For BRWD1, there is no documented biological connection between the available functional annotations with the gene and fertility. On BTA 24, the lead SNP was located at 29,556,826 bp (rs380439408). The closest gene is a novel gene, ENSBTAG00000044212. On BTA10, we identified SNP at 68,534,665 bp (rs211204488) as the lead SNP. The annotated element closest to this lead SNP is U7 snRNA and SLC39A12 is the closest gene. The lead SNP on BTA17 was located at 71,393,345 bp (rs110812733). This lead SNP is located close to LIF gene which is involved in many biological processes including “maternal process involved in female pregnancy”, thereby making it a good candidate gene for female fertility.

BTA	BP	effect	-log10(P)	Interval	Gene (distance in bp)	Annotation
1	69,742,415	3.68	11.03	68,745,058~69,992,433	UMPS	intron
1a	140,785,028	1.91	9.53	139,838,183~141,035,632	BRWD1	intron
6	95,867,927	2.15	11.80	95,377,650~96,118,034	GK2 (14,991)	intergenic
10	68,534,665	2.08	9.16	67,917,166~68,784,933	U7 snRNA (29,099)	intergenic
13	32,852,133	−3.09	16.48	31,866,274~33,102,422	SLC39A12 (100,247)	intergenic
17	71,393,345	1.83	8.54	70,501,607~71,656,270	LIF (20,510)	intergenic
24	29,556,826	2.32	10.18	28,794,256~29,807,379	ENSBTAG00000044212 (65,604)	intergenic

Base positions are given as positions in UMD 3.1.1 [49]. Lead SNPs identified in the second round are marked by^a. The method to define QTL intervals was described in Methods

Variants annotation of the variants in LD with lead SNP

Lead SNP may or may not be the causal variants [12]. To refine our search, we annotated all significant SNP in LD with the lead SNP (r² > 0.2). This included 2283 variants. Half of these variants (1153, Fig. 2A) were intron variants, followed by intergenic variants (1007, Fig. 2A). Among variants located in coding sequences, most (13, Fig. 2B) were synonymous variants, followed by missense variants (6, Fig. 2B). Some variants had predicted functional consequences: KALRN contained a splice acceptor variant (rs109819533, BTA1:69652189). Two missense variants (BTA13:32679690, rs464438710 and BTA1:69673871, rs209885271) were reported as “deleterious” by SIFT [20]. They were located in the genes SLC39A12 and KALRN.

Overlap of identified QTL in other cattle populations

Cow fertility in dairy cattle is a complex trait with many biological pathway involved and has a very low heritability. Non-genetic factors like body condition, farm management, feeding also affect fertility in cattle. Therefore, results show little agreement between studies [7, 10, 11, 13]. Besides, WGS has not been used extensively to study the genetic basis of variation in cow fertility in dairy cattle [21, 22]. The current study has identified seven QTL located on six different autosomes (Table 1). The QTL with its lead SNP at BTA17: 71,393,345 (rs110812733) was close to a QTL (70–73 Mb) for fertility reported in Brown Swiss cattle [21]. The lead SNPs for the QTL on BTA10, BTA13, and BTA24 were close to significant SNPs for female fertility traits previously reported in Danish and Swedish Holstein cattle [23]. The QTL on BTA1 and BTA6 are novel findings. Previous reported QTL were located outside QTL intervals obtained in the present study.

Search for candidate genes

We proposed candidate genes for seven QTL identified in this study. In this list, some candidate genes are supported by literature reports of relevant biological functions related to fertility in other species. In mouse, it has been shown that BRWD1 is essential for female fertility by epigenetic control of meiotic chromosome stability [24]. KLF6 is a member of the Krüppel-like factor, and have a similar function to other Krüppel-like factors. These are indispensable for normal implantation [25]. Leukemia inhibitory factor (LIF) is a cytokine. It is required for blastocyst implantation in mice. Several studies have shown that LIF is important for the establishment of pregnancy [26]. Previous studies have shown that SPIRE1 encodes a protein that drives two critical steps in asymmetric oocyte division in mice [27].

A similar strategy to prioritize candidate genes have helped identify several candidate genes for mastitis resistance [28]. In the present study, we faced two challenges like in our previous study [28]. The first challenge is choosing appropriate sources of annotations to prioritize candidate genes. Even though many biological processes are involved in fertility, we could only use the direct evidence that would affect conception and carrying the fetus since we cannot be sure what other pathway also involved. The second challenge is the long list of candidate genes provided by gene-based analysis. In the present study, we incorporated LD information in two ways: 1) using the linreg model from MAGMA [15]; 2) only considering genes close to significant SNPs in LD with lead SNP. We observed the cross-validation of DEGs from different studies or different animals in one study could provide better support. In this study, we utilized multiple RNA-seq datasets to validate our candidate gene list. One of the candidate genes, FRAS1, was found in two dataset with DEGs. In both datasets, differences in expression were in the same direction when comparing pregnant with non-pregnant animals. Three other candidate genes, ITGB5, ADCY5, and SEMA5B were DE in the third dataset in two different fertility-classified animals in the same directions when comparing pregnant with non-pregnant animals [19]. This indicated that, with the emergence of more functional study datasets permits the identification of better candidate genes with higher confidence based on biological support. Moreover, the successful combination of information from different sources to prioritize candidate genes suggested we would gain higher power with the improvement of the understanding of functional genomics.

Limitation of the present study

The GWAS for cow fertility using imputed WGS marker could identify only a small number of QTL in our study. Both limited sample size and small contribution of individual QTL to the total phenotypic variance for fertility have contributed to low power of detection. Even though the genome-wide markers can explain a sizable proportion of the trait heritability, the genome-wide significant SNPs often explain only a small proportion of the trait heritability. Meta-analysis of GWAS summary statistics from several populations can improve the power [29] and map the QTL location precisely [29]. We used the fertility index as phenotype in our study. The breeding value for fertility index is calculated based on several component traits, which bring many biological processes together. GWAS for individual component traits of fertility index will help to interpret the results biologically. The imputed WGS marker set helped us to precisely identify the location of QTL intervals. The fertility index included in the Nordic breeding goal and the identified whole genome sequence variants, if included in the EuroGenomics custom SNP chip [30], can help in routine genomic prediction.

Phenotype and genotype data

Phenotypic records of fertility in this study were obtained from the Nordic Cattle Genetic Evaluation database (NAV, http://www.nordicebv.info/). The phenotype records used for association were de-regressed breeding values [31, 32] from the routine genetic evaluation by NAV and were available for 5038 progeny tested Holstein bulls. The fertility index includes breeding values for interval from first to last insemination (heifers and cows), interval from calving to first insemination (cows) and number of inseminations (heifers and cows).

We used two-step method previously described by Iso-Touru et al. [33] and Wu et al. [34] to impute WGS data. All bulls were genotyped with the Illumina BovineSNP50 BeadChip (54 k) ver. 1 or 2 (Illumina, San Diego, CA, USA). In the first step, we imputed 54 k genotypes to high-density (HD) by IMPUTE2 v2.3.1 [35]. The reference population for imputation included 1222 Holsteins, 1326 Nordic Red Dairy Cattle, and 835 Danish Jerseys genotyped by Illumina BovineHD BeadChip. In the second step, these imputed HD genotypes were imputed to WGS by Minimac2 [36] with a multi-breed reference of 1228 animals from Run4 of the 1000 Bull Genomes Project [37] and additional data from Aarhus University (80 individuals, including 23 Holsteins, 30 Nordic Red Dairy Cattle, and 27 Danish Jersey) [38]. In this step of imputation, we performed in 5-Mb chunks with a buffer region of 0.25 Mb on either side. At the end, 22,751,039 bi-allelic variants were present in the imputed sequence data. After excluding SNPs with a minor allele frequency below 0.5% or with a large deviation from Hardy–Weinberg proportions (P < 1.0^− 6), 15,551,021 SNPs on 29 autosomes in Nordic Holstein cattle were retained for association analyses. The detailed of this WGS dataset was published previously [34].

Methodology of multiple QTL detection and estimation of genetic variants explained by QTL

We used our previous proposed method [12] to perform the association analysis. In brief, we literately ran GWAS using GCTA [39] for each chromosome by fitting the dose of lead SNP from previous literation as covariates. To reduce the false positive, we defined the lead SNP as the significant ones (experiment-wise 0.05 type I error rate after Bonferroni correction for 15,551,021 simultaneous tests corresponds to a threshold of –log₁₀(P) ≈ 8.5) with the largest –log₁₀(P) value in each literation and significant in first literation. To minimize the impact of random errors and imputation inaccuracy, we also checked whether the lead SNP is solo SNP. The solo SNPs are SNPs with no other significant SNP within a 1 Mb region. These SNPs were skipped in our analysis. The lead SNP in each round were collected to build a lead SNP list. The boundaries of each QTL region were defined as followed. If the SNP –log₁₀(P) value of the flanking 1 Mb region around the lead SNP decreased by more than three units compared to the value of the lead SNP and the region was larger than 0.25 Mb, then we set this SNP as the boundary; otherwise, we set ±0.25 Mb from the lead SNP as the QTL boundary.

LD calculation and variant annotation

The procedure for LD calculation and variant annotation were described in our previous study [28]. We used PLINK [40] to calculate the pairwise r² between the lead SNP and all other SNPs on the same chromosome. All significant SNPs with r² > 0.2 with the lead SNP were extracted for variants annotation. These SNPs were annotated by VEP (version 92) [41].

Candidate genes identification and confirmation with RNA-seq data

We included nearest genes and gene-based analysis as a list of potential candidate genes list. For the nearest genes, we used bedtools [42] closest function to find the nearest genes (or function annotated feature) to the lead SNPs. For gene-based analysis, we used MAGMA [15]. In order to run MAGMA [15] for cattle, we should provide customer gene annotation file and reference population to MAGMA. For cattle genome annotation, we downloaded the gene information file from Ensembl gene build 92 [43] and converted them to MAGMA style location file. For reference population, 455 Holstein animals from the 1000 Bull Genome Project (Run 6) [37, 38] were used for this purpose. We performed MAGMA gene analysis with the GWAS result using the model linreg. The number of genes (including 5′- and 3′-UTRs) with at least one SNP was 20,356; thus, the P-value threshold for genome-wide significance for the gene-analysis was 2.46 × 10^− 6. Significant SNPs in LD with lead SNPs were used to extract the closest or overlapping genes by bedtools [42] closest function from the list of significant genes from MAGMA [15]. These two sources of potential candidate genes were analyzed using DAVID [44] to retrieve the KEGG pathway annotations [45]. GO terms [46] associated with these genes were retrieved from Uniprot [47]. The MPD [16] was searched for mutations in these genes with known phenotypic effects related to fertility. At the meanwhile, the list of DEGs from three previous studies [17–19] were used to provide biological evidence for genes. All the potential candidate genes with any biological support (GO, KEGG, MPD and DE) were proposed as final candidate genes.