*Contributed equally
Asthenozoospermia (AZS) is characterized by reduced sperm motility and its pathogenesis remains poorly understood. Piwi-interacting RNAs (piRNAs) have been indicated to serve important roles in spermatogenesis. However, little is known about the correlation of piRNA expression with AZS. In the present study, small RNA sequencing (small RNA-seq) was performed on sperm samples from AZS patients and fertile controls. Reverse transcription-quantitative (RT-q) PCR was used to validate the small RNA-seq results. Bioinformatics analyses were performed to predict the functions of differentially expressed piRNAs (DEpiRNAs). 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. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses showed that the DEpiRNAs were mainly associated with transcription, signal transduction, cell differentiation, metal ion binding and focal adhesion. These results were verified by RT-qPCR analysis of eight selected piRNAs. The PCR results were consistent with the sequencing results in patients with AZS compared with controls in the first cohort. The expression of piR-hsa-32694, piR-hsa-26591, piR-hsa-18725 and piR-hsa-18586 was significantly upregulated in patients with AZS. The diagnostic power of the four piRNAs was further analyzed 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 achieved good diagnostic efficacy (AUC: 0.935).
According to a report from the World Health Organization (WHO), 25% of couples in developing countries experience infertility. Infertility often causes great psychological strain, which can eventually lead to mental illness and ~50% of cases of infertility are caused by male factors (
piRNAs are a unique class of small non-coding (nc)RNAs that are 26-30 nucleotides in length and guide Piwi proteins (
Therefore, the present study aimed to determine the expression patterns and functions of piRNAs in AZS patients. Bioinformatics analysis was conducted to elucidate the potential biological functions of piRNAs. These results may provide new biomarkers for the development of novel diagnostic and therapeutic strategies for AZS.
The protocol of the present study was approved by the Medical Ethics Committee of Nanchang University, the Second Affiliated Hospital and Jiangxi Maternal and Child Health Hospital (China; approval no. 2018089). All donors provided written informed consent for the collection and use of their samples for this protocol.
The present study selected 194 patients who were diagnosed with male sterility due to AZS and 143 normal healthy males who visited the Maternal and Child Health Hospital in Jiangxi Province between January 2019 and December 2020. A total of four patient specimens and three normal specimens were randomly selected for RNA expression analysis and 250 samples (150 AZS patient and 100 normal individuals) and 80 samples (40 AZS patient and 40 normal individuals) were randomly selected for the clinical verification experiment by reverse transcription-quantitative (RT-q) PCR assays. Semen samples (6-8 ml) were collected in a dry and sealed aseptic container. Sperm morphology, concentration, motility and quantity was used to evaluate sperm quality. According to the fifth edition of the AZS diagnostic criteria in 2010(
Total RNA was extracted from semen samples by TRIzol® (Thermo Fisher Scientific, Inc.; cat. no. 15596-018) reagent according to the manufacturer's specifications, the yield was determined by a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Inc.) and the integrity was assessed by agarose gel electrophoresis with ethidium bromide staining.
Quantification was performed with a two-step reaction process: Reverse transcription (RT) and PCR. Each RT reaction consisted of 0.5 µg RNA, 5 µl 2XTS miRNA Reaction Mix and 0.5 µl of TransScrip miRNA RT Enzyme Mix (TransGen Biotech, AT351), in a total volume of 10 µl. Reactions were performed in a GeneAmp PCR System 9700 (Applied Biosystems; Thermo Fisher Scientific, Inc.) for 60 min at 37˚C, followed by heat inactivation of RT for 5 sec at 85˚C. The 10 µl RT reaction mix was then diluted x10 in nuclease-free water and held at -20˚C.
RT-PCR was performed with a LightCycler 480 II Real-time PCR Instrument (Roche Diagnostics GmbH) in a 10-µl reaction system that included 1 µl of cDNA, 5 µl of 2x PerfectStart Green qPCR SuperMix (TransGen Biotech, AQ601), 0.2 µl of 10 µM universal primers, 0.2 µl of 10 µM miRNA-specific primer and 3.6 µl of nuclease-free water. The reactions were incubated in a 96-well plate (Beijing Jiaxinheng Biotechnology Co., Ltd.) at 94˚C for 30 sec, followed by 40 cycles of 95˚C for 10 sec and 65˚C for 35 sec. Each sample was run in triplicate. At the end of the reaction cycle, melting curve analysis was performed to detect the product specificity. The tailing method is used in the reverse transcription of piRNA, so the quantitative primer only needs to design the F primer and the R primer is the general primer sequence of the kit (TransGen Biotech Co., Ltd.; cat. no. AT351-01), GATCGCCCTTCTACGTCGTAT (TM=58˚C). The piRNA-specific primer sequences were synthesized by Beijing Qingke Xinye Biotechnology Co., Ltd. based on the piRNA sequences obtained from the piRNABase database (
Briefly, small RNA sequencing and analysis were conducted by Shanghai Oe Biotech. Co., Ltd. Data from Illumina HiSeq sequencing are called raw reads. First, the splice sequence was removed from the raw reads and then small RNA sequences of different lengths were obtained. Cutadapt (version 1.14) (
Differentially expressed genes (DEGs) between the AZS patients and the fertile controls were identified using DEG-seq software (version 1.18.0). Genes with q-values (FDRs) <0.05 and fold-changes >2.0 were considered differentially expressed (
Differences in demographic variables and expression levels of piRNA in the AZS group and the fertile control group were evaluated via the χ2 test. Fisher's exact test and hypergeometric distribution test were performed to assess the enrichment significance of GO items or KEGG pathways in the differentially expressed genes.
To evaluate the specificity and sensitivity of each piRNA in the diagnostic value of AZS, the ROC curves and AUC were analyzed. A high 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 the diagnostic efficiency. All the statistical analyses were performed in SPSS (Version 23, SPSS Inc.) and R software (v 3.5.0;
A total of seven semen samples were collected from AZS patients and male healthy controls according to the WHO guidelines. The basic information and sperm parameters of the AZS patients and healthy controls are presented in
Small RNAs were extracted from the semen specimens of healthy controls and AZS patients and analyzed 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 (nt) were obtained. Length distribution analysis showed that all the samples contained many RNAs shorter than 45 nt, which was consistent with the lengths of miRNAs (18-224 nt) and piRNAs (26-330 nt;
DEG-seq software was used to analyze the differences in piRNA expression levels between the AZS and normal groups. A total of 283 differentially expressed piRNAs were obtained with a significance threshold of |log2FC|≥1 and P-value ≤0.05, including 114 upregulated DEGs and 169 downregulated DEGs. Volcano plots and heatmaps of the DEGs are shown in
GO enrichment analysis was used to explore the potential functions of these DEpiRNA-target genes. The top 30 most frequent GO terms are shown in
KEGG pathway analysis was used to identify all the pathways and the top 20 most highly enriched pathways related to the DEpiRNA-target genes (
To validate the small RNA-seq results, RT-qPCR analysis was performed on 8 piRNAs among the DEGs. The present study selected four upregulated piRNAs (hsa-32694, hsa-26591, hsa-18725 and hsa-18586) and four downregulated piRNAs (hsa-6148, hsa-32713, hsa-2912 and hsa-15551) to quantify their expression in the AZS and normal group semen samples. The primer sequences used for investigating piRNA expression are listed in
Illumina high-throughput sequencing analysis indicated that the piR-hsa-32694, piR-hsa-26591, piR-hsa-18725 and piR-hsa-18586 expression levels were obviously different between AZS patients and controls. Therefore, the diagnostic value of these four piRNAs as potential biomarkers for AZS was assessed and evaluated. ROC curve analysis was used 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
The present study identified the comprehensive piRNA expression pattern in the semen samples of AZS patients. The small RNA sequencing data revealed that a total of 283 significant DEpiRNAs (114 upregulated and 169 downregulated) were obtained at a threshold of |log2FC|≥1 and P≤0.05. Further analyses, including DEpiRNA analysis and GO and KEGG pathway analyses, were performed to predict the functions of DEpiRNAs, which revealed the critical roles of piRNAs in regulating spermatogenesis. These sequencing results were validated by RT-qPCR analysis of eight piRNAs, specifically, four upregulated piRNAs (piR-hsa-32694, piR-hsa-26591, piR-hsa-18725 and piR-hsa-18586) and four downregulated piRNAs (piR-hsa-6148, piR-hsa-32713, piR-hsa-2912 and piR-hsa-15551). Of these, piR-hsa-32694, piR-hsa-26591, piR-hsa-18725, piR-hsa-18586 and piR-hsa-2912 produced RT-qPCR results consistent with the sequencing results, whereas piR-hsa-6148, piR-hsa-32713 and piR-hsa-15551 did not show significant differences in expression between the AZS and control samples. The expression of piR-hsa-32694, piR-hsa-26591, piR-hsa-18725 and piR-hsa-18586 was significantly upregulated in patients with AZS. The diagnostic power of these four piRNAs was further analyzed using ROC analysis and piR-hsa-26591 was found to be the piRNA with the most potential for diagnosing 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 achieved 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 (
Studies have shown that piRNAs affect gene silencing mainly by binding to Piwi subfamily proteins of the Ago/Piwi protein family. Human piwi family proteins, PIWIL1, PIWIL2, PIWIL3 and PIWIL4 are expressed in testis tissues (
Previous reports identified mRNAs that are related to sperm motility, particularly transcripts encoding flagella-associated protein 45 (CFAP45) and bromodomain-containing 2 (BRD2). The content of CFAP45 in low motile power sperm is significantly higher compared with that in high motile power sperm, while the content of BRD2 mRNA in sperm from AZS patients is significantly lower (
In summary, the present study identified a piRNA expression profile and defined the potential function of DEpiRNAs in AZS. Moreover, the findings also provide some information regarding the pathological mechanisms of male infertility.
Not applicable.
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
All the authors were involved in the conception and design of the study and data interpretation. LH and HC performed the experiments, XW and RW analyzed the data and PG, WH and WS helped perform the experiments. LH prepared the manuscript. HC oversaw the project and proofread the manuscript. LH and HC confirm the authenticity of all the raw data. All authors read and approved the final manuscript.
The protocol of the present study was approved by the Medical Ethics Committee of Nanchang University, the Second Affiliated Hospital and Jiangxi Maternal and Child Health Hospital (China; approval no. 2018089). All donors provided informed consent for the collection and use of their samples for this protocol.
Not applicable.
The authors declare that they have no competing interests.
Small RNA length distribution map. RNAs 18-24 nucleotides in length are miRNAs and RNAs 26-31 nucleotides in length are piRNAs. Patient: asthenozoospermia vs. control.
DEpiRNAs between patients with AZS and normal controls. (A and B) Visual images of DEGs between the AZS and normal groups. (C and D) Visual images of the top 20 DEGs between the AZS and fertile control groups. The green and red dots indicate the upregulated and downregulated genes, respectively (threshold of |log2fold change |≥1.0 and P≤0.05). The grey dots indicate that there were no differentially expressed genes. Rows and columns are used to represent DEGs and samples. Case: AZS vs. normal. Different degrees of color represent the expression level of DEGs. DEpiRNAs, differentially expressed piwi-interacting RNAs; AZS, asthenozoospermia; DEGs, differentially expressed genes.
Enrichment analyses of DEpiRNA. (A-C) Top 30 GO enriched terms among all differentially expressed genes and among upregulated and downregulated genes. (D-F) Top 20 KEGG pathways enriched among the total differentially expressed, upregulated and downregulated genes. DEpiRNAs, differentially expressed piwi-interacting RNAs; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.
Expression levels of piRNAs in the semen from AZS patients and controls. (A) piR-has-32694, (B) piR-has-26591, (C) piR-has-18725, (D) piR-has-18586, (E) piR-has-6148, (F) piR-has-3217, (G) piR-has-15551 and (H) piR-has-2912. ***P<0.01 control vs. AZS). piRNA/piR, piwi-interacting RNA; AZS, asthenozoospermia.
Diagnostic power of the four piRNA potential biomarkers, as determined using receiver operating characteristic curve analysis. (A) piR-has-18586, (B) piR-has-18725, (C) piR-has-26591, (D) piR-has-32694. [The optimal diagnostic point was assessed at the cut-off values with the highest Youden's indices (sensitivity and specificity-1)] piRNA/piR, piwi-interacting RNA; AUC, area under curve.
Diagnostic model of AZS using four piRNAs (expression of piR-has-18586, piR-has-18725, piR-has-26591 and piR-has-32694 used as a training set). piRNA/piR, piwi-interacting RNA; AUC, area under curve.
Basic information and parameters of AZS and fertile controls.
Characteristic | AZS (194.0) | Fertile controls (143) | P-value |
---|---|---|---|
Age (years) | 34.00±5.3 | 33.5±5.32 | 0.927 |
BMI | 23.50±3.69 | 24.6±2.79 | 0.923 |
Smoking | 0.87 | ||
Yes | 146 (75.00%) | 107 (75%) | |
No | 48 (25%) | 36 (25%) | |
Drinking | 0.96 | ||
Yes | 139 (71%) | 1033 (71%) | |
No | 55 (29%) | 40.00 (29%) | |
Sperm volume (ml) | 4±0.68 | 3±0.63 | 0.87 |
pH | 7.6±0.15 | 7.7±0.18 | 0.56 |
Sperm density (x106/ml) | 58.59±6.34 | 34.12±5.43 | <0.001 |
PR | 16.18±3.2 | 62.55±3.8 | <0.001 |
Sperm activity rate | 22.12±5.64 | 77.34±3.2 | <0.001 |
AZS, asthenozoospermia; PR, progressive motility.
aDifference was statistically significant (P<0.05).
Differentially expressed genes between AZS and controls (top 10).
Upregulated | Downregulated | ||||
---|---|---|---|---|---|
Gene | logFC | P-value | Gene | logFC | P-value |
piR-26591 | 8.490985 | 0.004693 | piR-15551 | -6.52024 | 0.009999 |
piR-32694 | 7.712583 | 0.000754 | piR-29529 | -6.49304 | 0.009822 |
piR-18725 | 7.221809 | 0.009867 | piR-25106 | -6.17803 | 0.006124 |
piR-16782 | 6.333706 | 0.006759 | piR-2645 | -5.63236 | 0.006583 |
piR-25491 | 6.224344 | 0.004153 | piR-9170 | -5.57597 | 0.002927 |
piR-5952 | 5.67597 | 0.006451 | piR-2912 | -5.32921 | 0.003 |
piR-18611 | 5.669888 | 0.008762 | piR-23203 | -4.82064 | 0.005891 |
piR-18586 | 5.527972 | 0.009949 | piR-6148 | -4.6561 | 0.007541 |
piR-15811 | 5.406727 | 0.007682 | piR-32713 | -4.10666 | 0.006346 |
piR-17725 | 4.956704 | 0.004188 | piR-29155 | -4.07376 | 0.009545 |
114 | 169 |
AZS, asthenozoospermia; piR, piwi-interacting RNA.
Top 30 GO terms of DEpiRNA-target genes.
Biological process | Cellular component | Molecular function | |
---|---|---|---|
Top 30 GO terms of total DEpiRNA-target genes | Transcription, DNA-templated | Nucleus | Metal ion binding |
Regulation of transcription, DNA-templated | Cytosol | ATP binding | |
Signal transduction | Cytoplasm | DNA binding | |
Positive regulation of transcription from RNA polymerase II promoter | Plasma membrane | RNA binding | |
Negative regulation of transcription from RNA polymerase II promoter | Nucleoplasm | Identical protein binding | |
Intracellular signal transduction | Integral component of membrane | DNA binding transcription factor activity | |
Multicellular organism development | Extracellular exosome | Zinc ion binding | |
Cell differentiation | Membrane | Calcium ion binding | |
Positive regulation of transcription, DNA-templated | Integral component of plasma membrane | ||
Protein homodimerization activity | |||
Negative regulation of transcription, DNA-templated | Mitochondrion | Protein serine/threonine kinase activity | |
Top 30 GO terms of upregulated | Transcription, DNA-templated | Nucleus | Metal ion binding |
Regulation of transcription, DNA-templated | Cytosol | ||
DEpiRNA-target genes | ATP binding | ||
Positive regulation of transcription from RNA polymerase II promoter | Cytoplasm | ||
DNA binding | |||
Signal transduction | Plasma membrane | RNA binding | |
Multicellular organism development | Integral component of membrane | Identical protein binding | |
Negative regulation of transcription from RNA polymerase II promoter | Nucleoplasm | Zinc ion binding | |
Intracellular signal transduction | Extracellular exosome | DNA binding transcription factor activity | |
G-protein coupled receptor signaling pathway | Membrane | Calcium ion binding | |
Cell adhesion | Integral component of plasma membrane | Actin binding | |
Cell differentiation | Extracellular region | Protein homodimerization activity | |
The top 30 GO terms of downregulated | Transcription, DNA-templated | Cytosol | Metal ion binding |
Regulation of transcription, DNA-templated | Nucleus | RNA binding | |
DEpiRNA-target gene | Positive regulation of transcription from RNA polymerase II promoter | Cytoplasm | ATP binding |
G-protein coupled receptor signaling pathway | Plasma membrane | DNA binding | |
Intracellular protein transport | Nucleoplasm | Identical protein binding | |
Intracellular signal transduction | Integral component of membrane | ||
Calcium ion binding | |||
Cell adhesion | Extracellular exosome | Ion channel binding | |
Signal transduction | Membrane | DNA binding transcription factor activity | |
Multicellular organism development | Golgi apparatus | Rab GTPase binding | |
Negative regulation of transcription from RNA polymerase II promoter | Cell junction | Protein serine;threonine kinase activity |
GO, Gene Ontology; DEpiRNAs, differentially expressed piwi-interacting RNA.
Top 20 KEGG pathways of DEpiRNA-target genes (upregulated and downregulated).
Gene item | Top 20 KEGG pathways |
---|---|
Total DEpiRNA-target genes | MicroRNAs in cancer |
Proteoglycans in cancer | |
Bacterial invasion of epithelial cells | |
Regulation of actin cytoskeleton | |
Signaling pathways regulating pluripotency of stem cells | |
Tight junction | |
Adherens junction | |
Cell adhesion molecules (CAMs) | |
Focal adhesion | |
Axon guidance | |
Sulfur relay system | |
ErbB signaling pathway | |
Sulfur metabolism | |
Glycosphingolipid biosynthesis-ganglio series | |
Inositol phosphate metabolism | |
Glycosaminoglycan biosynthesis-keratan sulfate | |
Neomycin, kanamycin and gentamicin biosynthesis | |
Other types of O-glycan biosynthesis | |
D-Arginine and D-ornithine metabolism | |
Valine, leucine and isoleucine biosynthesis | |
Upregulated DEpiRNA-target genes | MicroRNAs in cancer |
Protein digestion and absorption | |
Cushing syndrome | |
Aldosterone synthesis and secretion | |
Protein digestion and absorption | |
Insulin signaling pathway | |
Cholinergic synapse | |
Neurotrophin signaling pathway | |
Circadian entrainment | |
ECM-receptor interaction | |
Focal adhesion | |
Axon guidance | |
Phospholipase D signaling pathway | |
cAMP signaling pathway | |
Calcium signaling pathway | |
Ras signaling pathway | |
ErbB signaling pathway | |
Glycerophospholipid metabolism | |
Acute myeloid leukemia | |
Glioma | |
Downregulated DEpiRNA-target genes | Central carbon metabolism in cancer |
Proteoglycans in cancer | |
Salmonella infection | |
Mineral absorption | |
Protein digestion and absorption | |
Carbohydrate digestion and absorption | |
Endocrine and other factor-regulated calcium reabsorption | |
Cortisol synthesis and secretion | |
Thyroid hormone signaling pathway | |
Synaptic vesicle cycle | |
Circadian rhythm | |
Tight junction | |
Hippo signaling pathway | |
Ubiquitin mediated proteolysis | |
Calcium signaling pathway | |
MAPK signaling pathway | |
EGFR tyrosine kinase inhibitor resistance | |
Nicotinate and nicotinamide metabolism | |
Valine, leucine and isoleucine biosynthesis | |
Purine metabolism |
KEGG, Kyoto Encyclopedia of Genes and Genomes; DEpiRNAs, differentially expressed piwi-interacting RNA.
Primers sequences of the piRNA.
Gene symbol (piRNA) | Primer 5'→3' |
---|---|
piR-hsa-26591 | GACGAGGTGGCCGAGTGGTT |
piR-hsa-32694 | GAACAAGAAGACACTCGTGG |
piR-hsa-18725 | GACACTCGTGGAGGCGTCA |
piR-hsa-18586 | AGACACTCGTGGAGACGTC |
piR-hsa-32713 | CCTAGAAGACTGACTCTTTGC |
piR-hsa-6148 | TAGTGGTTATCACTTTCGCCT |
piR-hsa-2912 | TCCAACTGCTCTACGACACA |
piR-hsa-15551 | TGCGACGAGATGCTGCTTC |
hsa-5S | GGAGACCGCCTGGGAATA |
piR, piwi-interacting RNA.