Genome-Wide Microarray Analysis of Long Non-Coding RNAs in Eutopic Secretory Endometrium with Endometriosis

Endometriosis is a common, estrogen-dependent and refractory gynecological disorder able to induce multiple manifestations, including dysmenorrhea, pelvic pain, pelvic mass and infertility, and malignant-like behaviors [1]. The pathogenesis of endometriosis is poorly understood and the increased morbidity is seriously detrimental to the patient's quality of life. The dominant theory of the pathogenesis is retrograde menstruation of endometrial fragments sloughed through fallopian tubes into the peritoneal cavity [2]. The retrograde menstruation can keep these endometrial fragments alive in the ectopic site through adhesion, aggression and angiogenesis process and then form the ectopic lesions; therefore, the inherent biological characteristics of eutopic endometrium might have a decisive role in the pathogenesis of endometriosis. More and more pieces of evidence indicate the eutopic endometrium of endometriosis, especially secretory endometrium [3,4], has greater ability in the formation of new blood vessels, cell metastasis and invasion [5], and the level of gene expression is different compared to normal endometrium [1,6]. Therefore, abnormal alterations of eutopic endometrium might be the source of the pathogenesis of endometriosis.

In the whole genome, < 3% was found responsible for protein coding and > 80% of our genomic DNA is transcribed into noncoding RNAs (ncRNAs) [7]. MicroRNA (miRNA) is a kind of non-coding RNA, which is associated with the pathogenesis of many diseases [8]. MiRNA can bind its complementary sequence within the 3' untranslated region of its target mRNA transcript in varying degrees and can regulate the target mRNA at the post-transcription level. Earlier, we reported miRNAs are important in the development of endometriosis [9]. Long non-coding RNAs (lncRNAs) are another kind of ncRNAs, they are novel members of RNA transcripts >200 nucleotides and with no protein coding capacity [10]. Many studies have shown they can influence many biological processes, including cell differentiation, proliferation and apoptosis, material metabolism, immune regulation etc. [11]. Dysregulated lncRNAs are often involved in the pathogenesis and progression of many diseases. For example, HOTAIR is an oncogenic lncRNA up-regulated in pancreatic cancer [12], colorectal cancer [13] and breast cancer [14], and the level of HOTAIR expression is correlated closely with metastasis and poor prognosis. Many scholars focus on lncRNAs because they can act as regulatory molecules though a variety of mechanisms [15]. They can participate in X chromosome silence, genomic imprinting and chromatin modification as well as other important biological processes, including transcription activation, transcription interference. Because of the lack of specific biomarkers and different degrees of limitations of therapeutic targets and drugs, the pathogenesis of endometriosis is the key point of clinical research. The expression patterns of lncRNAs can help us to find differential expression of transcripts in endometriosis and provide information for exploring the pathogenesis.

This study reveals for the first time the different expression pattern of lncRNAs in eutopic compared to normal endometrium and were validated using qRT-PCR. We discuss the function of AC002454.1 and CDK6 in the pathogenesis and progression of endometriosis and verify their preliminary histologic correlations.

The patients in this study were enrolled at Shengjing Hospital, China Medical University from 2013-2015. We obtained eutopic endometrium from 40 patients (44.50 ± 3.12 years old) with endometriosis who were undergoing total hysterectomy as the experimental group. We collected normal endometrium from 28 patients (42.67 ± 4.02 years old) without estrogen-dependent disease at the same time. All patients had regular menstrual cycles and none had received gonadotrophin-releasing hormone (GnRH) analogues, other hormonal medications or antibiotic therapy in the 3 months before their surgery. The tissue samples were all in late secretory phase, which was confirmed pathologically. These samples were frozen in liquid nitrogen immediately and then stored at -80°C. We chose three samples of eutopic endometrium with endometriosis and three samples of normal endometrium using lncRNA microarray at random. This study was approved by the China Medical University Research Ethics Committee in accord with the Helsinki Declaration. Written informed consent was obtained from each patient before surgical procedures.

