Genetic effects of PDGFRB and MARCH1 identified in GWAS revealing strong associations with semen production traits in Chinese Holstein bulls

Using a genome-wide association study strategy, our previous study discovered 19 significant single-nucleotide polymorphisms (SNPs) related to semen production traits in Chinese Holstein bulls. Among them, three SNPs were within or close to the phosphodiesterase 3A (PDE3A), membrane associated ring-CH-type finger 1 (MARCH1) and platelet derived growth factor receptor beta (PDGFRB) genes. The present study was designed with the objectives of identifying genetic polymorphism of the PDE3A, PDGFRB and MARCH1 genes and their effects on semen production traits in a Holstein bull population. A total of 20 SNPs were detected and genotyped in 730 bulls. Association analyses using de-regressed estimated breeding values of each semen production trait revealed four statistically significant SNPs for one or more semen production traits (P < 0.05): one SNP was located downstream of PDGFRB and three SNPs were located in the promoter of MARCH1. Interestingly, for MARCH1, haplotype-based analysis revealed significant associations of haplotypes with semen volume per ejaculate. Furthermore, high expression of the MARCH1 gene was observed in sperm cells. One SNP (rs43445726) in the regulatory region of MARCH1 had a significant effect on gene expression. Our study demonstrated the significant associations of genetic variants of the PDGFRB and MARCH1 genes with semen production traits. The identified SNPs may serve as genetic markers to optimize breeding programs for semen production traits in Holstein bull populations.


Background
In livestock breeding, the diagnosis of male fertility is very important because about half of pregnancy failures can be attributed to decreased male fertility or male factor infertility [1]. Sires with subfertility problems lead to larger economic losses than infertile ones because the latter can be detected early while the former require a long period of observation. The significant economic importance of male fertility is also relevant in dairy cattle, especially in the situation in which artificial insemination is widely used. The quality and quantity of semen can be measured by semen production traits, such as semen volume, sperm motility and sperm concentration, as well as observations of abnormal spermatozoa. Semen production traits are complicated, being affected by many nongenetic factors such as age, season, interval between ejaculations and bull handlers, as well as genetic factors [2][3][4]. Semen volume, sperm concentration and the number of spermatozoa have been estimated to have moderate heritability (0.15-0.30), while sperm motility has been found to be highly heritable (close to 0.6) [4]. Similar results were obtained by Karoui et al., namely, that heritability estimates for semen production traits were moderate (0. 16-0.22) [5]. Therefore, genetic improvement of these traits via selection is possible.
Many candidate genes for semen production traits have been revealed using candidate association analysis and genome-wide association study (GWAS). Fortes et al. detected the most significant SNPs in the X chromosome associated with the percentage of progressive motile spermatozoa at 18 months of age and the percentage of normal spermatozoa at 24 months of age [6]. Hering et al. also highlighted several candidate genes associated with sperm concentration, sperm motility and sperm volume in Holstein-Friesian populations [7][8][9]. In addition, some genes such as FSHR, INHA, TNP1, TNP2, CAPN1 and SPAG11 have been widely studied as candidate genes for semen production traits of bulls [10][11][12][13][14]. Using Illumina Bovine SNP50 Beadchip (Illumina Inc., San Diego, CA, USA), our previous GWAS detected 19 significant SNPs for five semen production traits in a population of 692 Chinese Holstein bulls [15]. Of those, three SNPs located within or close to the phosphodiesterase 3A (PDE3A), platelet derived growth factor receptor beta (PDGFRB) and membrane associated ring-CH-type finger 1 (MARCH1) genes were significantly associated with initial sperm motility (P = 3.31E −05 ), sperm volume per ejaculate (P = 3.75E −05 ) and sperm volume per ejaculate (P = 6.00E −05 ), respectively.
In this study, we aim to investigate genetic variants potentially related to semen production traits in an independent dairy cattle population. We also explore the potential impact of SNP variation in regulatory regions of the above genes on gene expression.

