Two Mutations in the Caprine MTHFR 3'UTR Regulated by MicroRNAs Are Associated with Milk Production Traits

Abstract

Background

5,10-Methylenetetrahydrofolate reductase (MTHFR) plays a central role in folate metabolism by irreversibly converting 5,10-methylenetetrahydrofolate to 5-methylenetetrahydrofolate, a predominant circulating form of folate. Folate is reportedly important for milk protein synthesis, and MTHFR may be a key regulatory point of folate metabolism for milk protein synthesis in mammary epithelial cells. Prior to this study, polymorphisms of the MTHFR gene were not associated with milk production traits from a breeding perspective. Single nucleotide polymorphisms (SNPs) at microRNA (miRNA) binding sites (miR-SNPs) can affect gene expression. This study aimed to identify the effects of miR-SNPs (g.2244A>G and g.2264A>G) in the caprine MTHFR 3' UTR on the milk production traits of dairy goats.

Results

Guanzhong dairy (GD, n = 325) goats were used to detect SNPs in the caprine MTHFR 3' UTR by DNA sequencing. Two novel SNPs (g.2244A>G and g.2264A>G) were identified in the said region. The homozygous haplotype A-G of the SNPs g.2244A>G and g.2264A>G was significantly associated with milk yield and milk protein levels in GD goats (P < 0.05). Functional assays indicated that the MTHFR 2244 A → G substitution could increase the binding activity of hsa-miR-1266 with the MTHFR 3' UTR. The MTHFR 2264 A → G substitution could decrease the binding activity of hsa-miR-616 with the MTHFR 3' UTR. In addition, we observed a significant increase in the MTHFR mRNA levels of homozygous haplotype A-G carriers relative to those of homozygous haplotype G-A carriers. These results indicated that both SNPs altered the MTHFR mRNA levels. These altered levels of MTHFR mRNA may account for the association of SNPs with milk production traits.

Conclusions

This study is the first to report that the g.2244A>G and g.2264A>G polymorphisms were associated with milk production traits in GD goats. Further investigations should explore the underlying miRNA-mediated mechanisms that are modified by the g.2244A>G and g.2264A>G SNPs. The current study evaluated these SNPs as potential genetic markers in goats, with potential applications in breeding programs.

About the Authors:

Jinxing Hou

Contributed equally to this work with: Jinxing Hou, Xiaopeng An

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Xiaopeng An

Contributed equally to this work with: Jinxing Hou, Xiaopeng An

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Yuxuan Song

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Teyang Gao

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Yingnan Lei

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Binyun Cao

* E-mail: [email protected]

Affiliation: College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, 712100, P.R. China

Introduction

5,10-Methylenetetrahydrofolate reductase (MTHFR) is an essential enzyme for homocysteine and folate metabolism; this enzyme is involved in one-carbon metabolism and responsible for the final step in the conversion of dietary forms of folate to 5-methyltetrahydrofolate [1–3]. The caprine MTHFR gene is located on chromosome 16 and has 13 exons spanning 14.639 kb [4]. Analysis of the genetic variability of MTHFR has revealed its metabolic and functional aspects that are relevant to understanding the biology of folate metabolic pathways. The folate metabolism genes (such as MTHFR) play critical roles in regulating milk protein synthesis [5]. To the best of our knowledge, the association of MTHFR gene polymorphisms with milk production traits has not been reported from an animal breeding perspective.

MicroRNAs (miRNAs) are small RNAs of approximately 21 nucleotides. MiRNA can bind to the 3'untranslated region (3' UTR) of its target mRNA to post-transcriptionally regulate gene expression. Studies have shown that miRNAs play important roles in several biological processes, including embryonic development, cell proliferation and differentiation, apoptosis, fat metabolism, atherosclerosis and oncogenesis [6]. In previous studies, several miRNAs (e.g., miR-200a and miR-21) were found to be involved in the cellular differentiation of mammary glands [7,8].

