Seminal plasma piRNA array analysis identifies possible biomarker piRNAs for the diagnosis of asthenozoospermia

Asthenozoospermia (AZS) is characterized by reduced sperm motility, and its pathogenesis remains poorly understood. Piwi-interacting RNAs (piRNAs) have been recognized to play important roles in spermatogenesis. However, little is known about the correlation of piRNAs with AZS. In this study, small RNA sequencing was performed on samples from AZS patients and fertile controls. Real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR) was used to validate the small RNA-seq results. Bioinformatics analyses were performed to predict the functions of differentially expressed piRNAs. Logistic regression models were constructed, and receiver operating characteristic curve (ROC) analysis was used to evaluate their diagnostic performance. A total of 114 upregulated and 169 downregulated piRNAs were detected in AZS patients. GO and KEGG analyses showed that the differentially expressed piRNAs were mainly associated with transcription, signal transduction, cell differentiation, metal ion binding and focal adhesion. These results were verified by qRT-PCR of 8 selected piRNAs. Among the differentially expressed piRNAs tested, piR-hsa-32694, piR-hsa-26591, and piR-hsa-2912 produced qRT-PCR results consistent with the sequencing results in AZS compared with controls in the first cohort, whereas the other three genes did not show significant differences in expression. power piwi Hiwi2 PIWIL3, clear the piwi male in addition its main in transposon silencing, the piwi/piRNA mechanism the degradation of a large number of mRNA transcripts in various germ cells through miRNA- or/and siRNA-like . In this study, we predicted the target genes of piRNAs that were differentially expressed between the control group and the AZS group to explore how the differentially expressed piRNAs might change sperm function. In addition, enrichment analysis based on these predicted target genes showed that cell macromolecule biosynthesis, spermatogenesis, male gametogenesis, cytoplasmic transport and RNA metabolism were changed in the AZS group. KEGG pathway analysis showed that the differentially expressed piRNA targets were involved in focal adhesion kinase (FAK), actin cytoskeleton, axon guidance, the ErbB signalling pathway and other signal transduction pathways. An interesting study showed that FAK coordinates the transport of spermatozoa and prealbumin spermatids through the epithelium along the functional axis of the acroplasmic thinning (ES)-blood testicular barrier basement membrane (BM) 24 . A study on pig sperm found that many differentially expressed mRNAs related to infertility are involved in metabolic processes 25 . Specifically, during the sperm deformation stage, macromolecule metabolism and transportation are well controlled, and the production of fertile sperm is ensured through 26S proteasome exchange, digestion and absorption of nuclear protein, and elimination of RNA through phagocytosis of cytoplasmic droplets, residues and Sertoli cells . Our results also support the view that focal adhesion kinase and metabolism of fatty acids and other macromolecules key

piRNAs was further analysed using ROC analysis; piR-hsa-26591 exhibited an area under the ROC curve (AUC) of 0.913 (95% CI: 0.795-0.994). Logistic regression modelling and subsequent ROC analysis indicated that the combination of the 4 piRNAs reached good diagnostic efficacy (AUC: 0.935).