Resource population and analysis of de-regressed EBVs
A total of 730 Chinese Holstein bulls were selected without overlapping with the population of our previous GWAS to construct a single population in this study.
where y ijklmnopq is the phenotypic value of each trait of bulls; μ is the overall mean; F i represents the fixed effect of farm; H ij represents handlers of semen collection, which is nested in the farm effect; A k represents the fixed effect of age; S l represents the fixed effect of the season when frozen semen samples were collected; T m represents the number of collections on one day; I n represents the interval (in days) between collections; α o is the random polygenic effect, distributed as N (0, Aσ 2 a ) with the polygenic relationship matrix A and the additive genetic variance σ 2 a ; PE p is the permanent environment effect; and ε ijklmnopq is the random residual, distributed as N (0, I σ 2 e ) with identity matrix I and residual error variance σ 2 e . The EBVs were de-regressed, and the weights were calculated using the method proposed by Garrick et al. (2009) [16]. The descriptive statistics of the de-regressed and original EBVs for five semen production traits in the 730 bulls are listed in Table 1.
The levels of heritability of SVPE, SMOT, SCPE, NSPE and NMSPE were estimated to be 0.15, 0.12, 0.22, 0.16 and 0.12, respectively. Positive genetic correlations were observed among all traits, and the highest correlation was observed between NSPE and NMSPE ( Table 2).

SNP identification and genotyping
Genomic DNA was isolated from the frozen semen of 730 bulls using a standard phenol-chloroform method. The quality and quantity of extracted genomic DNA were measured with a NanoDrop™ Spectrophotometer (ND-2000c) (Thermo Scientific, Chelmsford, MA, USA) and gel electrophoresis. Then, each DNA sample was diluted to 50 ng/μL and stored at −20°C for subsequent use. A DNA pool was constructed from 50 randomly selected samples with equal amounts of DNA (50 ng/μL). A total of 83 pairs of primers were designed to amplify entire coding regions, partial introns, and 5′ upstream (3000 bp) and 3′ downstream regions (3000 bp) based on the genomic sequences of the bovine PDE3A, PDGFRB and MARCH1 genes (NCBI accession no. AC_000162.1, AC_000164.1 and AC_000163.1). PCR amplifications for pooled DNA were performed in a reaction volume of 20 μL comprising 2 μL of 50 ng/μL DNA, 1 μL of each primer, 10 μL of premix (containing dNTPs and DNA polymerase) (Tiangen, Beijing, China) and 6 μL of ddH 2 O. The amplification procedures were

Linkage disequilibrium (LD) analysis and haplotype construction
Hardy-Weinberg equilibrium was tested on each identified SNP using the chi-squared test at a P-value cut-off of 0.01. To estimate the extent of LD for the three genes, pairwise LD was measured among the SNPs of each gene based on the criterion of D' using the software Haploview 4.2 (Broad Institute of MIT and Harvard, Cambridge, MA, USA) [17]. Accordingly, haplotype blocks where SNPs were in high LD (D' > 0.90) were also determined based on confidence interval methods [18]. A haplotype with a frequency > 5% was treated as a distinguishable haplotype, and those haplotypes with relative frequency < 5% were pooled into a single group. Haplotype blocks within these SNPs were later employed to test their associations with the semen production traits in subsequent analyses.

