Circulating microRNAs in seminal plasma as predictors of sperm retrieval in microdissection testicular sperm extraction.

Background: Because of focal spermatogenesis in some nonobstructive azoospermia (NOA) patients, testicular spermatozoa can be retrieved by microdissection testicular sperm extraction (micro-TESE) for intracytoplasmic sperm injection (ICSI) to achieve successful fertilization. Currently, testicular biopsy is widely performed for the prognosis of micro-TESE; however, it might miss foci with active spermatogenesis because of the 'blind manner' of puncture, highlighting the needs for biomarkers that could indicate actual spermatogenesis conditions in the testis. Thus, we screened microRNAs in the seminal plasma for potential biomarkers to provide a non-invasive and reliable preoperative assessment for micro-TESE.
Methods: We screened the seminal plasma microRNAs from NOA patients with and without sperm retrieval (n=6 in each group) together with fertile men (n=6) by RNA sequencing, and the selected microRNAs were validated by quantitative polymerase chain reaction (qPCR). Next, a predictive model was established by performing ordered logistic regression using the qPCR data of 56 specimens, and the predictive accuracy of this model was evaluated using 40 more specimens in a blind manner.
Results: Four microRNAs (hsa-miR-34b-3p, hsa-miR-34c-3p, hsa-miR-3065-3p, and hsa-miR-4446-3p) were identified as biomarkers, and the predictive model Logit = 2.0881+ 0.13448 mir-34b-3p + 0.58679 mir-34c-3p + 0.15636 mir-3065-3p + 0.09523 mir-4446-3p was established by machine learning. The model provided a high predictive accuracy (AUC =0.927). Conclusions: We developed a predictive model with high accuracy for micro-TESE, with which NOA patients might obtain accurate assessment of spermatogenesis conditions in testes before surgery. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Approximately 10–15% of couples experience infertility, and approximately 50% of these cases are caused by male infertility, of which approximately 10% represent the most severe manifestation, named nonobstructive azoospermia (NOA) (1). NOA is characterized by the complete absence of sperm in semen due to impaired spermatogenesis in the testes. However, because of focal spermatogenesis, testicular spermatozoa can be retrieved by microdissection testicular sperm extraction (micro-TESE) for intracytoplasmic sperm injection (ICSI) to achieve successful fertilization in some NOA patients (2-5). The major challenge with micro-TESE is to find predictors of active spermatogenesis foci indicating the likelihood of recovering sperm (6). It is essential to evaluate the sperm retrieval likelihood to determine whether the patient is currently suitable for micro-TESE. In particular, for patients receiving certain medications to enhance spermatogenesis recovery before the operation, the optimal operation time must be estimated according to the prognosis (7). Currently, testicular biopsy is widely performed for the prognosis of micro-TESE; however, it might miss foci with active spermatogenesis because of the “blind manner” of puncture (8-11). Additionally, as an invasive method, testicular biopsy may lead to severe complications (12). Therefore, a more reliable and safer prognostic method is needed.

Seminal plasma is the liquid component of semen, which arises from secretions from the seminal vesicles (~65%), prostate (~25%), testes and epididymides (~10%) and periurethral glands (~1%) (13). Secretions from the testes and epididymides provide molecules related to spermatogenesis and spermatozoa maturation, which could be potential biomarkers prompting spermatogenesis conditions (13-15). Several studies have reported proteins, circulating cells and cell-free DNA/RNA as potential biomarkers of nonobstructive azoospermia subtypes (16-18). Among these factors, circulating microRNAs could be ideal candidates because of their smaller size and encapsulation within proteins as they are less susceptible to degradation than other RNA species. MicroRNAs are a family of noncoding RNAs of approximately 19–23 nucleotides (nt) long that function as regulators of messenger RNA (mRNA) stability and translation (19). Many studies have indicated that microRNAs are involved in spermatogenesis and that aberrant expression of some microRNAs is related to male infertility (20-25). MicroRNAs have been reported as reliable biomarkers in various biopsies due to their stable characteristics (26,27). As one of the most important regulators of spermatogenesis, they have also been detected in human spermatozoa and seminal plasma (28). Several studies have investigated microRNAs in seminal plasma using either microarray or exosome microRNA profiling for indicators of spermatogenesis (29-31). However, neither microarray nor exosome microRNA profiling represents the total microRNAs that exist in the seminal plasma. Furthermore, those studies were limited by one or more of the following reasons: a lack of model establishment, a limited number of patients, and a lack of clinical validation.

