Identification of LINC01279 as a cell cycle‑associated long non‑coding RNA in endometriosis with GBA analysis.

Endometriosis affects 6‑10% of women of reproductive age. Though a significant amount of research has explored the pathogenesis of endometriosis, little is clear. Elucidating the mechanisms is urgently required for improving the therapeutic efficiency of endometriosis treatment. Long non‑coding RNAs (lncRNAs) have recently acquired extensive attention as regulatory components in a variety of biological processes and diseases. However, the functions of many lncRNAs in endometriosis are poorly understood. Therefore, the exploration of the dysregulated genes in endometriosis, particularly lncRNAs, is of importance. In the present study, datasets for endometriosis, including GSE7305, GSE7846, GSE29981 and E‑MTAB‑694, were downloaded from Gene Expression Omnibus and ArrayExpress. Then, the limma and Affy packages were used to analyze the CEL file. The RankProd method was used to conduct meta‑analysis. Long intergenic non‑protein coding RNA 1279 (LINC01279) was significantly upregulated in the three datasets, and was the most upregulated lncRNA as determined by the RankProd method. Gene set enrichment and Gene Ontology analyses were conducted, which revealed that LINC01279 is likely to function as a cell cycle mediator in endometriosis. Finally, it was identified that LINC01279 is strongly associated with certain previously identified key factors in the development of endometriosis, including cyclin‑dependent kinase 14 and C‑X‑C motif chemokine ligand 12. Thus, it was demonstrated that LINC01279 may be associated with the pathogenesis of endometriosis. This may potentially represent a target in the therapy of endometriosis.

Introduction

Endometriosis, the presence of endometrial tissue outside of the uterus, affects 6–10% of women of reproductive age (1,2). Among patients with endometriosis, 40–50% of them experience problems with fertility, ~50% experience severe chronic pelvic pain, and some may have dyspareunia and irregular uterine bleeding. Endometriosis significantly affects life quality and work productivity (3–5). Though a significant amount of research has explored the pathogenesis of endometriosis, little is clear. Illuminating the mechanisms is urgent for improving the therapeutic efficiency of endometriosis treatment.

In recent years, the rapid advance of high-throughput sequencing-based gene expression profiling has facilitated the identification of more and more long non-coding RNAs (lncRNAs). Defined as non-coding RNAs over 200 nt in length, lncRNAs were initially considered to represent transcriptional noise, but have since received extensive attention as novel regulatory components in a variety of biological processes and diseases, including cancer, heart disease and celiac disease (6–8). In addition, studies have reported that a number of lncRNAs mediate biological processes in endometriosis. H19 regulated stromal cell growth via the IGF signaling pathway in the endometrium of patients with endometriosis (9). In addition, Lee et al (10) demonstrated an association between genetic polymorphisms of CDKN2B-AS and WNT4, and Korean patients with endometriosis. Nakaoka et al (11) identified an allelic imbalance in the regulation of ANRIL mediated by chromatin interaction at the 9p21 endometriosis risk locus. Although the ENCODE project and GENCODE annotation have identified thousands of lncRNAs, the effect of many lncRNAs to endometriosis is not yet understood (12,13).

In recent decades, microarrays have been widely used to identify candidate biomarkers and therapeutic targets through investigating the alteration of gene expression at a genome-wide level. Several studies on endometriosis have been performed to investigate the dysregulated genes and essential pathways involved in its pathogenesis. For example, Houshdaran et al (14) reported the aberrant endometrial DNA methylome and associated gene expression in patients with endometriosis; Wang et al (15) analyzed the serum microRNA profile by Solexa sequencing in patients with endometriosis; Wang et al (16) reported the differential expression of microRNAs in ectopic endometrial tissue during the implantation window for endometriosis patients with related infertility. Unfortunately, only a small number of significantly dysregulated genes between the normal endometrium and endometriosis have been identified. This motivated us to explore the dysregulated genes in endometriosis, specifically lncRNAs.

Materials and methods

Reannotation of Affymetrix human genome U133 plus 2.0

