Variation in Follicle-Stimulating Hormone Receptor Expression Is Associated with the Twinning Rate QTL Located on Bovine Chromosome 11 in Holstein Cattle

TL;DR

This study finds that higher expression of the follicle-stimulating hormone receptor gene in Holstein cows is linked to increased twin births.

Contribution

The study identifies differential gene expression, not coding variants, as the mechanism linking a chromosome 11 QTL to twinning rates in Holstein cattle.

Findings

FSHR gene expression is significantly associated with the twinning rate QTL genotype.

Increased FSHR expression correlates with higher twinning rates in Holstein cows.

LHCGR gene expression is not significantly linked to the twinning rate QTL.

Abstract

Multiple gene mapping studies have reported evidence for a gene contributing to twin births in Holstein dairy cattle on bovine chromosome 11. The identified region contains two genes which encode receptors for reproductive hormones that play critical roles in supporting development of follicles and eggs in the ovary. These gene mapping studies found no association with variants in the genes’ coding sequence and twinning rate, leading to the hypothesis that the effect is due to a difference in the expression of one of these genes. This study tested that hypothesis and found evidence that the follicle-stimulating hormone receptor gene has an elevated expression in conjunction with an increased twinning rate.

Twin births in dairy cattle present challenges for producers, resulting in increased prevalence of health issues for both cows and calves, thereby impacting profitability. Genome-wide association study (GWAS) analyses of the twinning rate in Holstein cattle have reported the most significant genomic association with twinning rate in a region containing two strong candidate genes: follicle-stimulating hormone receptor (FSHR) and luteinizing hormone/chorionic gonadotropin receptor (LHCGR). Coding-sequence variants of these genes were not associated with the twinning rate, suggesting that one of the two genes is differentially expressed in association with the twinning rate. Granulosa cells were collected from 38 Holstein cows that were selected to provide similar representation of genotypes for the twinning rate quantitative trait locus (QTL). RNA was extracted from granulosa cells and gene expression was assessed by quantitative PCR with data analyzed by the ΔΔCT method. Association of gene expression with QTL genotype was tested by the Kruskal–Wallis test with the QTL genotype based on the SNP most significantly associated with twinning rate. QTL genotype was significantly associated (p = 1.88 × 10−8) with the expression of FSHR but was not associated with LHCGR expression (p = 0.18). The increased FSHR expression was associated with an increasing copy number of the G allele and thus an increased twinning rate.

Figures

Full text

1. Introduction

Twin births in dairy cattle, while unavoidable, present significant challenges to the profitability of dairy farming. These challenges are the result of increased health risks for both the cows and their twin calves, including a higher likelihood of reproductive and metabolic disorders around the time of birth, as well as increased instances of abortion, stillbirth, neonatal calf death, and lower birth weights compared to single births [1,2].

Twinning is a complex trait influenced by both genetic and environmental factors. Factors that affect ovulation rate and the likelihood of twinning include genetics, season, parity, and milk production. The increased prevalence of twins in dairy cattle over time suggests that one or more of these factors have also changed during this period from the 1980s to the 1990s [1,2,3,4,5,6].

Despite the potential advantages of additional calves from twin births, the associated challenges and risks often result in net financial deficits for dairy farms. Consequently, most dairy operations strive to reduce the occurrence of twinning [7]. This issue has substantial repercussions for the dairy industry, yet dairy farmers and their advisors currently lack the necessary scientific data to make informed management decisions to mitigate the negative impacts of twinning [8]. Therefore, a comprehensive understanding of the factors contributing to twinning is crucial for devising future strategies to manage this issue in dairy farming [9].

To mitigate twin births in dairy cattle, several interventions have been proposed. These include estrus synchronization protocols to elevate progesterone levels and decrease the incidence of double ovulation, termination of twin pregnancies, selective reduction in twin pregnancies to singletons by manually rupturing one amnion, and genomic selection to reduce the genetic predisposition for twinning [10].

Genomic Analysis of Twinning in Holstein Cattle

The genomic analysis of twinning in Holstein cattle has provided significant insights into the genetic factors influencing this trait. Recent studies have identified key quantitative trait loci (QTL) and candidate genes that play crucial roles in reproductive processes, offering a deeper understanding of the genetic underpinnings of twinning.

