Abstract
Juan Ji1,2, 4,*, Yufeng Qin3,*, Ran Zhou1,2,*, Rujin Zang5, Zhenyao Huang1,2, Yan Zhang1,2, Minjian Chen1,2, Wei Wu1,2, Ling Song1,2, Xiufeng Ling4, Hongbing Shen1,6, Zhibin Hu1,6, Yankai Xia1,2, Chuncheng Lu1,2, Xinru Wang1,2
1State Key Laboratory of Reproductive Medicine, Institute of Toxicology, Nanjing Medical University, Nanjing 210029, China
2Key Laboratory of Modern Toxicology of Ministry of Education, School of Public Health, Nanjing Medical University, Nanjing 210029, China
3Epigenetics & Stem Cell Biology Laboratory, National Institute of Environmental Health Sciences, Research Triangle Park, NC 27709, USA
4Department of Children Health Care, Nanjing Maternity and Child Health Care Hospital Affiliated to Nanjing Medical University, Nanjing 210029, China
5Department of Pediatric Surgery, State Key Laboratory of Reproductive Medicine, Nanjing Children’s Hospital Affiliated Nanjing Medical University, Nanjing 210008, China
6Department of Epidemiology and Biostatistics and Key Laboratory of Modern Toxicology of Ministry of Education, School of Public Health, Nanjing Medical University, Nanjing 211166, China
*The authors have contributed equally to this study and they should be regarded as joint first authors
Correspondence to:
Chuncheng Lu, email: [email protected]
Xinru Wang, email: [email protected]
Keywords: X-linked miRNAs, non-obstructive azoospermia
Received: January 08, 2016 Accepted: March 28, 2016 Published: April 16, 2016
ABSTRACT
Human X chromosome has higher densities of microRNAs (miRNAs) compared to the average densities on autosomes. Given that numbers of X-linked miRNAs can escape from meiotic sex chromosome inactivation (MSCI) silencing, it is proposed that X-linked miRNAs may play critical roles in the process of spermatogenesis. To test the hypothesis, we performed DNA capture sequencing of human X-linked miRNAs, which was followed by a two-stage case-control study to identify the non-obstructive azoospermia (NOA) related single nucleotide variants (SNVs) in 1107 NOA cases and 1191 fertile healthy controls. Eventually, we found rs5951785, located near hsa-miRNA-506/507, increased the risk of NOA, while rs1447393, near hsa-miRNA-510, decreased the risk of NOA. Functional analysis revealed that rs5951785 significantly inhibited cell proliferation and induced cell apoptosis. Taken together, our results demonstrated that X-linked miRNAs played important roles in the pathogenesis of NOA.
INTRODUCTION
MicroRNAs (miRNAs), a novel class of small non-coding regulatory RNAs [1], can repress target gene expression through binding to the 3’ untranslated regions (3’ UTR) of mRNA [2, 3], and their dysregulation contributes to the aetiology of numerous diseases, such as cancer, metabolic diseases and viral pathogenesis [4]. Accumulating evidence demonstrates that miRNAs take part in the process of spermatogenesis [5].
Comparing to autosomes, human X chromosome, containing genes that are uniquely expressed in male spermatogenesis, has higher density of miRNAs (referred as X-linked miRNAs) [6] and 38% of them are testis-specific or testis-enriched [7]. They participate in a number of biological processes, including cellular differentiation, proliferation and apoptosis [8]. Due to the existence of the meiotic sex chromosome inactivation (MSCI), X chromosome is silenced during spermatogenesis [9]. However, some X-linked miRNAs can escape from MSCI and regulate gene expression during late meiosis and early post-meiotic stage of spermatogenesis [7, 10]. Recently, Royo et al. reported that X-linked miRNAs are expressed before prophase I but absent during pachynema [11], and X-linked miRNAs have been reported as potential contributor to the regulation of cell cycle and spermatogonial stem cell renewal, which is crucial for the regulation of spermatogenesis, especially the mitosis of spermatogonia and meiosis of spermatocytes [12].
Single nucleotide variations (SNVs) (including mutation, single nucleotide polymorphism (SNP) and insertion-deletion (INDEL)) in pre-miRNAs, flanking regions and miRNA binding regions have been demonstrated to affect miRNA-mediated regulatory function and related with human diseases (e.g. male infertility) [13]. Zhang et al. reported that SNPs in miR-1302-binding site of CGA increased the risk of idiopathic male infertility [13]. Our previous studies also demonstrated that genetic variants in miRNA biogenesis genes (DICER1 and DROSHA) increased the risk of idiopathic male infertility [14].
To systematically analyze X-linked miRNA-related SNVs and their effects on non-obstructive azoospermia (NOA), the most severe form of male infertility, we performed a chromosome-wide scan for SNVs near pre-miRNAs in X chromosome followed by SNPscan and SNaPshot in two independent validation phases. Eventually, two SNPs (rs5951785 and rs1447393), located near the regions of hsa-miR-506/507 and hsa-miR-510, were found to be associated with NOA, and in vitro analysis was performed to clarify their potential functions in spermatogenesis.
RESULTS
Identifying X-linked miRNAs’ SNVs by sequencing
Based on the miRNASNP (http://www.bioguo.org/miRNASNP/), customer designed arrays were used to capture the X-linked miRNAs regions including upstream (1KB) and downstream (1KB) in 96 NOA cases and 96 fertile controls. Totally, 91 regions were captured followed by high-throughput sequencing on Illumina HiSeq 2000 to generate 100 pair-end reads. On average, each sample was sequenced to an average depth of 115×, with nearly 90% of the targeted regions covered by ≥2×. Totally, we identified 139 SNVs (including one SNP with MAF<0.05 and 138 SNPs with MAF>0.05) in X-linked miRNAs regions, among which, 22 SNPs were associated with NOA risk (Supplementary Material, Table S1).
Two-stage validation in large cohort
For fast track replication, 22 signals were included in Stage I validation using an independent Chinese population (Supplementary Material, Table S2) by custom designed SNPscan. Unexpectedly, 18 markers identified in screening stage were not observed to be associated with NOA risk in Stage I with P values > 0.05 (Supplementary Material, Table S2). Only rs547043 near hsa-miR-4330, rs5951785 near hsa-mir-506/507, rs1447393 near hsa-mir-510, and rs5985440 near hsa-miR-652 were retained associated with NOA, among which rs547043 near hsa-miR-4330 was inconsistent with screening stage. To confirm the relationship between the other 3 SNPs and NOA risk, we carried out stage II validation in another large population (Supplementary Material, Table S3). Rs5951785 and rs1447393 were both observed to be still associated with NOA in the same direction as illustrated in the screening stage and validation stage I. Next, we performed a meta-analysis of the genotype data from stage I and II. In the combined analysis, we found that rs5951785 significantly increased the risk of NOA in the Han Chinese population (P meta = 1.01×10−3, OR = 1.45, 95% CI = 1.16-1.81). Rs1447393 acted as potential protective factor on NOA (P meta =1.31×10−4, OR=0.58, 95%CI=0.44-0.77) (Table 1). To further extend our analyses, we searched the GTEx database to see whether these two variants were quantitative trait loci (eQTL) variants, albert we did not find significant eQTLs in the available data sets [15].
Table 1: Two SNPs in human X-linked miRNAs were identified associated with NOA and validated in two independent cohorts
MiRNA | ACC | Location in X chr. (hg19) | SNPa | Study | Casesa | Controlsa | OR (95%CI)b | Pb | Qc |
|---|---|---|---|---|---|---|---|---|---|
hsa-mir-506/507 | MI0003193 | 146312238-146312361 | rs5951785 A > G | Validation I | 423/113 | 471/85 | 1.48(1.09-2.02) | 1.31×10−2 | |
Validation II | 367/99 | 448/85 | 1.42(1.03-1.96) | 3.16×10−2 | |||||
Combinedd | 790/212 | 919/170 | 1.45(1.16-1.81) | 1.01×10−3 | 0.856 | ||||
hsa-mir-510 | MI0003197 | 146353853-146353926 | rs1447393 C > G | Validation I | 462/47 | 469/89 | 0.54(0.37-0.78) | 1.15×10−3 | |
Validation II | 423/37 | 471/65 | 0.63(0.41-0.97) | 3.49×10−2 | |||||
Combinedd | 885/84 | 940/154 | 0.58(0.44-0.77) | 1.31×10−4 | 0.596 |
Chr, chromosome; ACC, Accession number in miRbase;NOA, non-obstructive azoospermia; OR, odds ratio. Combined P values with a fixed effect model are presented.
a Major allele/minor allele.
b OR and P values were calculated by additive model.
c P value for Cochrane Q statistics test.
d Combined P values with a fixed effect model are presented.
Effects of SNPs on MiRNAs and their targets
To understand the impacts of these two SNPs (rs5951785 near hsa-miR-506/507; rs1447393 near hsa-miR-510) on the miRNA expression, we transfer the wild-type pre-miRNAs and mutant pre-miRNAs into HEK-293T cells. Through qPCR, we found that the expression levels of these three miRNAs were all significantly down-regulated in mutant type (rs5951785 near hsa-miR-506, P=4.47×10−2, Figure 1A; rs5951785 near hsa-miR-507, P=4.80×10−3, Figure 2A; rs1447393 near hsa-miR-510, P=3.00×10−4, Figure 3A).
Figure 1: Functional analysis revealed that rs5951785 near miR-506 might contribute to the risk of NOA. A. Quantitative real-time PCR was applied to detected mature miR-506 expression levels in transfected HEK-293T cells. The miR-506 expression level was significant down-regulated in mutant type compared with the wild type. B. Dual-luciferase reporter assay was conducted to measure effects of rs5951785 near miR-506 on promoter-transcriptional activity, it was increased in HEK-293T cells with miR-506 mutant type for GLI3 and PIK3C2A. C. CCK8 was used to determine the influence of rs5951785 near miR-506 on cell growth. Rs5951785 near miR-506 significantly inhibited the HEK-293T cell growth at 48h and 72h. D, E. Flow cytometry assay was conducted to measure the effects of rs5951785 near miR-506 on cell apoptosis. (E) Histogram of cell apoptosis was presented to depict cell apoptotic percentages. The cell apoptosis population was markedly increased in HEK-293T cells transfected with miR-506 mutant type. F, G. Effects of rs5951785 near miR-506 on cell cycle were performed with flow cytometry. (G) Results quantitated in cell cycle were shown in histogram. There was no difference of cell cycle between miR-506 wild and miR-506 mutant type. Each data point represented the mean ± SE from three separate experiments in which treatments were performed in triplicate. *P<0.05, **P<0.01, ***P<0.001.
Figure 2: Functional analysis demonstrated that rs5951785 near miR-507 might contribute to the risk of NOA. A. RT-PCR was applied to detected mature miR-507 expression levels. The miR-507 expression level was significant down-regulated in mutant type compared with the wild type. B. Dual-luciferase reporter assay was conducted to measure the effects of rs5951785 near miR-507 on promoter-transcriptional activity. It was significantly increased in miR-507 mutant type for ADAM17. C. CCK8 was used to determine the influence of rs5951785 near miR-507 on cell growth. No significant difference was observed in cell proliferation of the rs5951785 near miR-507. D, E. Flow cytometry assay was conducted to measure the effects of rs5951785 near miR-507 on cell apoptosis. (E) Histogram of cell apoptosis was presented to depict cell apoptotic percentages. There was no significant difference of cell apoptosis in rs5951785 near miR-507. F, G. Effects of rs5951785 near miR-507 on cell cycle were analyzed with flow cytometry. (G) Results quantitated in cell cycle were shown in histogram. There was no difference of cell cycle between miR-507 wild and miR-507 mutant type. Each data point represented the mean ± SE from three separate experiments in which treatments were performed in triplicate. *P<0.05, **P<0.01, ***P<0.001.
Figure 3: Functional analysis indicated that rs1447393 near miR-510 might protect from the risk of NOA. A. RT-PCR was applied to detected mature miR-510 expression levels. The miR-510 expression level was significant down-regulated in mutant type compared with the wild type. B. Dual-luciferase reporter assay was conducted to measure the effects of rs1447393 near miR-510 on promoter-transcriptional activity. It was significantly increased in miR-510 mutant type for PRDX1. C. CCK8 was used to determine the influence of rs1447393 near miR-510 on cell growth. The cell growth was markedly increased with miR-510 mutant type at 48h. D, E. Flow cytometry assay was performed to analyze the effects of rs1447393 near miR-510 on cell apoptosis. (E) Histogram of cell apoptosis was presented to depict cell apoptotic percentages. There was no significant difference of cell apoptosis in rs1447393 near miR-510. F, G. Effects of rs1447393 near miR-510 on cell cycle were determined by flow cytometry. (G) Results quantitated in cell cycle were shown in histogram. There was no difference of cell cycle between miR-510 wild and miR-510 mutant type. Each data point represented the mean ± SE from three separate experiments in which treatments were performed in triplicate. *P<0.05, **P<0.01, ***P<0.001.
Next, we used TargetScan, miRanda and miRwalk to predict the putative targets of these three miRNAs (Supplementary Material, Tables S4-S5). Through conducting dual-luciferase reporter assay, we found that the luciferase activities of CDK4, GLI3, PIK3C2A, ADAM17, SFRP2 and PRDX1 were significantly decreased when compared to the vectors, suggesting that they were the potential targets of hsa-miR-506, hsa-miR-507 and hsa-miR-510, respectively (Figure 1B, 2B, 3B). Then, we wanted to know whether these SNPs changed miRNAs’ binding ability with its targets. As results shown, the relative luciferase activities were significantly increased with hsa-miR-506 mutant type for GLI3 (P=1.11×10−5) and PIK3C2A (P=4.00×10−4, Figure 1B), hsa-miR-507 mutant type for ADAM17 (P=2.00×10−4, Figure 2B) and hsa-miR-510 mutant type for PRDX1 (P=2.00×10−3, Figure 3B) compared with wild-type, respectively. However, to the CDK4 (Figure 1B) and SFRP2 (Figure 3B), no significant difference was observed between the wild-type and mutant-type.
Effects of MiRNAs and sNPs on cell functions
To investigate whether these X-linked miRNAs involved in cell function, HEK-293T cells were transfected with miRNAs containing wild type and mutant allele. As shown in Figure 1, cell growth was significantly decreased at 48h (P=1.00×10−2) and 72h (P=9.30×10−3) with hsa-miR-506 mutant type (Figure 1C). It suggested that rs5951785 near hsa-miR-506 exerted a growth-inhibiting function. The cell growth was markedly increased at 48h (P=4.16×10−2) and 72h (P=4.90×10−3) with hsa-miR-510 mutant allele (Figure 3C), indicating that rs1447393 near hsa-miR-510 promoted cell proliferation. No significant difference was observed in the rs5951785 near hsa-miR-507 (Figure 2C). Besides, we also found hsa-miR-506 mutant-type significantly increased cell apoptosis (P=1.69×10−2, Figure 1D, 1E), while there was no difference in hsa-miR-507 nor hsa-miR-510 between the wild and mutant alleles (Figure 2D, 2E; Figure 3D, 3E), implicating that growth inhibition was accompanied with increased apoptosis population. No significant difference was observed in cell cycle (Figure 1F, 1G; Figure 2F, 2G; Figure 3F, 3G).
DISCUSSION
X-linked miRNAs are important to spermatogenesis, but till now few studies have investigated the relationship between X-linked miRNAs and male infertility. To comprehensively investigate their relationship, we conducted X chromosome-wide miRNAs capture sequencing and large-scale follow-up genotyping. Our results demonstrated that i) rs5951785 near hsa-miR-506/507 was associated with significantly increased risk of NOA, while rs1447393 near hsa-miR-510 decreased the risk of NOA; ii) rs5951785 significantly decreased binding affinity of hsa-miR-506, inhibited cell proliferation and promoted cell apoptosis.
Several studies have reported that the SNPs in pre-miRNAs could reduce the production of mature miRNAs [16, 17]. Consistent with these related studies, the expression levels of these three miRNAs were all down-regulated with mutant alleles. Through dual luciferase report assays, we found rs5951785 and rs1447393 decreased their binding activity with mRNA 3’UTR of GLI3, PIK3C2A, ADAM17 and PRDX1, respectively.
GLI3, a potential target of miR-506, has been reported highly expressed in spermatogonia cells (both type A and type B) [18]. It could inhibit cell proliferation and promote cell death [19], which was consistent with our functional analysis of rs5951785 near hsa-miR-506, suggesting that GLI3 might be associated with NOA. ADAM17, a target of miR-507, is involved in germ cell apoptosis during spermatogenesis [20]. Environmental toxicants exposures can also active ADAM17 to induce the extrinsic pathway of apoptosis in the germ cell [21]. However, we didn’t observe apoptosis changes after overexpress hsa-miR-507. This might due to the low expression level of ADAM17 in HEK-293T cell. PRDX1, a target gene of hsa-miR-510, was reported to participate in removing the ROS [22]. Spermatozoa is highly susceptible to ROS-induced damage and dysregulation of PRDX1 might induce sperm DNA damage.
In conclusion, rs5951785 near hsa-miR-506/507 and rs1447393 near hsa-miR-510 were identified to be potential modifier of NOA. Our findings further highlighted that genetic variations in miRNAs might play important roles in the pathogenesis of NOA and might serve as therapeutic targets for NOA.
MATERIALS AND METHODS
Study design and participants
This study was approved by the institutional review board of Nanjing Medical University, China (FWA00001501) and conducted according to the Declaration of Helsinki. All participants were voluntary and completed the informed consent in written before taking part in this research.
A three-stage case-control study was conducted: we captured and sequenced all X-linked miRNAs by customer designed array (Agilent SureSelect DNA Capture Array) in 96 NOA affected individuals and 96 fertile controls in screening stage, and a subsequent 2-stage replication with large-scale follow-up genotyping in large cohorts of NOA-affected subjects and fertile controls was performed. The subjects included in screening stage were recruited from the Nanjing Center of Reproductive Medicine between March 2010 and January 2013. The first replication stage (stage I) included 548 NOA cases recruited from the infertility clinic at the Affiliated Hospitals of Nanjing Medical University at Jiangsu (NJMU Infertile study) and 560 healthy male controls from the same hospitals during the same period. And the second replication stage (stage II) included 465 NOA cases sampled from Renji Hospital, Shanghai, and 537 healthy male controls from the same hospital. Some cohorts within the sample sets have been reported in previously published data [23, 24].
X-linkd miRNAs sequencing
Through miRNASNP (http://www.bioguo.org/miRNASNP/), we included 91 X-linked miRNAs (Supplementary Material, Table S6) in the study [25]. Customer designed Agilent SureSelect Arrays, which target pre-miRNAs regions and enlonger +/- 1KB, were used to capture DNA by hybridization. The hybridization and enrichment was carried out on AB 2720 Thermal Cycler (Life Technologies, USA). Captured DNA was enriched by SureSelect Target Enrichment Kit (Agilent technologies, USA) and used to construct library by Truseq DNA Sample preparation Kit following the standard protocol (Illumina, USA). Sequencing was performed on the Illumina HiSeq 2000 to generate paired-end 100bp reads. Sequencing data was deposited in NCBI SRA with accession number SRP072055.
Identifying X-linked miRNA SNVs from sequencing data
After checking quality of sequencing reads by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), low quality, short reads and 3’ end of reads with a quality less than 15 were trimmed by FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html). For calling of variants, good quality reads were mapped to the human reference genome (UCSC, hg19) using Burrows-Wheeler Aligner (BWA) [26]. Genome Analysis Toolkit (GATK) was used for realignment, base quality score recalibration and duplicate remove [27]. In brief, VarScan and SAMtools were used to identify variants with in-house parameters [28, 29]. SNPs, In/Dels and mutations were annotated using dbSNP132 or the 1000 Genomes Project.
Follow-up genotyping of screened genetic variants by SNPscan and SNaPshot
Candidate SNVs were genotyped in additional cases and controls by a custom-by-design 48-Plex SNPscan™ Kit (Cat#:G0104; Genesky Biotechnologies Inc., Shanghai, China). This kit was developed according to patented SNP genotyping technology by Genesky Biotechnologies Inc., which was based on double ligation and multiplex fluorescence PCR [30]. In order to validate the genotyping accuracy using SNPscan™ Kit, 5% duplicate samples were analyzed by single nucleotide extension using the Multiplex SNaPshot Kit (Applied Biosystems Inc., Foster City, CA, USA), and the concordance rates were more than 99%.
Statistical analysis
Stata package (Stata Corp, LP) and P-link packages were used for data analyses [31]. Multivariate logistic regression under additive genetic model was used to calculate the odds ratios (ORs) and 95% confidence intervals (CIs) of the associations between SNPs and NOA risks. The Cochrane Q statistics test was used for the assessment of heterogeneity. Mann–Whitney test was used to compare the gene expression levels between different groups.
Target gene prediction and luciferase reporter assays
TargerScan (http://www.targetscan.org/), miRanda (http://www.microrna.org/microrna/home.do) and miRwalk (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/), were applied to identify the potential mRNA targets of the three X-linked miRNAs (Supplementary Material, Tables S4-S5). Luciferase Reporter Assayswas descripted previously [32]. All the sequences of pre-miRNAs were confirmed by Sanger sequencing and shown in the Supplementary Material, Table S7.
Cell proliferation, cell cycle and cell apoptosis assay
Cells were seeded in 96-well plates. 24, 48 and 72 hours after transfection, CCK8 solution (Vazyme, China) was add to each well, then continued incubating for 1 hour at 37°C. The absorbance of each well was detected using TECAN Infinite M200 multimode microplate reader (Tecan, Belgium). As to cell cycle and apoptosis, cells were collected 24 hours after transfections and measured by FACS (BD, USA). All experiments were performed in triplicate independently.
ACKNOWLEDGMENTS
We thank all the research staff and students who took part in this work. Funding was provided by grants from National Basic Research Program of China (973 Program, 2012CBA01306), the National Natural Science Foundation of China (81471500, 81330067 and 81322039), Distinguished Young Scholars of Jiangsu Province (BK20130041) and the Priority Academic Program for the Development of Jiangsu Higher Education Institutions (Public Health and Preventive Medicine).
CONFLICT OF INTEREST
None declared.
Author Contributions
CL and XW directed the study, obtained financial support and were responsible for study design. JJ performed overall project management with YQ and RZ, performed statistical analysis with ZH, MC and WW, and drafted the initial manuscript. RZ, YZ and LS were responsible for animal care and sample preparation. XL, HS, ZH and YX conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
REFERENCES
1. Ha M, Kim VN. Regulation of microRNA biogenesis. Nat Rev Mol Cell Biol. 2014; 15:509-524.
2. Wightman B, Ha I, Ruvkun G. Posttranscriptional regulation of the heterochronic gene lin-14 by lin-4 mediates temporal pattern formation in C. elegans. Cell. 1993; 75:855-862.
3. Lee RC, Feinbaum RL, Ambros V. The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell. 1993; 75:843-854.
4. Hammond SM. An overview of microRNAs. Adv Drug Deliv Rev. 2015; 87:3-14.
5. Hayashi K, Chuva de Sousa Lopes SM, Kaneda M, Tang F, Hajkova P, Lao K, O’Carroll D, Das PP, Tarakhovsky A, Miska EA, Surani MA. MicroRNA biogenesis is required for mouse primordial germ cell development and spermatogenesis. PloS ONE. 2008; 3:e1738.
6. Guo X, Su B, Zhou Z, Sha J. Rapid evolution of mammalian X-linked testis microRNAs. BMC genomics. 2009; 10:97.
7. Song R, Ro S, Michaels JD, Park C, McCarrey JR, Yan W. Many X-linked microRNAs escape meiotic sex chromosome inactivation. Nat Genet. 2009; 41:488-493.
8. Croce CM, Calin GA. miRNAs, cancer, and stem cell division. Cell 2005; 122:6-7.
9. Turner JM. Meiotic sex chromosome inactivation. Development. 2007; 134:1823-1831.
10. Papaioannou MD, Nef S. microRNAs in the testis: building up male fertility. J Androl. 2010; 31:26-33.
11. Royo H, Seitz H, Ellnati E, Peters AH, Stadler MB, Turner JM. Silencing of X-Linked MicroRNAs by Meiotic Sex Chromosome Inactivation. PLoS Genet. 2015; 11: e1005461.
12. Stouffs K, Lissens W. X chromosomal mutations and spermatogenic failure. Biochim Biophys Acta. 2012; 1822:1864-1872.
13. Zhang H, Liu Y, Su D, Yang Y, Bai G, Tao D, Ma Y, Zhang S. A single nucleotide polymorphism in a miR-1302 binding site in CGA increases the risk of idiopathic male infertility. Fertil Steril. 2011; 96:34-39.
14. Qin Y, Xia Y, Wu W, Han X, Lu C, Ji G, Chen D, Wang H, Song L, Wang S, Wang X. Genetic variants in microRNA biogenesis pathway genes are associated with semen quality in a Han-Chinese population. Reprod Biomed Online. 2012; 24:454-461.
15. GTEx Consortium. Human genomics. The Genotype- Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015; 348: 648-660.
16. Hu Z, Chen J, Tian T, Zhou X, Gu H, Xu L, Zeng Y, Miao R, Jin G, Ma H, Chen Y, Shen H. Genetic variants of miRNA sequences and non-small cell lung cancer survival. J Clin Invest. 2008; 118:2600-2608.
17. Wang SH, Wang SH, Li H, Sun GR, Lyu SJ, Liu XJ, Li ZJ, Kang XT. SNP in pre-miR-1666 decreases mature miRNA expression and is associated with chicken performance. Genome. 2015; 58:81-90.
18. Persengiev SP, Kondova, II, Millette CF, Kilpatrick DL. Gli family members are differentially expressed during the mitotic phase of spermatogenesis. Oncogene. 1997; 14:2259-2264.
19. Trnski D, Sabol M, Gojevic A, Martinic M, Ozretic P, Musani V, Ramic S, Levanat S. GSK3beta and Gli3 play a role in activation of Hedgehog-Gli pathway in human colon cancer - Targeting GSK3beta downregulates the signaling pathway and reduces cell proliferation. Biochim Biophysica Acta. 2015; 1852:2574-2584.
20. Urriola-Munoz P, Lizama C, Lagos-Cabre R, Reyes JG, Moreno RD. Differential expression and localization of ADAM10 and ADAM17 during rat spermatogenesis suggest a role in germ cell differentiation. Biological Res. 2014; 47:31.
21. Moreno RD, Urriola-Munoz P, Lagos-Cabre R. The emerging role of matrix metalloproteases of the ADAM family in male germ cell apoptosis. Spermatogenesis. 2011; 1:195-208.
22. Kampfer C, Spillner S, Spinnler K, Schwarzer JU, Terradas C, Ponzio R, Puigdomenech E, Levalle O, Kohn FM, Matzkin ME, Calandra RS, Frunqieri MB, Mayerhofer A. Evidence for an adaptation in ROS scavenging systems in human testicular peritubular cells from infertility patients. Int J Androl. 2012; 35:793-801.
23. Lu C, Zhang F, Yang H, Xu M, Du G, Wu W, An Y, Qin Y, Ji G, Han X, Gu A, Xia Y, Song L, et al. Additional genomic duplications in AZFc underlie the b2/b3 deletion-associated risk of spermatogenic impairment in Han Chinese population. Hum Mol Genet. 2011; 20:4411-4421.
24. Lu C, Wang Y, Zhang F, Wang Y, Zhang F, Lu F, Xu M, Qin Y, Wu W, Li S, Song L, Yang S, Wu D, et al. DAZ duplications confer the predisposition of Y chromosome haplogroup K* to non-obstructive azoospermia in Han Chinese populations. Hum Reprod. 2013; 28:2440-2449.
25. Gong J, Tong Y, Zhang HM, Wang K, Hu T, Shan G, Sun J, Guo AY. Genome-wide identification of SNPs in microRNA genes and the SNP effects on microRNA target binding and biogenesis. Hum Mutat. 2012; 33:254-263.
26. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009; 25:1754-1760.
27. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010; 20:1297-1303.
28. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012; 22:568-576.
29. Li Q, Wang Y, Bryant SH. A novel method for mining highly imbalanced high-throughput screening data in PubChem. Bioinformatics. 2009; 25:3310-3316.
30. Chen X, Li S, Yang Y, Yang X, Liu Y, Hu W, Jin L, Wang X. Genome-wide association study validation identifies novel loci for atherosclerotic cardiovascular disease. J Thromb Haemost. 2012; 10:1508-1514.
31. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, Sham PC. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007; 81:559-575.
32. Ji J, Qin Y, Ren J, Lu C, Wang R, Dai X, Zhou R, Huang Z, Xu M, Chen M, Wu Wei, Shen H, Hu Z, et al. Mitochondria-related miR-141-3p contributes to mitochondrial dysfunction in HFD-induced obesity by inhibiting PTEN. Sci Rep. 2015; 5:16262.
