A Post-GWAS Functional Analysis Confirming Effects of Three BTA13 Genes CACNB2, SLC39A12, and ZEB1 on Dairy Cattle Reproduction

In our previous GWAS of Chinese and Nordic dairy cattle, genes CACNB2, SLC39A12, and ZEB1 locating on BTA 13 were suggested as candidate genes for reproduction. In this study, validation of these associations was performed in an independent population with records of nine reproductive traits. More importantly, functions of these genes in the reproductive process were verified by employing the expression data of ovarian follicles. The potential variants within the three genes were firstly detected in 68 Chinese Holstein bulls, and then screened in 1,588 Chinese Holstein cows using the KASP (Kompetitive allele-specific PCR) method. There were nine variants with polymorphisms in CACNB2, five in SLC39A12, and four in ZEB1, respectively, of which one SNP was in the upstream regulatory region, two in exon region, four in downstream regulatory region, and 11 SNPs in intronic regions. Amongst the 18 variants, g.33267056T/G in CACNB2 explained the largest phenotypic variance for age at first calving (0.011%), interval from first to last insemination (0.004%), and calving ease (0.002%), while g.32751518G/A in SLC39A12 contributed the most to stillbirth in heifers (0.038%). Two haplotype blocks were constructed for CACNB2 while one each for SLC39A12 and ZEB1, which were significantly associated with five reproductive traits, including age at the first service, age at the first calving, calving ease in heifers and cows, and the interval from calving to the first insemination. We then studied the profile of gene expression in granulosa cells isolated from four developmental stages of ovarian follicles from eight dairy cows. All three genes were differentially expressed between ovarian follicles with different sizes (p < 0.05), indicating their potential roles in the reproductive process of dairy cows. This study successfully demonstrated the associations of three BTA 13 genes CACNB2, SLC39A12, and ZEB1 with reproduction and further examined their expression levels in ovarian follicles directly. These findings can be beneficial for the ongoing genomic selection program for reproductive traits which have long been considered as traits that are difficult to achieve genetic improvement due to the lack of efficient genetic markers.

1 Introduction

Poor reproductive performances of dairy heifers and cows can largely reduce the efficiency and overall profitability of dairy herds due to extra costs on herd management, veterinary care, and involuntary culling (Dobson et al., 2007; Guo et al., 2014; Liu et al., 2017a). This makes reproduction among the most important functional traits in dairy industry globally. However, due to the highly polygenic character and its negative genetic correlations with intensively selected milk production traits, genetic improvement of reproductive performance is relatively slow (Dobson et al., 2007; Berry et al., 2014). By integrating causal mutations or variants that are in tight LD (linkage disequilibrium) with causal mutations into genomic prediction, the prediction reliability of breeding values can be improved (Schultz et al., 2020), and therefore, accelerating the dairy cattle genetic improvement program.

Genome-wide association study (GWAS) has been served as the primary strategy to detect genes associated with complex traits in the past 15 years. In our previous GWAS of Chinese and Nordic Holsteins (Liu et al., 2017b), we were able to detect QTLs (quantitative trait loci) for reproductive traits on BTA 13 and concluded genes CACNB2 (calcium voltage-gated channel auxiliary subunit beta 2, for interval from calving to the first insemination (ICF) and interval from first to last insemination in cows (IFLc)), SLC39A12 (solute carrier family 39 member 12, for IFLc), and ZEB1 (zinc finger E-box binding homeobox 1, for ICF and IFLc) as candidate genes. The gene CACNB2 encodes the β subunit of voltage-gated calcium channels (Ca+), which is a ubiquitous secondary messenger that regulates most cellular processes and also a prerequisite for gonadotropin-releasing hormone receptors (GnRH)-induced activation of extracellular signal-regulated kinase that medicates the secretion of FSH (Clapham, 2007; Parekh, 2011). Another GWAS in Japanese Holsteins also identified CACNB2 as a candidate gene for conception rate and addressed its role in releasing follicle-stimulating hormone from the anterior pituitary gland (Sugimoto et al., 2013). However, the causative SNPs within CACNB2 for reproduction is still not fine confirmed; moreover, in addition to conception rate, its effects on other reproductive traits remain to be investigated. The gene SLC39A12 is a protein-coding gene and involves in the transport of zinc into the cytoplasm. The protein coded by SLC39A12 belongs to the subgroup of proteins that shows the structural characteristics of zinc transporters and serves as an important component of the internal homeostasis of the cellular zinc (Taylor and Nicholson, 2003; Lisle et al., 2013). The gene ZEB1, encoding zinc finger and homologous domain transcription factors, is one of the key transition factors of epithelial mesenchymal transition in epithelial membrane cells (Bilyk et al., 2017). However, there is a lack of functional analyses to confirm the effect of gene SLC39A12 and ZEB1 on reproduction in cattle.

Our previous GWAS in Chinese and Nordic Holsteins used the 54K SNP array that was designed to cover common SNPs in main cattle breeds to achieve acceptable performances of genomic prediction (Matukumalli et al., 2009). The variants with the most significant effect reported in GWAS are therefore not necessary to be causal for having a direct impact on reproduction, given the limitation in genotyping resolution. Although the use of multiple Holstein populations can narrow down QTL intervals to some extent, the QTL reported in GWAS was still extensively large because LD structures are generally consisted within the Holstein breed (Zhou et al., 2013). In addition, there could be multiple independent loci associated with reproductive traits within a gene. Although some significant associations (e.g., CACNB2) were able to be replicated in populations other than Chinese and Nordic Holsteins, the post-GWAS analyses are still necessary to narrow down QTLs and even fine-map causative mutations.