Introduction
According to a report evaluating reproductive system diseases that was published by the World Health Organization (WHO) in 2004, one in four couples in developing countries experience infertility. Infertility often causes great psychological strain, which can eventually lead to mental illness. Approximately 50% of cases of infertility are caused by male factors 1 . AZS is one of the major causes of human male infertility and is characterized by impaired motility of sperm, which causes infertility by preventing the sperm from reaching the egg 2 . However, the pathogenesis of AZS has not been absolutely elucidated. Recently, studies have shown that noncoding RNAs, especially piRNAs, play important roles in the pathogenesis of AZS.
Piwi-interacting RNAs (piRNAs) are a unique class of small noncoding RNAs that are 26-30 nt in length and bind to Piwi proteins. These short RNAs were originally discovered in mouse testicular cells and have been shown to play key roles in germ cell maintenance in mammals. Mammalian spermatogenesis is a highly regulated developmental process that can be divided into three steps, in chronological order: meiosis, mitosis and differentiation. This process produces haploid gametes. Abnormal expression of piRNAs in testicular cells may lead to the failure of spermatogenesis. Heyn et al found that five piRNAs (DQ589977, DQ591415, DQ598918, DQ601291 and DQ601609) were downregulated in cases of spermatogenic failure and that this downregulation was negatively correlated with the degree of methylation of the piwi2 and Tudor domain protein 1 (tdrd1) gene promoters. Bioinformatics analysis of these five piRNAs also revealed some other potential target genes, such as IL-16, kallikrein 1 (KLK1), G protein coupled receptor 156 (gpr156), histone cluster 1h2aa (hist1h2aa), rab24 (Ras oncogene family member), and sphingomyelin phosphodiesterase 3 (smpd3) 3 . Zeeba et al observed that the rs508485 mutation in the hiwi2 gene influences the stability of mRNA or changes the binding affinity of miRNA and also resulted in idiopathic nonobstructive azoospermia 4 . Hong and her colleagues identified 5 piRNAs as biomarkers for the diagnosis of AZS 5 . However, little is known about the expression patterns and functions of piRNAs in AZS.
Therefore, this study aimed to investigate the expression profile and functions of piRNAs in AZS patients. Bioinformatics analysis was also performed to describe the piRNA/mRNA interaction network, biological process and signal pathways. These results may provide potential targets for the development of novel diagnostic and therapeutic strategies against AZS. We also assessed the diagnostic value of the differentially expressed piRNAs in the seminal plasma as an additional tool for AZS.

Ethics statement
This study was approved by the Ethics Committee of The Second Affiliated Hospital of Nanchang University, China. All donors gave informed consent for the collection and use of their samples for this study.

Sample collection
We selected 4 patients diagnosed with male sterility due to AZS and 3 normal healthy males at the Maternal and Child Health Hospital in Jiangxi Province. Semen samples (6-8 ml) were collected in a sterile container. We used a sperm quality analyser to evaluate the basic semen parameters, including sperm morphology, concentration, motility, viability and quantity. According to the World Health Organization (WHO, 2010, 5th Edition) diagnostic criteria for AZS, the proportions of sperm with fast progressive motility (A < 25%) and total progressive motility (A + B < 50%) in fresh ejaculates were determined. Normal sperm samples were defined as those with the following parameters: semen concentration ≥ 15 × 10 6 /ml, semen volume ≥ 2 ml, pH ≥ 7.2, progressive motility (PR) ≥ 32%, and PR + nonprogressive motility (NP) ≥ 40%. The age range of both the AZS patients and the normal controls was 28-40 years old; the patients were also well matched to controls in terms of drinking and smoking status. Controls with underlying disease or abnormal sperm parameters were excluded from this study. Patients with gonadal dysplasia due to congenital or inherited disease were excluded from the study. The semen samples were liquefied at 37 ℃ for 30 min and separated by a discontinuous 45%/85% Percoll density gradient to eliminate the contamination from round cells, including germ cells and white blood cells. The spermatozoa in the lower (85% Percoll) layer were collected and washed 3 times with phosphate-buffered saline (PBS). The quality of the purified sperm was examined by microscopy.

RNA isolation and qRT-PCR assays
Total RNA was extracted from semen samples using TRIzol reagent according to the manufacturer's specifications. The yield of RNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and the RNA integrity was evaluated using electrophoresis on an agarose gel stained with ethidium bromide.  Table 2. The piRNA expression levels were normalized to that of (input the reference gene, e.g., RNU6B, 5S rRNA) and were calculated using the 2-ΔΔCt method (Livak and Schmittgen, 2001).

