Altered microRNA profiles of testicular biopsies from patients with nonobstructive azoospermia.

Many studies have shown that microRNAs (miRNAs) play vital roles during the spermatogenesis. However, little is known about the altered miRNA profiles of testicular tissues in nonobstructive azoospermia (NOA). Using microarray technology, the miRNA expression profiles of testicular biopsies from patients with NOA and of normal testicular tissues were determined. Bioinformatics analyses were conducted to predict the enriched biological processes and functions of identified miRNAs. The microarray data were validated by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR), the results of which were then validated with a larger sample size. Correlations between the miRNA expression levels and clinical characteristics were analyzed. Receiver operating characteristic (ROC) curve analysis was used to evaluate the diagnostic ability of miRNAs for azoospermia. Hierarchical clustering showed that 129 miRNAs were significantly differentially expressed between the NOA and control groups. Bioinformatics analysis indicated that the differentially expressed miRNAs were involved in spermatogenesis, cell cycle, and mitotic prometaphase. In the subsequent qRT-PCR assays, the selected miRNA expression levels were consistent with the microarray results, and similar validated results were obtained with a larger sample size. Some clinical characteristics were significantly associated with the expression of certain miRNAs. In particular, we identified a combination of two miRNAs (miR-10b-3p and miR-34b-5p) that could serve as a predictive biomarker of azoospermia. This study provides altered miRNA profiles of testicular biopsies from NOA patients and examines the roles of miRNAs in spermatogenesis. These profiles may be useful for predicting and diagnosing the presence of testicular sperm in individuals with azoospermia.

INTRODUCTION

According to the World Health Organization (WHO) survey, approximately 15% of couples suffer from infertility worldwide, reaching as high as 30% in some parts of the developing world, with the cause of infertility equally attributable to men and women.12 Azoospermia is a key contributor to male infertility, and it is caused by obstruction of the excurrent ductal system or testicular failure, with the latter occurring in approximately 1% of all adult men.3 Among the factors affecting male sterility, genetic factors are important causes of infertility, as 15%–30% of male infertility cases are idiopathic or thought to be caused by genetic defects. To date, studies have mainly investigated chromosomal abnormalities, Y chromosome microdeletions, and mutations and genetic polymorphisms in genes associated with infertility.4

Spermatogenesis is a complex process that occurs in the testicular seminiferous tubules and is dependent on many genetic and environmental factors. This process involves thousands of genes and includes three stages: mitosis, meiosis, and haploid cell differentiation.56 In addition to the classical posttranscriptional control of spermatogenesis by specific proteins that require a high degree of regulation, mRNA posttranscriptional regulation is also essential for the regulation of gene expression, growth, and development. In recent years, some studies have indicated that miRNAs, some short endogenous noncoding RNA molecules, play crucial roles in the posttranscriptional regulation of spermatogenesis in humans678910111213141516 and that miRNAs in seminal plasma might predict the presence of sperm in testicular tissue;1718 however, only a few miRNAs have been studied, and the altered miRNA profiles in nonobstructive azoospermia (NOA) need to be further elucidated. Although sperm can be found in some NOA patients through testicular aspiration or microdissection testicular sperm extraction, the majority of these patients are still unable to obtain sperms, which is a major problem for clinicians and has received wide attention from the society in recent years.

Therefore, in this study, we used miRNA microarray technology to compare and analyze differential miRNA expression in testicular tissues obtained from NOA and obstructive azoospermia (OA) males, respectively, the latter of which exhibit normal spermatogenesis. Bioinformatics studies, including biological functions, signal pathways, and disease enrichment analyses, were performed for the identified miRNAs, the results of which provide a solid foundation for the diagnosis and further study of the mechanisms of miRNAs in NOA.

PATIENTS AND METHODS

Patient selection