Reannotation of the probe IDs for Affymetrix Human Genome U133 Plus 2.0 (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) was required to obtain data regarding lncRNA expression. Initially, we downloaded the HG-U133 Plus 2 annotation file from Affymetrix (CVS format; edition no. 36; www.affymetrix.com/estore/index.jsp; Affymetrix; Thermo Fisher Scientific, Inc.), which included 54,674 probe IDs for coding RNAs, microRNAs, lncRNAs and other non-coding RNAs. The criteria for assigning lncRNA probe IDs were as follows: i) The Refseq ID started with NR_ or XR_, indicative of non-coding RNAs; ii) the Ensemble gene IDs' annotations were antisense, processed transcripts, sense overlapping, non-sense mediated decay, sense intronic or lincRNA; and iii) pseudogenes, rRNAs, microRNAs, and other small RNAs including tRNAs, snRNAs and snoRNAs were excluded. Finally, 3,809 probe IDs were reannotated as lncRNAs in Affymetrix Human Genome U133 Plus 2.0 (Affymetrix; Thermo Fisher Scientific, Inc.) representing 2,964 different lncRNAs (17).

In addition, mRNA probe IDs were identified by excluding non-coding RNAs (lncRNAs, microRNAs, pseudogenes, rRNAs, tRNAs, snRNAs and snoRNAs). Probe IDs encoding more than one mRNA and those without RefSeq Protein IDs were excluded. Finally, 38,429 probe IDs remained, representing 17,510 mRNAs.

Individual microarray data analysis

To identify co-expressed mRNAs that were differentially expressed lncRNAs and lncRNA in endometriosis, the raw CEL files for GSE7305 (10 normal samples and 10 endometriosis samples) (18), GSE7846 (5 normal samples and 5 endometriosis samples) (19) and GSE29981 (20 endometriosis samples) were downloaded from GEO (www.ncbi.nlm.nih.gov/gds/), and E-MTAB-694 (17 normal samples and 18 endometriosis samples) (20) from ArrayExpress (www.ebi.ac.uk/arrayexpress). The limma and Affy packages were used to read and process the raw CEL files. The datasets were normalized using the robust multichip averaging (RMA) method in the Affy package to return log2-transformed expression intensities (21). The normalized datasets were then subjected to analysis in limma to identify the differentially expressed genes. Genes significantly dysregulated in endometriosis compared with normal tissue were defined by the thresholds of P<0.05, and log2 fold change (FC)>1 (upregulated) or <-1 (downregulated). The pheatmap package was used to produce a heap map of the differentially expressed genes, with a green, black and red color scale from low to high expression. Finally, the differentially expressed genes were overlapped among the GSE7305, GSE7846 and E-MTAB-694, and the results were represented with a Venn diagram produced with the Venny online tool (bioinfogp.cnb.csic.es/tools/venny/index.html).

Meta-analysis of multiple microarray datasets

Meta-analysis was performed for the GSE7305, GSE7846 and E-MTAB-694 which contain normal samples and endometriosis using the RankProd package (22) to identify the upregulated and downregulated genes between normal endometrial tissue and endometriosis. Initially, the raw CEL files were normalized using RMA in the Affy package to return log2-transformed intensities (21). The normalized datasets were then merged using the inSilicoMerging package, and batch effects were corrected using COMBAT (23). To identify the most differentially expressed probe sets, the RP advance function within the RankProd package was used (22). False discovery rates (pfp) for differential expression were determined using 1,000 permutations. Upregulated or downregulated probe sets were identified based on false discovery rate (pfP<0.01) and FC value (FC>1, upregulated or FC<1, downregulated). Probes that mapped to multiple genes were discarded, to avoid the misinterpretation of the results and to increase specificity.

GO analysis and gene set enrichment

In order to investigation the function of the overlapping lncRNAs, guilt-by-association (GBA) analysis was used, in which the function of a poorly characterized lncRNA gene can be inferred on the basis of the functions of co-expressed protein-coding genes (PCGs) (24). Initially, the top 500 positively Pearson's correlated protein coding genes were subjected to GO analysis to identify the functional of the poorly characterized genes. In addition, gene set enrichment analysis (GSEA) was conducted to further investigate the functions of the lncRNAs. In GSEA, the median LINC01279 expression value was used to divide samples into high and low groups according to the associated LINC01279 expression level. As a metric for ranking genes in GSEA, the difference between the mean expression between samples with low and high LINC01279 expression was used; other parameters were used at their default values. Besides, the motifs c2.cp.kegg.v6.1.symbols.gmt and c5.bp.v6.1.symbols.gmt were enrolled to our analysis for GSEA. Finally, the Pearson correlations between LINC01279 and the candidate genes were computed by the Cor package.