SNPs are the most abundant form of DNA variation in the animal genome. Brodersen and Voinnet [9] showed that SNPs in miRNA binding sites can affect miRNA-induced genetic repression. SNPs that affect miRNA binding to target genes are called miR-SNPs [10]. Liu et al. [11] reported that the human SNP rs3735590 C → T influences miR-616 binding to the PON1 gene, thereby increasing the risk of ischemic stroke and carotid atherosclerosis. Clop et al. [12] demonstrated that the GDF8 allele of Texel sheep is characterised by a G → A transition in the 3' UTR that creates a target site for miR-1 and miR-206, which are highly expressed in skeletal muscle; consequently, translational inhibition of the myostatin gene occurs, thereby contributing to muscular hypertrophy in Texel sheep. Gao et al. [13] showed that the SNP g.1536 C>T in the TNP2 3' UTR alters the binding of TNP2 with bta-miR-154, and is associated with the semen quality of Chinese Holstein bulls. Given the regulatory role of miRNAs in gene expression, miR-SNPs may function as promising markers for milk production traits.

Based on the aforementioned considerations, we detected polymorphisms in the MTHFR gene of Guanzhong dairy (GD) goats by DNA sequencing, and investigated the associations between these genetic markers and milk production traits. After evaluating the associations between the candidate miR-SNPs and phenotypes of interest, we conducted reporter assays to confirm the effects of these miR-SNPs.

Materials and Methods

Ethics statement

The Animal Ethical Committee at the Northwest A & F University approved the experimental procedures.The animals were reared in GuanzhongDairy Goat Breeding Base in Zhouzhi county of Shaanxi province (34°14'N, 113 108°37'E and 1000 m altitude).

Animals and genomic DNA isolation

Blood samples were obtained from 325 GD goats. The goats were reared under the same standard conditions with dry-lot nutrition. All animals included in the study descended from 13 GD goat sires. Data on health, fertility and milk production were recorded by dairymen and veterinarians. The milk production traits were recorded from 2008 to 2013. For milk yield analysis, 325 GD goats were milked twice a day with a milking machine throughout the entire lactation period, and the milk weight was measured with an electronic scale. According to Luo [14], milk yield from the first to third lactation periods were standardised to 300 d in milk. For milk constituent analysis, a total of 325 milk samples were collected from the 325 GD goats once a month throughout the full lactation period. The first sampling occurred at least 20 d after parturition to exclude the risk of contamination with colostrum. The milk components (protein and fat) were determined with an ultrasonic S60SEC milk analyser (Milkotronic Company, Nova Zagora, Bulgaria). The consecutive milk records of the same goat were used in the study. For SNP analysis, 5 mL of blood was aseptically collected from the jugular vein of each individual (325 GD goats) and then individually stored in tubes containing the anticoagulant ACD (citric acid: sodium citrate: dextrose, 10:27:38).All samples were transported to the laboratory on ice. Genomic DNA was extracted from white blood cells by standard phenol‒chloroform extraction. All experiments were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.

SNP investigation and genotyping

The caprine MTHFR gene (GenBank Accession Nos. NC_022308 and XM_005690674) was used to design one pair of primers for amplifying the 3' UTR of caprine MTHFR (326bp). The exact position of amplicon in the 3' UTR of caprine MTHFR is in 2185bp-2510bp region (GenBank accession No. XM_005690674). Their optimal annealing temperatures are shown in Table A in S1 File. Subsequently, we screened the samples to identify SNPs of MTHFR by pooled DNA sequencing [15]. One DNA pool was created by mixing 5 μL of DNA (100 ng/μL) from each of the 325 samples. The pooled DNA as a template was used to amplify the 3' UTR of caprine MTHFR. A PCR Master Mix (2×) Kit (Thermo Scientific, NY, USA) was used for PCR amplification. The 25 μL PCR reactions contained 50 ng of pooled genomic DNA, 12.5 μL of 2× reaction mix (including 500 μM of each dNTP; 20 mMTris–HCl, pH 9; 100 mMKCl; 3 mM MgCl2), 0.5 μM of each primer and 0.5 U of Taq DNA polymerase (Thermo Scientific, NY, USA). The cycling protocol involved initial denaturation at 95°C for 5 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 50°C for 30 s and extension at 72°C for 35 s, with a final extension at 72°C for 10 min. Ten PCR products amplified with pooled DNA as a template were sent to the Beijing Genomics Institute (Beijing, China) for Sangersequencing in both directions.SNP discovery was conducted using the Chromas version 2.31 and DNAstar version 7.0 software programs. SNP nomenclature and numbering followed the recommendations of http://www.hgvs.org/mutnomen/recs-DNA.html#number. If SNPs were found in the 3' UTR of caprine MTHFR. The 325 samples were amplified independently, and then the 325 PCR products were sequenced independently to detect the different individual genotypes.