Total RNA was extracted from about 100 mg of tissue using TRIzol® reagent (TaKaRa Bio Inc., Japan) according to the manufacturer's instructions. The quantity of total RNA was assessed by spectrophotometric measurement (NanoDrop ND-1000; Thermo Fisher Scientific) of the ratio of absorbance at 260 nm to that at 280 nm (A₂₆₀/A₂₈₀) and RNA integrity was observed by standard denaturing agarose gel electrophoresis.

Human lncRNA Expression Microarray V3.0 (Arraystar Inc.) is designed for the global profiling of human lncRNAs and mRNAs, which can detect 30,586 lncRNAs and 26,109 coding transcripts. The lncRNAs were constructed using the authoritative public transcriptome databases (Refseq, UCSC knowngenes, Gencode, etc.). After hybridization, the microarrays were washed and then scanned by a DNA microarray scanner (Agilent Technology part number G2505C) and raw data were extracted using Feature Extraction software (Agilent Technology v11.0.0.1). Quantile normalization and subsequent data processing used the GeneSpring GX v12.1 software package (Agilent Technology). The microarray experiments and collection of expression data were done by KangChen Bio-tech (Shanghai, China) [16].

Total RNA was reverse-transcribed to cDNA using a PrimeScript™ RT Reagent Kit (TaKaRa Bio Inc., Japan) according to the manufacturer's instructions. The qRT-PCR primer sequences were designed by Primer 5.0 software and blasted for individual specificity in NCBI (Table 1). The qRT-PCR analysis was done with a LightCycler480 (Roche Applied Science) according to the manufacturer's protocol with SYBR Premix ExTaq (TaKaRa Bio Inc., Japan) and normalized using glyceraldehyde-3-phosphatedehydrogenase (GAPDH) as the endogenous control. The final qRT-PCR reaction mixture was as follows; 10 µl of 2 × SYBR Premix Ex Taq, 6 µl of double-distilled water, 1 µl of forward primers, 1 µl of reverse primers and 2 µl of cDNA, final volume 20 µl. The qRT-PCR reaction protocol was as follows; denaturation at 95°C for 30s, followed by 45 cycles at 95°C for 5 s, 60°C for 34 s then annealing and extension at 95°C for 15 s, 60°C for 60 s and 95°C for 15 s. Each sample was analyzed in triplicate and the expression levels of lncRNAs and mRNAs were expressed as fold change using the 2^-∆∆Ct method.

The GO project provides a controlled vocabulary to describe genes, gene products and descriptions of functions through the use of domain-specific ontologies. The ontology includes three domains: Biological Process, Cellular Component and Molecular Function. GO terms enrichment were considered statistically significant only if p<0.05 for Fisher's exact test.

Pathway analysis is a functional analysis mapping genes to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, which has become the first choice for gaining insight into the underlying biology of differentially expressed genes and proteins. The degree of enrichment of the genes in each pathway was judged by the EASE score. Fisher's exact test and the hypergeometric distribution and differences were considered statistically significant if p<0.05.

We constructed a CNC network that included the dysregulated lncRNAs and mRNAs. The Pearson correlation coefficient was calculated and the R value was used to compute the Pearson correlation coefficient between lncRNAs and their targeted mRNAs. The Pearson correlation coefficients ≥ 0.90 were selected as significantly correlated pairs to draw the CNC network using open source bioinformatics software Cytoscape (version 2.8.1). In this network, red diamond nodes represent lncRNAs and green circular nodes represent mRNAs. Pink unbroken lines represent positive correlation and blue broken lines represent negative correlation.

The relative expression levels of lncRNAs and mRNAs were presented as fold change using the 2^-∆∆Ct method. Fisher's exact test was used in GO analysis and pathway analysis. In the CNC network, Pearson correlation coefficient was used to test the correlation between expression of lncRNAs and co-expressed mRNAs. Student's t-test was used to analyze the difference between the expression levels of lncRNAs and mRNAs between eutopic and normal endometrium measured by qRT-PCR analysis. The correlation between AC002454.1 and CDK6 was analyzed using the linear correlation test. The SPSS 19.0 software package (SPSS Inc., Chicago, IL) was used for all analyses. Data were expressed as mean ± SEM and statistical significance was set at p<0.05.

