LncRNAs induce oxidative stress and spermatogenesis by regulating endoplasmic reticulum genes and pathways

Oligozoospermia or low sperm count is a leading cause of male infertility worldwide. Despite decades of work on non-coding RNAs (ncRNAs) as regulators of spermatogenesis, fertilization, and male fertility, the literature on the function of long non-coding RNAs (lncRNAs) in human oligozoospermia is scarce. We integrated lncRNA and mRNA sequencing data from 12 human normozoospermic and oligozoospermic samples and comprehensively analyzed the function of differentially expressed lncRNAs (DE lncRNAs) and mRNAs (DE mRNAs) in male infertility. The target genes of DE lncRNAs were identified using a Gaussian graphical model. Gene ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways were primarily enriched in protein transport and localization to the endoplasmic reticulum (ER). The lncRNA–mRNA co-expression network revealed cis- and trans-regulated target genes of lncRNAs. The transcriptome data implicated DE lncRNAs and DE mRNAs and their target genes in the accumulation of unfolded proteins in sperm ER, PERK-EIF2 pathway-induced ER stress, oxidative stress, and sperm cell apoptosis in individuals with oligozoospermia. These findings suggest that the identified lncRNAs and pathways could serve as effective therapeutic targets for male infertility.


INTRODUCTION
Infertility is a global problem affecting human reproductive health. According to the World Health Organization, global fertility and population growth rates have continuously declined over the past 50 years ( Figure 1). It is estimated that about 10 to 15% of couples in the reproductive ages are infertile, with approximately 40 to 50% of infertility cases attributed to "male factor" [1][2][3][4]. The most common causes of male infertility are low sperm concentration (oligozoospermia) and reduced sperm motility High-throughput screening of bull semen transcriptome identified the function of ncRNAs in sperm motility [37]. We compared the sperm transcriptome profiles of patients with obstructive azoospermia (OA) (n = 3) and individuals with normal spermatogenesis (including the control group) to systematically define the expression profiles of lncRNAs and mRNAs. We next used these expression data to identify oligozoospermia-related lncRNAs and functional genes and their involvement in male infertility.

RNA sequencing and identification of DE mRNAs and DE lncRNAs
We performed a case-control study involving six normozoospermic (NS) samples from healthy fertile AGING controls and six oligozoospermic (OS) samples from infertile patients to identify the major perturbations in oligozoosperms. Cases and controls were closely matched for age, body mass index (BMI), and sperm DFI ( Table 1). The levels of reproductive hormones, such as follicle-stimulating hormone (FSH), luteinizing hormone (LH), and testosterone (T), were analyzed. As summarized in Table 1, we did not observe any differences in the levels of these hormones between NS and OS samples.
We next assessed the association between semen characteristics and oligozoospermia. Five of the six examined semen characteristics (concentration, progressive rate, non-progressive rate, motility, and normal morphology) differed between the two groups (p = 0.0001-0.048; Table 1). Compared with individuals in the OS group, those in the NS group had a 9-fold higher concentration of collected sperms (NS: 81.78 ± 37.21; OS: 8.83 ± 3.55). Similar patterns were observed for other characteristics (progressive rate, non-progressive rate, motility, and normal morphology) using the t-test ( Table 1). The remaining characteristic, ejaculate volume, was not associated with the development of oligozoospermia (p = 0.939).
We performed a high-throughput whole transcriptome shotgun sequencing to compare the sperm RNA expression patterns between NS and OS groups. We extracted high-quality RNA from each sperm sample and constructed cDNA libraries. Next, we calculated the levels of transcripts of lncRNAs and mRNAs on an Illumina HiSeq X-Ten Platform, which supported a read length of 2 × 150 bp with quality scores of ≥ 75% of bases above Q30. After filtering the raw data, 24.4 to 79.8 million clear reads and 2,501 lncRNAs and mRNAs were identified distribution of base composition and distribution of quality are shown in Figure 2. Figure 2 shows the characteristics of identified lncRNAs and mRNAs. The lncRNA transcripts were shorter in length ( Figure 2D) and contained fewer exons than mRNAs (Figure 2A). The majority of identified lncRNAs contained one to two transcripts, whereas most mRNAs had one to four transcripts ( Figure 2C)findings that are in agreement with those of previous studies [30][31][32]. The length of open reading frames (ORFs) of mRNAs was shorter than that of lncRNA ORFs ( Figure 2B).