In the present study, we screened whole microRNAs in seminal plasma from different groups of patients by improved absolute quantitative high-throughput small RNA sequencing (sRNA-seq) for more reliable biomarkers, and 4 microRNAs were identified as potential predictors through multiplex reverse transcription-quantitative polymerase chain reaction (RT-qPCR). A microRNA predictive model was established by ordered logistic regression through machine learning with qPCR data. The model was used to assess testicular spermatogenesis ratings to accurately estimate the likelihood of sperm retrieval in micro-TESE. Independent clinical validation of the model was performed to evaluate the predictive accuracy of the model (32). We present the following article in accordance with the STARD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-5100/rc).

Methods

Patients and study design

All seminal plasma samples in this study were collected from Guangdong Provincial Reproductive Science Institute (Guangdong Provincial Fertility Hospital). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Research Ethics Committee of Guangdong Provincial Reproductive Science Institute (Guangdong Provincial Fertility Hospital) (No. 2019025), and individual consent for this retrospective analysis was waived.

The samples were divided into 3 sets for different purposes in the study, including the screening set, training set and test set (). Each set contained 3 groups that were divided according to the following testicular spermatogenesis conditions: (I) the “good” (G) group comprised fertile semen from the sperm bank, representing normal spermatogenesis; (II) the “fair” (F) group comprised NOA patients with 100% sperm retrieval in micro-TESE, representing hypospermatogenesis (existence of focal spermatogenesis); and (III) the “poor” (P) group comprised NOA patients with no sperm retrieval in micro-TESE, representing no active spermatogenesis. The screening set was used to identify significantly differentially expressed microRNAs as potential biomarkers for spermatogenesis. The training set was used for machine learning to construct ordered logistic regression modeling and establish a predictive model. The test set evaluated the model in which the RT-qPCR data of this set were applied to the model for prognosis, and the output result was compared with the clinical micro-TESE record to test the accuracy.

Figure 1

Study design. Flow diagram of the study. RT-qPCR, reverse transcription-quantitative polymerase chain reaction; UMI, unique molecular index; sRNA-seq, small RNA sequencing.

Seminal plasma small RNA purification

Polysaccharides make seminal plasma highly viscous, and routine small RNA (sRNA) purification methods for liquid samples, such as the solid-phase reversible immobilization (SPRI) bead-based method and TRIzol-LS method, produce significantly low amounts of sRNAs with large amounts of polysaccharide retention. To improve the quantity and quality of seminal plasma sRNAs, we combined the TRIzol-LS precipitation and column purification methods (Appendix 1). The concentration and purity of the enriched sRNAs were checked with NanoDrop One (Thermo Fisher Scientific, Waltham, MA, USA). The size distribution of the sRNAs was detected on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) (Table S1).

Absolute quantitative high-throughput sRNA sequencing (HTS) using a unique molecular index (UMI) and bioinformatic analysis

MicroRNA library construction for HTS includes a PCR amplification step, which can introduce a vast amount of bias because of the complicated structures of microRNAs. For seminal plasma microRNA, which has a low input of the starting materials, more cycles of PCR amplification are needed to obtain a sufficient amount for sequencing. Increased PCR amplification cycles produce more bias than microRNAs from other samples. To reduce PCR-derived artifacts, we adopted absolute quantitative microRNA sequencing with UMIs (Figure S1A) (33). This method could distinguish original copies of microRNAs from their amplification products by barcoding a random hexamer sequence as the UMI to each microRNA before PCR. After amplification, although multiple counts of each UMI could be observed, we precisely quantified the original copy number of each microRNA by simply counting each UMI. This UMI-based microRNA sequencing provided us with more reliable differentially expressed genes (DEGs) between groups and increased the efficiency of the validation progress by reducing the number of false DEGs. For normalization, we synthesized cel-miR-39 (TaKaRa Bio., Dalian, China) as a spike-in reference microRNA. Exactly 250 fmol of cel-miR-39 was added to each seminal plasma sample before sRNA purification. Approximately 100 ng of sRNAs per sample was used in library construction for HTS. During library construction, UMI-labeled primers were added (Figure S1A). All the samples were sequenced on a BGISEQ-500 system (Beijing Genomics Institution, Shenzhen, China), and bioinformatic analysis was performed following the standard pipeline (Figure S1B).

