The difference in expression of long noncoding RNAs in rat semen induced by high-fat diet was associated with metabolic pathways

Animal models and sperm collection

All protocols in this study were approved by the Animal Care Committee of the Institute of Zoology, Chinese Academy of Sciences, and all animals were treated according to the guidelines of the Animal Care Committee. Male SD rats (6-week old; Hua Fu Kang Company, Beijing, China) were used in this study. The field experiments were approved by the Research Council of Chinese Academy of Sciences (certification number SCXK [Jing] 2011-0024). After one week of acclimation, the rats were randomly divided into DIO and normal groups. The rats in the DIO group were fed on a high-fat diet (20% sucrose, 10% lard, 2.5% cholesterol, 0.2% sodium cholic acid, and 67.3% (w/w) standard chow) for 12 weeks to induce obesity (average body weight: DIO >(1 + 20%) normal, n = 7). The rats in the normal group (n = 7) were given the standard diet. Sperm was collected from the epididymis of all the rats. Sperm collection was performed according to the method described by Jiang et al. (2016). Briefly, the rats were sacrificed by cervical dislocation; the spermatozoa were extracted and placed in preheated human tubal cultures, which were then centrifuged (1, 500 × g, 4 °C, 8 min) to collect the supernatant fluid.

Sperm analysis and testis histomorphology evaluation

Semen was obtained from the tail of the epididymis and quickly placed in a Centrifuge tube with 1ml of Ham’s F10 medium (Bioway Biotechnology Co., Ltd., Beijing, China) for analysis of sperm motility and density using computer–assisted semen analyzer (CASA—TOX IVOS; Hamilton Thorne, Beverly, MA, USA). Testis tissue were collected from 19-week-old DIO and normal rats, and first fixed in 4% neutral formaldehyde fixative, then embded in paraffin. Sections of 4 µm thicknesses were used for hematoxylin and eosin (HE) staining, which was conducted according to convention methods. Finally, the morphological changes in testis tissue section was observed using an optical microscope (Olympus, Tokyo, Japan).

RNA extraction

The supernatant fluid samples from three rats each from the DIO and normal groups were selected. We extracted and purified total RNA from these samples using the miRNeasy Mini Kit (QIAGEN, GmBH, Germany) according to the manufacturer’s instructions and calculated the RIN number to assess the integration of RNA using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).

LncRNA microarray experiments