Enrichment analysis of DE mRNAs and DE lncRNAs
We next performed GO and KEGG pathway enrichment analyses to study the functions of DE genes corresponding to DE mRNAs and DE lncRNAs. Overrepresentation (ora) GO analysis revealed 777, 196, and 188 terms enriched in biological processes (BP), cell components (CC), and molecular functions (MF), respectively. Similar results were obtained in the AGING pathway-level analysis (plage). These GO terms were primarily related to intracellular protein transport and localization. GO terms with the highest differential expression were enriched in SRP-dependent cotranslational protein targeting to the membrane, cotranslational protein targeting to the membrane, and protein targeting to the ER ( Figure 5A). Table 4 lists the DE genes associated with these processes, including CHMP4B, SEC61A1, SEC61A2, SEC61G, SEC62, SGTB, SPCS1, SRP9, SRP19, RYR2, ANK2, DDRGK1, INSIG1, KDELR2, RER1, and RTN4. Furthermore, differential GO terms were enriched in ER stress, oxidative stress, protein unfolding, and cell apoptosis ( Figure 5B, Supplementary Table 7). Twenty-six overrepresented GO terms were directly related to negative regulation of response to ER stress (Supplementary Table  7). In the case of protein unfolding, unfold protein binding had a significant difference (Supplementary Table 7). Five GO terms were related to oxidative stress, including neuron death in response to aerobic stress, negative regulation of aerobic stress-induced intrinsic acoustic signaling pathway, positive regulation of aerobic-stress induced intrinsic acoustic pathway, and regulation of oxidative stress-induced neuron death (Supplementary Table 7). Other apoptosis GO terms were those involved in the negative regulation of apoptosis (Supplementary Table 7).
The majority of DE genes were also enriched in the identified KEGG pathways, such as Alzheimer's disease pathway, Parkinson's disease pathway, ribosome and fatty acid degradation pathway (Supplementary Table 8).
Alzheimer's disease and Parkinson's disease pathways were associated with ER stress or oxidative stress-induced cell apoptosis, consistent with the previous GO terms analysis results (Supplementary Table 7). However, certain genes involved in ER stress and oxidative stress initiation in these pathways, such as APP and PRKN, were not differentially expressed in oligozoospermic individuals (Supplementary Table 2).

Co-expression analysis of DE lncRNAs
The co-expression network analysis of lncRNAs and mRNAs in oligozoospermic individuals revealed AGING potential internal adjustment mechanisms. Several lncRNAs correlated with a single mRNA and vice versa. For example, 21,468 lncRNA and mRNA partial correlation pairs were found for the top 20 DE lncRNAs (Supplementary Table 9). Of these, 1,542 pairs showed a strong correlation coefficient, 8,421 pairs showed a moderate correlation coefficient, and 11,505 pairs had a weak correlation coefficient (Supplementary Table 9).
To study the function of lncRNAs in gene regulation and pathogenesis of oligozoospermia, cis-and transpredictions were performed based on co-expression network analysis. All top 20 DE lncRNAs were predicted to have trans-regulated target genes. Lnc-TICAM1-1 and lnc-PHLDB1-regulated the expression of KDM4B and TRAPPC4 with cis-interaction effect (Supplementary Table 14). Further, lnc-TICAM1-1 may be regulated by transcription factors LF-A1, p300, and EllaE-A (Supplementary Table 14).
To validate the reliability of the lncRNA microarray data, we selected five upregulated lncRNAs (NFKB1, XBP1, SRP9, EIF2AK1, and ATF4) that were abundantly and differentially expressed (FCs > 6.0). qRT-PCR was used to analyze the differences in their expression. The results of the qRT-PCR analysis were largely consistent with those of the microarray data ( Figure 8).