Statistical analysis

Polymorphism information content (PIC) represents the probability that a given offspring gets a certain markeralle from the same markeralle of its father or mother. PIC used to estimate the polymorphism of marker locus.PIC is calculated as follows:where Pi and Pj are the allelic frequencies of the ith and jth alleles, respectively, and n is the number of alleles [16].The allele frequencies, heterozygosity (He) and PIC were calculated using POPGENE (version 1.31; http://www.ualberta.ca/~fyeh/popgene_download.html). Association analysis between the mutations of the MTHFR genotypes and milk yield, as well as the fat and protein contents, was performed using repeated measures in the general linear model procedure of SPSS version 16.0 statistical software (http://en.softonic.com/s/spss-16-software). The mixed linear model was described by: Yikm = μ + Gi + Nk + Eikm, where Yikm is the trait measured in each of the ikmth female goat, μ is the overall population mean, Gi is the fixed effect associated with the ith genotype, Nk is the fixed effect associated with the kth female goat’s litter size and Eikm is the random error. Multiple comparisons of the mean milk yield, as well as mean fat and protein contents, for different genotypes were performed via Fisher’s least significant difference (LSD) test by SPSS version 16.0. Fisher’s LSD test is initially similar to the Bonferroni multiple comparison test. It takes the square root of the Residual Mean Square from the analysis of variance and considers that to be the pooled standard deviation. Taking into account the sample sizes of the two groups being compared, it computes a standard error of the difference between those two means. Then it computes a t ratio by dividing the difference between means by the standard error of that difference. The effects associated with the farm, birth year and season of birth were not included in the linear model because preliminary statistical analyses indicated that these effects did not significantly influence the variability of traits in the analysed populations.

Bioinformatics analysis of MTHFR 3'UTR

Bioinformatics analysis was used to predict the effects of mutations in the 3'UTR of MTHFR gene on the miRNA binding sites. The analysis was performed with five bioinformatics tools, namely, TargetScan(http://www.targetscan.org/), MicroInspector(http://bioinfo1.uni-plovdiv.bg/cgi-bin/microinspector/), miRWalk(http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/index.html), RNAhybrid(http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/submission.html) and Segal Lab(http://genie.weizmann.ac.il/pubs/mir07/mir07_prediction.html). The combination of these approaches was considered to reduce the possibility of false positives.

RNA isolation and real-time quantitative PCR (RT-qPCR)

To detect the correlation between expression levels of the MTHFR mRNA and SNPs (g.2244A>G and g.2264A>G) in vivo, mammary alveolar was stripped from 31 mammary gland tissues with different haplotypes (G-G, n = 12; G-A, n = 5; A-G, n = 4;A-A, n = 10) for RNA extraction. Total RNA was extracted from these samples using TRIzol Reagent (Invitrogen, Carlsbad, USA). The extracted RNA was homogenised for RT-qPCR. The quantity and integrity of each RNA sample were assessed with an Agilent 2100 Bioanalyser (Agilent Technologies, USA). A total of 31 RNA samples, which included the homozygous haplotypes G-G, G-A, A-G and A-A, were subjected to reverse transcription using a cDNA High Capacity Kit (Invitrogen, Carlsbad, USA). RT-qPCR analysis with the SYBR Green PCR Master Mix (Takara, Dalian, China) was performed on a CFX Connect Real-time PCR Detection System (Bio-Rad, CA, USA). The relative expression levels of objective mRNAs were calculated using the ∆∆Ct method [17]. The primers of the MTHFR, UXT, RPS9 and RPS15 genes are described in Table B in S1 File. The fold changes were normalised by the expression levels of UXT, RPS9 and RPS15, and each assay was performed in triplicate [18, 19].

Cell culture

Goat mammary epithelial cells (GMECs) were obtained from mammary gland biopsy of four 2-year-old (50 days post-parturition) Guanzhong dairy goats. Four genotype cells (GG-GG, AA-GG, GG-AA and AA-AA) were directly isolated from Guanzhong dairy goats with these genotypes. After mammary tissues were surgically removed from the dairy goats, they were placed in a sterile ice-cold D-Hank’ssolution supplemented with 300 U/mL penicillin and were immediately transported to the laboratory. The mammary tissues were trimmed of visible fat and connective tissues and washed with D-Hank’s solution several times until the solution became pellucid and devoid of milk. Tissue sections were then minced into about 1 mm3cubes and rinsed again with D-Hank’s solution and kept up for 10 min at room temperature. The smaller pieces of tissues were put onto the cell dishes coated with collagen and were reversely incubated at 37°C with saturated humidity and 5% CO2. After 30 min,2 mL culture medium was added to the culture dish. The culture medium consisted of DMEM-F12 supplemented with 10% FBS, insulin (5 μg/ml), hydrocortisone (0.25 μmol/L), penicillin (100 U/mL), streptomycin (100μg/mL) and epidermal growth factor 1 (10 ng/ml). The medium was replaced with fresh medium every 48 h until the cells were migrated out of the tissue and visibly spread across the bottom of the dish. Being cultured for about five passages, GMECs were obtained according to Lin et al. (2013) [20].

Luciferase reporter assays

To construct the luciferase reporter plasmids containing the MTHFR 3'UTR, fragments of the MTHFR 3' UTR with different haplotypes (A-G, G-A, A-A and G-G) were amplified by PCR. The PCR products of the different haplotypes (A-G, G-A, A-A and G-G) were isolated and separated by agarose gel electrophoresis. These products were then linked to a pMD19-T vector with a TA Cloning Kit (Invitrogen, CA, USA). The recombinant pMD19-T vectors with different haplotypes were digested by the XhoI and NotI endonuclease enzymes. Finally, the digested products with different haplotypes were inserted between the Renilla and firefly luciferase genes in a psiCHECK-2 vector (Promega, WI, USA).The plasmids containing the A-G, G-A, A-A and G-G haplotypes were obtained and then confirmed by sequencing. For the luciferase reporter assay, GMECs were placed in 12-well plates (5 × 105 cells per well) and then cotransfected with psiCHECK-2 vectors containing the 3'UTR-haplotypes A-G, G-A, A-A and G-G. The mimics of miRNAs (hsa-miR-1266, hsa-miR-1289, hsa-miR-505 and hsa-miR-616) and their negative controls (GenePharma, Shanghai, China) were cotransfected with the reporter plasmids to a final concentration of 20 nmol/μL. At 36 h after transfection in GMECs, Renilla luciferase activity in lysates was measured with the Dual-Luciferase Reporter Assay System (Promega, WI, USA).The results were normalised against firefly luciferase activity. Assays were performed according to the manufacturer’s suggestions. Each experiment was independently performed three times, and each sample was evaluated in triplicate. The Mann–Whitney test was used to compare data, and a two-sided P value of 0.05 was considered significant.

Results

SNP identification and genotyping

Two SNPs (g.2244A>G and g.2264A>G) were identified in the caprine MTHFR 3'UTRby DNA sequencing (Figs 1 and 2). The mutant sequence of MTHFR 3'UTR was submitted to NCBI (GenBank Accession No. KM243021). Four alleles of the g.2244A>G and g.2264A>G SNPs introduced several different miRNA sites (S1 Fig). The PIC was 0.37 at the g.2244A>G and g.2264A>G loci, respectively (Table C in S1 File). The genotypic distribution and allele frequencies of both SNPs are shown in Table C in S1 File. The g.2244A>G and g.2264A>G loci were in Hardy–Weinberg disequilibrium (P < 0.05) (Table C in S1 File). To reveal the linkage relationships between both SNPs, the linkage disequilibrium was estimated in the GD goats (Table C in S1 File). When r2 > 0.33, the linkage disequilibrium was considered strong [21]. Thus, the g.2244A>G and g.2264A>G loci were closely linked in the GD goats.

[Figure omitted. See PDF.]

Fig 1. Sequencing chromatograms of different genotypes at g.2244A>G locus.

https://doi.org/10.1371/journal.pone.0133015.g001

[Figure omitted. See PDF.]

Fig 2. Sequencing chromatograms of different genotypes at g.2264A>G locus.

https://doi.org/10.1371/journal.pone.0133015.g002

Association analysis of SNPs with milk production traits

Individual GD goats with AA and AG genotypes had greater milk yield than those with the GG genotype at the g.2244A>G locus (P < 0.05; Table D in S1 File). In addition, individuals with the AA genotype at the g.2244A>G locus had higher milk protein levels than those with the AG and GG genotypes (P < 0.05). The combined effects of the two SNPs (g.2244A>G and g.2264A>G) showed that individuals with the AA-GG (homozygous haplotype A-G) genotype had higher milk yield than individuals with the AA-AA (homozygous haplotype A-A), GG-AA (homozygous haplotype G-A), GG-GG (homozygous haplotype G-G) and AG-AG genotypes (P < 0.05; Table 1). In addition, individuals with the AA-GG (homozygous haplotype A-G) genotype had higher milk protein content than individuals with the AG-AG and GG-GG (homozygous haplotype G-G) genotypes (P < 0.05).

[Figure omitted. See PDF.]

Table 1. Combined effect of two loci on milk production traits (means ± standard errors) in goats.

https://doi.org/10.1371/journal.pone.0133015.t001

Effects of SNPs on miRNA binding ability

The miRNA target site predictions showed that hsa-miR-1266, hsa-miR-1289, hsa-miR-505 and hsa-miR-616 could bind to the 3'UTR region of the MTHFR gene (S2 Fig). To understand the functional significance of the g.2244A>G and g.2264A>G substitutions, we used a luciferase reporter assay (Fig 3A). Compared with the haplotype 2244A-2264G and negative control (NC), hsa-miR-1266 suppressed luciferase expression in the presence of the haplotype 2244G-2264A (P < 0.05; Fig 3B). Additionally, hsa-miR-1289 significantly suppressed luciferase expression in the presence of the haplotype 2244G-2264A compared with NC (P < 0.05), but the difference with the haplotype 2244A-2264G was not significant (P > 0.05; Table E in S1 File). The presence of hsa-miR-616 efficiently reduced luciferase expression in the 2264A-2244A-containing plasmid compared with that in the 2264G-2244G-containing plasmid (P < 0.05). However, hsa-miR-505 did not significantly affect luciferase expression (P > 0.05; Table F in S1 File). These results suggested that the haplotypes 2244G-2264A and 2264A-2244A had lower MTHFR protein expression levels than the haplotype 2244A-2264G and 2264G-2244G.

[Figure omitted. See PDF.]

Fig 3. Characterization and functional analysis of MTHFR 3′UTR.

(A) Double fluorescein enzyme expression vector linked with MTHFR 3′UTR. (B) GMECs cells seeded on 12-well plates were transiently cotransfected with psiCHECK-2-MTHFR 3′UTR and miRNA mimics (hsa-miR-1266, hsa-miR-1289, hsa-miR-505 and hsa-miR-616) or stable negative control (NC). Luciferase activity was measured after 36 h of transfection. Renilla luciferase was normalized with firefly luciferase activity, and the mean activities ± standard deviation from 3 independent experiments are shown. *P<0.05.

https://doi.org/10.1371/journal.pone.0133015.g003

In this study, we performed RT-qPCR assays to explore whether hsa-miR-1266 and hsa-miR-616 affect MTHFR expression by degrading mRNA or suppressing mRNA post-translational translation (Fig 4). A total of 31 mammary gland samples with different haplotypes for the MTHFRg.2244A>G and g.2264A>G polymorphisms were used to assess the expression of MTHFR mRNA. Compared with individuals with homozygous haplotype G-G, A-G and A-A, the individuals with homozygous haplotype G-A had the lowest mRNA expression. By contrast, individuals with homozygous haplotype A-G had the highest mRNA expression among all the samples (Table G in S1 File). Therefore, hsa-miR-1266 and hsa-miR-616 may affect MTHFR mRNA expression by degrading mRNA.

[Figure omitted. See PDF.]

Fig 4. Comparison of mRNA expression levels of caprine MTHFR in mammary glands among four haplotypes.

The number of animals of homozygous haplotypes G-G, G-A, A-G and A-A were 12, 5, 4 and 10, respectively. The fold change was normalized against UXT, RPS9and RPS15genes. The different small letters represent significant difference at 5% level.

https://doi.org/10.1371/journal.pone.0133015.g004

Discussion

PIC is related to the use efficiency and selection potential of genetic markers; the greater polymorphism information content values and heterozygous ratio are, the greater the potentials of genetic markers are, and their effects are better for animal genetic breeding [16]. PIC values are classified as low polymorphism when PIC < 0.25, moderate polymorphism when 0.25 < PIC < 0.50 and high polymorphism when PIC > 0.50). The GD goats have moderate genetic diversity at the g.2244A>G and g.2264A>G loci, so the both loci have a better potential as genetic markers.

