Testicular biopsies microarray analysis reveals circRNAs are involved in the pathogenesis of non-obstructive azoospermia

Circular RNAs (circRNAs) have been reported to be involved in many diseases. But there is no report on circRNAs in non-obstructive azoospermia (NOA). The purpose of this paper is to explore the circular RNA expression profile and potential functions of circRNAs in NOA patients. We first preformed circRNA expression profiling analysis using a circRNA microarray in testicular samples from NOA and obstructive azoospermia (OA) patients. CircRNAs were validated by qRT-PCR. Bioinformatics analysis were used to construct the ceRNA network. GO and KEGG enrichment analysis were performed by using DAVID. Microarray analysis identified 82 differentially expressed circRNAs in NOA specimens. The differential expression of hsa_circRNA_402130, hsa_circRNA_072697, hsa_circRNA_030050, hsa_circRNA_100812 and hsa_circRNA_406168 was confirmed by qRT-PCR. Enrichment analysis revealed the association of hsa_circRNA_402130 and hsa_circRNA_072697 with multiple signaling pathways. The data indicated that circRNAs were significantly dysregulated in NOA specimens and might involve in the pathogenesis of NOA.


INTRODUCTION
Circular RNAs (circRNAs) are a class of closed circular RNA molecules that rarely encode proteins and have no 5 'caps and 3' poly(A) tails. Their unique characteristics, such as better stability and resistance to nucleases make circRNAs potential candidates for clinical diagnostic/prognostic biomarkers [1]. The molecular mechanisms of circRNAs are diverse. In addition to acting as miRNA sponges, they may also play a role by binding protein or DNA sequences [2][3][4]. Some of them can also regulate gene expression by affecting RNA splicing or translating into short peptides [5]. They are conserved among different species, and also exhibit tissue-specific and developmental-stage specific traits [6]. Emerging evidences suggest that circRNAs are likely to play important function in different physiological or pathological processes. A multitude of studies have demonstrated that circRNAs are associated with various diseases, including cancer, Alzheimer's disease, ovarian endometriosis, preeclampsia [7][8][9][10][11]. Moreover, a previous study has reported that there are 15,996 circRNAs in normal human testis and more than 20 testis-derived circRNAs are stably expressed in seminal plasma [5]. All the above evidences indicate that circRNAs might play a role in spermatogenesis and maturation. However, the global expression and function of circRNAs in nonobstructive azoospermia (NOA) are largely unknown.
NOA is considered a major cause of male infertility [12,13], which can't be completely cured by assisted reproductive technology (ART). Several genetic and AGING epigenetic factors, especially Y chromosome microdeletion, are associated with the pathogenesis NOA [14][15][16][17][18]. In recent years, genome-wide association analysis (GWAS), has led to the identification of a few key NOA risk loci including PRMT6, PEX10, HLA-DRA [14,16]. In addition to these genetic factors, epigenetic factors also play an important role in the development of NOA. For instance, Chencheng Yao et al. found that there were many differentially expressed miRNAs in testicular spermatogonia, pachytene spermatocytes, and round spermatids between NOA and OA patients [13]. Further, another group also found that miRNA-188-3p could regulate spermatogenic cell apoptosis by targeting MLH1 [19]. However, the epigenetic regulatory mechanisms, especially the role of circRNAs, in NOA are yet unknown.
Hence, the purpose of this study was to elucidate the expression profile of circRNAs and unravel the possible functions and interactions of circRNAs in NOA.

Identification of the dysregulated circRNAs between NOA and OA groups
We performed microarray-based profiling analyses to screen the differentially expressed circRNAs in 6 NOA and 3 OA tissues. The box plot depicts similar distributions of the intensities from all the normalized datasets ( Figure 1A). Raw variation of circRNA expression between the two groups is shown in the scatter plot ( Figure 1B). Volcano plots and hierarchical clustering exhibit the differentially expressed circRNAs between the NOA and OA samples (Figure 2A and 2B). Using a threshold of fold change (FC) ≥ 1.5 and a pvalue < 0.05, we identified 82 differentially expressed circRNAs, of which 16 were up-regulated and 66 were down-regulated in NOA tissues compared with OA tissues (Table 1 and Table 2). According to the annotation of human circRNAs, we found 81.71% of the dysregulated circRNAs were exonic circRNAs, 9.76% were intronic circRNAs and 8.54% were sense overlapping circRNAs ( Figure 2C). Further analysis revealed that these dysregulated circRNAs were mostly distributed on chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr16, chr17, chr19, chr20, chr21 and chr22, but not on chr15, chr18, chrX and chrY ( Figure 2D).