DISCUSSION
Non-coding RNAs participate in the pathogenesis and development of several reproductive diseases. For example, certain miRNAs, such as mir-290-295 and mir-34c, are involved in spermatogenesis and oogenesis [27,38,39], and their dysregulation results in primordial germ cell migration disorders or downregulated expression of Nanos2 in spermatogonial stem cells [27,38,39]. We first screened the whole genome expression patterns of lncRNAs and mRNAs in spermatozoa of patients with oligospermia and normal individuals. We next compared the expression profiles of lncRNAs and mRNAs between the two groups using bioinformatics and systematically analyzed the characteristics of oligozoospermia-related lncRNAs and mRNAs [40,41].
Our results revealed that 2,364 lncRNAs and 2,599 mRNAs were abnormally expressed in patients with oligozoospermia. The proportion of downregulated lncRNAs and mRNAs was less than that of upregulated ones. Furthermore, we did not find any difference in the expression of lncRNAs by qRT-PCR and RNAsequencing [40,41], which could be attributed to different normalization processes used in the two methods [41][42][43].
The chromosomal distribution of DE lncRNAs and DE mRNAs was uneven, with the majority of them located on chromosomes 1, 17, and 19. The number of differential genes on the X chromosome was higher than that on the Y chromosome [44]. Bache conducted a cohort study of male infertility and reported 5 out of 10 disease-related autosomal bands on chromosome 1 [45], suggesting that chromosome 1 contains a domain   AGING regulating the male reproductive capacity. Further, Bache reported that men with pericentric inversion of chromosome 1 were at the risk of developing infertility due to spermatogenesis dysfunction [46]. Similar findings were reported by Balasar who found a large pericentric inversion of chromosome 1 (46, XY, inv (1) (p22q32)) during a routine chromosomal analysis [46].