The 3'UTR polymorphisms may be within or at the vicinity of the miRNA binding site; these polymorphisms may interfere with miRNA function and lead to differential gene expression, thereby affecting the animal phenotype [22]. Previous reports showed that such polymorphisms influence MTHFR gene expression [23, 24]. The human MTHFR rs4846049 (G>T) is a potential candidate SNP that modulates miRNAs; the MTHFR mRNA complex has the greatest change in the binding free energy for hsa-miR-149; based on luciferase analysis, hsa-miR-149 inhibited the activity of the reporter vector carrying the T allele, but not the G allele [25]. The effects of SNPs in the 3' UTR of MTHFR on animal production traits have remained largely unknown. In the present study, the results of our integrative analysis showed that the g.2244A>G and g.2264A>G SNPs may change the binding activity of miRNAs (hsa-miR-1266 and hsa-miR-616) with the MTHFR 3'UTR. We used a bioinformatics analysis method similar to that of Sun et al. [26], who reported that the transcribed target SNPs associated with schizophrenia alter the binding activity of miRNAs with the 3'UTR of mRNA. Therefore, we speculated that the SNPs g.2244A>G and g.2264A>G accounted for differences in the milk production traits of GD goats because hsa-miR-1266 and hsa-miR-616 could change the expression levels of the MTHFR gene.