MicroRNA response elements analysis of validated circRNAs
Several studies have demonstrated that circRNAs regulate the expression of their target genes by competing with their target miRNAs for binding [8,20,21]. The Arraystar's homemade miRNA target prediction software based on TargetScan and miRanda was utilized to predict the circRNA-microRNA interactions of two confirmed circRNAs (hsa_circRNA_402130 and hsa_circRNA_ 072697). The top 5 targeted miRNAs were then identified based on the miRNA support vector regression (mirSVR) scores. Studies have shown that some NOArelated miRNAs are dysregulated in the testicular tissues of patients with NOA [22,23]. So we selected the the potential miRNAs shared by the software and the reports to construct circRNA-miRNA interactions. The 2D structure was created from the sequence analysis of microRNA response elements (MREs) ( Figure 4A and 4B). For hsa_circRNA_402130, the potential miRNA targets included hsa-let-7b-5p, hsa-let-7c-5p and hsa-let-7i-5p, and for hsa_circRNA_072697 the miRNA included hsa-miR-23c and hsa-miR-182-5p.

CircRNA-miRNA-mRNA interaction analysis
To dissect out the functions of the two validated circRNAs, starBase v2.0 and miRDB identified down-stream target genes of the abovementioned candidate miRNAs. The two circRNA-miRNA-mRNA regulatory networks (built by OmicShare tools) included 3 miRNAs and 212 mRNAs in the ceRNA network of hsa_circRNA_402130 ( Figure 5A, Supplementary  Table 1), and 2 miRNAs and 411 mRNAs in the ceRNA network of hsa_circRNA_072697 ( Figure 5B, Supplementary Table 2). As shown in Figure 5A, we found that this network contains multiple transcription factors such as HMGA2, SOX13, E2F5, LIN28B and so on. Similarly, multiple transcription factors such as FOXO3, FOXP2 and SP3 are included in the regulatory network of Figure 5B. These regulatory networks revealed that circRNAs may play a crucial part in the process of NOA.

GO and KEGG enrichment analysis of the predicted network genes for hsa_circRNA_402130 and hsa_circRNA_072697
The functions of all the target genes were analyzed by online biological classification tool DAVID. The top 10 GO terms and KEGG terms are displayed as bubble diagram. The GO biological process (BP) analysis ( Figure 6A) showed that the target genes were significantly enriched in positive regulations including transcription from RNA polymerase II promoter, transcription, cell migration, protein phosphorylation, AGING palate development, pathway-restricted SMAD protein phosphorylation, circadian regulation of gene expression. KEGG pathway analysis ( Figure 6B), showed target genes were enriched for microRNAs in cancer, MAPK signaling pathway, signaling pathways regulating pluripotency of stem cells, FoxO signaling pathway, ubiquitin mediated proteolysis, axon guidance, transcriptional misregulation in cancer, proteoglycans  in cancer and TNF signaling pathway. These findings imply that circRNAs might participate in the development of NOA through multiple signaling pathways.