In a study involving approximately 2500 Swiss Holsteins and utilizing around 600,000 imputed SNPs, researchers identified a major QTL on chromosome 11 (BTA11) for multiple births (p < 5 × 10^−9^), explaining about 16% of the total genetic variance [11]. This QTL was fine-mapped to a 70 kb window between 31.00 and 31.07 Mb on BTA11. This window was near two candidate genes, luteinizing hormone/chorionic gonadotropin receptor (LHCGR) and follicle-stimulating hormone receptor (FSHR), including the 5′ end of LHCGR and adjacent area and being approximately 180 kb from the 3′ end of FSHR. These genes could be strong candidates as they encode receptors for hormones involved in folliculogenesis. Analysis of whole-genome sequence data revealed a potential regulatory variant in the 5′ region of LHCGR as a possible candidate causal variant for the identified major QTL. Furthermore, the identified haplotype showed significant associations with stillbirth and days to first service, suggesting a direct influence on the trait of twinning.

A similar study involving the North American Holstein population explored genetic factors and pathways affecting twinning rates while also assessing the feasibility of genomic selection [12]. By employing single-step GWAS with whole-genome sequence data, the study identified significant SNPs linked to twinning rates. As in the Swiss Holstein population, the 31 Mb region of BTA11 containing LHCGR and FSHR was identified as the most significantly associated genomic region (p < 6 × 10^−12^). However, despite the greater SNP density employed, there was no evidence of coding-sequence changes in candidate genes being associated with twinning rate. Specifically, a coding-sequence variant in LHCGR reported to be associated with superovulatory response had no significant association with the twinning rate [12]. This, together with the strongest associations in this study and the Swiss study being outside the candidate genes’ coding regions, suggests the underlying causative variant is regulatory. This raises the question of which of these two candidate genes is differentially expressed in association with the twinning-rate QTL. Additionally, the study evaluated the potential of genomic selection to enhance twinning rates, suggesting that the integration of genomic data could improve selection accuracy and efficiency.

The LHCGR/FSHR genomic region’s association with twinning rates has been additionally validated in a recent GWAS of the Irish Holstein Friesian population [13]. This study included 2656 Holstein Friesian sires, using approximately 50,000 SNPs imputed to a full genome sequence. The most significant SNP association was observed on BTA11 near the candidate genes FSHR and LHCGR, a finding consistent with previous studies in Holstein Friesian populations in the US and Switzerland. The adjacent windows flanking the most significant SNP collectively accounted for approximately 0.5% of the variance in de-regressed breeding value estimates, indicating a modest but highly significant genetic association (p < 1 × 10^−9^). This report, combined with the previous results from the Swiss and North American Holstein populations, provides extremely strong support for the contribution of this genomic region to twinning rate. Further investigation is needed as to which of the two candidate genes accounts for this association.

As described above, results from the GWAS studies have suggested the contribution of either LHCGR or FSHR to a twinning rate variation is likely due to differences in gene expression. Consequently, a preliminary question is in what cell types should gene expression be evaluated? Gene expression analysis in cattle has indicated that elevated expression of FSHR and LHCGR occurs in the gonads, both ovary and testis [14]. Within the ovary, gene expression in granulosa cells is closely related to ovulation rate due to the critical roles these cells play in folliculogenesis and hormone regulation [15,16]. Granulosa cells surround the developing oocyte and are responsible for producing hormones, such as estrogen and inhibins, which are essential for follicle maturation and ovulation. The expression of FSHR and LHCGR in granulosa cells directly influences the granulosa cells’ ability to respond to hormonal signals, affecting granulosa cell proliferation, differentiation, and follicle development [17]. Understanding these gene expression patterns can provide valuable insights into reproductive processes and potential genetic factors influencing ovulation rates [18]. Consequently, granulosa cells were chosen as the focus for a gene expression analysis of FSHR and LHCGR. Expression of LHCGR is a marker for follicular deviation, i.e., a follicle becoming a dominant follicle with the potential to ovulate. As a result, examination of FSHR and LHCGR expression would require examination of granulosa cells from post-deviation follicles.

The use of granulosa cells in gene expression analysis for FSHR and LHCGR is justified by their crucial roles in follicular development. Dynamic changes in the expression of LHCGR and FSHR in granulosa cells occur during varying stages of follicular development in cattle. The expression of LHCGR is significantly higher in granulosa cells from growing follicles compared to those in atretic stages, and its expression is a marker for acquisition of follicular dominance [19]. Similarly, the expression of FSHR is higher in growing follicles, emphasizing its importance in the initial stages of follicular growth [19,20]. These findings underscore the roles of LHCGR and FSHR in folliculogenesis [16]. Consequently, the focus of this study was examination of LHCGR and FSHR expression in granulosa cells near the time of deviation and the association with genotype for the BTA11 twinning rate QTL.