Post-GWAS screening and whole gene sequencing are capable to increase the resolution of variant associations discovered by GWAS (Visscher et al., 2017). One widely used strategy for post-GWAS is to highlight the role of coding variants (e.g., missense variants) via exome sequencing (Zhao et al., 2015). Additionally, genetic variations in non-coding regions (e.g., intron) can be missed by GWAS when using commercial SNP arrays although they could still explain part of the phenotypic variations due to the functional relevance to causal variants (Cooper, 2010; Lomelin et al., 2010). The primary objective of post-GWAS analysis is to better understand the physiological pathway of the target trait represented by the candidate genes at the molecular level (Emig and Albrecht, 2011). The mRNA expression analyses have already been used to validate GWAS findings for several complex traits in dairy cattle such as milk production traits (Han et al., 2019), mastitis (Mi et al., 2021), and heat stress resistance (Fang et al., 2021). However, post-GWAS studies on reproduction traits in dairy cattle is relatively scarce. One main reason could be the difficulty to collect relevant tissue samples (Struck et al., 2018). For studying reproductive performances, ovarian follicles can be perfect samples due to its indispensable roles in reproduction process. Furthermore, a strong genotyping cohort and powerful population data used in the current study is suggestive of reliable results (Visscher et al., 2017).

Collectively, the objectives of this study were to explore the effects of three BTA 13 candidate genes CACNB2, SLC39A12, and ZEB1 on dairy cattle reproduction, by employing haplotype analyses and gene expression data. The findings from this study could positively contribute to the genetic improvement of reproduction performances in dairy cattle.

2 Methods

In this study, to explore the roles of three GWAS candidate genes (CACNB2, SLC39A12, and ZEB1) on dairy cattle reproduction, we first screened for polymorphisms and performed genotype- and haplotype-based association analyses in an independent Holstein population, and then studied the profile of gene expression in granulosa cells isolated from 8 cows. The designed technology route of this study is shown in Figure 1.

FIGURE 1

FIGURE 1. The designed technology route of post-GWAS functional analysis for gene CACNB2, SLC39A12 and ZEB1 in Holstein cattle.

2.1 Genotype and Haplotype-Based Association Analysis

2.1.1 Screen for Polymorphisms in 68 Bulls

Premier 5.0 and Oligo 7.0 (Sangon Biotech Co., Ltd. Shanghai) were used to design primes for all exons and flanking regions (within 2 kb distance of 5′ and 3’ end) in three GWAS candidate genes. In total, 50 pairs of primers (Supplementary Table S1) were designed, including 16 for CACNB2, 18 for SLC39A12, and 16 for ZEB1. These primers were then synthesized at Sangon Biotech (Shanghai, China).