The results of the association analysis showed that individuals with the homozygous haplotype A-G had higher milk yield and protein content than those with other haplotypes. Functional assays indicated that the MTHFR 2244 A → G substitution could increase the binding activity of hsa-miR-1266 with the MTHFR 3'UTR, whereas the MTHFR 2264 A → G substitution could decrease the binding activity of hsa-miR-616 with the MTHFR3'UTR. These data were consistent with results of in silico analysis, thereby suggesting the functional interaction between hsa-miR-1266 or hsa-miR-616 and the MTHFR mRNA 3' UTR. We further examined the MTHFR mRNA levels in GD goats with different haplotypes. A significant increase in the MTHFR mRNA levels of the homozygous haplotype A-G carriers was observed relative to those of the homozygous haplotype G-A carriers. These findings consistently suggested that the g.2244A>G and g.2264A>G SNPs altered the MTHFR expression level, and they were associated with milk production traits. However, this study had several limitations. Firstly, the implied regulation of MTHFR expression by hsa-miR-1266 and hsa-miR-616 is based on preliminary results. Further functional experiments are required to analyse the effects of this polymorphism on MTHFR function. Secondly, the influence of other miRNAs on mutant forms of the MTHFR gene must be considered.

In conclusion, our data suggested that the g.2244A>G and g.2264A>G polymorphisms were associated with milk production traits in GD goats. Future investigations should explore the underlying miRNA-mediated mechanisms that are modified by the SNPs g.2244A>G and g.2264A>G. This study evaluated these SNPs as genetic markers for goat genetics and breeding, with potential applications in breeding programs.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Effects of SNPs g.2244A>G (A) and g.2264A>G (B) on miRNA binding sites in the caprine MTHFR 3MTHFR.

