Single-nucleotide polymorphisms g.151435C>T and g.173057T>C in PRLR gene regulated by bta-miR-302a are associated with litter size in goats

Theriogenology 83 (2015) 1477–1483

Contents lists available at ScienceDirect

Theriogenology journal homepage: www.theriojournal.com

Single-nucleotide polymorphisms g.151435C>T and g.173057T>C in PRLR gene regulated by bta-miR-302a are associated with litter size in goats Xiaopeng An, Jinxing Hou, Teyang Gao, Yingnan Lei, Guang Li, Yuxuan Song, Jiangang Wang, Binyun Cao* College of Animal Science and Technology, Northwest A&F University, Yangling, Shaanxi, P.R. China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 May 2014 Received in revised form 19 January 2015 Accepted 24 January 2015

Single-nucleotide polymorphisms (SNPs) located at microRNA-binding sites (miR-SNPs) can affect the expression of genes. This study aimed to identify the miR-SNPs associated with litter size. Guanzhong (n ¼ 321) and Boer (n ¼ 191) goat breeds were used to detect SNPs in the caprine prolactin receptor (PRLR) gene by DNA sequencing, primer-introduced restriction analysis-polymerase chain reaction, and polymerase chain reaction-restriction fragment length polymorphism. Three novel SNPs (g.151435C>T, g.151454A>G, and g.173057T>C) were identified in the caprine PRLR gene. Statistical results indicated that the g.151435C>T and g.173057T>C SNPs were significantly associated with litter size in Guanzhong and Boer goat breeds. Further analysis revealed that combinative genotype C6 (TTAACC) was better than the others for litter size in both goat breeds. Furthermore, the PRLR g.173057T>C polymorphism was predicted to regulate the binding activity of btamiR-302a. Luciferase reporter gene assay confirmed that 173057C to T substitution disrupted the binding site for bta-miR-302a, resulting in the reduced levels of luciferase. Taken together, these findings suggested that bta-miR-302a can influence the expression of PRLR protein by binding with 30 untranslated region, resulting in that the g.173057T>C SNP had significant effects on litter size. Ó 2015 Elsevier Inc. All rights reserved.

Keywords: Single-nucleotide polymorphism microRNA Primer-introduced restriction polymerase chain reaction PRLR gene

analysis-

1. Introduction Prolactin (PRL) is an anterior pituitary polypeptide hormone involved in many reproductive pathways and is essential for reproductive performance [1,2]. Biological effects of PRL are mediated by its interaction with the PRL receptor (PRLR). The prolactin receptor belongs to the same family as the growth hormone receptor and is part of the cytokine receptor superfamily characterized by their ability to activate JAK2 and three members of the Stat family, STATl, STAT3, and STAT5 [3]. PRLR gene was mapped on bovine and goat chromosome 20 [4]. There are two distinct * Corresponding author. Tel.: þ86 29 87092102; fax: þ86 29 87092164. E-mail address: [email protected] (B. Cao). 0093-691X/$ – see front matter Ó 2015 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.theriogenology.2015.01.030

PRLR isoforms: long isoform with 557 amino acids and short isoform with 272 amino acids produced by alternative splicing of the primary transcript [5]. Heterodimerization of different PRLR isoforms produces inactive complexes which might also be of importance in the physiological context because PRL target cells usually express more than a single PRLR isoform [6,7]. The long isoform of PRLR binds PRL and contributes to activation of JAK2 kinases and subsequent phosphorylation of STAT5 transcription factors which bind to recognition sequences located in promoters of milk protein genes [8,9]. Analysis of the genetic variability of the PRL and PRLR loci has shed light on metabolic and functional aspects that are relevant to understanding the biology of the PRL endocrine axis. From an animal breeding perspective, polymorphisms of

1478

X. An et al. / Theriogenology 83 (2015) 1477–1483