Bioinformatic analysis
Briefly, small RNA sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China). The base reads were converted into sequence data (also called raw data/reads) by base calling. Low-quality reads were filtered out, and reads with 5' primer or poly(A) contaminants were removed. The reads without 3' adapter and insert tags and reads shorter than 15 nt and longer than 41 nt from the raw data were filtered, and clean reads were obtained. The clean reads were matched to the reference human (Homo sapiens) genome 6 . Sequences corresponding to known piRNAs were identified by perfect sequence matching to the piRBase database (piRBase release v2.0) 7 . miRNA sequences were determined by matching to the miRNA database (miRBase 19.0), and other small RNAs (rRNA, tRNA and snoRNA) were annotated by BLAST search against the Rfram database [8][9][10] .
Differentially expressed piRNAs between the AZS group and the normal control group were identified using DEG-seq software. Genes with q-values (FDRs) lower than 0.05 and fold-changes higher than 2.0 were considered differentially expressed 11,12 . The targets of differentially expressed piRNAs were predicted by using the miRNA database software 9,13,14 with the following parameters: S ≥ 150 ΔG, ≤ − 30 kcal/mol, and strict 5' seed pairing.
Functional annotation of predicted target genes was performed by Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes pathway analysis was conducted 15 . Corrected P values < 0.05 were considered to indicate significantly enriched differentially expressed piRNA-target-genes.

Statistical analysis
Differences in demographic variables and expression levels of piRNA between the AZS group and the control group were evaluated with the chi-squared test. Fisher's exact test was used to assess the significance of GO terms or KEGG pathway identifiers in the differentially expressed genes.
To evaluate the specificity and sensitivity of each piRNA in the diagnostic value of AZS. we analysed ROC curves and AUC. A large Youden's index (sensitivity and specificity-1) was used to calculate the cut-off value of the optimal diagnostic point of the ROC curve. Binary logistic regression and ROC curves were used to improve diagnostic efficiency. All statistical analyses were performed in SPSS (Version 23, SPSS Inc.) and R software (version 3.5.0). A value of p < 0.05 was considered statistically significant.

Patient data analysis and seminal plasma piRNA profiling
Seven semen samples were collected from AZS and male normal healthy controls according to the WHO guidelines. The basic information and sperm parameters of AZS and healthy controls are presented in Table 1. There were no distinct differences in age, BMI, smoking or drinking status between the AZS and control groups (P>0.05). However, there were obvious differences in motility and sperm morphology (P <0.001) between the AZS and control groups. The semen samples were observed under a light microscope equipped with a 100 x oil objective and were found to be free of seminiferous epithelium cells, somatic cells and white blood cells.
Small RNAs were extracted from semen samples of healthy controls and AZS patients and analysed by Illumina high-throughput sequencing. First, the splice sequences were eliminated from the raw reads and small RNA sequences of different lengths were obtained. After quality control steps, small RNAs in the range of 10 to 45 nucleotides were obtained. Analysis of the length distribution showed that all three seminal plasma samples contained many small RNAs with sizes consistent with miRNAs (18-24 nucleotides) and piRNAs (26-31 nucleotides) (Figure 1). In general, piRNAs account for approximately 0.05% of the total small RNA content and 0.5% of the total small RNA reads in seminal plasma. Among the 77 million piRNAs annotated in the piRBase database, 30,348 and 53,932 piRNAs were detected in the seminal plasma of healthy controls and AZS patients, respectively. 4.2 Identification of differential piRNA expression We used DEG-seq software to analyse the differences in piRNA expression levels between the AZS and normal groups. A total of 283 DEGs were obtained with a significance threshold of |log2FC|≥1 and P-value ≤ 0.05, comprising 114 upregulated DEGs and 169 downregulated DEGs. We displayed volcano plots and heatmaps of DEGs in Figure 2A, B. Twenty differentially expressed piRNAs, including the top 10 most strongly upregulated and top 10 most strongly downregulated piRNAs, are presented in Figure 2C, D and Table 2. 4.3 GO and KEGG pathway analysis GO enrichment analysis was used to explore the potential functions of these DEGs. The top 30 most frequent GO terms are given in Figure 3. As shown in Figure 3A, B, and C, the top cellular component terms were nucleus, cytosol, cytoplasm and nucleoplasm; the top biological process terms were regulation of transcription-DNA templated, signal transduction and positive regulation of transcription from the RNA polymerase II promoter; and the top molecular function terms were metal ion binding, ATP binding, DNA binding and RNA binding.
KEGG pathway analysis was used to identify all pathways and the top 20 most highly enriched pathways related to the potential target genes of the differentially expressed piRNAs ( Figure 3C, D ,E). The target genes enriched among the AZS DEGs were mainly involved in focal adhesion, actin cytoskeleton, axon guidance, ErbB signalling pathway and other signal transduction pathways.