These 1277 dysregulaed lncRNAs were collected from authoritative public databases such as GENCODE, RefSeq, UCSC_knowngene, pseudogene, Noncoding RNA Expression Database (NRED), Ultraconserved Regions (UCR), lincRNA_Cabili [17] and lincRNA_Khalil [18] (Fig. 1A). Hierarchical clustering was analyzed according to the lncRNAs expression values in the microarray (Fig. 1B). A scatter plot was used for assessing the difference of expression levels between eutopic and normal endometrium (Fig. 1C). The statistical significance of differentially expressed lncRNAs between the two groups was identified through volcano plot (fold change cut-off 2.0) (Fig. 1D).

We found 1277 dysregulaed lncRNAs from the microarray. Among them, 488 lncRNAs were up-regulated and 789 were down-regulated (fold change > 2.0). The most significantly up-regulated lncRNA is AC068282.3 (fold change 31.3) and RP11-403H13.1 is the most significantly down-regulated lncRNA (fold change 44.3). The expression levels of 1216 mRNA transcripts were different between the two groups (fold change > 2.0). Of these, 578 mRNAs were up-regulated and 638 were down-regulated. KMO was the most significantly up-regulated mRNA (fold change 13.7) and the most significantly down-regulated was SCGB2A2 (fold change 178.5). We divided these differentially expressed lncRNAs and mRNAs into three groups according to the fold change data (Table 2). The ten most significantly differentially expressed lncRNAs and mRNAs are given in Table 3 and Table 4.

First, we investigated the general signatures of dysregulated lncRNAs with lncRNA classification and chromosome distribution (Fig. 2A, B). According to the relative position of lncRNAs and their nearby coding gene, dysregulated lncRNAs can divided into six categories, including bidirectional, intron sense-overlapping, exon sense-overlapping, intergenic, intronic antisense and natural antisense. LncRNA classification showed the majority of lncRNA in our microarray data was intergenic lncRNA.

By filtering the microarray data, our data include three subgroups lncRNAs and their nearby dysregulated mRNAs analysis, including antisense lncRNA, enhancer lncRNA and long intergenic non-coding RNA (lincRNA). Most current research is on antisense lncRNAs, which are reported to control the level of gene regulation in pre-transcriptional, transcriptional and post-transcriptional through DNA-lncRNA, RNA-lncRNA or protein-lncRNA inter-actions [19]. We detected 13 antisense lncRNAs and their overlapping coding genes and seven pairs were in a consistent direction (fold change >2) (Table 5). Ørom et al. [20] found > 3000 lncRNAs with enhancer-like function in the enhancer region of the human genome. We detected 25 dysregulated enhancer lncRNAs with nearby differentially expressed mRNAs (distance < 300 kb) by filtering the microarray data and we found 54 differentially expressed lincRNAs and their nearby differentially expressed mRNAs (distance <300 kb) [17].

GO analysis was based on enrichment of the differentially expressed mRNAs, which include three domains; biological processes, cellular components and molecular functions (Fig. 3A, B). We found the majority of the up-regulation transcripts are about aspects of the cell cycle, including cell cycle, mitotic cell cycle, cell cycle phase, cell cycle process etc. Part of the down-regulation was associated with immune factors, including the immune system process, the cellular response to epidermal growth factor stimulus and C-C chemokine receptor activity.

Pathway analysis based on the KEGG database showed 20 and 14 pathways corresponded to up- and down-regulated mRNA transcripts, respectively (Fig. 3C, D). The cell cycle pathway was the most enriched up-regulated pathway, which contained 12 differentially expressed mRNA transcripts, and the most significantly up-regulated mRNA transcript is CDK6. In addition, MAPK signaling pathway, Wnt signaling pathway, TNF signaling pathway and cytokine-cytokine receptor interaction etc. are also enrichment in various degrees.