PRLR gene have been associated with reproductive traits in pigs and sheep [10,11] as well as with milk production traits in dairy cattle [12]. These characteristics make PRLR a strong candidate gene for reproductive traits. MicroRNAs (miRNAs) are small RNAs of approximately 21 nucleotides. An miRNA can bind to the 30 untranslated region (UTR) of its target messenger RNA (mRNA) to posttranscriptionally regulate genes. Studies have shown that miRNAs play a role in several biological processes, including embryonic development, cell proliferation and differentiation, apoptosis, fat metabolism, atherosclerosis, and oncogenesis [13]. In previous studies, several miRNAs (e.g., miR-17-5p and miR-143) are found to be involved in the follicular development and ovulation [14,15]. Singlenucleotide polymorphisms (SNPs) are the most abundant form of DNA variation in the animal genome. Brodersen and Voinnet [16] showed that SNPs within the binding site of miRNA can affect miRNA-induced genetic repression. If an SNP can influence the binding of an miRNA to its target gene, this SNP is called an miR-SNP [17]. In human, Liu et al. [18] reported that this SNP rs3735590 (C>T) influences miR-616 binding to the PON1 gene, increasing the risk of ischemic stroke and carotid atherosclerosis. Gao et al. [19] showed that the SNP g.1536 C>T in the TNP2 30 UTR, which altered the binding of TNP2 with bta-miR-154, was found to be associated with the semen quality traits of Chinese Holstein bulls. Given the regulatory role of miRNAs in gene expression, miR-SNPs could provide promising markers for reproductive traits. On the basis of above considerations, we detected the polymorphisms of PRLR gene in Guanzhong (GZ) and Boer (BE) goat breeds by DNA sequencing, primer-introduced restriction analysis-polymerase chain reaction (PCR), and PCRrestriction fragment length polymorphism and investigated the associations between these genetic markers and litter size. After evaluating the associations between the candidate miR-SNPs and phenotypes of interest, we conducted reporter assays to confirm the effects of the miR-SNPs.

2.2. SNP investigation and genotyping According to the caprine PRLR gene (GenBank accession No. NC_022312), three pairs of primers were designed to amplify the exon 2 and 30 UTR of caprine PRLR gene. Their optimal annealing temperatures and name to each DNA fragment are shown in Table S1. Herein, we screened them for identifying SNPs of this gene by a DNA pooling sequencing assay [20]. Five microliters of DNA (100 ng/mL) per sample was collected to create a DNA pool for each goat breed. Polymerase chain reaction products of different primer pairs were sent to Beijing Genomics Institute (Beijing, China) for sequencing in both directions. Discovery of SNPs was conducted using the Chromas 2.31 and DNAstar 7.0 software. Single-nucleotide polymorphism nomenclature and numbering followed http://www.hgvs.org/ mutnomen/recs-DNA.html#number. The 25-mL volume contained 50-ng genomic DNA, 12.5 mL of 2 reaction mix (including 500-mM dNTP each, 20-mM Tris–HCl, pH 9, 100-mM KCl; 3-mM MgCl2), 0.5 mM of each primer, and 0.5 units of Taq DNA polymerase. The cycling protocol was 5 minutes at 95  C, 35 cycles of denaturing at 94  C for 30 seconds, annealing at 50  C to 55  C (Table S1) for 30 seconds, extending at 72  C for 35 seconds, and with a final extension at 72  C for 10 minutes. DNA fragments I to III (5 mL) were mixed with 0.7-mL 10 buffer, 2.5-U restriction enzyme (NEB, Ipswich, Britain), and 3.8-mL sterilized ddH2O and then incubated for 1.5 hours at 37  C (Table S1). The restriction enzymes are shown in Table 1. DNA fragments I and II were used for the primer-introduced restriction analysis analysis. The SNP in DNA fragment III was genotyped using restriction fragment length polymorphism. Digestion products were subjected to 3.5% horizontal agarose gel electrophoresis or 12% polyacrylamide gel electrophoresis. The agarose and polyacrylamide gels were stained with ethidium bromide and 0.1% silver nitrate, respectively, and then the genotypes were observed.

2. Materials and methods

2.3. Statistical analysis

2.1. Animals and genomic DNA isolation