AGING
LncRNAs are categorized into five subgroups depending on their different transcriptional forms: sense lncRNAs, antisense lncRNAs, intronic lncRNAs, intergenic lncRNAs, and divergent lncRNAs [47]. Intergenic and antisense lncRNAs constituted 23.2% and 48.1%, respectively, i.e., almost three-quarters of the identified DE lncRNAs. Antisense lncRNAs are transcribed from the antisense strand of the protein-coding genes, whereas intergenic lncRNAs are transcribed from both strands [47,48]. Several abnormally expressed intergenic and antisense lncRNAs identified in patients with oligozoospermia indicated that lncRNAs regulated protein-coding genes during oligozoospermia. Moreover, intergenic lncRNAs have higher tissue specificity than proteincoding genes and thus are excellent indicators of different cell subsets in the auxiliary diagnosis of oligozoospermia [47,49].
GO and KEGG analyses revealed several genes and pathways to be enriched in protein transport and localization, such as protein targeting to ER, SRPdependent co-translational protein targeting to the membrane, and protein localization, in the OS group. Thus, we speculated these DE genes to be associated with spermatogenesis dysfunction and used them as markers for predicting oligozoospermia and studying its molecular mechanism. However, we ignored the fact that the degree of spermatozoa apoptosis was considerably higher in patients with oligozoospermia than in normal individuals [50], whereas hyperactive regulation of protein translation, transcription, and transport could stimulate cell apoptosis. Therefore, there existed differences in protein transport and localization between patients and normal individuals.
Endoplasmic reticulum stress occurs due to the accumulation of unfolded proteins caused by impaired protein folding or mutations [51]. Abnormal activation of protein targeting processes transport excess peptides to ER, thus increasing the probability of protein unfolding [52]. Three branches of ER stress response have been identified, namely IRE1 (inositol-requiring-1), PERK (protein kinase RNA-like ER kinase), and ATF6 (activating transcription factor-6) [53][54][55][56].  Limitations of RNA-sequencing could result in the inaccurate measurement of underexpressed genes, thus eliminating or underestimating certain ER-related underexpressed genes in differential expression analysis, and consequently interfering with GO and KEGG pathway analyses [57,58]. Moreover, the current methods of GO and KEGG pathway analyses focus on the number of significant genes but ignore the participation of genes in different or opposite pathways.
PERK is a type-I transmembrane kinase located in the ER membrane that phosphorylates EIF2α [59,60]. Phosphorylated EIF2a initiates the translation of poorly translated mRNAs such as ATF4 [59,61]. ATF4, a transcription factor, controls redox homeostasis in the ER, and its translation results not only in amino acid biosynthesis, autophagy, protein folding, but also apoptosis via activation of the proapoptotic factor CEBPZ (growth arrest and DNA damage/CEBP homology protein, the major proapoptotic transcription factor induced by ER stress) [59,62]. We found that EIF2A, ATF4, EIF2AK1, and CEBPZ were upregulated, indicating activation of the PERK-EIF2α axis and DNA damage-induced cell death.  Genes involved in the PERK-EIF2α pathway are also well-known activators of NF-κB signaling [51,63], which is involved in the host's innate immunity responses and prevention of apoptosis by repressing CEBPZ expression [51,64]. However, we observed that both CEBPZ and NF-κB were upregulated, suggesting a complex mechanism of NF-κB-modulated gene expression and protein synthesis of CEBPZ. The IRE1-TRAF2 pathway can also activate NF-κB signaling, UPR-induced gene XBP1, and apoptosis-promoting AGING factor p53 suppressor gene JNK1 [65,66]. We observed that JNK1 was downregulated and XBP1 was upregulated; however, ERN1, the gene encoding TRAF2 and IRE1, was not expressed.
We also identified several pathways related to oxidative stress and neuronal apoptosis. Oxidative stress has been implicated in several neurodegenerative diseases, such as Alzheimer's disease and Parkinson's disease. Moreover, ER stress is known to trigger oxidative stress [67]. Thus, ER stress-triggered oxidative stress and consequent apoptosis in sperm cells could be responsible for the manifestations of oligozoospermia. Studies have reported a direct link between protein unfolding and the production of reactive oxygen species (ROS) [67][68][69][70]. Protein disulfide isomerase (PDI), endoplasmic reticulum oxidoreductin (ERO1), nuclear factor erythroid 2-related factors (Nrfs), and NADPH oxidases (NOXs) are major enzymatic components of ROS production during UPR [67][68][69][70]. Furthermore, CEBPZ induces ERO1 and PDI gene family to activate NOXs, consequently enhancing cytosolic ROS production [71,72] and depleting antioxidants such as glutathione (GSH). We found that both ERO1-α and ERO1-β, and PDI genes, such as ERP29, PDIA2, PDIA5, TMX2, and TXNDC12, were upregulated in patients with oligozoospermia. Although we did not observe the expression of NOX family genes, the downstream genes of oxidative stress, including MAPK14, SP1, GPX1, and NOX5, were activated. These findings indicated that the PERK-EIF2α pathway triggered ER stress and interfered with ROS generation via the PDI-ERO1 cycle. In addition, PERK phosphorylated and activated transcription factor Nrf2, translocating it to the nucleus to activate the antioxidant response elements and maintain GSH levels [73]. The Nrf2 gene was not differentially expressed in patients with oligozoospermia, suggesting the sperm antioxidant stress response system could be blocked by unfolded proteins-induced pro-survival and pro-apoptotic signaling events.
The ER stress-induced pre-emptive quality control (ERpQC) can selectively degrade ER-targeting proteins by inhibiting their translocation via rerouting or degradation, thus limiting further protein loading into the stressed ER [51,[74][75][76]. During ER stress, the DERLIN family proteins, the ribosome nascent chain, and signal recognition particle complex (RNC-SRP) are rerouted from the ER translocation pathway to the cytosolic degradation pathway [52,[74][75][76]. Furthermore, the E3 ligase ubiquitinates the fully translated proteins and effectively transports them to the proteasome by p97 and Bag6 chaperones for degradation by the ER-associated degradation system (ERAD) [77][78][79]. Several proteins involved in the organization of the retro-translocon channel, such as Sec61, Derlin family proteins, p97, and HRD1 E3 ligase, have been described [74][75][76][77][78][79]. However, genes encoding these proteins were not expressed in patients with oligozoospermia under PERK-EIF2α pathway-induced ER stress. This finding suggests that  the ERpQC system could become dysfunctional to protect the sperms under ER stress, resulting in conformational diseases due to accumulation of unfolded proteins.

AGING
GO and KEGG pathway analyses indicated that DE genes, such as SRP9, EIF2A, CEBPZ, EIF2AK3, and NFKB1, were mainly enriched in protein targeting, ER stress, and oxidative stress. LncRNAs are transcribed along with the co-expressed genes they regulate [47]. Therefore, the functions of lncRNAs in oligozoospermia could be studied by the function of related target genes. Their cis or trans regulation by lncRNAs, such as EIF2A trans-regulator lnc-FAM198B-1, could regulate spermatogenesis and sperm apoptosis-associated pathways and can be considered as potential diagnostic or therapeutic targets. AGING Some of the previously reported genes related to spermatogenesis disorders, such as the SOX gene family, SPATA gene family, and TAM receptor genes, were differentially expressed in our study. Consistent with the findings of a study by Jiang who reported that the Sox gene family maintains adult male fertility [80], our study found that reduced expression of SOX genes blocked the initiation or progression of spermatogenesis and affected cell maturation. Furthermore, SPATA family genes function in spermatogenesis, sperm maturation, and fertilization. Hypermethylation in the promoter regions of SPATA family genes correlates with oligozoospermia and male infertility [81]. Although SPATS2 and SPATA21 were downregulated in patients with oligozoospermia in our study, altered methylation patterns were not reported. Lu reported that mice knockout for genes (TYRO3, AXL, and MER) were deficient in sperm maturation [82].

CONCLUSIONS
We identified key lncRNAs, pathways, and gene subnetworks using RNA sequencing datasets that are associated with the pathogenesis of oligozoospermia ( Figure 9). The identified lncRNAs and pathways could serve as effective therapeutic targets for male infertility.

Sperm collection and experimental design
All 12 human sperm samples were obtained from Peking University International Hospital from March 2016 to December 2016. This study was approved by the Ethics Committees of Peking University International Hospital and was conducted in accordance with the Declaration of Helsinki. Participants provided their informed consent in writing. Participants were selected by excluding those with abnormal karyotype and Y chromosome microdeletion. The clinical data of all participants are shown in Table 1.

Collection, analysis, and purification of sperm samples
All human sperm samples were collected and analyzed according to the WHO guidelines using a Makler R chamber (Sefi Laboratories, Tel Aviv, Israel) [17]. Normozoospermic samples were defined as those with concentrations ≥ 15 × 10 6 /mL. Oligozoospermic samples had a concentration of < 15 × 10 6 /mL. All sperm samples were cryopreserved at 1:1 dilution (Irvine Scientific, Santa Ana, CA, USA) in 2 mL cryogenic vials to check the yolk buffer and equilibrated to room temperature for 15 min. Next, all sperm samples were suspended 10 cm above liquid nitrogen for 1 to 2 h and finally transferred to liquid nitrogen. All samples were thawed as follows: cryotubes were transferred to a 37° C water bath and incubated for 5 min. The cell suspension was centrifuged at 500 g for 5 min and the supernatant was discarded. Sperms were purified and the RNA was extracted immediately. The samples were stored at −80° C until further processing.

RNA extraction and quality control
The RNA was extracted from frozen human sperm samples using an RNA Extraction Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. The RNA was quantified by NanoDrop ND-1000, and the RNA quality was assessed on Agilent 2100 Bioanalyzer (Aligent Technologies, Santa Clara, CA, USA). The RNA integrity was assessed by running the sample on a standard denaturing agarose gel.

Library preparation and lncRNA sequencing
Total RNA was isolated using the Trizol reagent (Invitrogen, CA, USA). The RNA quantity and purity were detected by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent Technologies). High-quality RNA (10 μg; RNA integrity number > 7) per sample was obtained by removing ribosomal RNA using the Ribo-Zero Gold rRNA removal kit (Illumina, San Diego, USA) according to manufacturer's instructions. Approximately 1.5 μg of purified RNA per sample was used to construct the sequencing libraries using the NEBNext Ultra Directional RNA Library Prep Kit from Illumina (NEB, Ipswich, MA, USA) and index codes were added to the corresponding samples. Briefly, the purified RNA was fragmented into small pieces using divalent cations at elevated temperatures. The firststrand cDNA was synthesized using cleaved RNA fragments as templates, random hexamer primers, and M-MuLV reverse transcriptase. Subsequently, the second-strand cDNA was synthesized using Polymerase I and RNase H with reaction buffer (dTTPs were replaced by dUTPs). Exonuclease/polymerase was used to degrade the remaining overhangs and create blunt ends. After adenylation of the 3' ends of DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated for hybridization. The 300-bp DNA fragments were selected and purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). Paired-end sequencing was performed on an Illumina HiSeq X-Ten (Novogene Corporation, Beijing, China).