The green boxes represent potential miRNA target sites.

https://doi.org/10.1371/journal.pone.0133015.s001

(TIF)

S2 Fig. Sequence of the heptamer motif for each miRNA target site.

Green represents miRNA. Red represents MTHFR 3THFR.

https://doi.org/10.1371/journal.pone.0133015.s002

(TIF)

S1 File. Analysis results of primer information, genotypic distribution, association of SNPs with milk production traits, luciferase activity, MTHFR mRNA expression levels and miRNA target sites.

Primer information of MTHFR gene (Table A). Primer information for real-time quantitative PCR (Table B). Genotypic distribution and allelic frequencies of two SNP loci in the caprine MTHFR gene (Table C). Association analysis of g.2244A>G and g.2264A>G loci with milk production traits (means ± standard errors) in Guanzhong dairy goats (Table D). Luciferase activity (means ± standard deviation) of the psiCHECK-2-MTHFR 3′UTR vector that include haplotypes 2244G-2264A and 2244A-2264G (Table E). Luciferase activity (means ± standard deviation) of the psiCHECK-2-MTHFR 3′UTR vector that include haplotypes 2264G-2244G and 2264A-2244A (Table F). Comparison of expression levels of caprine MTHFR mRNA in mammary glands among four haplotypes (Table G).