To confirm our microarray data, we selected five up-regulated lncRNAs at random (RP11-369C8.1|ENST00000554711, RP11-432J24.5|ENST00000455414, AC068282.3|ENST00000433673, GBP1P1|ENST00000513638, and SNHG1|NR_003098) and two down-regulated lncRNAs (AC007246.3|ENST00000451547 and FTX|ENST00000455395). We focused also on CDK6 (fold change 7.69), which is the most significantly up-regulated mRNA of the Cell Cycle pathway and its neighboring down-regulation antisense lncRNA is AC002454.1. We used qRT-PCR to measure each candidate lncRNA and CDK6 expression levels in 40 eutopic and 28 normal endometrium. The results of qRT-PCR showed RP11-369C8.1, RP11-432J24.5, AC068282.3, GBP1P1, SNHG1 and CDK6 were consistent with the microarray data and showed the same trend of up- or down-regulation (p<0.05). However, qRT-PCR showed antisense lncRNA AC002454.1 was up-regulated, whereas the microarray data suggested it was down-regulated (p 0.002) (Fig. 4). FTX was not statistically significantly different between the two groups.

The CNC network was built for six differentially expressed lncRNAs and hundreds of differentially expressed mRNAs based on correlation analysis in the eutopic endometrium in women with versus without endometriosis (Fig. 5). There are 549 nodes in the CNC network, including six lncRNAs (ENST00000554711, ENST00000455414, ENST00000433673, ENST00000513638, ENST00000451547 and NR_003098) nodes and 543 mRNA nodes. These nodes made 629 network pairs (398 positive and 2 31 negative correlations). The CNC network indicates one mRNA can be correlated with one or multiple lncRNAs. For example, ISG20 is a down-regulated mRNA in the microarray data and is associated with both RP11-369C8.1|ENST00000554711 and RP11-432J24.5|ENST00000455414 negative correlations.

AC002454.1 is an antisense lncRNA located on human chromosome 7:92465802-92546437 next to the CDK6 mRNA, which is located on human chromosome 7: 92234234-92465941 (Fig. 6A). Aberrant expression levels of AC002454.1 and CDK6 were obtained from the microarray assay by qRT-PCR and the expression levels of them in tissues were significantly correlated in the positive direction. The Pearson correlation coefficient was 0.747, p<0.0001 (Fig. 6B).

The pathogenesis of endometriosis is complex and multifactorial. To date, it has not been possible to explain the variety of clinical manifestations. Study of eutopic endometrium with endometriosis is important for exploring source molecular mechanisms of this disease in the field of gynecology and the disordered characteristics of endometrial epithelial cells also cause many diseases in gynecology [3,4,21]. A growing number of candidate miRNAs of eutopic endometrium have been found. With deeper understanding of miRNAs, lncRNAs in the pathogenesis and progression of many diseases and the regulated relationship between lncRNAs and miRNAs are cause for concern.

Recent studies have demonstrated the dysregulation of lncRNAs in the human genome can induce many diseases [12]. LncRNAs can regulate gene expression by diverse mechanisms, including transcription [22] and post-transcriptional [23] chromatin modification processing [24]. Most studies about lncRNAs are in a very early stage at present and a large number of them need extensive research. Microarray expression profiling and bioinformatics prediction can help us to find dysregulated lncRNAs and speculate about their molecular mechanisms by examining the relative position of their nearby coding gene. In recent years, the study of lncRNAs in gynecology has made rapid progress and ongoing research is providing a new direction for studying the pathogenesis of many diseases. For example, HOST2 [25] and LOC134466 [26] are known to be associated with human ovarian cancer. H19 [27] and PTENP1 [28] are expressed more highly in endometrial carcinoma. MALAT1 is involved in cervical cancer cell invasion and cell cycle progression [29]. But reports of lncRNAs as candidates for endometriosis are rare. This study revealed for the first time the different expression pattern of lncRNAs in eutopic endometrium of endometriosis and normal endometrium. Microarray expression profiling suggests 1277 lncRNAs and 1216 mRNAs were expressed differentially. We selected five up- and two down-regulated lncRNAs at random and measured their relative levels of expression by qRT-PCR. We focused on CDK6, which is the most significantly up-regulated mRNA of the cell cycle pathway and its neighboring antisense lncRNA is AC002454.1. In our results, most of them were consistent with the microarray data. There are smaller differences of expression of lncRNAs detected by qRT-PCR compared to the microarray data. Microarray analysis is a very good technique for exploring target genes but, owing to the low specificity and sensitivity of this method, its credibility is limited. The low number of two groups of samples used for the microarray analysis might limit the microarray results and these patients could have individual differences. Candidate target transcripts need to be validated with more samples. Differences of data normalization between microarray and qRT-PCR can influence the validated results [30]. In addition, our results confirmed the presence of dysregulated mRNAs in endometriosis, including S100P [31], TNFSF10 [32], IL-11 [33], MAPK8 [34] etc. Sun et al. [35] used lncRNA microarray profiles in ovarian ectopic endometrium compared with eutopic endometrium and predicted the function of dysregulated lncRNAs using co-expressed protein-coding genes. Our microarray data and the results presented by Sun et al. have some overlap, for example, LOC255167, LOC100505776, KLKP1, HNF1A-AS1 etc. These lncRNAs are expressed differently in ectopic, eutopic and normal endometrium. Therefore, study of these lncRNAs in the pathogenesis of endometriosis will be the target of more research. We also found some new mRNAs with strongly different expression in eutopic endometrium and these, as yet unidentified, mRNAs might provide new ideas and direction for future research into the mechanism underlying endometriosis.