Analyses of associations with semen production traits
Pedigree information of the resource population was traced back for three generations to construct the numerator relationship matrix. The associations of SNPs and haplotypes with the five semen production traits were evaluated using the mixed procedure in SAS 9.3 (SAS Institute Inc., Cary, NC). The model was performed as follows: where y ijk is the de-regressed EBVs; μ is the overall mean of de-regressed EBVs; G i is the fixed effect corresponding to the genotype of polymorphisms or haplotypes; α j is the random familial polygenic effect, distributed as N (0, A σ 2 a ), with the polygenic relationship matrix A and the additive genetic variance σ 2 a ; and e ijk is the random residual, distributed as N (0, Iσ 2 e ), with identity matrix I and residual error variance σ 2 e . Values at P < 0.05 were considered significant while values at P < 0.01 were regarded as highly significant. The differences among the effects of three genotypes on each SNP or haplotype were compared using multiple t-test with Bonferroni correction. In addition, the Bonferroni-corrected significance levels of 0.05/3 = 0.0167 and 0.01/3 = 0.0033 were used for comparison of the three genotypes. For the haplotypes, Bonferroni-corrected significance levels of 0.05/N and 0.01/N were used, where N represents the number of formed haplotypes in a block. Moreover, the additive (a), dominance (d) and allele substitution (α) were calculated according to the equation proposed by Falconer & Mackay [19], namely, a = (AA − BB)/2; d = AB -(AA + BB)/2; and α = a + d *(p − q). Here, AA and BB are the genotype frequencies of the two homozygotes; AB is the heterozygous genotype frequency; and p and q are the allele frequencies at the corresponding locus.
The percentage of genetic variance accounted for by the significant i-th SNP was estimated according to the formula below [20]: where p i and q i are the allele frequencies for the significant i-th SNP, a 2 i is the estimated additive effect of the significant i-th SNP on the trait under analysis and σ 2 a is the additive genetic variance for the trait.
Gene expression assays of PDE3A, PDGFRB and MARCH1 genes To further confirm the potential functions of the PDE3A, PDGFRB and MARCH1 genes, we conducted gene expression analyses of different genotypes. Fresh semen samples were collected from ten fully genotyped bulls.
Fresh semen samples were carefully laid on a monolayer of 40% Percoll. Somatic cell contamination of the The last stage used for the dissociation curve was as follows: 95°C for 15 s, 65°C for 10 s and 97°C for 60 s. Quantitative realtime PCR analysis of each gene was performed in triplicate and the relative gene expression was normalized using GAPDH by the 2 −ΔΔCt method, as described previously [21].
To further detect the effect of variants of significantly associated genes, the mRNA expression of sperm cells with different genotypes at sites of functionally important mutations was analyzed. The results of mRNA expression were analyzed by the GLM procedure in SAS 9.3 software.

Identification of SNPs in PDE3A, PDGFRB and MARCH1
Sequence analysis revealed that a total of 20 SNPs were detected using the pooled DNA of 50 bulls. Of those, eight SNPs were located in PDE3A, being distributed in exons (n = 2), introns (n = 3) and the 3′ untranslated region (n = 3). In addition, five SNPs were located in MARCH1, being distributed in the promoter (n = 3), an exon (n = 1) and an intron (n = 1). Furthermore, seven SNPs were located in PDGFRB, being distributed in exons (n = 2), introns (n = 2) and downstream of the gene (n = 3). The identified SNPs were then subjected to genotyping in 730 bulls. However, not all individuals were successfully genotyped at all SNPs. The number of remaining individuals as well as genotype frequencies, allele frequencies, primers for amplification and results of the chi-squared tests for each SNP are shown in Additional file 1: Table S1. Only one SNP, rs109116577, was a missense mutation, causing an amino acid change of Pro/Ser in the PDE3A protein. Three SNPs (rs456212302, rs42393923 and rs378918630) that were not in Hardy-Weinberg equilibrium were excluded from subsequent analysis (P < 0.01).

Association between SNPs and semen production traits in Chinese Holstein bulls
Association studies showed that four statistically significant SNPs associated with at least one semen production trait. The estimated effects of the four significant SNPs on semen production traits are shown in Table 4.
Three SNPs (rs211260176, rs208093284 and rs43445726) located in the promoter of MARCH1 were significantly associated with SVPE (P = 0.0246); SVPE (P = 0.0341); and all three of SVPE (P = 0.0001), NSPE (P = 0.0203) and NMSPE (P = 0.0189), respectively. In addition, the results showed that homozygous genotypes of all the significant SNPs were dominant for semen production traits. The dominant, additive and allele substitution effects of the significant SNPs on the target semen production traits are presented in Table 5.