Results

LINC01279 and MSC-AS1 were upregulated in endometriosis

To identify the potentially dysregulated lncRNAs in endometriosis, we performed an integrative analysis of microarray expression profiles from GEO datasets and ArrayExpress datasets created with the Affymetrix Human Genome U133 Plus 2.0 platform (Affymetrix; Thermo Fisher Scientific, Inc.). To obtain the dysregulated lncRNAs, the GSE7305, GSE7846 and E-MTAB-694 were enrolled in our study. As GSE29981 do not include normal samples and we do not enrolled this profile for dys-regulated lncRNAs analysis. First, the upregulated and downregulated lncRNAs were identified with the limma package and presented as a heatmap (Fig. 1A). As depicted in Fig. 1A, based on the cutoffs of P<0.05 and log2 FC>1, there were 32 upregulated lncRNAs in E-MTAB-694, 38 in GSE7305 and 55 in GSE7846. In addition, there were 20 downregulated lncRNAs in E-MTAB-694, 32 in GSE7305 and 2 in GSE7846. The overlapping upregulated and downregulated lncRNAs were represented with a Venn diagram; as shown in Fig. 1B and C, LINC01279 and MSC-AS1 were consistently upregulated in all three datasets, while there were no downregulated genes common to all three datasets. As no previous studies have reported the functions of these two lncRNAs to the best of our knowledge, the exploration of the functions of LINC01279 and MSC-AS1 in endometriosis may be of importance.

Figure 1.

Differentially expressed lncRNAs in endometriosis. (A) Heatmaps of differently expressed lncRNAs in GSE7305 (n=20), GSE7846 (n=10) and E-MTAB-694 (n=35). (B) Overlapping upregulated lncRNAs in GSE7305, GSE7846 and E-MTAB-694. (C) Overlapping downregulated lncRNAs in GSE7305, GSE7846 and E-MTAB-694. lncRNAs, long non-coding RNAs.

LINC01279 was the most upregulated lncRNA in endometriosis, as determined with the RankProd package

As described above, LINC01279 and MSC-AS1 were upregulated in three independent datasets. In order to identify which was the most upregulated in endometriosis, multiple microarray datasets were analyzed by the RankProd package (Fig. 2A). As shown in Table I, LINC01279 was the most upregulated lncRNA as identified in RankProd. These results suggested that LINC01279 may play an important role in the development of endometriosis.

Figure 2.

LINC01279 is significantly upregulated in endometriosis. (A) (top panel) The differentially expressed long non-coding RNAs in GSE7305 (n=20), GSE7846 (n=10) and E-MTAB-694 (n=35) as identified with the RankProd method; (bottom panel) the heatmap of LINC01279 pre- and post-merged by RankProd. (B and C) GO analysis was performed for the top 500 genes positively correlated with LINC01279 in (B) GSE29981 and (C) E-MTAB-694. The significantly enriched terms included GO:0000086 (G2/M transition of mitotic cell cycle) and GO:0051301 (cell division). GO, Gene Ontology; LINC01279, long intergenic non-protein coding RNA 1279.

Table I.

Top 20 upregulated long non-coding RNAs in endometriosis by Rankprod analysis.

Probe ID	ENTREZ_GENE_ID	Symbol	FC_class1/class2	pfp	P-value
227061_at	100506621	LINC01279	10.7411	0	<0.0001
226582_at	400043	LOC400043	5.3191	0	<0.0001
229669_at	440416	LOC339260	4.1442	0	<0.0001
229669_at	440416	CCDC144NL-AS1	4.1442	0	<0.0001
228564_at	375295	LINC01116	4.0437	0	<0.0001
228438_at	100132891	MSC-AS1	3.1726	0	<0.0001
236304_at	100652988	LINC00702	2.6954	0	<0.0001
232298_at	401093	MBNL1-AS1	2.4114	0	<0.0001
230595_at	572558	PGM5-AS1	2.3164	0	<0.0001
226210_s_at	55384	MEG3	2.1848	0	<0.0001
235606_at	344595	DUBR	2.0825	0	<0.0001
235362_at	729970	LOC729970	2.0117	0	<0.0001
236118_at	100128893	GATA6-AS1	1.9728	0	<0.0001
235173_at	401093	MBNL1-AS1	1.9670	0	<0.0001
226211_at	55384	MEG3	1.9429	0	<0.0001
238081_at	404201	WDFY3-AS2	1.8997	0	<0.0001
236704_at	101929787	LOC101929787	1.8943	0	<0.0001
228601_at	401022	HAGLR	1.8205	0	<0.0001
238360_s_at	100505576	LINC00672	1.8155	0	<0.0001
239823_at	101927841	LOC101927841	1.8093	0	<0.0001