In our research, a total of 77 testicular tissue samples were obtained from men undergoing testicular sperm extraction (TESE) for assisted reproductive and/or diagnostic biopsy at Peking University Third Hospital (Beijing, China) between October 2017 and May 2018, whose demographic characteristics are shown in . All patients underwent a routine semen test (SuiJia Software, Beijing, China), according to the WHO fifth edition guidelines (2010),19 to show the absence of ejaculated sperm, and were first diagnosed with azoospermia. NOA was defined as follows: (1) no chromosomal disease or history of cryptorchidism or mumps orchitis, (2) no sperm obtained by TESE, and (3) testicular maturation arrest observed in spermatogenic cells via biopsy. The samples were classified into two groups according to the pathologic alterations of spermatogenesis observed (): the NOA group included samples demonstrating maturation arrest in spermatogenic cells (n = 39), whereas the control group samples exhibited normal spermatogenesis (n = 38). Preoperative serum follicle-stimulating hormone (FSH), luteinizing hormone (LH), and total testosterone (TT) levels as well as the Johnson's scoring system were used to evaluate testicular spermatogenesis. Our research program was approved by the Ethics Committee at Peking University Third Hospital (NO. 2017SZ-035), and each participant signed written informed consent.

RNA extraction, labeling, and hybridization

Total RNA containing small RNA was extracted from tissue using the Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified with a mirVana™ miRNA Isolation Kit (Ambion, Austin, TX, USA) according to the manufacturer's protocol. The purity and concentration of RNA were determined from OD260/280 readings using a spectrophotometer (NanoDrop ND-1000, Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was determined by 1% formaldehyde denaturing gel electrophoresis.

miRNA profiling was performed using an Agilent Human miRNA 8 × 60K Microarray Release 16.0 (Agilent Technologies, Santa Clara, CA, USA). The Agilent array was designed with eight identical arrays per slide (8 × 60K format), with each array containing probes interrogating 2549 human mature miRNAs from miRBase R21.0 (www.mirbase.org/). Each miRNA was detected by 30 replicate probes. The array also contained 2164 Agilent control probes.

The microarray experiment was performed using 13 samples from the NOA group (n = 7) and control group (n = 6), according to the manufacturer's instructions. Briefly, the miRNAs were labeled using the Agilent miRNA labeling reagent. Total RNA (200 ng) was dephosphorylated and ligated with pCp-Cy3, and the labeled RNA was purified and hybridized to the miRNA arrays. Images were scanned and gridded using an Agilent microarray scanner (Agilent Technologies), and the differentially expressed genes were analyzed using fold change threshold values of ≥2 and ≤−2 with Agilent feature extraction software version 10.10 (Agilent Technologies).

Quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) validation of differentially expressed miRNAs

After reviewing the literature and assessing the microarray fold changes of the miRNAs, four upregulated and four downregulated miRNAs were selected for PCR verification in tissue samples from a large number of individuals. A Prime Script RT kit was used to synthesize cDNA (Tiangen, Beijing, China) according to the manufacturer's specifications. miRcute Plus miRNA qPCR Kit (Tiangen) was then used to evaluate the abundance of each single miRNA using an Applied Biosystems 7500 PCR System (Thermo Fisher Scientific). All miRNA-specific reverse transcriptase and PCR primers were synthesized by Tiangen, and the process was conducted as agreed upon with the manufacturer. In short, each qRT-PCR was performed in triplicate at 95°C for 10 min, followed by 40 amplification cycles of 15 s at 95°C and 60 s at 60°C. The RNU6B snRNA primer assay (Tiangen) was used as the internal standardization control. The dynamic changes in specific selected miRNAs were measured using the 2−ΔΔCt relative quantitative method.

Bioinformatics analyses for differentially expressed miRNAs