Validation of differentially expressed piRNAs with qRT-PCR
To validate the small RNA-seq results, qRT-PCR analysis was performed on 8 randomly selected differentially expressed piRNAs. We selected 4 upregulated piRNAs (piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, piR-hsa-18586) and 4 downregulated piRNAs (piR-hsa-6148, piR-hsa-32713, piR-hsa-2912, piR-hsa-15551) and performed qRT-PCR to quantify their expression in AZS and normal group semen samples. As shown in Figure 4, the expression of piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, and piR-hsa-18586 were upregulated in the AZS group compared with the normal group, and piR-hsa-2912 indicated a downregulated expression pattern, which was consistent with the small RNA-seq results. However, the other three genes did not show significant differences in expression between the two samples.

Evaluate piRNAs as a diagnostic marker for AZS.
Illumina high-throughput sequencing analysis indicated that piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, and piR-hsa-18586 expression levels were obviously different between AZS and controls. Therefore, the diagnostic value these four piRNAs as potential biomarkers for AZS was assessed and evaluated. We used ROC curve analysis to calculate whether piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, and piR-hsa-18586 could distinguish AZS patients from controls. As shown in Figure 5, the AUCs of piR-has-18586, piR-has-18725, piR-has-26591, and piR-has-32694 were 0.850 (P<0.007), 0.875 (P<0.004), 0.913 (P<0.001), and 0.883 (P<0.003), respectively. To improve the diagnostic power, we combined these four piRNAs in 80 patients with AZS and 100 normal controls as a training cohort to build a diagnostic model using binary logistic regression. An AUC of 0.935 was found ( Figure 6) when the model was applied to an independent testing cohort (40 controls and 40 AZS patients) to estimate the diagnostic accuracy of the joint diagnosis model. These results showed that the combined diagnostic model can distinguish AZS patients from controls (specificity = 95%, sensitivity = 90%).