pfp, false discovery rates; LINC01279, long intergenic non-protein coding RNA 1279.

LINC01279 acted as a cell cycle mediator in endometriosis

To explore the function of LINC01279, we conducted a GBA analysis in the E-MTAB-694 and GSE29981 microarray datasets, which contained normal endometrial tissue and endometriosis samples. Firstly, the Pearson correlations between LINC01279 and all protein coding genes were computed. Then we select the top 500 positively correlated genes were enrolled in GO analyses. The result showed that the functions of these genes were enriched in many biological processes; in particular, the overlapping genes from all three datasets were consistently enriched in cell cycle and cell division (Fig. 2B and C). Next, GSEA was conducted for E-MTAB-694, GSE29981 and GSE7305 between the LINC01279 low and high expression groups. Overlapping the results showed that both the high and low LINC01279 expression groups were enriched in the cell cycle pathway category, including E-MTAB-694 [enrichment score (ES)=0.67], GSE7305 (ES=0.61) and GSE29981 (ES=0.59) (Fig. 3A-C) with c2.cp.kegg.v6.1.symbols.gmt motif. With c5.bp.v6.1.symbols.gmt motif, the result were also enriched in cell cycle (Fig. 3D-F). These results indicated that the lncRNA LINC01279 may act as a cell cycle mediator in endometriosis.

Figure 3.

LINC01279 acts as a cell cycle mediator in endometriosis. (A-C) The GSEA with c2.cp.kegg.v6.1.symbols.gmt motif of the high and low LINC01279 expression groups identified an enrichment in the cell cycle pathway category for the (A) E-MTAB-694 (ES=0.67), (B) GSE7305 (ES=0.61) and (C) GSE29981 (ES=0.59). (D-F) The GSEA with c5.bp.v6.1.symbols.gmt motif of the high and low LINC01279 expression groups identified an enrichment in the cell cycle biological process category for the (D) E-MTAB-694 (ES=0.49), (E) GSE7305 (ES=0.53) and (F) GSE29981 (ES=0.52). LINC01279, long intergenic non-protein coding RNA 1279; GSEA, gene set enrichment analysis; ES, enrichment score.

Possible mechanism for LINC01279 in promoting endometriosis development

In order to identify the mechanism for the promotion of the development of endometriosis, we calculated the correlation between LINC01279 and all protein coding genes in a merged dataset, including the data from GSE7305, GSE7846, GSE29981 and E-MTAB-694, with the inSilicoMerging method. The result showed that previously identified regulators of endometriosis, such as cyclin-dependent kinase 14 (CDK14) (Fig. 4A), CSCL12 (Fig. 4B), cytochrome P450 family 1 subfamily B member 1 (CYP1B1) (Fig. 4C), E26 transformation specific (ETS1) (Fig. 4D), insulin-like growth factor 1 (IGF1) (Fig. 4E), nitric oxide synthase 3 (NOS3) (Fig. 4F), secreted frizzled related protein 2 (SFRP2) (Fig. 4G) and slit guidance ligand 3 (SLIT3) (Fig. 4H), were positively associated with LINC01279. This demonstrated that LINC01279 may regulate the expression of these genes to promote the development of endometriosis.

Figure 4.

Correlation between the expression of LINC01279 and key factors in endometriosis. LINC01279 expression is positively correlated with (A) CDK14, (B) CXCL12, (C) CYP1B1, (D) EST1, (E) IGF1, (F) NOS3, (G) SFRP2 and (H) SLIT3. LINC01279, long intergenic non-protein coding RNA 1279; CDK14, cyclin-dependent kinase 14; CXCL12, C-X-C motif chemokine ligand 12; CYP1B1, cytochrome P450 family 1 subfamily B member 1; ETS1, E26 transformation specific; IGF1, insulin-like growth factor 1; NOS3, nitric oxide synthase 3; SFRP2, secreted frizzled related protein 2; SLIT3, slit guidance ligand 3.

Discussion