2. Materials and Methods

2.1. Sample Collection

The Institutional Animal Care and Use Committee of the College of Agricultural and Life Sciences, University of Wisconsin–Madison, approved this research. Granulosa cells were obtained from 41 Holstein cows from The University of Wisconsin–Madison’s dairy herd. Cows were chosen to provide similar representation of genotypes for the twinning rate quantitative trait locus (QTL) with 14, 14, and 13 animals of CC, CG, and GG genotypes, respectively. The cows sampled underwent estrous synchronization for the first insemination following calving using a Double-Ovsynch protocol [21]. Six days after timed artificial insemination (TAI), a dose of gonadotropin-releasing hormone (GnRH, 50 µg gonadorelin diacetate tetrahydrate, Cystorelin, Boehringer Ingelheim, Duluth, GA) was administered to the cows to induce the luteinization of dominant follicle(s) from the first follicular wave following breeding and initiate a new follicular wave by release of follicle-stimulating hormone (FSH). Consequently, dominant follicle(s) from the new follicular wave developed in a hormonal background of higher progesterone levels, mimicking follicular development during the ovulatory follicular wave (Figure 1).

The new dominant follicle was aspirated seven days after GnRH administration to collect granulosa cells for RNA extraction and subsequent gene expression analyses. The aspirate was examined under a low-power microscope (Nikon SMZ-1B) to identify and separate the oocyte before centrifugation of the sample at 500× g using an Eppendorf 5415C centrifuge to separate the supernatant (follicular fluid) from the pelleted granulosa cells. The supernatant was placed on an ice pack and 90 µL of RNAlater (Ambion, Austin, TX, USA) was added to the cell pellet (granulosa cells) to preserve RNA. Similarly, 30 µL of RNAlater was combined with the oocyte for RNA preservation. All sample types were subsequently stored at −80 °C. The concentration of estradiol (E2) in the follicular fluid was measured following an Elisa assay protocol [22] and used to assess the dominance status of the follicle (E2 > 200 ng/mL) [23].

2.2. Gene Expression Analysis Preparation

RNA was extracted using the NucleoSpin^®^RNA Plus kit (MACHEREY-NAGEL, Bethlehem, PA, USA) per the manufacturer’s protocol, and the integrity of the RNA samples was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Subsequently, cDNA was synthesized using the SuperScript IV Single Cell/Low-Input cDNA PreAmp kit (Thermo Fisher Scientific, Waltham, MA, USA) to facilitate qPCR reactions. The expression of mRNA was quantified using the SsoAdvanced Universal Probes Supermix (Bio-Rad, Hercules, CA, USA) on a Bio-Rad CFX Connect instrument (Bio-Rad, Hercules, CA, USA). qPCR amplification of the candidate genes FSHR and LHCGR, as well as the reference gene splicing factor 3a subunit 1 (SF3A1), was performed using TaqMan gene expression assays (Thermo Fisher Scientific, Waltham, MA, USA) specific for Bos taurus (Table 1). SF3A1 was chosen as the reference gene based on previous work indicating an expression level comparable to target genes and no genotype-related differences in expression in granulosa cells from Trio allele carriers and non-carriers [18]. To validate primer efficiency in qPCR, a series of cDNA dilutions were created and qPCR was run. C_T_ values were plotted against the dilution factors to generate a standard curve and calculate efficiency, with an acceptable efficiency being in a range of 90–110% [24].

2.3. Genotype Data and Statistical Analyses

Single nucleotide polymorphism genotype data (79k) were available for all animals, but the SNP most strongly associated with the twinning rate in the most recent analysis (11: 31065831; chromosome:basepair in bovine genome assembly ARS-UCD1.2) [13] was not included. Consequently, high-density genotypes on BTA11 (n = 274,462), including 11:31065831, were imputed. Imputation was done directly from BTA11 79k data (2029 SNPs) to high density using Beagle version 5.4 [25] and genotype data from run7 of the 1000 Bull Genomes Project as a reference [26]. Only animals with a dairy breed composition and genome coverage of 8× or greater were included in the reference data, as previous work had found that the inclusion of lower coverage samples lessened imputation accuracy [13]. Imputation accuracy was examined by masking four reference samples to 79k and determining concordance between imputed and original genotypes, the resulting concordance of which was 97%.

