RNA-sequencing and bioinformatics analysis of long noncoding RNAs and mRNAs in the asthenozoospermia

Infertility is a worldwide health problem, affecting approximately 15% of couples with child-bearing age, and male factor has been deemed to 50% causes for infertility [1]. The main cause of male infertility is sperm quality loss displaying azoospermia, oligozoospermia, asthenozoospermia and teratospermia. In these cases, asthenozoospermia is one of the major causes of human male infertility, characterized by reduced forward motility of spermatozoa (grade A + B sperm motility < 50% or A < 25%), which prevents the sperm from moving to the egg and penetrating it, eventually leading to infertility [2]. However, the cause and pathogenesis of asthenozoospermia are not completely understood.

Exosomes are nanosized extracellular vesicles (40–180 nm) that can be released by almost all types of cells and exist stably in various biological fluids (e.g. blood, urine, saliva, seminal and follicular) [3]. The importance of exosomes lies in their role in intercellular communication by transmitting bioactive molecules, including proteins, DNA, mRNAs and non-coding RNAs. Exosomes could transfer information to recipient cells and thereby influencing the cell functions [4]. These unique properties make exosomes a potential candidate as reliable biomarkers. The recent studies have shown that exosomes are secreted along the male reproductive tract and are thought to be involved in spermatozoa maturation and function [5,6]. Murdica et al. [7] have reported that seminal plasma of men with severe asthenozoospermia contains exosomes that affect spermatozoa motility, and could promptsperm capacitation through an increased induction of tyrosine phosphorylation, therefore, inducing the acrosome reaction. Meanwhile, Murdica et al. [8] revealed the negative modulator of sperm function glycodelin as over-represented in semen exosomes isolated from asthenozoospermic patients by proteomics analysis.

Long noncoding RNAs (lncRNAs) are a group of nonprotein-coding RNAs with a length longer than 200 nucleotides that distributed in the genome broadly [9]. Accumulating evidence suggests that lncRNAs regulate gene expression in the form of RNA in epigenetic regulation, regulation of transcription and post-transcriptional regulation levels [10]. The regulatory and structural functions of lncRNAs have been reported to play pivotal roles in a broad scope of biological processes, and disorderly regulation of lncRNAs leads to diverse human diseases [11]. Recent discovery of lncRNAs as critical regulators in normal and disease development provides new clues for delineating the molecular regulation in male germ-cell development [12]. Necsulea et al. [13] performed a large-scale evolutionary study of lncRNA repertoires and expression patterns, in 11 tetrapod species and found the potential functions for lncRNAs in spermatogenesis processes. Wen et al. [14] revealed that rapidly evolving testis-specific lncRNAs play critical roles in late Drosophila spermatogenesis. Wichman et al. [15] found that dysregulation of specific lncRNAs led to sperm counts or infertility reduced in mice. Recently, Zhang et al. [16] explored the expression profiles and characteristics of human lncRNAs in normal and asthenozoospermia sperm indicating an association between lncRNAs expression and sperm motility. However, the expression and characteristics of lncRNAs in seminal plasma exosomes of normal and asthenozoospermia have not been addressed.

Therefore, to investigate the intricate regulatory network associated with asthenozoospermia, we performed RNA-sequencing and analyzed the differentially expressed genes (DEGs), differentially expressed lncRNAs (DElncRNAs) and differentially expressed messenger RNAs (DEmRNAs) and new lncRNAs between normal and asthenozoospermia groups with bioinformatics in the present study. Subsequently, we predicted the biological function of the DEGs using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and DO (Disease Ontology) enrichment analyses. In addition, we constructed co-expression networks of the lncRNA–miRNA–mRNA with specific bioinformatics approaches to search the hub lncRNAs in the development of asthenozoospermia. To confirm sequencing results with five of the identified DElncRNAs, we performed quantitative reverse-transcription polymerase chain reaction (qRT-PCR) between normal (n=20) and asthenozoospermia (n=20) groups.