LD among identified SNPs and haplotype association results
Pairwise D' measures between genotyped SNPs of the three genes were investigated and the inferred haplotype blocks are shown in Fig. 1. For PDE3A, one block consisting of three SNPs (rs110167512, rs42393928 and rs42393903) was inferred and three haplotypes were formed in the studied population. Association analysis  P-value is the significance level from analyses of the association of SNPs with semen production traits. **: P < 0.01; *: P < 0.05. Different superscript letters (lower-case letters: P < 0.05; upper-case letters: P < 0.01; Bonferroni-adjusted value after multiple testing) refer to significant differences among the genotypes. %Var indicates the percentage of genetic variance explained by the significant SNPs for traits   P-value is the significance level from analyses of the association of haplotypes with semen production traits. **: P < 0.01. Superscript letters (P < 0.05; Bonferroniadjusted value after multiple testing) refer to a significant difference among the genotypes revealed that haplotypes of PDE3A did not reach significance for five semen production traits. For PDGFRB, four haplotypes in one block with seven SNPs did not significantly associate with semen production traits. For MARCH1, three SNPs (rs211260176, rs208093284 and rs43445726) constituted a block in the studied population. The main haplotypes of TTT (H1), CCT (H2) and CCC (H3) accounted for frequencies of 57.2%, 20.7% and 19.5% of the total, respectively. Haplotype association study of MARCH1 demonstrated a significant association with SVPE (P = 0.0013) ( Table 6).
Functional prediction of the allele-dependent transcription factor binding sites The mutations in the regulatory regions of a gene can affect the transcription rate by changing the transcription factor binding sites [22]. Therefore, three significant SNPs located in the promoter of MARCH1 may be involved in altered transcription factor binding sites and may subsequently lead to gene expression differences. Sequences including the significant SNPs (21 bp) were subjected to a comparison with the reference transcription factor binding sites in the JASPAR CORE Vertebrata database (http://jaspar.genereg.net/cgi-bin/jaspar_db.pl), using a relative profile score threshold of 85%. As a result, the three regulatory SNPs were predicted to create some new transcription binding sites via the substitution of C to T. The details of this are shown in Fig. 2.

Expression regulation of the mutations in the MARCH1 Gene
The mRNA expression of the three genes was determined by quantitative real-time PCR and normalized using internal GAPDH expression in sperm cells. MARCH1 had a higher expression than the other two genes in sperm cells (Fig. 3a). To investigate the potential regulatory role of SNPs in the 5′ regulatory region, association analyses between different genotypes and MARCH1 expression level were conducted. We found that the two genotypes (CC was not observed) of rs43445726 were associated with a significant difference in gene expression of MARCH1 (P = 0.0035) (Fig. 3d). Similarly, another two SNPs of MARCH1 showed a tendency for genotypespecific gene expression (both P = 0.1918) (Fig. 3b, c).

Discussion
In support of our previous GWAS, we provide further evidence for the significant genetic effects of the PDGFRB and MARCH1 genes on semen production traits in another population of Chinese Holstein bulls. The PDE3A enzyme, mainly located in the postacrosomal segment of the sperm head [23], has been reported to exhibit activity in catalyzing cAMP into 5′ AMP in spermatozoa. cAMP is an important secondary messenger in the control of sperm functions, encompassing activation of motility, the acrosome reaction, hyperpolarization of sperm plasma membrane and ATP analysis [24]. In normal spermatozoa, PDE3A activity is inhibited by cGMP, maintaining a high cAMP level [25]. PDE3A is only 0.1 Mb away from another SNP found to be significantly associated with sperm motility in our previous GWAS. It is also close to an SNP (1.8 Mb away) shown to be significantly associated with sperm motility in a Polish Holstein bull population [8], and about 2.3 Mb away from an SNP significant for sire conception rate in another GWAS [26]. However, we did not observe significant associations between mutations in PDE3A and five semen production traits. We suspected that PDE3A might have genetic effects on semen production traits, but the causal mutation has not been detected in the studied population. Further studies should be conducted to confirm the genetic effects of PDE3A on semen production traits.
One SNP situated downstream of PDGFRB was found to have significant associations with SVPE, NSPE and NMSPE, which had already been declared to be significantly associated with SVPE in our previous GWAS. The mRNA of PDGFRB has been detected in gonocytes [27], Leydig and Sertoli cells, but not round spermatids Fig. 2 The changed transcription factor binding sites due to allele substitution of significant SNPs in the regulatory region of MARCH1. The significant SNPs in sequences are highlighted in red. The red dotted lines indicate the predicted transcription factor binding sites or primary spermatocytes [28]. Similarly, we also did not detect the expression of PDGFRB in sperm cells. Gonocytes, the precursors of spermatogonial stem cells, are located in the center of the seminiferous tubules. At a defined species-specific period of time (cattle: prepuberty), their proliferation and migration to the basement membrane give rise to spermatogonial stem cells, which maintain spermatogenesis in the mature testis [29]. The PDGFRB protein has been reported to play a leading role during the proliferation and migration of gonocytes [27]. The inhibition of PDGFRB tyrosine kinase activity was shown to reduce testis size, delay the initiation of spermatogenesis and thus provoke a drastic reduction of epididymal sperm count [27].
MARCH proteins are ubiquitin ligases and target glycoproteins for lysosomal destruction via ubiquitination of the cytoplasmic tail. Unlike the above two candidate genes, the functions of MARCH1 are seldom analyzed in relation to spermatogenesis. However, previous studies identified that three MARCH family members, MARCH-XI, MARCH10 and MARCH7, are highly expressed in developing spermatids [30][31][32]. The MARCH-XI protein is postulated to be a ubiquitin ligase that mediates transmembrane glycoproteins in the trans-Golgi network and multivesicular body transport pathway, which is associated with acrosomal formation in developing spermatids [30]. MARCH10 is abundantly expressed in elongated spermatids. Furthermore, immunohistochemical analysis of MARCH10 proteins revealed that they are predominantly located in the cytoplasmic lobes and the principal piece of the flagella. It is supposed that MARCH10 proteins are synthesized in the cytoplasm and then transported to the developing flagella [31]. Similarly, MARCH7 proteins that are localized to the acroplaxome and flagella mediate K48-linked ubiquitination in the acrosome/acroplaxome region and may be related to the regulation of head shaping and flagellar formation in developing spermatids [32]. In our study, the high expression of MARCH1 in spermatozoa and significant effects of MARCH1 on SVPE, SNPE and SNMPE support the asumption that MARCH1 functions in spermatogenesis, as MARCH-XI, MARCH10 and MARCH7 do.
In the present study, alleles involving a substitution of C to T in the regulatory region of MARCH1 were predicted to add a series of transcription factor binding sites. Specifically, one of the added transcription factors, RHOXF1, is encoded by an X-linked reproductive homeobox gene and has been observed to be specifically expressed in testis, especially in pachytene spermatocytes and round spermatids [33]. RHOXF1 critically upregulates many genes in male reproduction and may also modulate the transcription of MARCH1 [34]. Furthermore, one of the significant SNPs, rs43445726, resulted in a marked difference in expression level, with expression associated Fig. 3 Normalized mRNA expression of PDE3A, PDGFRB and MARCH1 genes and associations between significant SNPs and the expression level of MARCH1 in semen samples. a Normalized mRNA expression of PDE3A, PDGFRB and MARCH1 genes in semen samples. b SNP of MARCH1, rs211260176: no significant difference of expression levels was detected among CC (n = 2), CT (n = 5) and TT (n = 3), P-value: 0.1918. c SNP of MARCH1, rs208093284: no significant difference of expression levels was detected among CC (n = 2), CT (n = 5) and TT (n = 3), P-value: 0.1918. d SNP of MARCH1, rs43445726: extremely significant difference of expression was detected between TC (n = 2) and TT (n = 8), P-value: 0.0035 with the CT genotype being seven times that for the TT genotype. As for the phenotypes, CT was associated with lower sperm volume and sperm number than TT, which reflected the negative effect of MARCH1 on SVPE, NSPE and NMSPE. Furthermore, it was predicted that the SNP rs43445726 explained 12.22% of the genetic variance of SVPE, implying significant genetic effects of this mutation. However, the detected effects may be limited to the specific population studied here, so further analysis should be conducted to reveal the function of MARCH1 in spermatogenesis and verify the functional implications of its mutations.

Conclusion
Our findings demonstrated that PDGFRB and MARCH1 were significantly associated with semen production traits and presented the high expression of MARCH1 in mature sperm, which were consistent with previous GWAS and functional analyses. Our results not only provide new insight into the functions of the PDGFRB and MARCH1 genes, but also contribute useful information for marker-assisted selection or genome selection strategies of genetic improvement programs for semen production traits.

Additional file
Additional file 1: Table S1.