Identification of lncRNAs and mRNAs
To obtain high-quality lncRNA and mRNA reads, the adaptor sequences and low-quality sequences (N bases > 10% or base quality < 10; greater than 50%) were removed. The remaining reads were aligned to the UCSC (human genome, hg19) using HISAT (v0.1.6). The mapped reads were assembled by StringTie (v.1.0.4) separately. Next, all assembled transcripts from samples were merged to reconstruct a comprehensive transcriptome with Cuffcompare (v.2.1.1). To obtain highly reliable novel lncRNAs, we first removed the transcripts with a length of less than 200 bp. Transcripts with peak expression of less than 2.0 across all samples and present only in one sample were considered background and removed [18]. The transcripts overlapping with known RNAs (lncRNAs and mRNAs) on the same strand were removed. The remaining transcripts encoding any conserved protein domains were excluded by Pfam database alignment with HMMER. Finally, the transcripts predicted by the lcRNA prediction software (CPC v0.9-r2 and CNCI v2.1) were used for further analysis.

Analysis of DE lncRNAs and DE mRNAs
HTseq (v0.6.1) feature extraction software was used to calculate the number of readings of each gene (and lncRNA). The R software Limma package was used to normalize the expression data. Genes with low expression were eliminated by filterByExpr function in edgeR package with parameters: min.count = 0, min.total.count = 15, large.n = 10, min.prop = 0.7. Raw counts were converted to log2 count per million (CPM) reads. Next, we checked the outliers of the standardized data, calculated the count outside the threshold using the mean and standard deviation method (× ± 2 SD) as the outliers, and replaced the count with the predicted value of the trimmed average of all samples. This method is similar to replacing outliers in the DESeq2 package using Cook's distances. Subsequently, the limFit function of the Gaussian regression model in the Limma package was used for model fitting, and DE genes (and lncRNAs) in RNA sequencing data were identified using the empirical Bayes estimation. The regression model was adjusted for age and BMI that could interfere with the analysis. The corrected p-value (corrp) was calculated based on the Benjamini-Hochberg controlled error detection rate (FDR) [1]. The thresholds of DE genes and lncRNAs were corrected as p-value < 0.05 and fold change ≥ 2, respectively.