In recent years, with the development of microarray technology for the determination of SNPs, copy number, methylation status and mRNA expression, more and more large-scale genotyping analyses using microarrays have been performed to identify the genomic regions associated with endometriosis (25–28). Computational models have been widely used to compare the genotypes of healthy subjects and patients, thus identifying variants associated with a certain disease. lncRNAs, which may act as biological mediators in various diseases, have recently gained significant attention. This inspired us to compare the lncRNA expression profiles of normal endometrial tissue and endometriosis.

In our study, we found that LINC01279 and MSC-AS1 were consistently upregulated in three independent datasets. However, there was no overlap in downregulated lncRNAs in the same datasets. The RankProd method identified that LINC01279 was the most up-regulated lncRNA in endometriosis. To infer the function of this lncRNA in endometriosis, we conducted a GBA analysis, in which the co-expressed mRNAs were identified, and GO and GSEA analyses were performed to predict the function of LINC01279. This method may allow us to explore the functions of previously uncharacterized lncRNA.

As described above, from the GBA analysis, it was identified that LINC01279 may act as a cell cycle mediator. In order to further predict the possible mechanisms of LINC01279 in endometriosis, we conducted the Pearson correlation coefficient for all protein coding genes with LINC01279. Some previously identified key factors in the development of endometriosis were found to be strongly correlated with LINC01279, including the cell cycle regulating factor CDK14 (29). Endometriosis is conventionally viewed as a proliferative disorder characterized by the growth of endometrial cells in ectopic locations, and previous data have shown proliferation defects in the eutopic endometrium of women with endometriosis compared with women without endometriosis. For example, CSCL12 may promote the proliferation, migration and invasion of endometriosis (30). In addition, CYP1B1 (31), ETS1 (32), IGF1 (33), NOS3 (34), SFRP2 (35) and SLIT3 (36), which were previously identified to serve important roles in endometriosis, were positively correlated with LINC01279.

LINC01279 is 5,176 bp in length. In recent years, many studies have reported the mechanisms of how lncRNAs may regulate biological processes. For example, some lncRNAs may act as competing endogenous RNAs (ceRNAs) to ‘sponge’ microRNAs, thus affecting the expression of their target mRNAs. Yang et al (37) found that UCA1 could function as an endogenous sponge by directly binding miR-485-5p to regulate the expression of its target gene, matrix metallopeptidase 14, thus suppressing epithelial ovarian cancer metastasis. Sun et al (38) reported that HOXA11-AS acted as a ceRNA for miR-1297, thereby depressing EZH2 expression and imposing an additional level of post-transcriptional regulation in gastric cancer cells. Yu et al (39) reported that GAS5 functioned as a ceRNA for miR-222 to increase the level of p27 protein, thereby inhibiting the proliferation of hepatic stellate cells. Previous studies have also reported that some lncRNAs may scaffold the chromatin modification factors PRC2, LSD1 and DNMT1, or bind EZH2 to regulate cellular biological processes. For instance, the intronic lncRNA ANRASSF1 recruited PRC2 to the RASSF1A promoter, thus reducing the expression of RASSF1A and increasing cell proliferation. The lncRNA HOXA-AS2 has been reported to repress P21 and KLF2 expression transcription by binding with EZH2 and LSD1 in colorectal cancer (40). Yao et al (41) reported that the lncRNA ADAMTS9-AS2 was regulated by DNMT1, and inhibited migration in glioma cells. However, the mechanisms of how lncRNAs regulate biological processes are still unclear, and more research efforts are required.

To the best of our knowledge, this is the first study to date to report that the lncRNA LINC01279 acted as a cell cycle mediator in endometriosis, as determined by GBA analysis. Over recent decades, many CDK inhibitors have been used in treating cancer, such as the PARP inhibitor, olaparib, in ovarian cancer (42) and bortezomib in multiple myeloma (43). The results of this study indicate that LINC01279 may have potential as a new target for treating endometriosis. As described above, various mechanisms for the regulation of cellular functions by lncRNAs have been reported. Nevertheless, the exact mechanism of how LINC01279 affected the cell cycle in endometriosis is unclear and the in vivo and in vitro study should be done in further study.

Taken together, an lncRNA, LINC01279, was overexpressed in endometriosis. Furthermore, GO and GSEA analyses determined that LINC01279 potentially acted as a cell cycle mediator in endometriosis.