We selected 25 male sterility diagnosed with asthenozoospermia and 25 male normal healthy individuals at the First Affiliated Hospital of Hainan Medical College. The semen samples (8 ml) were collected by masturbation into a sterile container after 3–7 days of abstinence. We used a sperm quality analyzer to assess the semen basic parameters, including morphology, concentration, motility, viability and counts of spermatozoa. According to the World Health Organization (WHO, 2010, Fifth Edition) standards on the diagnostic criteria of asthenozoospermia, the rapid forward progressive motile sperm (grade A <25%) and total progressive motile sperm (grades A + B <50%) in fresh ejaculation [17]. The semen samples were considered normal with the following parameters: sperm concentration ≥ 15 × 10⁶/ml, semen volumes ≥ 1.5 ml, pH ≥ 7.2, progressive motility (PR) ≥ 32%, PR + Non-progressive motility (NP) ≥ 40%. The age range of asthenozoospermia patients and normal healthy controls was 24–48 years old. Control subjects with underlying diseases or abnormal sperm parameters were excluded from the study. All the participants are ethnic Han Chinese with no genetic relationship. Patients with congenital or hereditary disease that cause gonadal dysplasia were excluded from the study. Patients with genitourinary infections, endocrine diseases such as diabetes and hyperthyroidism, severe heart, liver, kidney or other functional abnormalities were also excluded.

The semen samples were allowed to liquefy for 30 min at 37°C, 5% CO₂ incubator. Then, semen samples were centrifuged at 1000 × g for 10 min to separate seminal plasma containing exosomes. Subsequently, the each seminal plasma sample was transferred to a 15 ml centrifuge tube, and was centrifuged at 300 × g for 10 min, 2000 × g for 20 min at 4°C, respectively. We discard the precipitate and remove the cells. Then, supernatants were in turn centrifuged at 10,000 × g for 30 min. We discard the precipitate and remove subcellular components. The remaining supernatants were centrifuged at 10,000 × g for 60 min. We discard the supernatant, and the resulting precipitate is the seminal plasma exosomes. Exosome-containing pellets were re-suspended with 30 ml PBS, followed by centrifugation at 10,000 × g for 60min. The Precipitator-purified exosomes were suspended with 1 ml PBS. Then, samples were loaded separately into the eppendorf tube and stored in a −80°C refrigerator for future use. The morphological characteristics were observed by transmission electron microscopy (H-7650, Hitachi Limited, Japan) according to the manufacturer’s instructions.