GO and KEGG enrichment analyses
To study the biological effects of DE lncRNAs in sperm samples, over-representation analysis method (ora) of R software gene ontology for RNA-seq package (goseq, v1.38.0) and the pathway level analysis method (plage) of gene set variation analysis package (gsva, v1.34.0) were used to determine GO terms and KEGG pathways. Based on GO and KEGG databases, these functions mapped the target genes to the corresponding pathways or biochemical processes. The p-value was calculated to indicate the significance of gene enrichment. For both GO and pathway analyses, a p-value < 0.05 was considered significant.

Co-expression network analysis
Compared with the traditional method of constructing a co-expression network using Pearson's correlation coefficient, the Gaussian graphical model (GGM) and weighted correlation network analysis (WGCNA) efficiently removed the pseudo correlation effect on the network structure. We used the R software highthroughput Gauss map model package (Statistical Inference of Large-Scale Gaussian Graphical Model [SILGGM]) to construct the co-expression network of mRNAs and lncRNAs. The identified co-expressed mRNAs and lncRNAs (p-value < 0.05) consisted of strong co-expression pairs with a partial correlation coefficient 1 ≥ r > 0.95 or -1 ≤ r < -0.95, moderate coexpression pairs with a partial correlation coefficient 0.95 ≥ r > 0.80 or -0.95 ≤ r < -0.80, and weak coexpression pairs with a partial correlation coefficient 0.8 ≥ r or -0.8 ≤ r. Further, we classified the co-expressed lncRNAs with a partial correlation coefficient > 0.95 or < -0.95 through cis-trans analysis. LncRNAs having flanked genes within 300 kb were considered as cisacting lncRNAs, and those with lncRNAs-gene pair sequence homology (blastn, E-value < 1.0E-10 and identity > 99 and matched length ≥ 20 bp) were considered as trans-acting lncRNAs.

QPCR analysis
Isolated RNA was reverse-transcribed to cDNA using a Reverse Transcription Kit (Takara, Dalian, China). The qRT-PCR analyses were performed using a StepOnePlus RT-PCR Instrument with Power SYBR Green (Takara, Dalian China). The qRT-PCR conditions were as follows: 95° C for 2 minutes, followed by 40 cycles of 95° C for 15 seconds and 60° C for 30 seconds. All experiments were performed and analyzed in triplicate. The primers used in this study were listed in Table 1. Then, lncRNA expression levels were normalized to GAPDH and calculated using the 2−ΔΔCt method.

Statistical analysis
Data are expressed as mean ± standard error of the mean (SEM). The p-values < 0.05 were considered significant. Graphs were plotted and analyzed using GraphPad Prism version 6.0 (GraphPad Software, La Jolla, CA, USA).

Availability of data and materials
Data generated during the study and presented in the manuscript are available from the corresponding author upon request.

Ethics approval and consent to participate
All experiments involving human sperm samples in this study were approved by the Ethics Committees of Peking University International Hospital. The participants provided their informed consent in writing.

Editorial note
& This corresponding author has a verified history of publications using a personal email address for correspondence.

AUTHOR CONTRIBUTIONS
Li Tian, Ya-Peng Wang, and Shou-Long Deng conceived the study and designed the experiments. Tie-Cheng Sun and Kun Yu performed the experiments. Yi Zhang and Tie-Cheng Sun analyzed the data. Tie-Cheng Sun, Shan-Jie Zhou, Hong Yu, and Yi Zhang contributed reagents/materials/analysis tools. Tie-Cheng Sun, Yi Zhang, and Kun Yu wrote the manuscript.