https://doi.org/10.1371/journal.pone.0133015.s003

(DOC)

Acknowledgments

This study was supported by the National Support Program of China (2011BAD28B05-3), the national sparking plan project (2013GA850003), the Science and Technology Innovation Project of Shaanxi Province (2011KTCL02-09), the PhD Research Funding of Northwest A&F University (2013BSJJ003) and China Province Postdoctoral Science Foundation (129869).

Author Contributions

Conceived and designed the experiments: BYC. Performed the experiments: JXH XPA. Analyzed the data: TYG YNL. Contributed reagents/materials/analysis tools: TYG YNL. Wrote the paper: JXH XPA YXS.

Citation: Hou J, An X, Song Y, Gao T, Lei Y, Cao B (2015) Two Mutations in the Caprine MTHFR 3'UTR Regulated by MicroRNAs Are Associated with Milk Production Traits. PLoS ONE 10(7): e0133015. https://doi.org/10.1371/journal.pone.0133015

References

1. Gaysina D, Cohen S, Craddock N, Farmer A, Hoda F, Korszun A, et al. No association with the 5,10-methylenetetrahydrofolate reductase gene and major depressive disorder: Results of the depression case control (DeCC) study and a meta-analysis. Am J Med Genet B. 2008; 147B(6):699–706. pmid:18165972.

2. Goyette P, Pai A, Milos R, Frosst P, Tran P, Chen ZT, et al. Gene structure of human and mouse methylenetetrahydrofolate reductase (MTHFR). Mamm Genome. 1998; 9(8):652–656. pmid:9680386

3. Senhaji N, Serbati N, Diakite B, Arazzakou S, Hamzi K, Badre W, et al. Methylenetetrahydrofolate reductase C677T variant in Moroccan patients with inflammatory bowel disease. Gene. 2013; 521(1):45–49. pmid:23542077.

4. Dong Y, Xie M, Jiang Y, Xiao NQ, Du XY, Zhang WG, et al. Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Capra hircus). Nat Biotechnol. 2013;31(2):135–141. pmid:23263233.

5. Menzies KK, Lefevre C, Macmillan KL, Nicholas KR. Insulin regulates milk protein synthesis at multiple levels in the bovine mammary gland. Funct Integr Genomic. 2009; 9(2):197–217. pmid:19107532.

6. Ambros V. The functions of animal microRNAs. Nature. 2004;431(7006):350–355. pmid:15372042.

7. Galio L, Droineau S, Yeboah P, Boudiaf H, Bouet S, Truchet S, et al. MicroRNA in the ovine mammary gland during early pregnancy: spatial and temporal expression of miR-21, miR-205, and miR-200. Physiol Genomics. 2013; 45(4):151–161. pmid:23269700.

8. Nagaoka K, Zhang HL, Watanabe G, Taya K. Epithelial cell differentiation degulated by microRNA-200a in mammary glands. PlOSONE. 2013;8(6): e65127. pmid:23750238.

9. Brodersen P, Voinnet O. Revisiting the principles of microRNA target recognition and mode of action. Nat Rev Mol Cell Bio. 2009;10(2):141–148. pmid:19145236.

10. Verbraak EJC, van 't Veld EM, Koerkamp MG, Roelen BAJ, van Haeften T, Stoorvogel W, et al. Identification of genes targeted by FSH and oocytes in porcine granulosa cells. Theriogenology. 2011; 75(2):362–376. pmid:20965553.