Western blot was used to detect exosome-specific antigen molecules CD63 and CD81. We used the reducing Laemmli buffer to dissolve the seminal plasma exosomes (5 μg) and boiled for 5 min at 95°C. Proteins were resolved in a 10% sodium dodecyl sulfate-polyacrylamide gel (SDS-PAGE), followed by being transferred to a PVDF membrane. The membranes were blocked in 5% skimmed milk in PBS containing 0.5% Tween-20 at room temperature for 1 h prior to being probed with the anti-CD63 (1:20,000, #556019; BD Pharmingen), anti-CD81 (1:5,000, #555675; BD Pharmingen) at 4°C overnight. After washing with T-TBS, membranes followed by being incubated with the appropriate horseradish peroxidase-labeled secondary antibody (goat anti-mouse (1: 10,000) or anti-rabbit (1: 2000) immunoglobulin G) for 45 min. The enhanced chemiluminescence (ECL) reagent (Advansta, Menlo Park, CA, U.S.A.) was used to detect the positive immunoreactive bands.

LncRNA-sequencing was performed by the Bioacme Biotechnology Co., Ltd. (Wuhan, China). Briefly, total RNA was extracted with Trizol reagent (Invitrogen, Carlsbad, CA, U.S.A.) from the exosomes of asthenozoospermia patients and normal healthy controls. Total RNA was quantified using a spectrophotometer (NanoDrop 2000; Thermo Fisher Scientific, Waltham, MA, U.S.A.) according to manufacturers’ instructions. We used the Agilent 2100 system and an RNA Nano 6000 Assay kit (Agilent Technologies, Santa Clara, CA, U.S.A.) to assess RNA integrity. The ribosomal RNA (rRNA) was removed using Ribo-Zero rRNA removal beads (Illumina, Inc., San Diego, CA, U.S.A.), fragmentation (the average segment length is approximately 200 nt), to enable accurate lncRNA analysis. The extracted total RNA was reverse transcribed into single-stranded complementary DNA (cDNA), then synthesized into a double-stranded cDNA, terminal repair, 3′ terminal addition, ligation, and addition of primers, PCR amplification and purification, and quality inspection of the library. The generated libraries were sequenced on an Illumina HiSeq 3000 (Illumina Inc, San Diego, CA) with PE150 according to the manufacturer’s protocol. Subsequently, data analyses were performed in silico.

The htseq-count software was used to calculate the gene expression, and the alignment of the mass below 10 in the bam file were excluded, and the obtained gene expression count matrix was used as the downstream analysis data. The raw count data were filtered based on the filter criteria that the counts per million (CPM) was greater than 1. The filtered data were standardized using the R software DESeq2 package. All samples were quality controlled based on standardized data to eliminate outlier samples. We used principal component analysis (PCA) to evaluate inter-sample relationships.

We used the R software DESeq2 package to analyze the differences in gene’s expression levels between the asthenozoospermia group and the normal control group. Genes with P value ≤ 0.05 and the log₂ fold change (FC) absolute value ≥ 1 were considered differentially expressed.

The GO and KEGG and DO enrichment analyses were applied to determine the functional roles and related pathways of DEGs using the clusterProfiler package in R [18]. The GO analysis (http://www.geneontology.org) is a functional analysis associating DEGs with the GO biological process (BP), cellular component (CC) and molecular function (MF) categories [19]. We make the GO function annotations for DEGs and explore the biological significance of each gene. The KEGG database (http://www.genome.jp/kegg/) includes biochemical reactions, signaling pathways, metabolic pathways and biological processes. KEGG pathway analysis to identify signaling pathways associated with DEGs. The DO (http://disease-ontology.org) analysis was used for the identification of DEG-disease associations [20]. P values < 0.05 were considered significantly enriched by the DEGs.

DElncRNAs and DEmRNA were identified through fold change filtering. The correlation between the expression levels of each DElncRNAs and DEmRNAs was calculated using the Pearson’s correlation test. The correlation coefficient was greater than 0.8 and P value less than 0.05 were defined as the DElncRNA–DEmRNA pairs coexpressed. In addition, we performed GO and KEGG enrichment analyses of DElncRNA co-expressing differential genes.

To identify the targeted genes of DElncRNAs, we used the bedtools software to search the nearby coding genes of DElncRNAs and excluded the gene with a distance of more than 100 kilobases upstream or downstream of DElncRNAs. Simultaneously, we took genes that intersect with the co-expressed DEGs with DElncRNAs. These obtained DEGs were defined as the cis-regulated target genes.

We used hypergeometric distribution to test whether a DEG target microRNA (miRNA) set was enriched on lncRNA target miRNA set, and P value less than 0.05 was considered statistically significant. Accordingly, a co-expression regulatory network representing the possible lncRNA–miRNA–mRNA interaction between key lncRNA and target mRNA/miRNA was established by Cytoscape software (version 3.7.1).

Each sample was assembled separately using Stringtie software, then all sample transcripts were combined using the Stringtie merges command, and then the combined results were compared with known gene models using gffcompare software to discover new transcript information. Moreover, new transcripts with lengths greater than 200 for coding potential identification were screened using the LGC software (http://bigd.big.ac.cn/lgc). Abnormal samples Z1, Z4, R2 were excluded, and we used the R software DESeq2 package for differential analysis to screen differential expression of new lncRNAs.

Trizol reagent (Invitrogen, Carlsbad, CA, U.S.A.) was used to extract RNA from the exosomes of asthenozoospermia patients (n=20) and normal healthy controls (n=20) according to the manufacturer’s instructions. RNA was reverse transcribed to cDNA using the Reverse Transcription Kit (Takara Co., Ltd., Dalian, China). The quantitative reverse-transcription polymerase chain reaction (qRT-PCR) was performed using the Applied Biosystems 7900HT system (Applied Biosystems) with an SYBR Premix Ex Taq™ kit (Takara Bio, Inc.) according to the manufacturer’s instructions. The primers were designed using Primer BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Supplementary Table S1) and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). The lncRNA expression level was verified by qRT-PCR using GAPDH as the reference gene with the 2^−ΔΔCT method. Each sample was run as three replicates and qPCR experiments were repeated at least three times (n=3).

Statistical analysis was carried out using the SPSS13.0 software (IBM Corp, Armonk, NY). The significant differences in expression levels between asthenozoospermia patients and normal healthy controls groups were tested using a two-tailed Student’s t test. The significance of the GO terms or enrichment of pathway identifiers in the DEGs was evaluated using the Fisher’s exact test. A value of P less than 0.05 was considered statistically significant.

The basic characteristics and sperm parameters of case and control are summarized in Table 1. A total of 50 subjects were recruited in the case–control study, including 25 patients (mean age ± standard deviation (SD): 32.8 ± 6.53 and 35.0 ± 6.22) and 25 controls (mean age ± SD: 32.4 ± 5.32 and 34.7 ± 5.38). There were no statistically significant differences between case and control group with regarded to age, BMI, smoking and drinking distribution (P>0.05). The sperm parameters, including sperm volume, sperm density, PR, IM, and sperm activity rate are significantly different between the case group and the control group (P>0.05).

The characterization of exosomes from human seminal plasma was shown in Figure 1. The translucent electron microscope observation results showed that the vesicles derived from normal and asthenozoospermia men displayed similar shape and size, and the ranged of vesicles from 50 to 150 nm, with most of them <100 nm in both groups (Figure 1A,B). To further characterize the vesicles, we used the Western blot to detect the exosome-specific antigen molecules CD63 and CD81 (Figure 1C). The analysis results indicated that all samples expressed these markers, and the total protein content has no significant differences between normal and asthenozoospermia groups.

We used the R software DESeq2 package to analyze the differences in gene's expression levels between asthenozoospermia and normal groups. A total of 4228 significantly DEGs, containing 2344 up-regulated DEGs and 1884 down-regulated DEGs, in asthenozoospermia group compared with normal group were obtained by |log₂FC| ≥ 1 and P-value ≤ 0.05 (Supplementary Table S2). Top 10 significantly up-regulated and top 10 significantly down-regulated DEGs (asthenozoospermia vs. normal group) were displayed in Table 2. Volcano plots and heatmaps of DEGs between asthenozoospermia and normal groups were displayed in Figure 2A,B.

To explore the potential functions of these DEGs, GO and KEGG pathway and DO enrichment analyses were applied with the down-regulated and up-regulated DEGs, separately. Through the GO analysis of DEGs, the up-regulated DEGs were found to be mostly enriched as follows: protein localization to endoplasmic reticulum and establishment of protein localization to membrane ribosome (BP), ribosome and ribosomal subunit (CC), channel activity and passive transmembrane transporter activity (MF) in the top 5 GO enriched terms (Figure 3A and Table 3); the down-regulated DEGs were found to be mostly enriched as follows: nuclear division and organelle fission (BP), lateral element and P granule (CC) and nuclease activity and vitamin binding (MF) in the top 5 GO enriched terms (Figure 3B and Table 3). KEGG pathway analysis revealed that the up-regulated DEGs were mainly associated with ribosome, transcriptional misregulation in cancer, alcoholism, cocaine addiction and ferroptosis in the top 10 KEGG pathways enriched (Figure 3C and Table 4). In addition, the down-regulated DEGs exhibited a strong association with inflammatory and metabolism pathways, according to KEGG pathway analysis (Figure 3D and Table 4). Through the DO analysis, we found that the up-regulated DEGs were significantly enriched among cancers (Figure 3E and Table 5), while the down-regulated DEGs were associated with reproductive system disease, retinal degeneration and mycosis in the top 10 DO enriched (Figure 3F and Table 5).

We compared the lncRNA and mRNA expression levels in the asthenozoospermia group with the normal group from the obtained RNA-sequencing data. In total, 995 DElncRNAs were obtained, including 656 up-regulated lncRNAs and 339 down-regulated lncRNAs in the asthenozoospermia compared with the normal control group (Supplementary Table S3). Top 10 significantly up-regulated and top 10 significantly down-regulated DElncRNAs were displayed in Table 5. The volcano plot and heatmaps for clustering analysis of DElncRNAs between asthenozoospermia and normal groups were displayed in Figure 4A,B, respectively. To investigate possible lncRNA–mRNA interactions, we carried out RNA-sequencing analysis and revealed 2338 DEmRNAs between asthenozoospermia and normal control. Among these, 1128 mRNAs were up-regulated in asthenozoospermia group, whereas 1210 mRNAs were down-regulated (Supplementary Table S4). Top 10 significantly up-regulated and top 10 significantly down-regulated DEmRNAs (asthenozoospermia vs. normal) were displayed in Table 5. The volcano plot and heatmaps for clustering analysis of DEmRNAs between asthenozoospermia and normal groups were displayed in Figure 4C,D, respectively. Moreover, we found 11,706 novel DElncRNAs (4321 up-regulated and 7385 down-regulated) in asthenozoospermia group compared with the normal group (Supplementary Table S5, Table 6 and Figure 4E,F).

To identify the correlation between DElncRNAs and DEmRNAs, we performed the Pearson’s correlation test to calculate the correlation coefficients between lncRNA and mRNAs. The correlation coefficient was greater than 0.8 and P value less than 0.05 were defined as the DElncRNA–DEmRNA pairs coexpressed. We found 254,981 positive correlations DElncRNA–DEmRNA pairs. Table 7 showed the top 10 positive correlations DElncRNA–DEmRNA pairs.

LncRNA functions by acting on protein-coding genes via cis-acting elements and trans-acting factors. In the present study, to explore how the DElncRNAs may alter sperm function. We predicted the cis- and trans-regulated target genes of DElncRNAs between asthenozoospermia and normal control groups. The closest coding genes to lncRNAs 100 kb upstream and downstream of them were screened, and their associations with lncRNA were analyzed using the Bedtools program. Seven cis-regulated correlation pairs lncRNA–mRNA (IGF2-AS/IGF2, AL139353.2/HEATR5A, AC003102.1/UBTF, AC010247.2/POU2F2, AC005034.3/MRPL19, AL031666.2/ZMYND8, and AFAP1-AS1/AFAP1), as showed in the Table 8 and Figure 5.

We used DElncRNA–DEmRNA co-expression data combined with miRNA target relationship information in the starBase database to perform hypergeometric distribution tests to construct the lncRNA–miRNA–mRNA regulatory network (Figure 6). The interaction network included 11 lncRNAs (LINC00893, AC005034.3, COX10-AS1, MIR497HG, LINC00894, AC015813.1, AP000424.1, MIR17HG, LINC00667, LINC00662 and SNHG3), 35 miRNAs and 59 mRNAs (Supplementary Table S6).

To confirm the sequencing and bioinformatics results, two up-regulated (LINC00667 and COX10-AS1) and three down-regulated (LINC00893, MIR497HG and IGF2-AS) lncRNAs in asthenozoospermia group were chosen for qRT-PCR. Consistent with sequencing results, all five lncRNA were found to be differentially expressed in the asthenozoospermia group compared with the normal (P<0.05; Figure 7).

In the present study, we performed an RNA-sequencing analysis between asthenozoospermia and normal group and identified 4228 significantly DEGs (2344 up-regulated and 1884 down-regulated) in seminal plasma exosomes RNA. About 995 known DElncRNAs (656 up-regulated and 339 down-regulated) and 11,706 novel DElncRNAs (4321 up-regulated and 7385 down-regulated) and 2338 DEmRNAs (1128 up-regulated and 1210 down-regulated) were identified in the present study. We used the GO and KEGG pathway and DO enrichment analyses to further understand the potential biological functions of these DEGs. Moreover, the lncRNA–miRNA–mRNA co-expression network was constructed, and the result indicated that some hub lncRNAs and mRNAs were involved in asthenozoospermia. From the co-expression analyses, we identified seven cis regulated correlation pairs lncRNA–mRNA. In addition, five of the DElncRNAs were verified with qRT-PCR to further validate the reliability of the RNA-sequencing results.

Exosomes are 40–180 nm double membrane microvesicles originated from cell membrane, which thought to be involved in intercellular communication by transferring information via small biological molecules (i.e. lipids, proteins, DNA, mRNA and noncoding RNA) [21]. On the membrane surface of exosomes display several markers, including CD9, CD63, CD81, LAMP1 and TSG101 [21,22]. We detect the antigen molecules CD63 and CD81 in the seminal plasma exosomes using the Western blot. The results indicated that CD63 and CD81 expressed in all samples, and the total protein content has no significant differences in asthenozoospermia group compared with the normal group. Previous studies indicated that semen exosomes affect spermatozoa motility and functions [7,8].

LncRNAs can regulate gene expression and are abundant in the genomes of various organisms [23]. LncRNAs are usually more than 200 nucleotides in length and exhibit greater species, tissue, and cell specificity than do shorter-length microRNAs (miRNAs) and messenger RNAs (mRNAs) due to their evolutionary unconserved characteristics [24]. The main functions of lncRNAs include roles in chromatin modification, transcriptional regulation and post-transcriptional regulation. In addition, lncRNAs affect the regulation of mRNA translation and stability largely based on the competitive endogenous RNA (ceRNA) regulation mechanism of binding to miRNA. Increasing evidence indicated that the ceRNA regulatory network of lncRNA–miRNA–mRNA was important in different diseases [25,26]. Accumulating evidence indicates that lncRNAs act as regulators in diverse biological processes and play important roles in many human diseases [27]. Moreover, several lncRNAs were reported to serve as disease biomarkers [28].

Recently, Zhang et al. [16] investigated the expression profiles and characteristics of lncRNAs in normal and asthenozoospermia sperm, and indicated the association between lncRNAs expression and sperm motility. In the present study, we performed an RNA-sequencing analysis of asthenozoospermia patients and normal controls seminal plasma exosomes and identified the DEGs, known DElncRNAs, DEmRNAs, the novel DElncRNAs, DEGs functions and involved pathway, and the ceRNA regulatory network of lncRNA–miRNA–mRNA. Seven cis regulated correlation pairs lncRNA–mRNA were identified in the present study.

Previous studies on the insulin-like growth factor II (IGF2) gene are the most widely reported. Spermatozoa is the only source of IGF2 mRNA since IGF2 is a paternally inherited gene. Poplinski et al. [29] conclude that low sperm counts were clearly associated with IGF2/H19 ICR1 hypomethylation, and idiopathic male infertility was strongly associated with imprinting defects at IGF2/H19 ICR1. Ni et al. [30] indicated that the aberrant methylation of IGF2 and KCNQ1 gene may be associated with sperm DNA damage. Tang et al. [31] revealed that miR-210 was involved in spermatogenesis by targeting IGF2 in male infertility. Cannarella et al. [32] found that IGF2 mRNA was found to be present in human spermatozoa, and their transcription levels were positively correlated with sperm concentration and total sperm count.

It has been reported that OCT2 (octamer binding protein 2, POU2F2) acts as a spermatogonial marker to distinguish the early stages of spermatogenesis [33]. TLE3, a transcriptional co-regulator that interacts with DNA-binding factors, plays a role in the development of somatic cells. Lee et al. [34] revealed that TLE3 is differentially expressed in Sertoli cells and plays a crucial role in regulating cell-specific genes involved in the differentiation and formation of Sertoli cells during testicular development. Huang et al. [35] found that MGAT1 was expressed at lower levels in mouse all germ-cell types, as well as somatic Sertoli cells. Previous studies showed that conditional deletion of the mouse Mgat1 gene (Mgat1 cKO) in spermatogonia causes a germ-cell autonomous defect leading to infertility, and MGAT1 regulateed ERK1/2 signaling during spermatogenesis via different mechanisms [36,37]. The loss of caseinolytic peptidase P (CLPP) leads to infertility in mice [38]. EPB41 is one of the pleiotropic monogenic disorder's genes that cause male infertility, and manifests as azoospermia [39]. Bao et al. demonstrated that LRRC8A-dependent VRAC activity is essential for male germ cell development and fertility [40]. LRRC8/VRAC anion channels are required for late stages of spermatid development in mice [41].

In the present study, our results indicated that the IGF2, POU2F2, TLE3, MGAT1, CLPP and LRRC8A were co-expressed with lncRNAs IGF2-AS, AC010247.2, LINC00893, MIR497HG, MIR497HG and LINC00893, respectively. These mRNAs expression levels were down-regulated in the asthenozoospermia group compared with the normal group. The EPB41 and ODF2L were co-expressed with lncRNAs LINC00667 and COX10-AS1, respectively. However, the two genes mRNAs expression levels were up-regulated in the asthenozoospermia group. Moreover, the qRT-PCR analysis confirmed our sequencing results. However, further studies are needed to evaluate the roles of DElncRNAs and DEmRNAs in the molecular mechanism of human asthenozoospermia.

There are some limitations that cannot be ignored in the present study. First, the study sample for RNAsequencing was relatively small, which limited the ability to determine the DElncRNA and DEmRNAs between asthenozoospermia and normal group. Further validation analysis with larger sample size and biological test was needed to confirm our conclusions. Second, the biological function study of significantly different lncRNA, mRNA and miRNA was not performed in our research. Therefore, the biological function analyses are required to provide improved understanding of the roles and mechanisms of DElncRNA and DEmRNAs in the pathogenesis of asthenozoospermia.

In conclusion, we identified key DElncRNA and DEmRNAs between asthenozoospermia and normal group. Furthermore, bioinformatics analysis identified the potential functions of DEgenes and the regulatory network was constructed. Our data provide a bioinformatics analysis of genes and pathways that may be involved in the pathological mechanisms of asthenozoospermia. However, further studies are still required to investigate their mechanisms in the occurrence and development of asthenozoospermia.

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

The authors declare that there are no competing interests associated with the manuscript.

This work is supported by the Hainan Provincial Natural Science Foundation of China [grant number 818QN315]; the Hainan Province Science and Technology Department [grant number ZDKJ2017007]; the Ministry of Science and Technology of China [grant number 2014DFA30180]; and the National Natural Science Foundation of China [grant number 81660433].

Yanlin Ma contributed to the design of the study. Hui Lu was responsible for manuscript preparation. Dongchuan Xu and Ping Wang performed the statistical analysis. Wenye Sun and Xinhuai Xue performed the experiments. Yuxin Hu and Chunli Xie contributed to collect the samples and data. All authors were responsible for drafting the manuscript, read and approved the final version.

This study was approved by the Institutional Ethical Committee of the First Affiliated Hospital of Hainan Medical College (Approval number: 2018 (Scientific Research) No. (2)). All the experiments were performed in accordance with the principles of the World Medical Association Declaration of Helsinki. Written informed consent was obtained from each subject prior to enrollment.

We are grateful to the patients and control individuals for their participation in the study. We also thank the clinicians and hospital staff that contributed to sample and data collection for this study.

DEG

differentially expressed gene

DO

Disease Ontology

FDR

false-discovery rate

GO

Gene Ontology

KEGG

Kyoto Encyclopedia of Genes and Genomes

lncRNA

long noncoding RNA

MF

molecular function

qRT-PCR

quantitative reverse-transcription polymerase chain reaction