qPCR data were analyzed using the ΔΔC_T_ method [27] to quantify relative expression of FSHR and LHCGR normalized to the reference gene SF3A1 and to calculate fold-change differences among QTL genotypes. ΔC_T_ values were linearly transformed (2^−ΔC^T) to generate gene expression values for statistical analyses [27]. Assumptions of normality and homogeneity of variance were evaluated in R (version 4.4.2) using the Shapiro–Wilk test [28] and Levene’s test [29], respectively. Expression values (2^−ΔC^T) violated both assumptions (p < 0.05), therefore differences among genotypes (CC, CG, GG) were assessed using the Kruskal–Wallis test [30] followed by Dunn’s post hoc pairwise comparisons with p-values adjusted for multiple testing [31]. Expression distributions were visualized using jittered dot plots with group means and confidence intervals. Fold-change values were computed as the ratio of mean expression between genotypes. The effect of the genotype on follicle diameter and follicular fluid estradiol level was evaluated by analysis of variance.

3. Results

No association was observed between genotypes and either follicle diameter or estradiol level in follicular fluid. The Kruskal–Wallis analysis indicated a significant association (p = 1.88 × 10^−8^) of the expression of FSHR with the twinning rate QTL genotype (Figure 2). Pairwise comparisons showed significant differences among all genotypes, with adjusted p-values of p = 5.95 × 10^−3^, p = 7.52 × 10^−9^, and p = 1.03 × 10^−2^ for CC vs. CG, CC vs. GG, and CG vs. GG, respectively. Specifically, FSHR expression increased with the increasing twinning rate QTL genotype. In contrast, there was no significant association (p = 0.18) of LHCGR expression with the twinning rate QTL genotype (Figure 3). The fold change in FSHR expression for the highest (genotype GG) versus lowest (CC) twinning rate QTL genotype was a 17.78-fold increase. FSHR expression for genotype GG vs. CG and CG vs. CC increased by 3.88-fold and 4.58-fold, respectively. These results demonstrate that the BTA11 twinning rate QTL is associated with FSHR gene expression, supporting the hypothesis of a regulatory mechanism underlying the QTL.

4. Discussion

The results of this study provide strong evidence supporting the hypothesis that the BTA11 twinning rate QTL acts through a regulatory variant that alters FSHR expression. Further work is needed to test genotype association with FSHR protein abundance. Twinning in Holstein cattle is primarily dizygotic [32], with an estimated frequency of approximately 95% of twin births. This implies that the association of FSHR expression with twinning is through the underlying trait of ovulation rate. FSHR is primarily expressed in granulosa cells, where it facilitates response to FSH, leading to follicle growth and development. Initiation of a wave of follicular development occurs in response to a transient increase in circulating FSH. As one follicle in the recruited cohort acquires responsiveness to LH by expressing LHCGR and producing LH receptors, it becomes dominant, producing estradiol and inhibin that reduces the FSH available to sustain development of subordinate follicles; consequently, the subordinate follicles become atretic. We speculate that elevated expression of FSHR could support continued follicular development as the FSH level declines, permitting a second follicle to become dominant with an ensuing double ovulation event. These findings, combined with the absence of coding-sequence associations with twinning rates in three GWAS studies, suggest that the causal variant is regulatory rather than structural. The increase in FSHR expression is dramatic, but an increase of this magnitude is not unprecedented. In previous work with another regulatory variant affecting ovulation rate, the Trio allele, an approximately seven- to nine-fold increase in SMAD6 expression was observed for Trio allele carriers versus wild-type homozygotes [18,23]. In the current study, the change between heterozygote and homozygote genotypes (CG vs. CC, CG vs. GG) was approximately four-fold.

Other studies have also suggested an association between FSHR and ovulation rate or multiple births in cattle. Two studies comparing FSHR expression in cattle selected for ovulation and twinning rate with cattle not selected for twinning have been reported with somewhat contradictory results. One study reported increased FSHR expression in granulosa cells from small pre-deviation follicles in twinner versus non-twinner cows, similar to the results reported herein [15]. However, a second study comparing the same populations reported more equivocal results. Follicles were collected at two different time points (days 3–4 vs. days 5–6) and were categorized as estradiol-active or estradiol-inactive based on an estradiol:progesterone ratio in follicular fluid, with FSHR expression assessed in the two largest follicles of each type on each day. For estradiol-active follicles, a difference was observed only for the largest follicle on days 3–4, with control cows having higher FSHR expression versus twinner cows. Conversely, for estrogen-inactive follicles, FSHR expression in the largest follicle on days 3–4 was higher in twinner cows versus non-twinners, while the result was reversed for the second largest follicle on days 3–4 [33]. A recent study examining a missense variant in FSHR and its relation to reproductive function and in vitro embryo production failed to find an effect with those traits, but an association with ovulation rate was observed [34]. Whether this reflects linkage disequilibrium between the missense mutation and the unknown FSHR regulatory variant suggested by the twinning GWAS studies and our current results is unknown.