DISCUSSION
Results from the Human Genome Project and the DNA Elements Encyclopedia Project indicated that proteincoding gene sequences account for only 1-3% of the human genome sequence. Also, most of the transcribed sequences in the human genome are non-coding RNA sequences [24]. Studies demonstrated that non-coding RNAs including miRNAs and lncRNAs participate in the occurrence and development of NOA [13,18,22,23,[25][26][27]. However, the role of circRNAs, a class of non-coding RNA, has not been reported in NOA.
Here, we constructed the expression profile of circRNAs in NOA by high-throughput circRNA microarray for the first time. 16 up-regulated and 66 down-regulated circRNAs in NOA were identified. As most of these circRNAs originate from exons or introns, these results imply that dysregulated expression of circRNAs may be involved in the process of NOA. Subsequently, differential expression of circRNAs was confirmed by qRT-PCR. As circRNAs usually function as RNA sponge to regulate gene expression [4,20,21], bioinformatics was used to predict target miRNAs of hsa_circRNA_402130 and hsa_circRNA_072697. The results showed that hsa-let-7b-5p, hsa-let-7c-5p and hsa-let-7i-5p had a binding site for hsa_circRNA_402130, and hsa-miR-23c and hsa-miR-182-5p had a binding site for hsa_circRNA_ 072697. These indicated that hsa-let-7c-5p, hsa-let-7b-5p and hsa-let-7i-5p were the targets for hsa_ circRNA_402130, and hsa-miR-23c and hsa-miR-182-5p were the targets for hsa_circRNA_072697. In addition, these predicted miRNAs were previously reported that hsa-let-7b-5p, hsa-let-7c-5p and hsa-let-7i-5p were up-regulated and hsa-miR-23c and hsa-miR-182-5p were down-regulated in NOA [22,23]. Among these target miRNAs, let-7 family is an important regulator of stemness by regulating LIN28A, a gene that plays pivotal roles in organogenesis, embryonic development and tumorigenesis [28][29][30]. Jixin Tang et al. [31] indicated that let-7 miRNA family (let-7c-5p and let-7i-5p) regulate spermatogenesis in mice through LIN28A. Another member of the let-7 miRNA family let-7g could significantly reduce the germ cell pool in mice [32]. Based on these findings, we speculate that hsa_circRNA_402130 may cause NOA by regulating the expression of let-7 miRNA family.
Moreover, based on the above two circRNAs, we constructed "circRNA-miRNA-mRNA" ceRNA regulatory networks. From these regulatory networks, we found many transcription factors that regulate spermatogenesis. Studies have shown that SP3 presents a unique expression pattern during spermatogenesis in mice and other mammalian [33,34]. It has been confirmed that HMGA2 was expressed during the process from spermatocyte to spermatid in mice, and the phenotype of HMGA2-null mice was few spermatids and complete absence of spermatozoa [35]. We also found some protein coding genes that are AGING affected by let-7 miRNAs involved in spermatogenesis. TET3 is expressed in round spermatids and sperm and its effect on sperm methylome might be essential in human [36]. As early as 2006, a group found that CDC25A was down-regulated in the testis tissue of male infertile patients [37]. Subsequent studies suggested that CDC25A might be regulated by BOLL and participate in human spermatogenesis [38]. These evidences all indirectly indicated that circRNAs were involved in spermatogenesis. However, experiments are required to confirm these interactions and, to elucidate the localization of circRNAs and target miRNAs.
Additionally, we also performed GO and KEGG enrichment analysis of miRNA target genes for the first time. GO enrichment analysis revealed miRNA target genes that positively regulate transcription from RNA polymerase II promoter, transcription, cell migration, protein phosphorylation, palate development, pathwayrestricted SMAD protein phosphorylation, circadian regulation of gene expression and so on. Given that phosphorylation of SMAD promotes differentiation of mouse SSC [39], and proliferation of germ cells in mice [40], it is possible that circRNAs may regulate spermatogenesis by affecting SMAD activity in NOA patients. Among the KEGG pathways associated with the two circRNAs, FoxO signaling pathway and MAPK signaling pathway are strongly associated with spermatogenesis in mice [41][42][43]. Based on our and previous data, we can easily say that humans and mice share many common ways in the regulation of spermatogenesis. More studies on NOA patients should be conducted to improve our understanding of the mechanisms of spermatogenesis and NOA to validate the results of animal experiments. In addition, we identified a signaling pathway involved in regulation pluripotency of stem cells within the network of the dysregulated circRNAs. This is significant as self-renewal and mitosis of spermatogonial stem cells are important steps in human spermatogenesis. Based on the abovementioned ceRNA networks and enrichment analysis results, we have discovered many molecules and pathways related to stemness. Thus, we suspect that these two circRNAs may be involved in the development of NOA by regulating the stemness of spermatogonial stem cells. In a word, our study provides new ideas and strategies for NOA research. But whether circRNAs are involved in these processes deserves to be confirmed by subsequent in vitro and in vivo experiments.
Our study also has two limitations. First, more testicular tissue samples could be included in this study. This might reduce the error caused by individual differences of patients. In addition, due to the difficulties in acquiring testicular samples, we were unable to use the testicular samples of volunteers who had already been known fertility and normal spermatogenesis as normal controls.
In summary, this is the first report that reveals the expression profile of differentially expressed circRNAs and, in conjunction with bioinformatics analysis, provides the first assessment of ceRNA networks and circRNAs associated pathways in NOA. Results indicated that circRNAs may play important functions and have potential to become therapeutic targets for NOA.

Testicular tissue samples collection
The study was approved by the Ethics Committee of the Institute of Human Reproduction and Stem Cell Engineering of Central South University. Testicular biopsies were collected from 6 NOA and 3 OA patients. All participants signed informed consent, and routine semen analysis based on the World Health Organization criteria showed no sperm. Most patients have been tested for follicle-stimulating hormone (FSH), luteinizing hormone (LH) and testosterone (T). All the patients underwent testicular fine needle aspiration for histological examination at Reproductive and Genetic Hospital of CITIC-Xiangya. Those patients without testicular sperm were defined as NOA patients. OA patients had normal spermatogenic function accompanied by blockage of the vas deferens, but no other congenital diseases. The results of the pathological examinations are in the supplementary figure (Supplementary Figure 1). None of the patients were exposed to radiation, high temperature or toxic substances prior to their hospitalization. NOA patients with chromosomal abnormalities, Y chromosome microdeletion, and varicocele were excluded from the study.

Total RNA extraction
Nine RNA samples of testicular biopsies from NOA and OA patients were extracted using TRIzol (Invitrogen, USA) following the manufacturer's protocol. The quantity and purity of total RNA was measured on the NanoDrop® ND-1000 spectrophotometer (Thermo Fisher, USA). RNA integrity was measured by denaturing agarose gel electrophoresis.

Microarray detection
All RNA samples were analyzed by Arraystar circRNA Microarray at Kangchen Bio-tech (Shanghai, China). Sample labeling and array hybridization were performed according to the manufacturer's protocol AGING (Arraystar Inc.). The hybridized arrays were washed, fixed and scanned on the Agilent Scanner G2505C. Agilent Feature Extraction software was used for raw data extraction. The data was then processed using the limma package of R software. Low intensity was filtered after quantile normalization of the raw data. Student's t-test was used to estimate the statistical significance between groups. CircRNAs with a fold change ≥ 1.5 and p < 0.05 were considered significant.

Annotation and functional prediction of hsa_circRNA_402130 and hsa_circRNA_072697
The circRNA-microRNA interactions of hsa_circRNA_ 402130 and hsa_circRNA_072697 were predicted with Arraystar's home-made microRNA target prediction software based on TargetScan [44] and miRanda [45]. The analysis process of the software is as follows: First, it performs target prediction by using the miRNA information in the latest miRBase on the target sequence. Then, candidate miRNAs are screened by Context+<=-0.05 and Context<=-0.05. Finally, according to the miRNA coverage of MREs (MicroRNA Recoginition Elements) on the target sequence, the required circRNA-microRNA interactions are obtained.

Statistical analysis
Student's t-test (two-tailed) was used to estimate statistical significance between groups. Data analysis was performed using GraphPad Prism 5.0 (GraphPad Software, La Jolla, CA). A p-value < 0.05 was considered significant.

AUTHOR CONTRIBUTIONS
Liqing Fan designed and guided the research. Hao Bo and Fang Zhu performed the experiments. Hao Bo, Han Zhang, Zhizhong Liu and Xingming Wang analyzed the data. Ruiling Tang, Fang Zhu, Dai Zhou and Guanghui Gong collected and prepared the samples. Hao Bo and Zhizhong Liu wrote the paper. Wenbin Zhu, Yueqiu Tan and Liqing Fan revised the manuscript, and all authors had approved the final manuscript.