Pathway analysis is the best choice for obtaining insight into the underlying function of dysregulated genes [36]. GO and pathway analysis can suggest what changes in gene function and which are changes of cellular pathways in the eutopic and normal endometrium. GO analysis found differently expressed transcripts associated with the cell cycle and immune regulation. In the mechanism of endometriosis, abnormal immune function has become a consensus [37,38]. The source of the molecular mechanism leading to immune escape has been one direction of exploration. This research undoubtedly provides a new clue. In addition, the role of cell cycle disorders in the occurrence and development of endometriosis is a new focus for study of the pathogenesis of endometriosis. Delayed and incomplete transition from the proliferative phase to secretory phase is a feature of eutopic endometrium of women with endometriosis. There is a greater ability for cell proliferation of the eutopic endometrium with endometriosis compared to the control eutopic endometrium [4]. Retrograde menstruation of late secretory endometrium with high proliferation ability makes it easier to form ectopic lesions. The eutopic endometrium can have malignant-like behaviors, including abnormal proliferation, cell metastasis or invasion ability because of some genes with dysregulated expression that can regulate the cell cycle. Some workers are seeking genes that can regulate the cell cycle to inhibit abnormal proliferation of endometrium cells as treatment of endometriosis. For example, AMH was able to decrease the percentage of cells in S phase and increase the percentage in G1 and G2 phase [39]. GPER1 induced endometriotic stromal cell cycle arrest and accumulation in the sub-G1 phase, which can suppress proliferation and induce apoptosis [40]. In the KEGG pathway analysis, the most significant pathway is the cell cycle, consistent with characteristics of the secretory eutopic endometrium cells in a state of continuous proliferation. In dysregulation of mRNAs of the cell cycle pathway, there is a high level of expression of CDC6 (fold change 2.16), which belongs to the AAA+ family of ATPases. The main function is to regulate the initiation of DNA replication and normal mitosis. Sideridou et al. [41] reported CDC6 can repress the transcription of E-cadherin, reduce cell connection and induce loss of cell polarity. CDC6 also serves as an important driver for epithelial-mesenchymal transition and cancer development. The expression of CDC6 was increased in the peritoneum of patients with endometriosis [34]. Therefore, CDC6 might be one of the candidate genes of endometriosis and its action might be related to the theory of epithelial-mesenchymal transition. A biology prediction method found CDC6 and the miR-199a-5p 3'UTR region have binding sites. MiR-199a-5p was down-regulated significantly in eutopic and ectopic endometrium of women with endometriosis [42]. MiRNAs can regulate the expression of lncRNA negatively or positively [43] LINC00313 is a kind of lncRNA. Its transcript-006 was down-regulated in our microarray (fold change 3.96). We used a biology prediction method based upon the CLIP-Seq experimental data (http://starbase.sysu.edu.cn/) and found LINC00313 transcript-006 and miR-199a-5p have binding sites (accession number GSE42701), which adds credibility to this study. MiR-199a-5p, LINC00313 and CDC6 might have important regulation relationships. They might be of importance to the occurrence and development of endometriosis, which warrants further in-depth study. In addition, several pathways have been proven associated with endometriosis, for example, the MAPK pathway, the cytokine-cytokine receptor interaction pathway and Wnt signaling pathway etc.