The allelic frequencies, heterozygosity and polymorphism information content (PIC) were calculated using PopGene (version 1.31) (http://www.ualberta.ca/wfyeh/ popgene_download.html). The linkage disequilibrium was performed by SHEsis online software (http://analysis.bio-x. cn/myAnalysis.php). Association analyses between PRLR genotypes for three SNPs and litter size were performed with SPSS 16.0 software (http://en.softonic.com/s/spss-16software). Multiple comparisons of the means are performed using the least significant difference method. Data were analyzed with the following mixed linear model for SNPs and traits: Yilm ¼ m þ Gi þ Sl þ Eilm, where Yilm is the trait measured on each of the ilmth animal, m is the overall population mean, Gi is the fixed effect associated with ith genotype, Sl is the random effect associated with the lth sire, and Eilm is the random error. Effects associated with farm, birth year, and season of birth are not matched in the linear model, as the preliminary statistical analyses indicated that these effects did not have a significant influence on variability of traits in the analyzed populations.

Blood samples were obtained from 512 goats belonging to both of the following breeds: GZ (n ¼ 321) and BE (n ¼ 191) goats. They were reared in Zhouzhi and Linyou county of Shaanxi province, respectively. All diets were based on alfalfa, corn silage, and a combination of concentrates including corn, soya meal, and bone meal. Health, fertility, and production records were maintained by the dairymen and veterinarians. The litter size from the first to fourth parity was obtained from production records. First, 5 mL of blood per individual was collected aseptically from the jugular vein and kept in a tube containing anticoagulant ACD (citric acid:sodium citrate:dextrose, 10: 27: 38), and then all samples were delivered back to the laboratory in an ice box. The genomic DNA was extracted from white blood cells using a standard phenol–chloroform extraction protocol. All experiments were performed in accordance with the National Institute of Health Guide for the Care and Use of Laboratory Animals.

X. An et al. / Theriogenology 83 (2015) 1477–1483 Table 1 Genotypic distribution and allelic frequencies of three single-nucleotide polymorphism loci in the caprine PRLR gene. Locus

Restriction Genetic structure enzyme

g.151435C>T

MwoI

g.151454A>G BslI

g.173057T>C

Genotype CC TC TT Allele C T He PIC Equilibrium c2 test Genotype AA GA GG Allele A G He PIC Equilibrium c2 test

Breed GZ

BE

33 185 103

25 106 60

0.39 0.61 0.58 0.36 P < 0.01

0.41 0.59 0.56 0.37 P ¼ 0.04

107 117 37

70 104 17

0.61 0.39 0.55 0.36 P < 0.01

0.64 0.36 0.54 0.36 P ¼ 0.01

RsaI

Genotype CC 245 CT 76 Allele C 0.88 T 0.12 He 0.24 PIC 0.19 Equilibrium c2 test P ¼ 0.02 LD of g.151435C>T and g.151454A>G r2 ¼ 0.34 LD of g.151435C>T and g.173057T>C r2 < 0.01 LD of g.151454A>G and g.173057T>C r2 < 0.01

133 58 0.85 0.15 0.30 0.22 P ¼ 0.01 r2 ¼ 0.36 r2 ¼ 0.01 r2 ¼ 0.02

Abbreviations: BE, Boer; GZ, Guanzhong; He, heterozygosity; LD, linkage disequilibrium; PIC, polymorphism information content.

2.4. Bioinformatics analysis of PRLR 30 UTR On the basis of bioinformatics analysis, we predicted the effects of mutation in the 30 UTR region of PRLR gene on miRNA-binding sites by using four common websites (MicroInspector: http://bioinfo1.uni-plovdiv.bg/ cgi-bin/microinspector/, miRWalk: http://www.umm. uni-heidelberg.de/apps/zmf/mirwalk/index.html, RNA hybrid: http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/ submission.html, and Segal Lab: http://genie.weizmann. ac.il/pubs/mir07/mir07_prediction.html; Fig. S1). We considered that the combination of these approaches would greatly reduce the possibility of false positive.

2.5. Ovarian RNA isolation and real-time quantitative PCR For detection of the correlation between the expression levels of PRLR mRNA and SNP g.173057T>C in vivo, a total of 53 goats with different genotypes (38 for CC and 15 for CT genotypes) were subjected to extraction of the total RNA using the Trizol Reagent (Invitrogen, Carlsbad, USA). The total RNAs were homogenized and pooled for real-time quantitative PCR (RT-qPCR). The quantity and integrity of the total RNAs were assessed using an Agilent 2100

1479

Bioanalyzer (Agilent Technologies, USA). Two total RNA libraries were subjected to reverse transcription using a complementary DNA High Capacity Kit (Invitrogen). The RT-qPCR analysis using SYBR Green PCR Master Mix (Takara, Dalian, China) was performed on a CFX Connect Real-Time PCR Detection System (Bio-Rad, CA, USA) to explore whether miRNA affected PRLR expression by degrading mRNA or suppressing mRNA posttranslational translation. Three housekeeping genes (glyceraldehyde-3phosphate [GAPDH], tyrosine 3-monooxygenase [YWHAZ], and b-actin) were used for normalization in RT-qPCR [21]. Relative expression levels of PRLR mRNA were calculated using the DDCt method [17]. The primer information of PRLR, b-actin, GAPDH, and YWHAZ genes are showed in Table S2. Fold changes were normalized by the expression levels of b-actin, GAPDH, and YWHAZ genes, and each assay was performed in triplicate.

2.6. Cell culture Human renal epithelial cell line 293T (HRE293T) was recovered from frozen stocks and cultured in Dulbecco’s modified Eagle’s medium with 10% fetal bovine serum, 100 units/mL penicillin, and 100 mg/mL streptomycin. All cells were cultured in a 5% CO2 humidified incubator at 37  C. Cell transfection was performed when cells reached a confluence of 70% to 80%.

2.7. Luciferase reporter assays To construct the luciferase reporter plasmids of PRLR 30 UTR, the PRLR 30 UTR fragments (443 bp) carrying the alleles, 173057T and 173057C, were amplified by PCR. The PCR products containing the alleles 173057T and 173057C were extracted and separated by agarose gel, and then, they were linked to a pMD19-T vector with TA cloning Kit (Invitrogen). The recombined pMD19-T vectors carrying different alleles were digested with XhoI and NotI endonuclease enzyme. Finally, the digested products carrying different alleles were inserted between renilla and firefly luciferase genes in a psiCHECK-2 Vector (Promega, WI, USA), and then the plasmids containing alleles 173057T or 173057C were conducted (Fig. 1A), which were confirmed by sequencing. For luciferase reporter assay, HRE293T cells were placed in 24-well plates (1  105 cells per well) and then cotransfected with psiCHECK-2 Vector 30 UTR-allele T or psiCHECK-2 Vector 30 UTR-allele C. The mimics of btamiR-302a and their negative controls (GenePharma, Shanghai, China) were cotransfected with the reporter plasmids at a final concentration of 20 nmol/mL. Thirty-six hours after transfection in HRE293T cells, renilla luciferase activity in lysates was measured with a Dual-Luciferase Reporter Assay System (Promega) and normalized against the activity of firefly luciferase. Assays were followed by the manufacture’s suggestions. Each experiment was independently repeated three times, and each sample was evaluated in triplicate. Using the Mann– Whitney test, a two-sided P value of 0.05 was considered significant.

1480

X. An et al. / Theriogenology 83 (2015) 1477–1483

Fig. 1. Characterization and functional analysis of the caprine PRLR 30 untranslated region (30 UTR). (A) Potential bta-miR-302a target sites at the caprine PRLR 30 UTR for allele C. PRLR 30 UTR was cloned into the psiCHECK plasmid at the XhoI enzyme site. The triangles show the 173057 T/C polymorphism. (B) HRE293T cells seeded on 24-well plates were transiently cotransfected with psiCHECK-PRLR 30 UTR and bta-miR-302a mimics or stable negative control (NC). Luciferase activity was measured after 36 hours of transfection. Renilla luciferase activity was normalized with firefly luciferase, and the mean activities  standard deviation from three independent experiments are shown. *P < 0.05.

3. Results 3.1. SNP identification and genotypes The g.151435C>T and g.151454A>G mutations were in the intron 2 of the caprine PRLR gene (GenBank accession No. KJ572973). The g.173057T>C mutation was in the 30 UTR (GenBank accession No. KJ572972). Both alleles of g.173057T>C SNP introduced several different miRNA sites, whereas the C allele introduced bta-miR-302a binding site that was abrogated in the presence of the T allele (Fig. S1). In both goat breeds, the PIC was 0.36 and 0.37 at the g.151435C>T locus, 0.36 at the g.151454A>G locus, 0.19 to 0.22 at the g.173057T>C locus (Table 1). According to the classification of PIC (low polymorphism if PIC value < 0.25, moderate polymorphism if 0.25 < PIC value < 0.50, and high polymorphism if PIC > 0.50), the GZ and BE goat breeds at the g.151435C>T and g.151454A>G loci had moderate genetic diversity. Genotypic distribution and allelic frequencies of three SNPs are shown in Table 1. It was shown that both goat breeds were in Hardy–Weinberg disequilibrium at three SNP loci (P < 0.05; Table 1). To reveal the linkage relationships between the three SNPs, the linkage disequilibrium was estimated in these breeds (Table 1). If r2 is greater than 0.33, the linkage disequilibrium was considered strong [22]. Following the result, both g.151435C>T and g.151454A>G loci were closely linked in both goat breeds. 3.2. Association analysis of SNPs with litter size In GZ and BE goat breeds, individuals with the TT genotype had greater litter size than those with the CC genotype in the third and average parity at the g.151435C>T

locus (P < 0.05; Table 2). At the g.173057T>C locus in the GZ goat breed, individuals with the CC genotype had higher litter size than those with the CT genotype in the fourth and average parity (P < 0.05); in addition, in the BE goat breed, individuals with the CC genotype had higher litter size than those with the CT genotype in the third and average parity (P < 0.05). In the GZ goat breed, individuals with the C6 genotype (TTAACC) had higher litter size than those with the C2 genotype (CCAACT) in the third parity and C5 genotype (TCGACT) in the first parity (P < 0.05; Table 3). Individuals with the C6 genotype (TTAACC) had higher litter size than those with the C1 (CCAACC), C3 (TCAACC), and C5 (TCGACT) genotypes in the average parity (P < 0.05). In the BE goat breed, individuals with the C6 genotype (TTAACC) had higher litter size than those with the C2 genotype (CCAACT) in the first, second, and third parity (P < 0.05) and C8 genotype (TCGACT) in the first and second parity; in addition, individuals with the C6 genotype (TTAACC) had higher litter size than those with the C5 (TCAACC), C7 (TCGACC), and C8 (TCGACT) genotypes in the average parity (P < 0.05). Other results are showed in Table 3. 3.3. Effects of SNP g.173057T>C on bta-miR-302a binding ability The caprine PRLR 30 UTR containing the C allele was predicted a binding site for bta-miR-302a (Fig. S1). Luciferase reporter assays showed that the 173057C clone had a significantly lower expression level of luciferase than did the 173057T clone and negative control (P < 0.05; Fig. 1). These results suggested that the allele C had a lower PRLR protein expression level than did the allele T. No significant difference levels of PRLR mRNA were observed among individuals with the CC and CT genotypes

X. An et al. / Theriogenology 83 (2015) 1477–1483

1481

Table 2 Association analysis of g.151435C>T, g.151454A>G, and g.173057T>C loci with litter size (mean  standard error) in goats. Locus

Breed

Genotype

First parity litter size

g.151435C>T

GZ

CC TC TT CC TC TT AA GA GG AA GA GG CC CT CC CT

1.39 1.45 1.53 1.28 1.28 1.32 1.52 1.45 1.43 1.33 1.26 1.35 1.49 1.41 1.36 1.14

BE

g.151454A>G

GZ

BE

g.173057T>C

GZ BE

(33) (185) (103) (25) (106) (60) (107) (177) (37) (70) (104) (17) (245) (76) (133) (58)

               

0.09 0.04 0.05 0.09 0.05 0.06 0.05 0.04 0.08 0.06 0.05 0.11 0.03 0.06 0.04b 0.06a

Second parity litter size 1.55 1.65 1.69 1.60 1.71 1.82 1.64 1.67 1.65 1.70 1.76 1.65 1.64 1.70 1.77 1.64

               

0.09 0.04 0.05 0.10 0.05 0.07 0.05 0.04 0.08 0.06 0.05 0.13 0.03 0.06 0.05 0.07

Third parity litter size 1.67 1.75 1.86 1.68 1.83 1.92 1.78 1.76 1.87 1.80 1.89 1.71 1.80 1.70 1.90 1.69

               

Fourth parity litter size

0.09a 0.04 0.05b 0.10a 0.05 0.06b 0.05 0.04 0.08 0.06 0.05 0.12 0.03 0.05 0.04b 0.06a

1.91 1.95 1.91 1.88 1.95 1.97 1.90 1.96 1.92 1.94 1.95 1.94 1.98 1.78 1.99 1.86

               

0.09 0.04 0.05 0.10 0.05 0.06 0.05 0.04 0.09 0.06 0.05 0.12 0.03b 0.06a 0.04 0.06

Average litter size 1.63 1.70 1.75 1.61 1.69 1.75 1.71 1.71 1.72 1.69 1.71 1.66 1.73 1.64 1.75 1.58

 0.05a 0.02  0.03b  0.06a  0.03  0.04b  0.03  0.02  0.04  0.03  0.03  0.09  0.02b  0.03a  0.02b  0.04a

a,b

Values with different superscripts within the same column in particular population differ significantly at P < 0.05. The number in parentheses represent the number of samples. Abbreviations: BE, Boer; GZ, Guanzhong.

(P > 0.05; Fig. S2). We speculated that bta-miR-302a may affect the expression levels of PRLR mRNA by suppressing mRNA posttranslational translation.

c.1201G>A (exon 9), and c.1355C>T (exon 9). The three SNPs (g.151435C>T, g.151454A>G, and g.173057T>C) were first found in the caprine PRLR gene. Prolactin receptor seems to be especially promising because it not only affects reproduction and growth traits but also affects milk production traits [11,25]. Studies on the intron 1 and exon 10 of PRLR gene showed that it is either a major gene influencing the prolificacy or is in close linkage with such a gene in Small Tail Han sheep, ewes with the genotype BB or AB had 0.64 to 0.76 or 0.44 to 0.54 more lambs than those with the genotype AA, respectively [11]. In this study, the g.151435C>T SNP was located in the intron region, but the result showed that it was associated with litter size in GZ and BE goat breeds. The reason could be that introns could influence many other stages of mRNA metabolism,

4. Discussion The study identified three SNPs in the caprine PRLR gene. The g.151435C>T and g.151454A>G mutations were in the intron 2, and the g.173057T>C mutation was located in the 30 UTR. Di et al. [23] found two SNPs (G35A and A86G) in the intron 2 of caprine PRLR gene. Sequence alignment revealed that g.151435C>T and g.151454A>G mutations were not the same mutations as A86G and G35A, respectively. Zidi et al. [24] detected the following three SNPs in the caprine PRLR coding region: c.1131G>A (exon 9), Table 3 Combined effect of three loci on litter size (mean  standard error) in goats. Breed

Genotype

Number

First parity litter size

GZ

C1 (CCAACC) C2 (CCAACT) C3 (TCAACC) C4 (TCGACC) C5 (TCGACT) C6 (TTAACC) C7 (TTAACT) C8 (TTGACC) C9 (TTGACT) C10 (TTGGCC) C11 (TTGGCT) C1 (CCAACC) C2 (CCAACT) C3 (CCGGCC) C4 (CCGGCT) C5 (TCAACC) C6 (TTAACC) C7 (TCGACC) C8 (TCGACT) C9 (TTGACC)

26 5 17 123 40 40 15 9 5 28 6 6 5 6 7 15 29 60 28 11

1.42 1.40 1.53 1.49 1.35 1.63 1.53 1.44 1.40 1.43 1.67 1.33 1.00 1.50 1.14 1.27 1.48 1.35 1.11 1.27

BE

                   

0.10 0.22 0.12 0.05 0.08a 0.07b 0.13 0.17 0.22 0.10 0.21 0.19 0.20ac 0.19 0.17 0.12 0.08b 0.06bc 0.09a 0.14

Second parity litter size 1.42 2.00 1.53 1.67 1.68 1.78 1.67 1.67 1.60 1.57 1.83 1.83 1.00 2.00 1.57 1.60 1.79 1.75 1.71 1.91

                   

0.10a 0.22b 0.12 0.04b 0.08b 0.08b 0.13 0.16 0.22 0.09 0.20 0.21b 0.23a 0.21b 0.20 0.13b 0.09b 0.07b 0.10b 0.16b

Third parity litter size 1.73 1.40 1.82 1.76 1.73 1.88 1.73 2.00 1.80 1.89 1.83 1.83 1.40 1.67 1.71 1.69 2.03 1.93 1.75 1.91

                   

0.10 0.22a 0.12 0.04 0.08 0.07b 0.13 0.16b 0.22 0.09b 0.20 0.19 0.21a 0.19 0.18 0.12ac 0.09b 0.06bc 0.09ac 0.14bc

Fourth parity litter size 1.92 1.80 1.77 2.04 1.78 2.00 1.80 1.89 1.60 1.93 1.83 2.00 1.80 2.00 1.86 1.93 2.00 2.00 1.93 1.82

                   

0.10 0.23 0.13a 0.05b 0.08a 0.08 0.13 0.17 0.23 0.10 0.21 0.20 0.22 0.20 0.18 0.13 0.09 0.06 0.09 0.15

Average litter size 1.63 1.65 1.66 1.74 1.63 1.82 1.68 1.75 1.60 1.71 1.79 1.75 1.30 1.79 1.57 1.62 1.83 1.76 1.63 1.73

                   

0.05ac 0.12 0.06c 0.02bc 0.04ac 0.04b 0.07 0.08 0.12 0.05 0.11 0.11bde 0.12a 0.11bde 0.10 0.07de 0.05b 0.04d 0.05e 0.08bde

a,b,c,d,e Values with different superscripts within the same column in particular population differ significantly at P < 0.05. Combination genotypes less than five were not in table. Abbreviations: BE, Boer; GZ, Guanzhong.

1482

X. An et al. / Theriogenology 83 (2015) 1477–1483

including initial transcription of the gene, editing and polyadenylation of the pre-mRNA, nuclear export, and translation and decay of the mRNA product [26]. The SNP g.173057T>C was located in the 30 UTR that is complementary to the seed region of bta-miR-302a. The C to T substitution of this SNP leads to noncomplementarity between bta-miR-302a and PRLR mRNA pairing. The allele C can cause lower PRLR protein expression than the allele T. A recent study identified SNPs in the porcine PRLR (pPRLR) short isoforms 30 UTR that reduce protein expression and four haplotypes of pPRLR short isoforms that suppress pPRLR long isoforms signaling and may differentially impact the phenotypic effects of PRL in vivo [6]. Therefore, endogenous bta-miR-302a might regulate PRLR protein expression, resulting in that the g.173057T>C SNP had significant effects on litter size. Furthermore, no significant difference of the PRLR mRNA expression level among different genotypes was observed in ovarian tissues using the RT-qPCR assay. These results suggested that the SNP did not affect mRNA expression; however, given that miRNA binding to mRNAs does not always lead to transcript cleavage, and sometimes it leads to translation repression, it is possible that the SNP leads to a change in PRLR protein. However, we did not test this in our study, and therefore remains a possibility, but not proven. These data suggested that bta-miR-302a may affect PRLR expression by suppressing mRNA posttranslational translation. The reproductive traits are complex quantitative traits involving multiple genes, loci, and interactions, so it is important to analyze the combined effect of multiple genes or loci on litter size [27]. In this study, associations between multiple loci and litter size from the first to the fourth parity were analyzed. The litter size of goats tended to increase in later lactations. Single SNP trait association analyses revealed that g.151435C>T and g.173057T>C loci had effects on litter size in both goat breeds. In further study, combinative genotypes of three loci reported that C6 genotype (TTAACC) was more favorable for litter size in both goat breeds. Accumulating evidence further showed that PRL interacts with the conserved extracellular domain of PRLR in various target cells and activates a cascade of intracellular events, mainly the JAK2/STAT signaling pathway, via specific sites on the PRLR intracellular domain tail [28]. The biochemical and physiological functions, together with the results obtained in our study, indicate that the SNPs associated with litter size had potential applications in a marker-assisted selection program for goat breeding. In conclusions, endogenous bta-miR-302a might regulate PRLR protein expression, resulting in that the g.173057T>C SNP had significant effects on litter size. Further analysis revealed that the combinative genotype C6 (TTAACC) was better than the others for litter size. Therefore, the present study contributes to evaluate them as genetic markers in goat genetics and breeding and have potential applications in breeding programs. 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), and the PhD research funding of Northwest A&F University (2013BSJJ003) Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j. theriogenology.2015.01.030.

References [1] Krasnow JS, Hickey GJ, Richards JS. Regulation of aromatase mRNA and estradiol biosynthesis in rat ovarian granulosa and luteal cells by prolactin. Mol Endocrinol 1990;4:12–3. [2] van Rens BT, Evans GJ, van der Lende T. Components of litter size in gilts with different prolactin receptor genotypes. Theriogenology 2003;59:915–26. [3] Chan SR, Rickert CG, Vermi W, Sheehan KCF, Arthur C, Allen JA, et al. Dysregulated STAT1-SOCS1 control of JAK2 promotes mammary luminal progenitor cell survival and drives ERaþ tumorigenesis. Cell Death Differ 2014;21:234–46. [4] Hayes H, Le Chalony C, Goubin G, Mercier D, Payen E, Bignon C, et al. Localization of ZNF164, ZNF146, GGTA1, SOX2, PRLR and EEF2 on homoeologous cattle, sheep and goat chromosomes by fluorescent in situ hybridization and comparison with the human gene map. Cytogenet Cell Genet 1996;72:342–6. [5] Scott P, Kessler MA, Schuler LA. Molecular cloning of the bovine prolactin receptor and distribution of prolactin and growth hormone receptor transcripts in fetal and utero-placental tissues. Mol Cell Endocrinol 1992;89:47–58. [6] 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 2013;45:75–85. [7] Belugin S, Diogenes AR, Patil MJ, Ginsburg E, Henry MA, Akopian AN. Mechanisms of transient signaling via short and long prolactin receptor isoforms in female and male sensory neurons. J Biol Chem 2013;288:34943–55. [8] Molitch ME. Mutant prolactin receptor and familial hyperprolactinemia. N Engl J Med 2014;370:977. [9] Bole-Feysot C, Goffin V, Edery M, Binart N, Kelly PA. Prolactin (PRL) and its receptor: actions, signal transduction pathways and phenotypes observed in PRL receptor knockout mice. Endocr Rev 1998;19:225–68. [10] Terman A. Effect of the polymorphism of prolactin receptor (PRLR) and leptin (LEP) genes on litter size in Polish pigs. J Anim Breed Genet 2005;122:400–4. [11] Chu MX, Mu YL, Fang L, Ye SC, Sun SH. Prolactin receptor as a candidate gene for prolificacy of small tail han sheep. Anim Biotechnol 2007;18:65–73. [12] Lu AJ, Hu XC, Chen H, Dong YJ, Pang YH. Single nucleotide polymorphisms of the prolactin receptor (PRLR) gene and its association with growth traits in chinese cattle. Mol Biol Rep 2011;38:261–6. [13] Ambros V. The functions of animal microRNAs. Nature 2004;431: 350–5. [14] Otsuka M, Zheng M, Hayashi M, Lee JD, Yoshino O, Lin SC, et al. Impaired microRNA processing causes corpus luteum insufficiency and infertility in mice. J Clin Invest 2008;118:1944–54. [15] Yao N, Lu CL, Zhao JJ, Xia HF, Sun DG, Shi XQ, et al. A network of miRNAs expressed in the ovary are regulated by FSH. Front Biosci 2009;14:3239–44. [16] Brodersen P, Voinnet O. Revisiting the principles of microRNA target recognition and mode of action. Nat Rev Mol Cell Biol 2009;10:141–8. [17] 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:362–76. [18] 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:161–7. [19] Gao Q, Ju ZH, Zhang Y, Huang JM, Zhang XJ, 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.

X. An et al. / Theriogenology 83 (2015) 1477–1483 [20] 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:16871–4. [21] Zhang Y, Zhang XD, Liu X, Li YS, Ding JP, Zhang XR, et al. Reference gene screening for analyzing gene expression across goat tissue. Asian Australas J Anim 2013;26:1665–71. [22] Sun P, Zhang R, Jiang Y, Wang X, Li J, Lv H, et al. Assessing the patterns of linkage disequilibrium in genic regions of the human genome. FEBS J 2011;278:3748–55. [23] Di R, Yin J, Chu MX, Cao GL, Feng T, Fang L, et al. DNA polymorphism of introns 1 and 2 of prolactin receptor gene and its association with litter size in goats. Anim Sci Pap Rep 2011;29: 343–50. [24] Zidi A, Serradilla JM, Jordana J, Carrizosa J, Urrutia B, Polvillo O, et al. Pleiotropic effects of the goat prolactin receptor genotype

[25]

[26]

[27] [28]

1483

on milk fatty acid composition. Domest Anim Endocrinol 2010; 39:85–9. Lu A, Hu X, Chen H, Dong Y, Zhang Y, Wang X. Novel SNPs of the bovine PRLR gene associated with milk production traits. Biochem Genet 2011;49:177–89. Lemieux C, Marguerat S, Lafontaine J, Barbezier N, Bahler J, Bachand F. A pre-mRNA degradation pathway that selectively targets intron-containing genes requires the nuclear poly(A)-binding protein. Mol Cell 2011;44:108–19. An X, Ma T, Hou J, Fang F, Han P, Yan Y, et al. Association analysis between variants in KISS1 gene and litter size in goats. BMC Genet 2013;14. Neradugomma NK, Subramaniam D, Tawfik OW, Goffin V, Kumar TR, Jensen RA, et al. Prolactin signaling enhances colon cancer stemness by modulating Notch signaling in a Jak2-STAT3/ERK manner. Carcinogenesis 2014;35:795–806.

1483.e1

X. An et al. / Theriogenology 83 (2015) 1477–1483

Table S1 Primer information of PRLR gene for screening polymorphisms. 

Primer

Sequence (50 /30 )

Gene region

Name to each DNA fragment

Product size (bp)

Ta ( C)

F1 RAa FAa R2 F3 R3

F: ATATTCAGCAAGGAGCAAGA R: ATCCTCCTCTTGCGCAAATA F: ATTATTTGCACCAGAGGAGG R: AATGAGGATGGAAGTCAGAG F: AGTGAGAGTTATGGAAGGATG R: AAGGTTAAGCAACTGGTCTT

Exon 2 and partial intron 2

I

296

53

Exon 2 and partial intron 2

II

200

50

30 UTR

III

443

55

Abbreviations: Ta, annealing temperature; 30 UTR, 30 untranslated region. a RA and FA were used in primer-introduced restriction analysis-polymerase chain reaction. RA and FA were deliberately introduced into a point mutation (C) to create MwoI and BslI restriction sites, respectively.

Table S2 Primer information of PRLR and housekeeping genes for real-time quantitative polymerase chain reaction. 

Primer Sequence (50 /30 )

Gene region Product Ta ( C) size (bp)

RT

PRLR-CDS

b G Y

F: GCGAGGATTTGCTGATGGAA R: CGCCTTGCTCCATGTGTTCT F: CCACACCTTCTACAACGAGC R: ATCTGGGTCATCTTCTCACG F: AACCTGCCAAGTATGATGAG R: AGTGTCGCTGTTGAAGTC F: TCCATCGAACATCCAATCAT R: GCTCCTTGCTCAGTTACA

93

60

b-actin-CDS 105

60

GAPDH-CDS 118

60

YWHAZ-CDS 118

60

Abbreviation: Ta, annealing temperature.