Discussion
In this study, we identified the comprehensive piRNA expression pattern in semen samples of AZS patients. The small RNA sequencing data revealed that a total of 283 significant DEGs (114 upregulated and 169 downregulated) were obtained at a threshold of |log2FC|≥1 and P-value ≤ 0.05. Further in silico analyses, including piRNA/mRNA interaction network, GO and KEGG pathway analyses, were performed to predict the functions of differentially expressed piRNAs manifesting the critical roles of piRNAs in regulating spermatogenesis. These results were validated by qRT-PCR assays of 8 randomly selected piRNAs, specifically, 4 upregulated piRNAs (piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, piR-hsa-18586) and 4 downregulated piRNAs (piR-hsa-6148, piR-hsa-32713, piR-hsa-2912, piR-hsa-15551). Of these, piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, piR-hsa-18586 and piR-hsa-2912 produced qRT-PCR results consistent with the sequencing results, whereas the other three genes were did not show significant differences in expression between AZS samples and controls. The piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, and piR-hsa-18586 piRNAs were significantly upregulated in AZS. The diagnostic power of these four piRNAs was further analysed using ROC analysis, and piR-hsa-26591 was found to be the most diagnostic for AZS, with an AUC of 0.913 (95% CI: 0.795-0.994). A logistic regression model including all four piRNAs was constructed, and ROC analysis revealed that the combination model reached good diagnostic efficacy (AUC: 0.935). Therefore, these piRNAs were identified as hub genes that participate in AZS pathogenesis and could serve as biomarkers for this disease. Human spermatogenesis is a carefully coordinated and precisely regulated biological process that includes rigorous regulation of gene expression at specific phases 3 . Noncoding RNAs such as circular RNAs (circRNAs), long noncoding RNAs (lncRNAs) and microRNAs (miRNAs) are critical regulators of gene expression in multiple phases of spermatogenesis [16][17][18] . piRNAs are a recently discovered class of small RNAs that are mainly expressed in germ cells. piRNAs play important roles in maintaining the integrity of germline DNA, inhibiting transposon transcription, inhibiting translation, participating in the formation of heterochromatin, epigenetic regulation and germ cell genesis 19 . To date, there have been few reports describing the changes in piRNA profile in human AZS. In view of the importance of piRNAs in male reproduction, we studied the expression pattern of piRNAs in AZS sperm, analysing the differential expression of piRNAs in AZS and controls to explore the potential pathological mechanism of AZS.
Studies have shown that piRNAs regulate gene silencing mainly by binding with piwi subfamily proteins to form piRNA complexes (piRCs). The human genome encodes four piwi protein family members, Hiwi, Hili, Hiwi2 and PIWIL3, which are highly expressed in testicular cells 20 . Liu et al. showed that a Hiwi mutation is pathogenic, causing male infertility, and established a clear relationship between the human piwi gene and male infertility 21 . Gou and other research groups have demonstrated that in addition to its main functions in transposon silencing, the piwi/piRNA mechanism mediates the degradation of a large number of mRNA transcripts in various germ cells through miRNA-or/and siRNA-like mechanisms 22,23 . In this study, we predicted the target genes of piRNAs that were differentially expressed between the control group and the AZS group to explore how the differentially expressed piRNAs might change sperm function. In addition, enrichment analysis based on these predicted target genes showed that cell macromolecule biosynthesis, spermatogenesis, male gametogenesis, cytoplasmic transport and RNA metabolism were changed in the AZS group. KEGG pathway analysis showed that the differentially expressed piRNA targets were involved in focal adhesion kinase (FAK), actin cytoskeleton, axon guidance, the ErbB signalling pathway and other signal transduction pathways. An interesting study showed that FAK coordinates the transport of spermatozoa and prealbumin spermatids through the epithelium along the functional axis of the acroplasmic thinning (ES)-blood testicular barrier basement membrane (BM) 24 . A study on pig sperm found that many differentially expressed mRNAs related to infertility are involved in metabolic processes 25 . Specifically, during the sperm deformation stage, macromolecule metabolism and transportation are well controlled, and the production of fertile sperm is ensured through 26S proteasome exchange, digestion and absorption of nuclear protein, and elimination of RNA through phagocytosis of cytoplasmic droplets, residues and Sertoli cells 26 . Our results also support the view that focal adhesion kinase and metabolism of fatty acids and other macromolecules may play key roles in spermatogenesis and in AZS. Previous reports identified mRNAs that were related to sperm motility, particularly transcripts encoding flagella-associated protein 45 (CFAP45) and bromodomain-containing 2 (BRD2). The content of CFAP45 in low-motility sperm was significantly higher than that in high-motility sperm, while the content of BRD2 mRNA in AZS was significantly lower 27,28 . This indicates that low sperm motility could be attributed to the abnormal expression of CFAP45 and BRD2. Here, we found that the expression of four piRNAs, piR-hsa-32694, piR-hsa-26591, piR-hsa-18725 and piR-hsa-18586, was significantly upregulated in AZS patients. Furthermore, the ROC curve analysis revealed that the AUC of the combination of the 4 piRNAs was equal to 0.935 and demonstrated 93.8% sensitivity and 89.3% specificity. These results suggested that these four piRNAs were significantly and independently associated with sperm motility and male infertility. Bioinformatics analysis showed that CAFP45 and BRD2 were the potential target genes of piR-hsa-26591. Whether pir-hsa-26591 silences the expression of CAFP45 and BRD2 and leads to decreased sperm motility and infertility in AZS patients needs further study. Overall, this study may provide insights into the use of human sperm cell-specific piRNAs as biomarkers for AZS screening. Our study has some limitations. First, the small number of samples limited the statistical power of the microarray analysis. Second, we did not use a luciferase reporter assay to further verify the relationship between pir-hsa-26591 and CAFP45.
In summary, we have defined a distinctive seminal plasma piRNA signature in AZS. Moreover, our findings may provide additional information regarding the molecular mechanisms of male infertility.