Multiplex RT-qPCR for microRNAs

Five microliters of sRNA were used as the template for a 20-µL reaction. RT-qPCR was designed based on the stem-loop method (Figure S2) (34). All the primers and probes used in this study are provided in Table S2. A multiplex RT reaction with 13 microRNA-specific stem-loop primers (SLPs) was performed, and an optimized pulsed RT reaction condition was applied to reduce nonspecific reactions and increase the RT efficiency (Appendix 1) (35). To increase the sensitivity of the qPCR, preamplification (pre-PCR) of the complementary DNA (cDNA) templates was necessary. The total RT product was used as a template, together with a specific forward primer (FP) cocktail (mixture of 13 microRNA-specific primers) and universal reverse primer (URP) (Appendix 1). To remove additional primers and deoxy-ribonucleoside triphosphate (dNTPs) in the pre-PCR product, which might affect qPCR, we performed cleanup treatment with ExoSAP-IT reagent (Thermo Fisher Scientific, Waltham, MA, USA; #75001) (Appendix 1). The products were purified using a MinElute PCR Purification Kit (Qiagen, Hilden, Germany; #28004) to obtain template pools for qPCR. Fluorescent dye-based qPCR was performed for candidate biomarker selection, and hydrolysis probe-based qPCR was performed for the construction of the predictive model (Appendix 1). To test the sensitivity of the stem-loop RT-qPCR-based microRNA expression profiling method in the multiplex format, we constructed the standard curve and R2 of each microRNA qPCR (Figure S3).

Statistical analysis

The input data for the analysis were the ∆Cq values (the Cq values of specific microRNAs minus the Cq values of the reference microRNA cel-miR-39). The fold change (FC) of the expression level was calculated as 2−∆Cq or 2−∆∆Cq. All statistical analyses were performed using R script. We established a predictive microRNA model based on training data using the function Polr (proportional odds logistic regression) of the MASS package in R-studio. A model fitting test was performed to evaluate the model. Receiver operating characteristic (ROC) curves were constructed to test the predictive accuracy of the model.

Results

Patient characteristics

The baseline characteristics of the NOA patients in this study are presented in . Testicular sperm was successfully retrieved in the fair group (age 30±4), and no sperm was retrieved in the poor group (age 31±5). There was no significant difference in the distribution of age, female partner age or testicular volume between the fair group and poor group. The mean follicle-stimulating hormone (FSH) and luteinizing hormone (LH) levels in the fair group (FSH 15.33 IU/L, 95% CI: 12.22–18.43; LH 7.02 IU/L, 95% CI: 5.06–8.97) were not significantly different (P=0.84; P=0.55) from those of the poor group (FSH 15.75 IU/L, 95% CI: 12.60–18.90; LH 7.87 IU/L, 95% CI: 5.64–10.09).

Table 1

Baseline characteristics of enrolled NOA patients

Characteristics (mean ± SD)	Fair group (n=48)	Poor group (n=68)
Male age (years)	30±4	31±5
Female partner age (years)	28±5	28±4
Infertile duration (years)	4±2	4±2
Testicular volume (mL)	14±3*	12±5*
Serum FSH (IU/L)	15.33±8.01*	15.75±8.12*
Serum LH (IU/L)	7.02±5.05*	7.87±5.73*

*, no significant difference, Student’s t-test. NOA, nonobstructive azoospermia; FSH, follicle-stimulating hormone; LH, luteinizing hormone.

MicroRNA screening for candidate biomarkers

Our RNA-seq detected 786 known microRNAs from 18 samples, among which 158 were identified as differentially expressed microRNAs among groups based on the bioinformatic analysis with DEGseq, and 80 were defined as significantly differentially expressed microRNAs based on the expression FC (log2FC >2; P<0.01) (). We further narrowed the range of candidates to 12 microRNAs (hsa-miR-34b-3p, hsa-miR-676-3p, hsa-miR-376a-3p, hsa-miR-449b-3p, hsa-miR-202-5p, hsa-miR-942-5p, hsa-miR-518e-3p, hsa-miR-891b, hsa-miR-4446-3p, hsa-miR-3065-3p, hsa-miR-30b-3p, and hsa-miR-34c-3p) based on their expression levels and patterns. We first eliminated 41 low-abundance microRNAs (expression level <20) and then divided the other 39 microRNAs into 4 groups based on their expression patterns in the 3 groups (). We selected 2 groups of microRNAs whose expression patterns were consistent with the trend of spermatogenic situations from good to poor (). The selected 2 groups contained 12 microRNAs in total, among which 9 microRNAs showed a positive correlation with active spermatogenesis, and 3 showed an opposite trend (). In summary, 12 significantly differentially expressed microRNAs were identified as candidates for further validation via RT-qPCR.

Figure 2

Selection of differentially expressed microRNAs from sRNA-seq. (A) A total of 489 microRNAs were identified as differentially expressed microRNAs among 18 samples. (B) Intergroup analysis identified 80 differentially expressed microRNAs among the good, fair and poor groups; sRNA-seq, small RNA sequencing.

Figure 3

Expression patterns of the selected microRNAs in the three groups. Four patterns were identified based on the expression level changes from the “good” group to the “poor” group.

Differential expression profile of selected microRNAs

The 12 candidate microRNAs were tested by RT-qPCR with the 3 groups of sequencing samples to validate the FC of the expression level. We first tested the sensitivity of the stem-loop RT-qPCR-based microRNA expression profiling method in multiplex format. The standard curve and R2 of each microRNA qPCR were constructed. Among the 12 microRNAs, hsa-miR-202-5p, hsa-miR-518e-3p and hsa-miR-891b had R2 values of 0.9748, 0.9535 and 0.9312, respectively, and the other 9 microRNAs all had R2 values >0.99, indicating a sensitive and optimal amplification efficiency (Figure S3). Among these 12 microRNAs, only 4 (hsa-miR-34b-3p, hsa-miR-34c-3p, hsa-miR-3065-3p, and hsa-miR-4446-3p) displayed significantly differential expression patterns (). The expression levels of hsa-miR-34b-3p and hsa-miR-34c-3p were significantly decreased from the “G” to the “P”, indicating that they are potential biomarkers for active spermatogenesis. The FC of the expression level of hsa-miR-34b-3p for “F”/“G” was 0.12 (P<0.001) and that for “P”/“G” was 0.03 (P<0.001). The FC of the expression level of hsa-miR-34c-3p for “F”/“G” was 0.62 (P<0.05) and that for “P”/“G” was 0.09 (P<0.001). The expression levels of hsa-miR-3065-3p and hsa-miR-4446-3p were significantly increased from the “G” group to the “P” group, indicating the severity of impaired spermatogenesis. The FC of the expression level of hsa-miR-3065-3p for “F”/“G” was 2.94 (P<0.01) and that for “P”/“G” was 3.44 (P<0.05). The FC of the expression level of hsa-miR-4446-3p for “F”/“G” was 2.66 (P<0.001) and that for “P”/“G” was 3.36 (P<0.01). Therefore, hsa-miR-34b-3p, hsa-miR-34c-3p, hsa-miR-3065-3p and hsa-miR-4446-3p were identified as candidate biomarkers for model establishment with the training set.

Figure 4

Selection of potential predictors by RT-qPCR. Expression profile of the 12 selected microRNAs in the different groups. RT-qPCR data of the 12 selected microRNAs showed the relative expression levels in the “good”, “fair” and “poor” groups. Among these 12 microRNAs, miR-34b-3p and miR-34c-3p were significantly downregulated from the “good” group to the “poor” group, and miR-4446-3p and miR-3065-3p were significantly upregulated from the “good” group to the “poor” group. (***P<0.001, **P<0.01, *P<0.05). RT-qPCR, reverse transcription-quantitative polymerase chain reaction.

Establishment of the microRNA predictive model

Higher expression levels of hsa-miR-34b-3p and hsa-miR-34c-3p prompted more normal spermatogenesis and a higher possibility of retrieving sperm from micro-TESE. In contrast, higher expression levels of hsa-miR-3065-3p and hsa-miR-4446-3p implied more seriously impaired spermatogenesis and a lower likelihood of retrieving viable sperm. Thus, to examine how the expression levels of these 4 microRNAs predict the sperm retrieval outcome of micro-TESE, we performed an ordered logistic regression analysis. The input data for the ordered logistic regression analysis were the ∆Cq values of the 4 microRNAs from the training set (G group =6; F group =20; P group =30). We established the formula by estimating regression coefficients and the intercept in R using the Polr function in the MASS package (). The predicted probabilities of being estimated as “good”, “fair” and “poor” were denoted as “Pg”, “Pf” and “Pp”, respectively. According to the ordered logistic regression, we established the following formulas: Logit (Pg) = 2.9414 + 0.08932 hsa.miR.34b.3p + 0.59866 hsa.miR.34c.3p + 0.23658 hsa.miR.3065.3p + 0.19814 hsa.miR.4446.3p; Logit (Pg + Pf) = 3.9007 + 0.08932 hsa.miR.34b.3p + 0.59866 hsa.miR.34c.3p + 0.23658 hsa.miR.3065.3p + 0.19814 hsa.miR.4446.3p. Pp = 1 − (Pg + Pf).

Table 2

The predictive formula model for micro-TESE

Variable	Value	Std.Error	t value
Coefficients
hsa.miR.34b.3p	0.13448	0.3400	0.3956
hsa.miR.34c.3p	0.58679	0.3704	1.5843
hsa.miR.3065.3p	0.15636	0.4112	0.3802
hsa.miR.4446.3p	0.09523	0.3542	0.2688
Intercepts
G|F†	2.0881	2.7436	0.7611
F|P‡	3.0533	2.7656	1.1040

polr(formula = as.ordered(result) − hsa.miR.34b.3p + hsa.miR.34c.3p + hsa.miR.3065.3p + hsa.miR.4446.3p). †, Logit(Pg) = 2.0881 + 0.13448 hsa.miR.34b.3p + 0.58679 hsa.miR.34c.3p + 0.15636 hsa.miR.3065.3p + 0.09523 hsa.miR.4446.3p. ‡, Logit(Pg + Pf) = 3.0533 + 0.13448 has.miR.34b.3p + 0.58679 hsa.miR.34c.3p + 0.15636 hsa.miR.3065.3p + 0.09523 hsa.miR.4446.3p. Pp = 1 − Pg-Pf. micro-TESE, microdissection testicular sperm extraction.

Accuracy test of the predictive model

To evaluate the performance of the predictive model, we tested 40 samples with the model in a “blind” manner, meaning no clinical information was provided beforehand. We compared the results acquired from the predictive model to the clinical micro-TESE results. The clinical information revealed the following: 6 samples were from fertile men, denoted as “good”; 12 samples were from NOA patients with positive sperm retrieval, denoted as “fair”; and 22 samples were from NOA patients with negative sperm retrieval, denoted as “poor”. The prediction for the 40 samples from our model showed that all 6 fertile men were identified as in the “good” group, and 21 out of 22 patients with negative sperm retrieval were identified as in the “poor” group, indicating 100% and 95.45% accuracy rates in predicting patients who might show 100% and 0% sperm retrieval rates in micro-TESE, respectively (). For the 12 patients with positive sperm retrieval in micro-TESE, 8 were identified as in the “fair” group, 2 were identified as in the “good” group and 2 were identified as in the “poor” group. Because we claimed that both the “good” and “fair” groups represented the existence of active spermatogenesis in the testes, both samples appeared in the “good” and “fair” groups, implying a positive sperm retrieval likelihood in micro-TESE. Therefore, the accuracy rate for predicting sperm retrieval was 83.33%. We also calculated the sensitivity and specificity of the predictive model, which were 94.12% and 91.30%, respectively (). ROC curves were constructed based on the predictive performance of identifying positive sperm retrieval versus negative sperm retrieval (AUC =0.927), as well as for identifying “good” versus “poor” (AUC =0.955), “fair” versus “poor” (AUC =0.913) and “good” versus “fair” (AUC =0.707) ().

Table 3

Predictive model output

Predictive result	Good	Fair	Poor	Total
Good	6	2	0	8
Fair	0	8	1	9
Poor	0	2	21	23
Total	6	12	22	–
Accuracy	100%	66.67%	95.45%	–

Table 4

Predictive model evaluation

Predictive result	Clinical result
	Positive (+)	Negative (−)	Total
Positive (+)	16	2	18
Negative (−)	1	21	22
Total	17	23	–

Predictive result: FPR 8.70%; FNR 5.88%; TPR 88.89%; TNR 95.45%. Clinical result: sensitivity 94.12%; specificity 91.30%. FPR, false positive rate; FNR, false negative rate; TPR, true positive rate; TNR, true negative rate.

Figure 5

ROC curve analysis for the performance of the predictive model. (A) ROC curve for the predictive model discriminating the “success” group from the “fair” group, AUC 0.927; (B) ROC curve for the predictive model discriminating the “good” group from the “poor” group, AUC 0.955; (C) ROC curve for the predictive model discriminating the “fair” group from the “poor” group, AUC 0.913; (D) ROC curve for the predictive model discriminating the “good” group from the “fair” group, AUC 0.707. AUC, area under the curve. 0.5< AUC <1 indicates for good predictive value. ROC, receiver operating characteristic; AUC, area under the curve.

Discussion

As an effective method to retrieve testicular sperm from NOA patients, micro-TESE appears to be safe and has a higher yield of sperm than other methods. However, positive sperm retrieval was not possible from all patients. Thus, providing doctors and patients with a method for assessing the likelihood of retrieving sperm is essential, particularly when the patient requests another operation attempt following a failed sperm retrieval. Testicular biopsy, testicular ultrasound and hormone levels are assessment methods before micro-TESE; however, because of the heterogeneity of testicular tissue, these methods lack particular effectiveness (36). Several studies have been performed to identify molecular biomarkers of active spermatogenesis from seminal plasma, such as PIWI-interacting RNAs (piRNAs), mRNAs and microRNAs; however, these studies either lacked model construction or only established binary classification models (37).

In this study, we performed an ordered logistic regression with the 4 microRNAs identified from small RNA-seq and RT-qPCR for rating testicular spermatogenesis. With the predictive model we developed, patients could be sorted into 3 groups (good, fair and poor) based on their spermatogenesis conditions, allowing them to receive the appropriate treatment for a higher sperm retrieval rate. The model identified the “good” and “poor” samples with higher accuracy, and the more extreme the situation was, the more accurate the result was. However, for moderate-quality samples, the accuracy of the model must be improved. To improve the accuracy of micro-TESE outcome prediction, other factors should be considered together with the microRNA model to form a “predictive system”. For example, seminal plasma piRNAs, mRNA elements, long noncoding RNAs (lncRNAs), circulating DNAs, exfoliative cells, together with the histological diagnosis could be considered together for outcome prediction. We also summarized the histological results of the patients with different micro-TESE outcomes in this study (Table S3, Figure S4). Multiple panels should be established and combined by machine learning to make more accurate predictions. For the outcomes of assisted reproduction after micro-TESE, if testicular spermatozoa are retrieved by micro-TESE, NOA patients could achieve fatherhood via intracytoplasmic sperm injection (ICSI). The time of sperm retrieval and the time of egg retrieval from the partner might not be the same. There are 3 ways to select the timing of ICSI: egg and sperm retrieval in the same cycle, frozen egg recovery combined with sperm retrieval, and frozen sperm recovery combined with egg retrieval. Each method has certain limitations, such as the recovery rate of frozen eggs and sperm. In this study, according to the clinical records, among the 53 NOA patients who achieved successful sperm retrieval, 24 couples underwent ICSI, and 17 couples achieved a successful pregnancy. The other 29 patients froze their sperm for later use.

Conclusions

In this study, we established a predictive model to detect the possibility of sperm retrieval outcomes in micro-TESE. Our study has demonstrated that this model could provide considerable clinical value in identifying samples with a high likelihood of sperm retrieval and is valuable to direct the selection of the optimal time for operation planning.

The article’s supplementary files as