The miRNA and mRNA array data were analyzed for data summarization, normalization, and quality control using GeneSpring software V13 (Agilent). The default 90th percentile normalization method was performed for data preprocessing. To select the differentially expressed genes, we used fold change threshold values of ≥2 and ≤−2 and a Benjamini–Hochberg corrected P = 0.05. The data were Log2-transformed and median centered by genes using the Adjust Data function in CLUSTER 3.0 software (Stanford University, Palo Alto, CA, USA) and were then further analyzed via hierarchical clustering with average linkage. To better understand the roles of the differentially expressed mRNAs, Gene Ontology (GO) categories derived from the GO database (www.geneontology.org) were determined, and pathway analysis was performed. Three families of GO terms were identified: biological process, cellular component, and molecular function.20 Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) analysis allowed us to determine the biological pathways that were significantly enriched with differentially expressed mRNAs.21 Disease enrichment was analyzed by a web server, KOBAS 2.0 (http://kobas.cbi.pku.edu.cn/), which annotates an input set of genes with putative pathways and disease relationships based on mapping to genes with known annotations.22

Statistical analyses

All data were analyzed by SPSS 20.0 (IBM SPSS Inc., Chicago, IL, USA) and reported as the mean±standard deviation (s.d.). The variables were analyzed by unpaired two-sided Student's t-tests. Correlations between miRNA expression and demographic characteristic parameters were assessed by Spearman's correlation tests. Multivariate logistic regression and receiver operating characteristic (ROC) analyses were used to evaluate the diagnostic value of the miRNAs for azoospermia. P < 0.05 was used to indicate significant differences.

RESULTS

Distinctly different miRNA profiles in patients with NOA and control group

Hierarchical clustering of the systematic variations in the expression of miRNAs between the tissues from patients with NOA and controls is shown in (adjusted P < 0.05, multiple change >2 or <1/2). There were 129 human miRNAs that were differentially expressed in samples from patients with NOA compared with those in control samples. There were 83 miRNAs expressed significantly higher, and 46 miRNAs expressed significantly lower in the samples from patients with NOA compared with the corresponding controls (). The scattered plot intuitively reflects the number of genes with significant differences. The genes labeled red are upregulated, those marked green are downregulated, and those labeled black are nonsignificant difference genes. The vertical coordinates are the logarithmic values of the corresponding samples (). A circos diagram showed the location of these miRNAs on chromosomes. Most chromosomes had the differentially expressed miRNAs, especially the X chromosome and chromosomes 1, 5, and 9 (). Thus, our data suggested that a set of miRNAs is frequently aberrantly expressed in the tissues from the patients with NOA.

Microarray-generated miRNA expression profiles of testicular tissues from the NOA and control groups. (a) Hierarchical clustering of the miRNAs produced trees with significant differences between the patients with NOA and the control group. Rows: miRNAs; columns: samples. The legend at the top depicts how the relative miRNA expression levels are illustrated: blue indicates low expression, whereas yellow indicates high expression. (b) Volcano plot visualizing miRNA expression changes in the two groups. The vertical lines indicate fold changes of 2.0 and −2.0, respectively, and the horizontal line represents P = 0.05. (c) The scattered plot intuitively reflects the number of genes with significant differences. The genes labeled red are upregulated, those marked green are downregulated, and those labeled black are nonsignificant difference genes. The vertical coordinates are the logarithmic values of the corresponding samples. NOA: nonobstructive azoospermia; miRNAs: microRNAs.

GO, KEGG pathway, and disease enrichment analyses

To further investigate the functions of the differentially expressed miRNAs, GO was used to analyze these identified miRNAs. The differentially expressed miRNAs were showed to be involved in 20 biological processes or molecular functions, including spermatogenesis, male gamete generation, and sexual reproduction (). A total of twenty KEGG pathways were enriched, including the cell cycle, mitotic prometaphase, resolution of sister chromatid cohesion, reproduction, and fertilization pathways (). The diseases with top enrichment scores included primary ciliary dyskinesia, primary microcephaly, congenital disorders of development, episodic ataxias, and spermatogenic failure (). These analyses indicated the possibility that the differentially expressed miRNAs in the tissues from the patients with NOA are involved in the important biological processes including infertility.

(a) GO, (b) KEGG pathway, and (c) disease enrichment analyses for the target genes of the differentially regulated miRNAs. The abscissa indicates the ratio of input frequency/background frequency in the enrichment analysis, the bubble size indicates the number of genes annotated in each functional item, and the color corresponds to the P-value from the enrichment analysis. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; AURKA: aurora kinase A; TPX2: targeting protein for xenopus kinesin-like protein 2; PLK1: polo like kinase 1; G2: second gap; M: metaphase; RHO: ras homologous oncogenes; GTPases: guanosine triphosphatase; IMO: interphase microtubule organization; MR: mental retardation; miRNAs: microRNAs.

Comparison of qRT-PCR and microarray expression data

To further validate the results from the microarray data, we selected eight miRNAs among the differentially expressed miRNAs and tested their expression in testicular tissues from the patients with NOA (n = 7) and corresponding controls (n = 6) using real-time RT-PCR. Based on the data from real-time RT-PCR, the expression of the eight miRNAs showed the similar tendency with those from the microarray data (Figure and ). Meanwhile, a subset of miRNAs were selected for validation in additional 64 tissues from patients with NOA and controls by qRT-PCR. The differential expression of these miRNAs between the two groups was basically consistent with the microarray results (Figure and ).

Comparison of miRNA expression. (a) Comparison of qRT-PCR and microarray expression data for four upregulated miRNAs. (b) Comparison of qRT-PCR and microarray expression data for four downregulated miRNAs. (c) Comparison of the four upregulated miRNAs in NOA and control group samples. (d) Comparison of the four downregulated miRNAs between the two groups. NOA: nonobstructive azoospermia; qRT-PCR: quantitative reverse transcriptase polymerase chain reaction; miRNAs: microRNAs.

Correlations between miRNA expression and demographic characteristic parameters

To determine whether the expression level of the selected miRNAs in the tissues from the patients with NOA was associated with the development and progression of infertility, we analyzed their relationship with the patient's clinical features, including age; the levels of FSH, LH, and TT; the testicular volume; and Johnson's score determined from testicular biopsy. The upregulated miRNAs were positively correlated with the FSH level and negatively correlated with the Johnson's score, whereas the relationships between these parameters and the downregulated miRNAs were reversed, except for those between miR-31-5p and the FSH level ().

Diagnostic efficacy of the miRNAs for azoospermia

ROC analysis was used to evaluate the diagnostic value of the eight miRNAs for azoospermia using their expression in 77 samples. The results showed that four upregulated miRNAs (miR-370-3p, miR-10b-3p, miR-539-5p, and miR-22-5p) and three downregulated miRNAs (miR-34b-5p, miR-516b-5p, and miR-122-5p) could differentiate azoospermia. We selected three miRNAs (miR-10b-3p, miR-34b-5p, and miR-122-5p, all of which had an area under curve [AUC] >0.8 and P ≤ 0.001) and analyzed them via multivariate logistic regression analysis. Remarkably, this analysis yielded a model that included miR-10b-3p and miR-34b-5p, whose sensitivity and specificity for predicting azoospermia were 76.3% and 89.5%, and 82.1% and 76.9%, respectively (). When the miR-10b-3p and miR-34b-5p values were both included in the model, an increased AUC (0.962), sensitivity (97.4%), and specificity (87.2%) were obtained ( and ).

ROC analysis of miR-10b-3p, miR-34b-5p, and miR-10b-3p + miR-34b-5p for predicting azoospermia. ROC: receiver operating characteristic.

DISCUSSION

In recent decades, researchers have begun to pay attention to miRNA, an important new regulatory factor involved in posttranscriptional gene silencing. Although miRNA accounts for only 2% of the total number of human genes, researchers speculate that these small molecules regulate nearly 30% of human genes,23 many of which are related to the occurrence of diseases.242526 The miRNAs also regulate posttranscriptional gene expression, growth, and development, maintaining the normal physiological function of organisms.727 Studies have shown that miRNAs play crucial roles in spermatogenesis, but little research has been reported and the altered miRNA profiles in NOA need to be further elucidated. In this study, using a miRNA expression microarray, 83 upregulated and 46 downregulated miRNAs were identified in testicular tissue from patients with NOA compared with normozoospermic tissue. Furthermore, four upregulated (miR-370-3p, miR-539-5p, miR-10b-3p, and miR-22-5p) and four downregulated (miR-34b-5p, miR-31-5p, miR-516b-5p, and miR-122-5p) miRNAs were selected and verified via qRT-PCR in a total 77 samples (NOA, n = 39; control, n = 38). The results kept consistent between the two methods, with the fold changes in regulation occurring in the same direction. Moreover, we found that the upregulated genes were positively correlated with the FSH level and negatively correlated with the Johnson's score, while the relationships between these parameters and the downregulated genes were reversed, except for that between miR-31-5p and the FSH level. These results further indicate that miRNA expression is closely related to testicular spermatogenesis and may be a predictive indicator of the presence of sperm in the testes for patients with azoospermia.

Bioinformatics analysis revealed that molecular processes and signaling pathways (such as the cell cycle and mitosis) that control the fate of cells were significantly disrupted in NOA.28 In this study, GO and KEGG pathway analyses were performed to study the biological functions and molecular mechanisms of miRNAs during spermatogenesis. GO analysis revealed that the differential expression of miRNAs was involved in spermatogenesis, male gamete generation, and sexual reproduction. The KEGG pathways that were found to be enriched included the cell cycle, mitotic prometaphase, reproduction, and fertilization pathways. Furthermore, disease enrichment was also involved in spermatogenic failure. These results demonstrate that miRNAs are crucial for spermatogenesis and provide a reliable reference for further research.

In the last decade, the role of miRNAs in testicular spermatogenesis had received widespread attention. Studies have shown that the miR-34 family (miR-34b and miR-34c) are highly expressed from late meiosis to the sperm stage.1129 This miRNA family is also involved in the p53 tumor suppressor network, which may play a crucial role in apoptosis and p53-mediated cell death.3031 In addition, compared with prepubertal testes, miR-34b is highly expressed in adult testes, suggesting that miR-34b plays a potential role in the differentiation of male germ cells.3233 Abu-Halima et al.3334 examined the expression of hsa-miR-34b-5p and hsa-miR-122 in human spermatozoa, and the different pathological types of human testes were verified, indicating that hsa-miR-34b-5p and hsa-miR-122 are involved in apoptosis, cell proliferation, and differentiation. These studies are consistent with our conclusions. Meanwhile, we verified that miR-370-3p, miR-539-5p, miR-10b-3p, and miR-22-5p were overexpressed and that miR-31-5p and miR-516b-5p were downexpressed in the testis tissues of NOA patients, suggesting that they may lead to spermatogenesis failure. However, these miRNAs have not previously been reported to be associated with testicular spermatogenesis and further research is still needed.

Currently, the diagnosis of azoospermia is mainly based on FSH levels or testicular volume, but the diagnostic efficacy is not satisfactory35 and requires testicular biopsy. Exploring specific biomarkers with diagnostic and predictive functions is very necessary. The presence of high concentrations of miRNAs in seminal plasma or the exosomes of semen facilitates identification and quantification, and studies have found that some miRNAs may be potential biological markers.1836 Barceló et al.36 corroborated the preferential testicular expression of miR-122-5p, miR-34c-5p, and miR-449a, whereas miR-31-5p, miR-539-5p, and miR-941 are expressed in the testis, epididymis, and prostate. In addition, Abu-Halima et al.34 revealed that the combination of miR-34b and miR-122 with other conventional tests could improve the diagnostic accuracy for detecting different forms of NOA. These miRNAs may be produced in the testicles or epididymis and could serve as adjunct biomarkers for diagnosing male infertility. Similarly, our results showed that the three selected miRNAs (miR-10b-3p, miR-34b-5p, and miR-122-5p) had good diagnostic accuracy in testicular spermatozoa. Interestingly, in this study, multivariate logistic regression analysis yielded a model that included miR-10b-3p and miR-34b-5p expression. The combination of miR-10b-3p with miR-34b-5p obviously improved diagnostic accuracy for NOA patients. However, further researches are still needed to focus on exploring the potential functions of these miRNAs to clarify their association with spermatogenesis.

CONCLUSION

We systematically analyzed and verified the altered miRNA profiles of testicular tissue from NOA patients. Our bioinformatics and correlation analyses provided an important platform for future studies of miRNAs in spermatogenesis. In particular, we demonstrated that the combination of miR-10b-3p and miR-34b-5p has high diagnostic accuracy for NOA and may serve as a biomarker for diagnosing male infertility.

AUTHOR CONTRIBUTIONS

HTZ, ZZ, HCL, and HJ designed this experiment and analyzed the data. DFL and JMM provided samples. WHT and YZY collected and analyzed the data. HTZ and ZZ participated in statistical analysis and drafting of the manuscript. HCL, KH, and HJ reviewed the article. All authors read and approved the final manuscript.

COMPETING INTERESTS

All authors declared no competing interests.

Histopathological features of the NOA and control testicular tissues. HE staining elucidated the histopathological characteristics of the (a and b) NOA testis and (c and d) the normal testis. NOA: nonobstructive azoospermia; HE: hematoxylin and eosin.

The circos diagram shows the degree of difference between different genes based on gene location. In this image, red indicates upregulated differentially expressed miRNAs, whereas green indicates downregulated differentially expressed miRNAs; moreover, the length of each column indicates the magnitude of differential expression, with longer columns representing greater differences. miRNAs: microRNAs.

Supplementary Table 1

Clinical features of patients in the microarray and validation groups (mean±standard deviation)

Parameters	Microarray group (n=13)		Validation group (n=64)
	NOA (n=7)	Control (n=6)	NOA (n=32)	Control (n=32)
Age	29.429±2.507	31.000±3.162	31.000±4.621	32.531±5.858
FSH	15.909±4.517**	5.735±3.281	16.550±8.381**	6.132±3.207
LH	5.597±2.089	5.415±2.882	6.086±2.817*	4.049±2.060
TT	9.606±2.307	14.180±7.370	9.563±3.294	11.060±4.667
Testicular volume	9.000±2.517*	14.833±2.927	8.750±3.592**	14.500±4.166
Johnson’s scores	2.429±0.787**	8.667±1.033	2.813±1.693**	7.875±1.184

*P≤0.01, **P≤0.001. No significant between-group differences were found for age and TT (both P>0.05). Although the mean testosterone level of the control group was higher than that of the NOA group, the difference was not significant. However, the FSH, LH, and TT levels; the testicular volumes; and Johnson’s scores had significant between-group differences (P>0.05, ). NOA: nonobstructive azoospermia; FSH: follicle-stimulating hormone; LH: luteinizing hormone; TT: total testosterone; s.d.: standard deviation

Table 1

Correlations between microRNA expression and demographic characteristic parameters (n=77)

miRNAs	Correlation coefficients
	Age	FSH	LH	TT	Testicular volume	Johnson’s score
miR-370-3p	−0.040	0.329**	0.058	0.028	−0.266*	−0.518***
miR-10b-3p	0.211	0.230*	0.164	−0.068	−0.201	−0.237*
miR-539-5p	0.109	0.393***	0.313**	−0.084	−0.419***	−0.539***
miR-22-5p	0.176	0.277*	0.145	−0.008	−0.361**	−0.248*
miR-34b-5p	−0.023	−0.575**	−0.273*	0.209	0.361**	0.574***
miR-31-5p	0.062	−0.072	−0.144	−0.102	−0.025	0.233*
miR-516b-5p	0.073	−0.275*	−0.088	0.076	0.323**	0.316**
miR-122-5p	0.075	−0.379**	−0,147	−0.101	0.248*	0.468***

*P≤0.05, **P≤0.01, ***P≤0.001. Correlations among the selected miRNAs (2−ΔΔCt expression value) and age; FSH, LH, and TT levels; testicular volume; and Johnson’s score were assessed by Spearman’s correlation analysis. FSH: follicle-stimulating hormone; LH: luteinizing hormone; TT: total testosterone; miRNAs: microRNAs

Table 2

The predictive efficiency of differential microRNAs expression between the nonobstructive azoospermia (n=39) and control groups (n=38)

miRNAs	AUC	95% CI	P	Sensitivity (%)	Specificity (%)
miR-370-3p	0.771	0.668–0.874	<0.001	84.2	61.5
miR-539-5p	0.681	0.559–0.803	0.006	84.2	51.3
miR-10b-3p	0.827	0.736–0.917	<0.001	76.3	82.1
miR-22-5p	0.686	0.567–0.805	0.005	73.7	59.0
miR-34b-5p	0.891	0.822–0.960	<0.001	89.5	76.9
miR-31-5p	0.597	0.608–0.833	0.144	NS	NS
miR-516b-5p	0.721	0.608–0.833	0.001	55.3	76.9
miR-122-5p	0.841	0.756–0.926	<0.001	89.5	64.1
miR-10b-3p+miR-34b-5p	0.962	0.927–0.998	<0.001	97.4	87.2

ROC curve analysis showed the predictive efficiency of the miRNAs for identifying azoospermia. NOA: nonobstructive azoospermia; ROC: receiver operating characteristic; AUC: area under curve; CI: confidence interval; NS: not significant; miRNAs: microRNAs