11. Liu ME, Liao YC, Lin RT, Wang YS, Hsi E, Lin HF, et al. A functional polymorphism of PON1 interferes with microRNA binding to increase the risk of ischemic stroke and carotid atherosclerosis. Atherosclerosis. 2013; 228(1):161–167. pmid:23497787.

12. Clop A, Marcq F, Takeda H, Pirottin D, Tordoir X, Bibe B, et al. A mutation creating a potential illegitimate microRNA target site in the myostatin gene affects muscularity in sheep. Nat Genet. 2006; 38(7):813–818. pmid:16751773.

13. Gao Q, Ju Z, Zhang Y, Huang J, Zhang X, Qi C, et al. Association of TNP2 gene polymorphisms of the bta-miR-154 target site with the semen quality traits of Chinese Holstein bulls. PLOS ONE. 2014; 9(1):e84355. pmid:24416221.

14. Luo MF.Study on milk yield correction coefficient in dairy goat. Animal Husbandry and Veterinary Medicine. 1989; 1:5–6. (In Chinese).

15. Bansal A, van den Boom D, Kammerer S, Honisch C, Adam G, Cantor CR, et al. Association testing by DNA pooling: An effective initial screen. Proc Natl Acad Sci U S A. 2002; 99(26):16871–16874. pmid:12475937.

16. Botstein D, White RL, Skolnick M, Davis RW. Construction of a genetic linkage map in man using restriction fragment length polymorphisms. Am J Hum Genet. 1980; 32(3):314–331. pmid:6247908.

17. Ji ZB, Wang GZ, Xie ZJ, Wang JM, Zhang CL, Dong F, et al. Identification of novel and differentially expressed microRNAs of dairy goat mammary gland tissues using Solexa sequencing and bioinformatics. PlOS ONE. 2012;7(11):e49463 pmid:23166677.

18. Bionaz M, Loor JJ. Identification of reference genes for quantitative real-time PCR in the bovine mammary gland during the lactation cycle. Physiological Genomics. 2007; 29(3):312–319. pmid:17284669.

19. Bionaz M, Loor JJ. ACSL1, AGPAT6, FABP3, LPIN1, and SLC27A6 are the most abundant isoforms in bovine mammary tissue and their expression is affected by stage of lactation. J Nutr. 2008; 138(6):1019–1024. pmid:18492828.

20. Lin XZ, Luo J, Zhang LP, Zhu JJ. MicroRNAs synergistically regulate milk fat synthesis in mammary gland epithelial cells of dairy goats. Gene Expr. 2013;16(1):1–13. pmid:24397207.

21. Ardlie KG, Kruglyak L, Seielstad M. Patterns of linkage disequilibrium in the human genome. Nat Rev Genet. 2002; 3(4):299–309. pmid:11967554.

22. Trott JF, Freking BA, Hovey RC. Variation in the coding and 3' untranslated regions of the porcine prolactin receptor short form modifies protein expression and function. Anim Genet. 2014;45(1):74–86. pmid:24350599.

23. Yamada K, Chen ZT, Rozen R, Matthews RG. Effects of common polymorphisms on the properties of recombinant human methylenetetrahydrofolate reductase. Proc Natl Acad Sci U S A. 2001; 98(26):14853–14858. pmid:11742092.

24. Huang T, Tucker KI, Lee YC, Crott JW, Parnell LD, Shen JA, et al. Methylenetetrahydrofolate reductase variants associated with hypertension and cardiovascular disease interact with dietary polyunsaturated fatty acids to modulate plasma homocysteine in Puerto Rican adults. J Nutr. 2011;141(4):654–659. pmid:21270364.

25. Wu C, Gong Y, Sun A, Zhang Y, Zhang C, Zhang W, et al. The human MTHFR rs4846049 polymorphism increases coronary heart disease risk through modifying miRNA binding. Nutr Metab Cardiovasc Dis. 2013;23(7):693–698. pmid:22647417.

26. Sun XH, Zhang JY. Identification of putative pathogenic SNPs implied in schizophrenia-associated miRNAs.BMC Bioinformatics. 2014;15:194. pmid:24934851.