Unlike microRNA, it is difficult to predict the function of lncRNAs based solely on nucleotide sequences. We can use CNC network analysis to identify mRNAs and lncRNAs with the same expression patterns and can associate these lncRNAs with specific signaling pathways or diseases through these mRNAs and, thus, predict the function of lncRNAs. We constructed this network using six differentially expressed lncRNAs. Pearson correlation coefficients for mRNAs and lncRNAs ≥ 0.90 were selected as significantly correlated pairs. Among them, the level of SNHG1|NR_003098 expression is correlated with some members of the TNF pathway, including CREB1 and TAB2. SNHG1 is also over-expressed in gastric cancer [44]. You et al. [45] reported SNHG1 is up-regulated significantly in NSCLC cell lines resulting in enhanced proliferation. CREB1, a kind of transcription factor, is a member of the leucine zipper family of DNA-binding proteins and has an important role in cell proliferation and differentiation and the function of oncogenes [46]. The expression of CREB1 was increased in the peritoneum of patients with endometriosis [34]. Therefore, the CNC network results further verify the reliability of the microarray analysis.

In this study, we found CDK6 is the most significantly up-regulated mRNA transcript in the cell cycle pathway. CDK6, an oncogenic kinase that can regulate the cell cycle [47], is important for cell cycle G1 phase progression and G1/S transition. This kinase has been shown to phosphorylate downstream tumor suppressor protein pRb. Phosphorylation of pRb leads to the release of E2F transcription factors and induction of the expression of genes necessary for cell cycle progression. CDK6 was also over-expressed in erythroid cells [48] and in gastric cancer [49]. Otto and Sicinski [50] reported CDK6 can control the formation of blood vessels in lymphoid tumors. They first demonstrated over-expressed CDK6 induces transcription of VEGF-A, which can stimulate the formation of new blood vessels by endothelial cells. This function of CDK6 is unrelated to its kinase activity. The formation of new blood vessels was also key to endometriosis. AC002454.1 is an antisense lncRNA located on human chromosome 7:92465802-92546437 next to the CDK6 mRNA, which is located on human chromosome 7: 92234234-92465941, and they have 139 nucleotides complementarity. In our study, aberrant expression levels of AC002454.1 and CDK6 obtained by qRT-PCR and the expression levels of them in tissues were correlated significantly in a positive direction, the Pearson correlation coefficient was 0.747, p < 0.0001. Antisense lncRNAs have multiple functions and regulate protein coding genes positively or negatively. Sense RNAs and their antisense lncRNAs can form RNA duplexes by nucleotides complementarity. For example BACE1 and its antisense lncRNA BACE-AS have 104 nucleotides complementarity. The complementarity of BACE1 and BACE-AS increases the stability of BACE1 mRNA [19]. According to the relative position relations between antisense lncRNA and their sense partners in our genome, they can be classified as “head-to-head'',“nearby to tail” or “nearby to head”. The relative position relationship of AC002454.1 and CDK6 is “head-to-head” which means the 5' end of both genes align together [19]. Combination of the two complement parts might form an lncRNA-mRNA double chain involved in maintaining the stability of CDK6 mRNA and heighten its level of protein expression and enhance the proliferation and angiogenesis of endometrial cells in cell cycle disorders. We suggest AC002454.1 might affect CDK6 expression levels, leading to cell cycle disorder and they can participate together in the pathogenesis of endometriosis. This is an interesting line of study, however, the function of lncRNAs in endometriosis and related studies is not known. Attention has increased gradually in cell cycle disorder in the pathogenesis of endometriosis. CDK6 is important in regulation of the cell cycle network but there is little evidence available about the effect of CDK6 dysregulation on endometriosis. AC002454.1 regulation of CDK6 might be one of the important molecular mechanisms of higher cell proliferation rates of secretory phase eutopic endometrium with endometriosis. Hence, investigation of how AC002454.1 regulates the expression of CDK6 could be important. Our group has started to study the overall effect and mechanism of AC002454.1 and CDK6 to gain scientific information that will aid further work.

This project was supported by grants from the National Natural Science Foundation of China (NO.81270675), Science Project of Shenyang (NO.F14-231-1-55) and the Plan of Free Researchers in Shengjing Hospital (NO.201204).

All authors declare no conflict of interest.