Reports of association of FSHR variants with litter size differences in sheep also support the suggestion that FSHR contributes to genetic variation for multiple births. Variants in the 5′ regulatory region of prolific Chinese sheep breeds were tested for association with litter size in a candidate gene study with nominally significant (unadjusted p < 0.05) associations reported [35]. In a comparison of lowly prolific (Ile-de-France) and highly prolific (Romanov) sheep breeds, FSHR expression was significantly higher in pre-deviation (early follicular phase) granulosa cells from the more prolific breed; no significant difference between breeds was observed for FSHR expression post-deviation [36]. Furthermore, a recent report of a GWAS for litter size in two prolific Chinese sheep breeds found the strongest evidence for genomic association with twinning rates for SNPs in the FSHR locus [37].

Production of bovine embryos by superovulation, flushing, and embryo transfer is a well-known reproductive tool with more than 350,000 in vivo derived embryos produced annually worldwide [38]. Traditional superovulation protocols entail twice-daily administration of FSH over four days to increase the number of ovulatory follicles [39]. Response to this treatment is known to be highly variable [40,41]. Another way in which the BTA11 twinning rate QTL’s association with FSHR expression may have practical significance is that this QTL may be predictive of animal variation in response to superovulatory treatment. Two GWAS studies of responses to superovulatory treatment for the production of in vivo and in vitro embryos in dairy cattle have reported significant SNP associations on BTA11, but the regions of BTA11 identified are on the distal end of the chromosome, greater than 70 Mb from the locations of FSHR and LHCGR [42,43]. Other studies have taken a candidate gene approach in examining the association of FSHR genotype with outcomes from superovulatory treatment. Two SNPs located 320 and 270 bases upstream of the FSHR coding sequence were tested for association with superovulatory response in Chinese Holstein cattle, with significant differences in total number of ova and transferable embryos observed [44]. Non-synonymous coding-sequence variants in FSHR were examined for association with superovulatory response in Canadian Holstein cows, with all three missense variants showing association with some aspect of the response [45]. For c337C>G, the genotype was associated with percent viable embryos and percent unfertilized oocytes, the genotype for c871A>G was associated with embryo yield, and the C1973C>G genotype was associated with both embryo yield and percent unfertilized oocytes. Hirayama et al. [46] examined association of the same three variants with superovulatory response in Japanese Black cattle and found only the c337C>G variant to be significantly associated with number of transferable embryos. A recent study followed up on the reported association of the c337C>G missense variant with embryo viability by assessing association with in vitro blastocyst development, but no significant association was observed [34]. Whether the genetic variation in FSHR expression described herein impacts response to superovulatory treatment remains to be determined. If there is an association of this QTL with superovulatory response, the QTL genotype could aid in identifying females that are the best candidates for superovulation or eventually lead to the development of protocols most suitable for cows of a different genotype.

In the context of dairy cattle breeding, the favorable twinning rate QTL allele is the one associated with a lower twinning rate, given the negative consequences associated with twin births. As reported in one of the three studies identifying the BTA11 genomic association with twinning rate [13], no detrimental association of selection for the lower twinning rate QTL allele with either lactation or other reproductive traits would be expected. A remaining challenge is the identification of the actual causative variant. The information reported herein may help facilitate that effort as the differentially expressed gene is now strongly suggested to be FSHR, and evidence from the current analysis suggests that the underlying variant must be regulatory rather than coding, consistent with the three GWAS studies that mapped the QTL to non-coding regions. Regulatory variants are a challenge to identify but also provide the opportunity for in vitro validation. Candidate variants could be edited into granulosa cell lines containing the alternative allele, with effects of the gene edit on FSHR expression potentially confirming the causative nature of a variant. Once identified, the desired allele could be added to a growing list of known alleles that could be jointly edited into improved germplasm [47]. More immediately, a causative variant could be used to facilitate more effective genomic selection relative to twinning rates.

5. Conclusions

Results of the study reported herein provide strong evidence that the genetic variant underlying a genomic association on BTA11 with twinning rates in Holstein cattle is regulatory and affects expression of the follicle-stimulating hormone receptor gene. Further investigation is needed to identify the actual genetic variant and to test if the variant accounts for variation in response to superovulatory protocols as well as twinning rates.