The DNA samples were extracted from frozen semen of 68 Holstein bulls (Beijing Dairy Cattle Center, Beijing, China) using a standard phenol-chloroform method, and then, DNA concentration were diluted to 200 ng/μl by adding TB buffer. We randomly mixed DNA samples into three pools (one pool contained samples from 28 bulls and two pools contained 20 samples each) and then used pooled samples for the subsequent PCR amplification [ABI3730XL DNA analyzer (Applied Biosystems, Foster, CA, United States)]. Sequencing data were aligned to the reference genome (UMD 3.1) using Ensemble (https://asia.ensembl.org/) and examined for potential polymorphism using Chromas (http://technelysium.com.au/wp/chromas/) and DNAman (https://www.lynnon.com/dnaman.html/). In total, 37 variants were detected from genes CACNB2, SLC39A12, and ZEB1.

2.1.2 Genotype and Phenotype of 1,588 Cows

We considered 1,588 Holstein cows from eight dairy farms (Beijing Sunlon Livestock Development Co. Ltd, Beijing, China), of which complete pedigree and reproductive information have been recorded. All cows were elaborately selected to be independent of the GWAS discovery population (Liu et al., 2017b). DNA samples of cows were extracted from the whole blood samples using TIANamp Blood DNA Kit (TIANgen, Beijing, China). To shortlist variants included for KASP (Kompetitive allele-specific PCR) genotyping, a pilot analysis was performed in 46 (out of 1,588) randomly selected samples. Only one variant was kept if the distance between two variants was less than 200 bp. Ultimately, 18 variants (Table 1) were remained and used to genotype all 1,588 cows with the KASP method. A chi-square test was used to determine whether allelic frequencies of any variant were deviated from Hardy-Weinberg equilibrium. Bonferroni correction was applied to control for familywise false positive.

TABLE 1

TABLE 1. The detailed information of 18 SNPs within the gene CACNB2, SLC39A12 and ZEB1 in 1,588 Holstein cows.

We explored the relationship of 18 variants with nine reproductive traits which covered the broad mechanisms of cattle reproduction (Ettema and Santos, 2004; Nayeri et al., 2017), including age at the first service (AFS) and calving (AFC), the interval from calving to the first insemination (ICF), the interval from the first to last insemination in heifers (IFL_H) and cows (IFL_C), stillbirth in heifers (SB_H) and cows (SB_C), and calving ease in heifers (CE_H) and cows (CE_C). Estimated breeding values (EBVs), derived from single-trait animal models (Liu et al., 2017a) by national routine genetic evaluation (Independent Innovation League of Dairy Breeding, Beijing, China), were used as the response variable for association analyses. Descriptive statistics of EBVs for each reproductive trait were presented in Supplementary Table S2.

2.1.3 Construct Haplotypes Based on LD Structures

For each candidate gene, haplotype blocks were constructed based on LD structures of identified SNPs by using the Haploview4.0 software (Gabriel et al., 2002; Barrett et al., 2005). Haplotypes with relative frequency less than 5% were merged into one group (Table 2).

TABLE 2

TABLE 2. Association of haplotype blocks within gene CACNB2, SLC39A12 and ZEB1 with reproductive traits in Holstein cows (Least Square Mean ± Standard Error).

2.1.4 Association Analysis

Association analyses between every identified SNP or haplotype block and each reproductive trait were performed by using the GLM procedure of SAS 9.2. For both SNP- and haplotype block-based analyses, Bonferroni correction was used to control for false positives resulting from multiple testing. The model used for association analysis of each reproductive trait is as follows:

$\text{y}=\text{µ}+\text{G}+\text{e}$

Where y is individual’s EBV; µ is the overall mean; G is the fixed effect of either genotype or haplotype block; and e is the random residual effect. The proportion of the phenotypic variance $\left({\sigma }_P^2\right)$ explained by each identified SNP was calculated as $\left(\frac{{\sigma }_g^2}{{\sigma }_P^2}\right)$, where ${\sigma }_g^2$ is the variance explained by SNP (Hill and Mackay, 2004) and ${\sigma }_P^2$ is the phenotypic variance (Supplementary Table S2).

2.2 Gene Expression Assay of CACNB2, SLC39A12, and ZEB1

2.2.1 Ovary Collection, Follicle Selection, and Granulosa Cell Isolation

Bovine ovaries of eight dairy cows were collected at an abattoir and put into thermally insulated bottles (28–30°C) containing sterile physiological saline (with 100 U/ml Penicillin and 0.1 mg/ml Streptomycin) immediately. In less than 2 hours, ovaries were transported to the laboratory. After washing three times with warm (37.5°C) 0.9% NaCl solution and rinsing in 70% warm ethanol for 30 s, ovaries were washed with Dulbecco’s Phosphate-Buffered Saline (DPBS) for three times. Healthy ovarian follicles with amber-colored fluid were kept and classified into four developmental stages according to their diameters (d), including stage 1 (d ≤ 3 mm), 2 (3 mm < d ≤ 7 mm), 3 (7 mm < d ≤ 10 mm), and 4 (d > 10 mm), as shown by Figure 4A. The follicular fluid contained both cumulus-oocyte complexes (COCs) and granulosa cells (GCs). We used a filter (diameter = 70 μm) to filter out COCs and then centrifuged the filtrates (contained GCs) at 1,500 × g for 5 min. The supernatant of the follicular fluid was discarded by aspiration. The GC cells were washed three times in phosphate-buffered saline (pH 7.4) and placed in DMEM/F-12 (Gibco, Life Technologies Inc., Grand Island, NY, United States) which contained 1% penicillin-streptomycin and 10% fetal bovine serum (FBS, Gibco, Life Technologies Inc., Grand Island, NY, United States).

2.2.2 Quantitative Reverse Transcription PCR

The RT-qPCR was conducted to capture the relative expression level of genes CACNB2, SLC39A12, and ZEB1 in GCs collected from ovarian follicles at different developmental stages. Total RNA was extracted from three biological replicates of GCs from follicles with different stages. Reverse transcription was carried out through first-strand cDNA synthesis kit (Thermo Fisher Scientific, Germany) with oligo (dT) 18 primers according to the manufacturer protocols. Primer 3 (version 4.0, http://bioinfo.ut.ee/primer3/) and Primer blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used to design gene-specific primers spanning the exons (Supplementary Table S3).

The RT-qPCR was carried out through iTaq™ universal SYBR^® Green Super-mix (Bio-Rad Laboratories GmbH, Germany) in Applied Biosystem^® Step-OnePlus™ (Ap-plied Biosystems, CA, United States). A reaction volume of 20 µl was used, including 7.4 µl of double-distilled H2O, 0.3 µl of forwarding primer, 0.3 µl of reverse primer, 10 µl of 1x SYBR Green master mix (Bio-Rad Laboratories GmbH, Germany), and 2 µl of cDNA template. Light Cycler 480 instrument (Roche, Germany) was used to perform qPCR. The second derivative maximum method was employed for data acquiring. Although its reliability as a housekeeping gene in GCs has not been validated in this study, the gene GAPDH was used as the reference gene in RT-qPCR. The CT value adjusted by reference gene was used to measure the gene expression level at each developmental stage. The ANOVA procedure applied in the SAS software was used to test the differences in relative mRNA expression levels among four developmental stages. Bonferroni t test was used to perform multiple comparisons. The relative mRNA expression levels of each gene at different developmental stages were presented with least square means.

3 Results

3.1 Polymorphism Screening

In this study, by using 68 bulls, a total of 37 variants were identified from polymorphism screening, including 13 in CACNB2, 10 in SLC39A12, and 14 in ZEB1. Although most genetic variants of dairy cattle are from the gene pool of bull, some variants presenting in bulls may be absent or with poor polymorphism in cow population. To efficiently perform KASP genotyping in 1,588 cows for target SNPs with variants, we performed a pilot analysis conducting of 46 cows. Some loci without polymorphism or with poor polymorphism in cow population was firstly filtered out, and the variants that do not meet primer design criteria of KASP are also removed. Furthermore, considering the high LD, only one SNP within 200 bp was remained for association analysis. From the pilot analysis using 46 selected samples, 18 variants with polymorphism in gene CACNB2 (nine), SLC39A12 (five) and ZEB1 (four) were further genotyped in 1,558 Holstein cows. The allele frequencies and the corresponding p-values derived from the Hardy-Weinberg equilibrium test for 18 SNPs obtained in 1,588 cows are presented in Table 1. Except for g.34061171T/C locating on ZEB1 (p-value > 0.002), allele frequencies of all remaining SNPs were not significantly deviated from the Hardy-Weinberg equilibrium. The distributions of SNP clustering (Supplementary Figure S1) suggested that the genotyping qualities were generally good for all 18 SNPs.

From annotation, one SNP (g.34066997C/G, gene ZEB1) is in the upstream regulatory region, two (g.33267172C/T, CACNB2 and g.34063562C/G, ZEB1) in exon region, and four in downstream regulatory region (g.32664855G/A, g.32665313G/A and g.32751518G/A in SLC39A12, and g.34061171T/C in ZEB1). In addition, 11 SNPs were detected from intronic regions, including eight in CACNB2 (g.33253706A/C, g.33284577T/C, g.33258042G/T, g.33258138T/G, g.33258186G/A, g.33258354A/G, g.33267056T/G, and g.33267296G/C), two in SLC39A12 (g.32698687A/G and g.32668290G/A), and one in ZEB1 (g.34110507T/A). The locus g.33258042G/T located in the intron of CACNB2 is a newly discovered variant.

3.2 Association Analyses of SNPs With Reproductive Traits

Significant associations with reproductive traits (p < 0.05 after Bonferroni correction) were observed in all three candidate genes (Supplementary Table S4). The total number of significant associations detected for each examined SNP or haplotype block is shown in Figure 2A. Most number of significant associations were detected for CACNB2, in which nine SNPs were associated with eight reproductive traits. Furthermore, all nine SNPs used in association analysis had significant effects on multiple reproductive traits, with the number of associated reproductive traits ranged from two (g.33267296G/C for CE_C and ICF) to five (g.33258042G/T for AFC, CE_C, CE_H, ICF, and IFL_C, and g.33258186G/A for AFC, AFS, CE_C, IFL_H, and SB_H). When counting by traits, there were nine and eight SNPs in CACNB2 associated with CE_C and AFC, respectively. However, no SNP within CACNB2 had significant association with SB_C.

FIGURE 2

FIGURE 2. Number of SNPs (A) or haplotype blocks (B) that were significantly associated with reproductive traits in gene CACNB2, SLC39A12 and ZEB1. The legends for nine reproductive traits are shown on the right side. In graph “B”, the “*” represents the significant association on the level of p < 0.05, “**” represents p < 0.01, and “***” represents p < 0.001.

Regarding the gene SLC39A12, g.32664855G/A, g.32665313G/A, and g.32751518G/A were significantly associated with CE_C, while g.32664855G/A and g.32751518G/A within SLC39A12 were significantly associated with both ICF and SB_H. For ZEB1, there were four SNPs significantly associated with at least one reproductive trait, including g.34066997C/G and g.34063562C/G for CE_C and ICF, g.34061171T/C for CE_C, CE_H, ICF, and SB_C, and g.34110507T/A for AFC.

The proportions of phenotypic variance explained by each SNP for every reproductive trait are presented in Supplementary Table S5. Among these 18 SNPs, the SNP g.33267056T/G in gene CACNB2 explained the largest phenotypic variance for AFC (0.011%), CE_C (0.002%), and IFL_C (0.004%), while the SNP g.32751518G/A in gene SLC39A12 explained the largest variance for SB_H (0.038%). For gene ZEB1, g.34061171T/C explained a large phenotypic variance for ICF (0.002%), SB_C (0.004%), and SB_C (0.009%).

3.3 Association Analyses of Haplotype Blocks with Reproductive Traits

The LD between every two SNPs were estimated for each candidate gene. Two haplotype blocks (Figure 3) were constructed for CACNB2, including three haplotypes for block 1 (consisting of g.33253706A/C and g.33258042G/T) and four haplotypes for block 2 (consisting of g.33258186G/A, g.33258354A/G, g.33267056T/G, and g.33267172C/T). Only one haplotype block (Figure 3) was constructed in SLC39A12 (three haplotypes, consisting of g.32664855G/A, g.32665313G/A, and g.32668290G/A) and also one for ZEB1 (three haplotypes, consisting of g.34061171T/C, g.34063562C/G, and g.34066997C/G).

FIGURE 3

FIGURE 3. Haplotype blocks constructed based on Linkage disequilibrium (LD) in Chinese Holstein cows for gene CACNB2, SLC39A12 and ZEB1.

The significant associations detected between haplotype blocks and reproductive traits (p < 0.05) are shown in Figure 2B and Table 2. Consistent with the results from SNP-based analyses, the block 1 within CACNB2 was significantly associated with AFC (p < 0.05) and CE_C (p < 0.01). The haplotype block 2 in CACNB2 was significantly associated with AFC (p < 0.01), CE_C (p < 0.01), AFS (p < 0.05), and CE_H (p < 0.05). Besides, haplotype block in ZEB1 was significantly associated with ICF (p < 0.01) and CE_C (p < 0.05).

3.4 Gene Expression Analysis in GCs Isolated from Follicle of Different Developmental Stages

By using qRT-PCR, the relative expression level of gene CACNB2, SLC39A12, and ZEB1 was measured in GCs which were isolated from bovine ovarian follicular representing different development stages (Figure 4A). The relative mRNA expression of CACNB2 in the first stage was significantly higher than that observed in other stages (p < 0.05) (Figure 4B). The expression level of the gene SLC39A12 in GCs was first increased in the early stages and then decreased significantly in later stages (p < 0.05) (Figure 4C). There was no significant difference in the relative expression of ZEB1 gene at follicular development stage 1 and 2 (Figure 4D), but the gene expression level was significantly higher than that from stage 3 and 4. In general, the differential expression level of gene CACNB2, SLC39A12, and ZEB1 in GCs at different development stages implied that they might play important roles in reproductive process of cows.

FIGURE 4

FIGURE 4. The bovine follicles at different developmental stages (A), and the relative mRNA expression level of the gene CACNB2 (B), SLC39A12 (C) and ZEB1 (D) in isolated granulosa. Based on the diameter (D), follicles were divided into four stages, including stage 1 (d ≤ 3 mm), stage 2 (3 mm < d ≤ 7 mm), stage 3 (7 mm < d ≤ 10 mm) and stage 4 (d > 10 mm). Different letters are used when the expression levels are significantly different between two stages (p < 0.05).

4 Discussion

In our previous GWAS of Chinese and Nordic Holstein populations, CACNB2, SLC39A12, and ZEB1 were suggested as candidate genes for reproductive traits. This study confirmed the association of genes CACNB2, SLC39A12 and ZEB1 with reproductive performances in dairy cattle.

The SNP g.33258186G/A locating in the intron region of CACNB2 was observed to have significant associations with all heifer traits (AFC, AFS, IFL_H, SB_H, and CE_H), with the genotype GG as the advantageous genotype. Similarly, another CACNB2 SNP (g.33267172C/T) in the exon region was significantly associated with AFC, IFL_H, CE_C, and AFS and explained large phenotypic variations in heifer traits. Although not reaching significant threshold, the advantageous genotype CC of g.33267172C/T also has favorable effects on SB_H and CE_H. All these findings suggested the critical roles of g.33258186G/A and g.33267172C/T in heifer reproduction. The relationship between SLC39A12 and reproductive performance of dairy cattle have been less reported previously. Amongst nine traits studied in this study, SLC39A12 was mostly responsive for SB_H. In SLC39A12, the SNP g.32751518G/A located in the 3′UTR region (related to post-transcriptional gene expression and affecting the regulation of mRNA expression) was significant associated with both CE_C and SB_H. Further studies are necessary to understand its biological relevance (Edwards et al., 2013; Sammeth et al., 2008). Interestingly, for two SNPs (g.32665313G/A (significantly associated with CE_C) and g.32668290G/A) which were significantly deviated from HWE, only a few individuals (2 and 13, respectively) carrying the rare genotypes. This together with the reported association of SLC39A12 with bovine hereditary zinc-deficiency (BHZD) syndrome (Yuzbasiyan-Gurkan and Bartlett, 2006), suggesting that there could be lethal alleles located in SLC39A12. Previous GWAS also reported strong associations with fertility traits in a region containing ZEB1 along with other genes on BTA13 (Höglund et al., 2009; Höglund et al., 2014; Liu et al., 2017b). In this study, we did observe significant associations for ZEB1, including SNPs g.34066997C/G and g.34063562C/G (both located in the exon region) for ICF and CE_C.

The mammalian ovarian follicle consists of one oocyte that undergoes a series of biological events, including ovulation, fertilization, and formation of an embryo. Oocyte is surrounded by granulosa and theca cells, which can produce signals and hormones to enable the development of oocyte (Su et al., 2004). During the follicle development, GCs replicate, secrete hormones, and provide a critical microenvironment for follicular growth (Petro et al., 2012). Proliferation and differentiation of GCs are essential for normal follicular growth, oocyte development, ovulation, and luteinization (Eppig, 2001; Bettaieb and Averill-Bates, 2008). In this study, the expressions of CACNB2, SLC39A12, and ZEB1 were quantified in GCs from different ovarian follicle developmental stages in mature cows. For example, the relative expression of CACNB2 was at its highest level at the beginning of follicular development and then drastically attenuated at later stages.

In Sahiwal cattle (Vineeth et al., 2020), Japanese Holstein (Sugimoto et al., 2013), and Nordic Holstein (Höglund et al., 2014), the gene CACNB2 has been reported to be the candidate gene for reproductive traits. In a meta-analysis of heifer reproductive traits, serval genes involved in the calcium signaling pathway, including CACNB2, was identified as the candidate genes (Tahir et al., 2021). The CACNB2 is a member of the voltage-gated calcium channel superfamily. Calcium regulation, through voltage-gated channels, is a prerequisite for GnRH-induced activation of extracellular signal-regulated kinase, which in turn leads to the secretion of FSH (Bliss et al., 2012; Dang et al., 2014). Knockdown of CACNB2 can reduce FSH secretion, while the overexpression of CACNB2 can increase its release (Sugimoto et al., 2013). There is also evidence that showing CACNB2 can control for the secretion of FSH. Following our previous GWAS in Chinese and Nordic Holstein population, the present study further confirmed the genetic effect on reproductive traits in cattle, at both SNP and mRNA expression levels.

Zinc is important for oocyte meiotic progression (Kim et al., 2010), of which the deficiency of zinc during pregnancy can lead to many pregnancy complications in human, such as early embryo death, intrauterine growth retardation and teratogenesis, or miscarriage (Uriu-Adams et al., 2010). The SLC39 family, including SLC39A12, is mainly involved in the transport of zinc into the cytoplasm (Kong et al., 2014), via secreting a kind of protein with a structure characteristic of zinc transporter. In the cumulus-oocyte complex, the removal of cumulus granulosa cells leads to an increase in zinc content in oocytes (Lisle et al., 2013). Our results showed that the expression of SLC39A12 was different across different follicle stages. Based on these observations, our hypothesis is that in cattle, SLC39A12 could contribute to oocyte maturation by involving in the process of zinc homeostasis and have the same capacity as reported in humans. Further studies towards understanding zinc metabolism needs to be done in order to discover the mechanism of the action of SLC39A12 gene at the junction of GCs and oocyte.

The gene ZEB1 (an alternative name of AREB6) was reported to play a role in the expression of epithelial-mesenchymal transition (Yamakoshi et al., 2012), which is partly required during the process of vegetative ectoderm junction. It was found that miR-200b, miR-429, and their target gene (ZEB1), act in the regulation of mammalian reproduction. Eliminating these miRNAs inhibits luteinizing hormone (LH) synthesis by repressing transcription of its b-subunit gene, which can cause lowered serum LH concentration, an impaired LH surge, and failure to ovulate. The hypothalamo-pituitary-ovarian axis was shown to require miR-200b and miR-429 to support ovulation (Hasuwa et al., 2013). In this study, differential expression with varying follicular developmental stages indicated that ZEB1 was likely to be involved in the process of follicular development. In bovine, it has been reported that ZEB1 was a candidate gene for temperament in Brahman and Yunling cattle (Chen et al., 2020), residual feed intake in Angus (Weber et al., 2016), and growth and meat quality traits in Brazilian Nellore (Mudadu et al., 2016). To our knowledge, no literature ever linked it to reproductive traits except for our previous GWAS in Chinese and Nordic Holsteins (Liu et al., 2017b).

5 Conclusion

Our previous GWAS identified genes CACNB2, SLC39A12, and ZEB1 as promising candidates for reproduction traits in dairy cattle. In this study, we further validated effects of these genes on reproductive traits by association analyses at a population level and mRNA expression assay at a cell level. Numerous SNPs within CACNB2, SLC39A12, and ZEB1 were significantly associated with reproductive traits, such as g.33258186G/A and g.33267172C/T within CACNB2, g.32751518G/A within SLC39A12, and g.34066997C/G and g.34063562C/G within ZEB1. A novel SNP g.33258042G/T located in the intron of CACNB2 was discovered. There were significant differences on the relative expressions of CACNB2, SLC39A12, and ZEB1 in granulosa cell at the different developmental stage of ovarian follicle. Findings from this study would contribute towards genetic selection of reproductive traits by improving the reliability of genomic evaluation in dairy cattle.

Data Availability Statement

The data presented in the study are deposited in the figshare repository, accession number 10.6084/m9.figshare.19846765.

Ethics Statement

Ethical review and approval was not required for the animal study because Ethical review and approval were waived for this study. Healthy ovaries were collected from the slaughterhouse. The blood samples and frozen semen were collected along with the regular quar-antine inspection of the farms and breeding stations, so no ethical approval was required for this study. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

Conceptualization, YW. Methodology, YD and HL. Formal analysis, YD. Resources, GG and LL. Data curation, YD, ZC, and HZ. Writing—original draft preparation, AS and RS. Writing—review and editing, HZ, AL, and YW. Visualization, RS, AL, and HZ. Supervision, YW. Project administration, YW. All authors have read and agreed to the published version of the manuscript.

Funding

This research was kindly funded by China Agriculture Research System of MOF and MARA; the Program for Changjiang Scholar and Innovation Research Team in University, IRT-15R62 and Beijing Sanyuan Breeding Technology Ltd Co. funded project, SYZYZ20190005. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Conflict of Interest

Author GG is employed by Beijing Sunlon Livestock Development Co., Ltd. and author LL is employed by the Beijing Dairy Cattle Center. Both authors are involved in custody and analysis of the data of dairy herd improvement program.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the Dairy Association of China (Beijing, China) for providing the pedigree datasets.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.882951/full#supplementary-material

Abbreviations

AFS, age at the first service; AFC, age at the first calving; CACNB2, calcium voltage-gated channel auxiliary subunit beta 2; CE_C, calving ease in cows; CE_H, calving ease in heifers; ICF, the interval from calving to the first insemination; IFL_C, the interval from the first to last insemination in cows; IFL_H, the interval from the first to last insemination in heifers; SB_C, stillbirth in cows; SB_H, stillbirth in heifers; SLC39A12, solute carrier family 39 member 12; ZEB1, zinc finger E-box binding homeobox 1.

References

Barrett, J. C., Fry, B., Maller, J., and Daly, M. J. (2005). Haploview: Analysis and Visualization of LD and Haplotype Maps. Bioinformatics 21, 263–265. doi:10.1093/bioinformatics/bth457

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry, D. P., Wall, E., and Pryce, J. E. (2014). Genetics and Genomics of Reproductive Performance in Dairy and Beef Cattle. Animal 8, 105–121. doi:10.1017/S1751731114000743

PubMed Abstract | CrossRef Full Text | Google Scholar

Bettaieb, A., and Averill-Bates, D. A. (2008). Thermotolerance Induced at a Fever Temperature of 40°C Protects Cells against Hyperthermia-Induced Apoptosis Mediated by Death Receptor Signalling. Biochem. Cell Biol. 86, 521–538. doi:10.1139/O08-136

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilyk, O., Coatham, M., Jewer, M., and Postovit, L.-M. (2017). Epithelial-to-Mesenchymal Transition in the Female Reproductive Tract: From Normal Functioning to Disease Pathology. Front. Oncol. 7, 145. doi:10.3389/fonc.2017.00145

PubMed Abstract | CrossRef Full Text | Google Scholar

Bliss, S. P., Navratil, A. M., Xie, J., Miller, A., Baccarini, M., and Roberson, M. S. (2012). ERK Signaling, but not c-Raf, Is Required for Gonadotropin-Releasing Hormone (GnRH)-Induced Regulation of Nur77 in Pituitary Gonadotropes. Endocrinology 153, 700–711. doi:10.1210/en.2011-0247

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q., Zhang, F., Qu, K., Hanif, Q., Shen, J., Jia, P., et al. (2020). Genome‐Wide Association Study Identifies Genomic Loci Associated with Flight Reaction in Cattle. J. Anim. Breed. Genet. 137, 477–485. doi:10.1111/jbg.12461

PubMed Abstract | CrossRef Full Text | Google Scholar

Clapham, D. E. (2007). Calcium Signaling. Cell 131, 1047–1058. doi:10.1016/j.cell.2007.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Cooper, D. N. (2010). Functional Intronic Polymorphisms: Buried Treasure Awaiting Discovery within Our Genes. Hum. Genomics 4, 284–288. doi:10.1186/1479-7364-4-5-284

PubMed Abstract | CrossRef Full Text | Google Scholar

Dang, A. K., Murtazina, D. A., Magee, C., Navratil, A. M., Clay, C. M., and Amberg, G. C. (2014). GnRH Evokes Localized Subplasmalemmal Calcium Signaling in Gonadotropes. Mol. Endocrinol. 28, 2049–2059. doi:10.1210/me.2014-1208

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobson, H., Smith, R., Royal, M., Knight, C., and Sheldon, I. (2007). The High-Producing Dairy Cow and its Reproductive Performance. Reprod. Domest. Anim. 42, 17–23. doi:10.1111/j.1439-0531.2007.00906.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, S. L., Beesley, J., French, J. D., and Dunning, A. M. (2013). Beyond GWASs: Illuminating the Dark Road from Association to Function. Am. J. Hum. Genet. 93, 779–797. doi:10.1016/j.ajhg.2013.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Emig, D., and Albrecht, M. (2011). Tissue-Specific Proteins and Functional Implications. J. Proteome Res. 10, 1893–1903. doi:10.1021/pr101132h

PubMed Abstract | CrossRef Full Text | Google Scholar

Eppig, J. (2001). Oocyte Control of Ovarian Follicular Development and Function in Mammals. Reproduction 122, 829–838. doi:10.1530/rep.0.1220829

PubMed Abstract | CrossRef Full Text | Google Scholar

Ettema, J. F., and Santos, J. E. P. (2004). Impact of Age at Calving on Lactation, Reproduction, Health, and Income in First-Parity Holsteins on Commercial Farms. J. Dairy Sci. 87, 2730–2742. doi:10.3168/jds.S0022-0302(04)73400-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, H., Kang, L., Abbas, Z., Hu, L., Chen, Y., Tan, X., et al. (2021). Identification of Key Genes and Pathways Associated with Thermal Stress in Peripheral Blood Mononuclear Cells of Holstein Dairy Cattle. Front. Genet. 12, 662080. doi:10.3389/fgene.2021.662080

PubMed Abstract | CrossRef Full Text | Google Scholar

Gabriel, S. B., Schaffner, S. F., Nguyen, H., Moore, J. M., Roy, J., Blumenstiel, B., et al. (2002). The Structure of Haplotype Blocks in the Human Genome. Science 296, 2225–2229. doi:10.1126/science.1069424

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, G., Guo, X., Wang, Y., Zhang, X., Zhang, S., Li, X., et al. (2014). Estimation of Genetic Parameters of Fertility Traits in Chinese Holstein Cattle. Can. J. Anim. Sci. 94, 281–285. doi:10.4141/cjas2013-113

CrossRef Full Text | Google Scholar

Han, B., Yuan, Y., Li, Y., Liu, L., and Sun, D. (2019). Single Nucleotide Polymorphisms of NUCB2 and Their Genetic Associations with Milk Production Traits in Dairy Cows. Genes 10, 449. doi:10.3390/genes10060449

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasuwa, H., Ueda, J., Ikawa, M., and Okabe, M. (2013). MiR-200b and miR-429 Function in Mouse Ovulation and Are Essential for Female Fertility. Science 341, 71–73. doi:10.1126/science.1237999

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, W. G., and Mackay, T. F. C. (2004). D. S. Falconer and Introduction to Quantitative Genetics. Genetics 167, 1529–1536. doi:10.1093/genetics/167.4.1529

PubMed Abstract | CrossRef Full Text | Google Scholar

Höglund, J. K., Guldbrandtsen, B., Su, G., Thomsen, B., and Lund, M. S. (2009). Genome Scan Detects Quantitative Trait Loci Affecting Female Fertility Traits in Danish and Swedish Holstein Cattle. J. Dairy Sci. 92, 2136–2143. doi:10.3168/jds.2008-1104

PubMed Abstract | CrossRef Full Text | Google Scholar

Höglund, J. K., Sahana, G., Brøndum, R., Guldbrandtsen, B., Buitenhuis, B., and Lund, M. S. (2014). Fine Mapping QTL for Female Fertility on BTA04 and BTA13 in Dairy Cattle Using HD SNP and Sequence Data. BMC Genomics 15, 790. doi:10.1186/1471-2164-15-790

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, A. M., Vogt, S., O'Halloran, T. V., and Woodruff, T. K. (2010). Zinc Availability Regulates Exit from Meiosis in Maturing Mammalian Oocytes. Nat. Chem. Biol. 6, 674–681. doi:10.1038/nchembio.419

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, B. Y., Duncan, F. E., Que, E. L., Kim, A. M., O'Halloran, T. V., and Woodruff, T. K. (2014). Maternally-Derived Zinc Transporters ZIP6 and ZIP10 Drive the Mammalian Oocyte-To-Egg Transition. Mol. Hum. Reprod. 20, 1077–1089. doi:10.1093/molehr/gau066

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisle, R. S., Anthony, K., Randall, M. A., and Diaz, F. J. (2013). Oocyte-Cumulus Cell Interactions Regulate Free Intracellular Zinc in Mouse Oocytes. Reproduction 145, 381–390. doi:10.1530/REP-12-0338

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, A., Lund, M. S., Wang, Y., Guo, G., Dong, G., Madsen, P., et al. (2017). Variance Components and Correlations of Female Fertility Traits in Chinese Holstein Population. J. Anim. Sci. Biotechnol. 8, 56. doi:10.1186/s40104-017-0189-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, A., Wang, Y., Sahana, G., Zhang, Q., Liu, L., Lund, M. S., et al. (2017). Genome-Wide Association Studies for Female Fertility Traits in Chinese and Nordic Holsteins. Sci. Rep. 7, 8487. doi:10.1038/s41598-017-09170-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lomelin, D., Jorgenson, E., and Risch, N. (2010). Human Genetic Variation Recognizes Functional Elements in Noncoding Sequence. Genome Res. 20, 311–319. doi:10.1101/gr.094151.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Matukumalli, L. K., Lawley, C. T., Schnabel, R. D., Taylor, J. F., Allan, M. F., Heaton, M. P., et al. (2009). Development and Characterization of a High Density SNP Genotyping Assay for Cattle. Plos One 4, e5350. doi:10.1371/journal.pone.0005350

PubMed Abstract | CrossRef Full Text | Google Scholar

Mi, S., Tang, Y., Dari, G., Shi, Y., Zhang, J., Zhang, H., et al. (2021). Transcriptome Sequencing Analysis for the Identification of Stable lncRNAs Associated with Bovine Staphylococcus A Mastitis. J. Anim. Sci. Biotechnol. 12, 120. doi:10.1186/s40104-021-00639-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mudadu, M. A., Porto-Neto, L. R., Mokry, F. B., Tizioto, P. C., Oliveira, P. S. N., Tullio, R. R., et al. (2016). Genomic Structure and Marker-Derived Gene Networks for Growth and Meat Quality Traits of Brazilian Nelore Beef Cattle. BMC Genomics 17, 235. doi:10.1186/s12864-016-2535-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Nayeri, S., Sargolzaei, M., Abo-Ismail, M. K., Miller, S., Schenkel, F., Moore, S. S., et al. (2017). Genome-Wide Association Study for Lactation Persistency, Female Fertility, Longevity, and Lifetime Profit Index Traits in Holstein Dairy Cattle. J. Dairy Sci. 100, 1246–1258. doi:10.3168/jds.2016-11770

PubMed Abstract | CrossRef Full Text | Google Scholar

Parekh, A. B. (2011). Decoding Cytosolic Ca2+ Oscillations. Trends Biochem. Sci. 36, 78–87. doi:10.1016/j.tibs.2010.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Petro, E. M. L., Leroy, J. L. M. R., Van Cruchten, S. J. M., Covaci, A., Jorssen, E. P. A., and Bols, P. E. J. (2012). Endocrine Disruptors and Female Fertility: Focus on (Bovine) Ovarian Follicular Physiology. Theriogenology 78, 1887–1900. doi:10.1016/j.theriogenology.2012.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Sammeth, M., Foissac, S., and Guigó, R. (2008). A General Definition and Nomenclature for Alternative Splicing Events. PLoS Comput. Biol. 4, e1000147. doi:10.1371/journal.pcbi.1000147

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, B., Serão, N., and Ross, J. W. (2020). “Genetic Improvement of Livestock, from Conventional Breeding to Biotechnological Approaches,” in Animal Agriculture. Editors F. W. Bazer, G. C. Lamb, and G. Wu (Ames: Academic Press), 393–405. doi:10.1016/b978-0-12-817052-6.00023-9

CrossRef Full Text | Google Scholar

Struck, T. J., Mannakee, B. K., and Gutenkunst, R. N. (2018). The Impact of Genome-Wide Association Studies on Biomedical Research Publications. Hum. Genomics. 12, 38. doi:10.1186/s40246-018-0172-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, Y.-Q., Wu, X., O'Brien, M. J., Pendola, F. L., Denegre, J. N., Matzuk, M. M., et al. (2004). Synergistic Roles of BMP15 and GDF9 in the Development and Function of the Oocyte-Cumulus Cell Complex in Mice: Genetic Evidence for an Oocyte-Granulosa Cell Regulatory Loop. Dev. Biol. 276, 64–73. doi:10.1016/j.ydbio.2004.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugimoto, M., Sasaki, S., Gotoh, Y., Nakamura, Y., Aoyagi, Y., Kawahara, T., et al. (2013). Genetic Variants Related to Gap Junctions and Hormone Secretion Influence Conception Rates in Cows. Proc. Natl. Acad. Sci. U.S.A. 110, 19495–19500. doi:10.1073/pnas.1309307110

PubMed Abstract | CrossRef Full Text | Google Scholar

Tahir, M. S., Porto-Neto, L. R., Gondro, C., Shittu, O. B., Wockner, K., Tan, A. W. L., et al. (2021). Meta-Analysis of Heifer Traits Identified Reproductive Pathways in Bos I Cattle. Genes 12, 768. doi:10.3390/genes12050768

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, K. M., and Nicholson, R. I. (2003). The LZT Proteins; the LIV-1 Subfamily of Zinc Transporters. Biochimica Biophysica Acta (BBA) - Biomembr. 1611, 16–30. doi:10.1016/s0005-2736(03)00048-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Uriu-Adams, J. Y., and Keen, C. L. (2010). Zinc and Reproduction: Effects of Zinc Deficiency on Prenatal and Early Postnatal Development. Birth Defects Res. Part B Dev. Reproductive Toxicol. 89, 313–325. doi:10.1002/bdrb.20264

CrossRef Full Text | Google Scholar

Vineeth, M. R., Surya, T., Sivalingam, J., Kumar, A., Niranjan, S. K., Dixit, S. P., et al. (2020). Genome-Wide Discovery of SNPs in Candidate Genes Related to Production and Fertility Traits in Sahiwal Cattle. Trop. Anim. Health. Prod. 52, 1707–1715. doi:10.1007/s11250-019-02180-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., Mccarthy, M. I., Brown, M. A., et al. (2017). 10 Years of GWAS Discovery: Biology, Function, and Translation. Am. J. Hum. Genet. 101, 5–22. doi:10.1016/j.ajhg.2017.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, K. L., Welly, B. T., Van Eenennaam, A. L., Young, A. E., Porto-Neto, L. R., Reverter, A., et al. (2016). Identification of Gene Networks for Residual Feed Intake in Angus Cattle Using Genomic Prediction and RNA-Seq. Plos One 11, e0152274. doi:10.1371/journal.pone.0152274

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamakoshi, S., Bai, R., Chaen, T., Ideta, A., Aoyagi, Y., Sakurai, T., et al. (2012). Expression of Mesenchymal-Related Genes by the Bovine Trophectoderm Following Conceptus Attachment to the Endometrial Epithelium. Reproduction 143, 377–387. doi:10.1530/REP-11-0364

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuzbasiyan-Gurkan, V., and Bartlett, E. (2006). Identification of a Unique Splice Site Variant in SLC39A4 in Bovine Hereditary Zinc Deficiency, Lethal Trait A46: An Animal Model of Acrodermatitis Enteropathica. Genomics 88, 521–526. doi:10.1016/j.ygeno.2006.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., Peng, F., Zhou, Y., Hu, G., He, H., He, F., et al. (2015). Exon Sequencing Identifies a novelCHRNA3-CHRNA5-CHRNB4variant that Increases the Risk for Chronic Obstructive Pulmonary Disease. Respirology 20, 790–798. doi:10.1111/resp.12539

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Ding, X., Zhang, Q., Wang, Y., Lund, M. S., and Su, G. (2013). Consistency of Linkage Disequilibrium between Chinese and Nordic Holsteins and Genomic Prediction for Chinese Holsteins Using a Joint Reference Population. Genet. Sel. Evol. 45, 7. doi:10.1186/1297-9686-45-7

PubMed Abstract | CrossRef Full Text | Google Scholar