We performed the microarray analysis using the Rat Genome Oligo nucleotide 4,644 k Microarrays (Agilent, Santa Clara, CA, USA) at the Shanghai Biotechnology Corporation (SBC, Shanghai, China). We amplified and labeled the total RNA using the Low Input Quick Amp WT Labeling Kit (Cat.# 5190-2943, Agilent Technologies, Santa Clara, CA, USA), according to the manufacturer’s instructions; the labeled cRNAs were purified using the RNeasy mini kit (Cat.# 74106, QIAGEN, GmBH, Germany). Based on the instructions in the Agilent microarray supporting kit for the Hybridization Oven (Cat.# G2545A; Agilent technologies, Santa Clara, CA, USA), the conditions for hybridization were set as 65 °C at 10 rpm for 17 h, and the volume of the cRNA sample for hybridization was 1.65 µg. The slides were then washed in staining dishes (Cat.# 121, Thermo Shandon, Waltham, MA, USA) using the Gene Expression Wash Buffer Kit (Cat.# 5188–5327, Agilent Technologies, Santa Clara, CA, USA), according to the manufacturer’s instructions. The information obtained from the scanner was loaded into the image analysis program, Feature Extraction software 10.7 (Agilent Technologies, Santa Clara, CA, USA), and the data were normalized using the Quantile algorithm from GeneSpring Software 12.6.1 (Agilent technologies, Santa Clara, CA, USA).

Bioinformatics data analysis

After the original data was normalized using the GeneSpring Software (Agilent Technologies, Santa Clara, CA, USA), we screened high-quality probes for further data analysis. We analyzed fold-change (fold-differences in expression) and used t-tests (Student’s t-test) for investigating the differentially expressed genes. After the raw data from the microarray was standardized and converted to log2 values, a scatter plot was generated in a two-dimensional coordinate system. Using the online analysis software DAVID (https://david.ncifcrf.gov/), we analyzed the Gene Ontology (GO) enrichment of the differentially expressed mRNAs and their functions, based on three aspects: biological processes (BP), cellular components (CC), and molecular functions (MF). The log10 values (p-value) denote enrichment scores and represent the significance of the GO term enrichment among the differentially expressed genes (DEGs). We also performed KEGG pathway analysis to reveal pathway clusters covering the DEGs; here, the log10 values (p-value) denote the enrichment score and represent the significance of the pathway correlations.

Quantitative reverse transcription-polymerase chain reaction (RT-PCR) analysis

RT-PCR was used to confirm the lncRNA expression profile data obtained from the microarray. Total RNA was isolated using the Trizol reagent (Life Technologies). Single-stranded cDNA was prepared from 2 μg of total RNA, according to the manufacturer’s instructions (Promega, USA), and the lncRNA expression was measured using quantitative PCR using SYBR Premix ExTaq. Two microliters of each cDNA was subjected to PCR amplification using primers specific for CUST_2117_PI428311958 (uc008nvu.1), CUST_4640_PI428311958 (AY621350), CUST_9613_PI428311958 (XR_009220.3), CUST_5105_PI428311958 (FQ225056), and CUST_6638_PI428311958 (BC058491). The lncRNA primers used in this study are shown in Table 1.

LncRNA target prediction and lncRNA-mRNA co-expression network

Differentially expressed lncRNAs were selected for target prediction, as previously described (Han et al., 2012). We used two independent algorithms to identify the target genes. The first algorithm searched for those acting in cis. The University of California Santa Cruz (UCSC) gene annotations (http://genome.ucsc.edu/) were used to pair and visualize the lncRNAs in the UCSC genome browser. All genes transcribed within a 10-kbp window up- or downstream of the lncRNA were considered cis-target genes. The second algorithm was based on mRNA sequence complementarity and RNA duplex energy prediction; it evaluated the impact of lncRNA binding on complete mRNA molecules using the BLAST software for first-round screening. The RNAplex software was used to screen target genes in trans (Tafer & Hofacker, 2008), with the RNAplex parameter set as e ≤  − 30. Because the majority of identified LncRNAs functions were not clear, we established an lncRNA-mRNA co-expression network that comprised differentially expressed lncRNAs for cis- and trans-targeted mRNAs from the re-annotated Affymetrix Rats Genome Array data to reveal the connection between lncRNAs and mRNAs.

Statistical analysis

Quantile normalization and subsequent processing of the raw data were performed using the GeneSpring Software GX 12.6.1 (Agilent technologies, Santa Clara, CA, USA). The results were reported as mean ± SEM from three independent tests. Student’s t-test was performed using SPSS (13.0) to estimate the statistical significance of differences between groups. P < 0.05 was considered statistically significant.

Effects of high-fat diet on glycolipid metabolism in SD rats

We compared the fasting levels of blood glucose, low-density lipoprotein (LDL), high-density lipoprotein (HDL), and triglycerides (TG), as well as the body weight of DIO and normal rats. The results showed that the levels of HDL, LDL, and TG in the DIO group were not significantly higher than those in the control group. However, the fasting blood glucose level and body weight of the DIO rats were significantly higher than those of the control group (Fig. 1, Data S1).

Effects of high-fat diet on sperm motility and testicular morphological structure in SD rats

Semen analysis show that the sperm concentration and motility in the DIO group rats (53.63 ± 13.82∗10⁶/ml, 44.33 ± 7.81%, respectively) was significantly lower than that in the normal control group (92.18 ± 6.99∗10⁶/ml, 69.14 ± 2.46%, respectively) (p < 0.05). In the normal control group, the testicular tissue showed that the sperm cells of the seminiferous tubules were normal and the sperm cells were arranged closely. The cell structure was clear at all stages, and a large number of mature spermatozoa were found in the lumen. At the same time, the testicular interstitial cells between seminiferous tubules can be clearly observed round and distributed in clusters; compared with the control group, the number of sperm cells of obese rats induced by high fat diet was decreased and arranged disorderly in the testes. The number of mature spermatozoa decreased significantly, and the sperm cells were found fallen in clusters in the lumen, sperms appeared deformity, and substantial cells were significantly reduced (Fig. 2).

Expression levels of lncRNA and mRNA in the spermatozoa of DIO and normal rats

Rat lncRNA Microarray (V6.0) is capable of detecting 23,260 lncRNAs and 26,623 mRNAs. In the present study, 9,843 lncRNAs and 23,183 mRNAs were detected in the sperm samples. After microarray scanning and normalization, 973 lncRNAs and 2,994 mRNAs were found to be differentially expressed, with fold change ≥2.0 and P < 0.05. Among these, 457 and 516 lncRNAs were up-regulated and down-regulated, respectively, while, 1,316 and 1,678 mRNAs were up-regulated and down-regulated (fold change ≥2.0 and P < 0.05), respectively, in the sperms of the DIO rats, compared with controls. Thirty-three lncRNAs displayed fold change >10, among which 16 were up-regulated and 17 were down-regulated (Data S2; Table 2). CUST_6253_PI428311958 (fold change: 29.3) was the most up-regulated lncRNA, while CUST_3471_PI428311958 (fold change: 36.8) was the most down-regulated lncRNA in the sperms of DIO rats, compared with controls. Twenty-five mRNAs displayed fold change >15, among which 9 were up-regulated and 16 were down-regulated (Data S3; Table 3). Visualization using scatter plots showed significant variations in the expression levels of lncRNAs and mRNAs (Fig. 3).

GO and pathway analysis

LncRNAs are known to be involved in the function of the corresponding mRNA gene, and the mRNAs found to be significantly differentially expressed based on the GO enrichment analysis could reveal differences in the regulation of lncRNAs. In this study, prediction terms with p-value <0.01 were selected and ranked based on the enrichment factor ([Count/ Pop. Hits]/[List. Total/Pop. Total]) or enrichment score (−log10 [p-value]). From our data, 263 BP terms, 39 CC terms, and 40 MF terms were found (p < 0.01) in the differentially expressed mRNAs (Data S4). Here, we showed that the GO terms with the top 10 enrichment scores and those with the top 30 enrichment factors for differentially expressed mRNAs were associated with biological processes. Further, the cellular components were most relevant to an extracellular matrix component, the extracellular matrix. In addition, the GO terms for molecular function were correlated with growth factor binding, which is most important for sperm formation and semen quality (Figs. 4A– 4C).

KEGG pathway analysis of mRNAs that were significantly differentially expressed was performed to detect the pathways and molecular interactions associated with these genes. A total of 38 important KEGG pathways were found with P value <0.05, and they were ranked based on their enrichment scores (−log10 [p-value]) (Data S5). Our data showed that the pathways with the top 11 enrichment scores were associated with mRNAs. The metabolic pathway was the top pathway in protein-coding genes such as “mucin type O-glycan biosynthesis,” “tyrosine metabolism,” “protein digestion and absorption,” “complement and coagulation cascades,” “drug metabolism-cytochrome P450,” “peroxisome,” and “carbon metabolism.” The result suggested that these pathways might contribute significantly to the pathogenesis and development of DIO-associated male infertility (Fig. 5).

Verification of the microarray data by RT-PCR

We randomly selected five dysregulated lncRNAs, including both up-regulated (CUST_2117_PI428311958, CUST_4640_PI428311958, CUST_9613_PI428311958, CUST_5105_PI428311958) and down-regulated (CUST_6638_PI428311958) ones, for verification with sperm samples from three other rats, using GAPDH and RPL19 as the internal standards. The dissolution curve analysis showed a single peak, indicating that the specificity of PCR amplification and sample triplet repeat was satisfactory. The results from the RT-PCR and microarray were consistent with each other. Thus, the results of qRT-PCR verified the accuracy of the microarray data, providing valid evidence that lncRNAs might play an important role in the pathogenesis of male infertility caused by DIO (Fig. 6).

Coding-non-coding gene network

We established a form that included the differential expression lncRNAs for cis- (Data S6) and trans- (Data S7) targeted coding genes from the re-annotated Affymetrix Rat Genome Array data. In addition, we selected four up-regulated lncRNAs, CUST_2117_PI428311958 (uc008nvu.1), CUST_4640_PI428311958 (AY621350), CUST_5105_PI428311958 (FQ225056), and CUST_6805_PI428311958 (FQ212903), and two down-regulated lncRNAs, CUST_1425_PI428311958 (uc008cdl.2) and CUST_6637_PI428311958 (AF139830), for cis- and trans-targeted gene prediction. We then established an lncRNA-mRNA network. Through target gene prediction, target genes of the 26 aforementioned mRNAs were detected (Fig. 7).