Exploring the Potential Key IncRNAs with Endometriosis by Construction of a ceRNA Network

Endometriosis is characterized by endometrial tissue (endometrial glandular and stromal) abnormally growing outside the uterine cavity. 1 The incidence rate in fertile women is up to 10%, and may be up to 50% in women who are experiencing infertility and/or pelvic pain. 2 , 3 Due to the high prevalence and the resulting physical and mental pain in fertile women, the etiology and pathophysiology of endometriosis needs further research. Long non-coding RNA (lncRNA) is a non-coding RNA varying in length from 200 nucleotides to 100 kilobases, which has been involved in a wide range of biological processes from regulation transcription to immune responses. 4–7 Many researchers have discovered that lncRNA may play an important role in regulated genes expression and the progression of endometriosis. 8–10 Huang et al. reported that lncRNA UCA1 was downregulated in endometriosis, 9 and Zhu et al. showed an effect of lncRNA BANCR on endometriosis through the ERK/MAPK pathway. 11 Therefore, more attention should be paid to lncRNAs in the study of progression of endometriosis. Current reports have demonstrated an intrinsic interaction among diverse RNA, referring to protein-coding messenger RNAs and non-coding RNAs such as long non-coding RNAs. 4 These RNA transcripts could behave as competing endogenous RNAs (ceRNAs) or natural microRNA sponges - they co-regulate each other by competing for binding to shared microRNAs, a family of small non-coding RNAs that are important post-transcriptional regulators of gene expression. 4 , 12 , 13 By studying this novel RNA crosstalk, new insight about gene regulatory networks and the association with human disease will be clear. Besides, our previous study has identified other regulatory networks involved in microRNA, transcription factors, and genes in order to try to understand the abnormal moleculer changes in the occurrence and development of endometriosis. 14 This study aims to construct the lncRNA-related ceRNA network and further survey the potential key lncRNAs for understanding the molecular mechanism of endometriosis by performing high-throughput sequencing techniques to detect lncRNA and mRNA expression profiling in eight paired ectopic endometria (EC) samples and eutopic endometria (EU) samples from eight patients with ovarian endometriosis.

A total of 1437 differential expression lncRNAs (DElncRNAs), 4225 differential expression mRNAs (DEmRNAs), 207 differential expression miRNAs (DEmiRNAs) were obtained by DESeq2. In total 742 lncRNAs were up-regulated and 695 lncRNAs were down-regulated; 2345 mRNAs were up-regulated and 1880 mRNAs were down-regulated; and 105 miRNAs were up-regulated and 102 miRNAs were down-regulated (see Tables S1–S6 ). The volcano plot of DEG (lncRNA, miRNA, mRNA) is shown in Figures S2–S4 . The heat map of DEG (lncRNA, miRNA, mRNA) to classify tissue and stage is shown in Figures S5–S7 . A total of 56,935 DElncRNA-DEmiRNA pairs and 178,004 DEmRNA-DEmiRNA pairs were attained by calculating the negative co-expression regulatory patterns of these DElncRNA-DEmiRNA interactions and DEmRNA-DEmiRNA interactions. Then, a sum of 832 interactions involving 353 lncRNAs and 94 miRNAs, and 21,533 interactions involving 2912 mRNAs and 195 miRNAs were attained by intersecting DElncRNA-DEmiRNA and DEmRNA-DEmiRNA pairs from negative co-expression regulatory patterns with pairs from five databases introduced in ‘Data collection’ section. Finally, the interaction with 94 in common miRNA pairs of lncRNA and mRNA were retained to construct the lncRNA-associated ceRNA network, which included 2811 DEmRNAs, 353 DElncRNAs, and 94 DEmiRNAs. As shown in Figure 1 , the lncRNA-associated ceRNA network was composed of two separate networks naturally. Especially, network 1 including 14,137 dysregulated lncRNA–mRNA interactions, referred to 242 lncRNAs, 55 miRNA and 1600 mRNAs, and network 2 including 4459 dysregulated lncRNA–mRNA interactions, referred to 111 lncRNAs, 39 miRNA and 1151 mRNAs (seen in Table 1 ). Table 1 Details of the lncRNA-Associated ceRNA Networks and Significant Ternary Relationships Network Type Number in Network 1 Number in Network 2 LncRNA-associated ceRNA network lncRNA–mRNA interactions 14,137 4459 lncRNA 242 111 miRNA 55 39 mRNA 1660 1151 Significant ternary relationship lncRNA–mRNA interactions 666 197 lncRNA 107 22 miRNA 40 22 mRNA 274 136 Figure 1 Constructed lncRNA-associated ceRNA network. Orange represents lncRNA, green represents miRNA, blue represents mRNA. Details of the lncRNA-Associated ceRNA Networks and Significant Ternary Relationships Constructed lncRNA-associated ceRNA network. Orange represents lncRNA, green represents miRNA, blue represents mRNA. Functional enrichment analysis was performed on the mRNAs in these two networks. Notably, mRNAs in network 1 were enriched in biological processes including myeloid leukocyte activation, lymphocyte activation, positive regulation of cytokine production and activation of immune response (seen in Figure 2A ), network 2 was mainly enriched in biological processes including tissue morphogenesis, embryonic organ development, epithelial cell development, and regulation of hormone levels (seen in Figure 2B ). The network and functional analysis based on co-dysregulated pairs whose CC were ranked in the top 0.05 and 0.1 percentile and p-value < 0.01 is shown in Figures S8–S13 . The functional enrichment analysis was consistent. Figure 2 Functional enrichment analysis of mRNA for two networks. ( A ) Functional enrichment analysis for network 1. ( B ) Functional enrichment analysis for network 2. Functional enrichment analysis of mRNA for two networks. ( A ) Functional enrichment analysis for network 1. ( B ) Functional enrichment analysis for network 2. Based-on the lncRNA-associated ceRNA network, we further screened the significant ternary relationship (lncRNA-miRNA-mRNA), which also involved two independent modules (shown in Figure 3A and B , Table 1 ). As shown in Figure 3C and D , the degree of lncRNA is significantly larger than mRNA in two modules, which means lncRNA plays a major role in regulation expression of mRNA, so according to the degree of lncRNA in the significant ternary relationship, the top six lncRNAs were selected to research their function by these enriched lncRNA related mRNAs. Figure 3 Network visualization of significant ternary relationships and degree of lncRNA and mRNA in significant ternary relationships. ( A ) Network visualization of significant ternary relationship for module 1. ( B ) Network visualization of significant ternary relationship for module 2. ( C ) Degree of lncRNA and mRNA in significant ternary relationship for module 1. ( D ) Degree of lncRNA and mRNA in significant ternary relationship for module 2. Network visualization of significant ternary relationships and degree of lncRNA and mRNA in significant ternary relationships. ( A ) Network visualization of significant ternary relationship for module 1. ( B ) Network visualization of significant ternary relationship for module 2. ( C ) Degree of lncRNA and mRNA in significant ternary relationship for module 1. ( D ) Degree of lncRNA and mRNA in significant ternary relationship for module 2. The top six degrees of lncRNAs in module 1 were LINC01140, MSC-AS1, HAGLR, CKMT2-AS1, JAKMIP2-AS1, and AL365361.1 (shown in Figure 4A ). Functional enrichment analysis of these lncRNA related mRNAs revealed that these lncRNAs take part in positive regulation of phagocytosis and myeloid leukocyte activation, cell adhesion mediator activity and tissue remodeling, etc. (shown in Figure 4B ). Figure 4 Subnetwork of top six degrees of lncRNA in significant ternary relationships and functional enrichment analysis with top six degrees of lncRNA related mRNAs. ( A ) Subnetwork of top six degrees of lncRNA in significant ternary relationship for module 1 ( B ) Functional enrichment analysis with top six degrees of lncRNA related mRNAs for module 1. ( C ) Subnetwork of top six degrees of lncRNA in significant ternary relationship for module 2. ( D ) Functional enrichment analysis with top six degrees of lncRNA related mRNAs for module 2. Subnetwork of top six degrees of lncRNA in significant ternary relationships and functional enrichment analysis with top six degrees of lncRNA related mRNAs. ( A ) Subnetwork of top six degrees of lncRNA in significant ternary relationship for module 1 ( B ) Functional enrichment analysis with top six degrees of lncRNA related mRNAs for module 1. ( C ) Subnetwork of top six degrees of lncRNA in significant ternary relationship for module 2. ( D ) Functional enrichment analysis with top six degrees of lncRNA related mRNAs for module 2. As for module 2, the top six degrees of lncRNAs were PAX8-AS1, MIR17HC, PART1, HOXA-AS3, PLAC4, and LINC00511 (shown in Figure 4C ). Functional enrichment analysis of these lncRNA related mRNAs revealed that these lncRNAs take part in negative regulation of cell proliferation, blood vessel development and regulation of epithelial cell differentiation, etc. This may account for the endometrium migration ( Figure 4D ).

This study was approved by the Chinese PLA General Hospital (No. S2016-123-01), and was conducted in accordance with the Declaration of Helsinki. All patients were informed about the purpose of this study, and each patient signed informed consent. Eight paired EC and EU samples were selected and prepared for RNA sequencing. The detailed information about patients can be seen in our previous work, 14 which were the same samples with miRNA sequencing. Endometrial tissue samples were obtained by pipelle sampling with inner membrane in the secretory phase including two endometrial gland cells and stromal cells. 15 All tissue samples were preserved in RNA-later solution (Sigma Aldrich, Poland) at 4°C for 24 h and then kept at −80  °C until utilized. RNA extraction, purification and RNA integrity analysis were performed as previously described in our previous work. 14 Three μg of RNA in each sample was prepared for RNA sequencing. Library construction and rRNA removal were performed by the Ribo-Zero Gold Kit (Epicentre, USA) as well as NEB Next Ultra RNA Library Prep Kit (NEB, Ipswich, USA) according to the manufacturer’s protocols. Then, an Illumina HiSeq4000 system (Illumina, USA) was used to sequence paired-end 150-bp of the cDNAs conducted by Annoroad Genome (Beijing, China). The raw data were processed to filter out adaptor-polluted reads, low-quality reads and reads with the number of N bases accounting for more than 5% with Perl scripts. Bowtie2 (v2.2.3) 16 was used as gapped read aligner, and TopHat (v2.0.12) 17 was used to align RNA-seq reads to the human genome (GRCh37) ( http://grch37.ensembl.org/index.html ). HTSeq (v0.6.0) was used to count the read counts of each gene, 18 and genes expression in each sample were estimated by the reads per kilobase of one gene per million reads (RPKM). The miRNA expression profiles were obtained from GSE105765 , which involved eight paired ectopic endometria (EC) samples and eutopic endometria (EU) samples from eight patients with ovarian endometriosis in our previous study. 14 miRNA-mRNA and miRNA-lncRNA interaction data were obtained from miRcode (291,978 miRNA-lncRNA interactions and 1,588,069 miRNA-mRNA interactions), 19 miRDB (5332 miRNA-lncRNA interactions and 3,224,437 miRNA-mRNA interactions), 20 miRTarBase (968 miRNA-lncRNA interactions and 477,771 miRNA-mRNA interactions), 21 RNAInter (24,142 miRNA-lncRNA interactions and 2,123,235 miRNA-mRNA interactions) 22 and StarBase (71,952 miRNA-lncRNA interactions and 1,075,989 miRNA-mRNA interactions). 23 In sum, a total of 377,739 miRNA-lncRNA interactions and 4,757,629 miRNA-mRNA interactions were collected. A total of 31,078 genes were involved in the study after experimental and data preprocessing, referring to 18,244 mRNAs, 11,325 lncRNAs, and 1509 miRNAs. DESeq2 24 was used to analyze differentially expressed genes (DEGs) with padj 1. Metascape ( http://metascape.org/ ) was performed for GO and KEGG pathway enrichment analysis. 25 Cytoscape software was used to visualize the lncRNA-related ceRNA networks and the lncRNA–mRNA interactions. 26 First, collecting all the miRNA-lncRNA and miRNA-mRNA interactions from five databases mentioned in Data collection. Second, calculating the negative co-expression regulatory pattern of miRNA-lncRNA and miRNA-mRNA by differential expression miRNA, lncRNA, and mRNA with Spearman correlation coefficient < −0.5 and p-value < 0.01. Also, the network and metascape results of the co-dysregulated pairs whose CC were ranked in the top 0.05 and 0.1 percentile and p-value < 0.01 were calculated to compare. Third, intersecting the miRNA-lncRNA and miRNA-mRNA interactions mentioned in step one and two above, then all the interactions with common miRNA make up the lncRNA-associated ceRNA network. We further selected the significant ternary relationship of lncRNA-miRNA-mRNA by calculating the score of any pair of the i th lncRNA and the j th mRNA, 27 which was evaluated by the Jaccard similarity coefficient, the formula is below. (1) \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document} $$scor{e_{i,j}} = {{{L_i}I{M_j}} \over {{L_i}U{M_j}}}$$ \end{document} Where L i represents the number of miRNA which interactwith i th lncRNA and M j represents the number of miRNA which interaction with the j th mRNA. And p-value of any pair of the i th lncRNA and the j th mRNA calculating by hypergeometric test as follows: (2) \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document} $$p = \mathop \sum \limits_{i = {m_c}}^{\min ({m_p},{m_n})} {{\left({\matrix{ {{m_n}} \cr i\cr } } \right)\left({\matrix{ {{M_T} - {m_n}} \cr {{m_p} - i} \cr } } \right)} \over {\left({\matrix{ {{M_T}} \cr {{m_p}} \cr } } \right)}}$$ \end{document} where M T represents all the number of miRNA, m p represents miRNA number which regulated the j th mRNA, m n represents miRNA number which regulated the i th lncRNA, m c represents miRNA number which regulated both the j th mRNA and the i th lncRNA. And significant lncRNA-miRNA-mRNA ternary relationships were determined by score > 0.5 and P < 0.01. The Flowchart of the whole work can be seen in Figure S1.

In summary, we provide some potential key lncRNAs in endometriosis by constructing a lncRNA-associated ceRNA network, and further screened the significant ternary relationship (DEmiRNAs-DEmRNAs-DElncRNAs) to select the top degrees of lncRNAs. Functional enrichment analysis on these lncRNA related mRNAs indicated that these lncRNAs take part in positive regulation of epithelial cell differentiation, cell adhesion mediator activity, myeloid leukocyte activation, negative regulation of cell proliferation, and blood vessel development. While our previous research identified global transcriptome abnormalities in endometriosis and defined some crucial miRNAs, TFs and genes involved in the pathogenesis of endometriosis, the key lncRNAs identified in this work will provide a comprehensive insight in the study of the molecular mechanisms involved in the development and proliferation of endometriosis, which could also have implications in therapeutic strategies for endometriosis.

The competitive endogenous RNA (ceRNA) hypothesis proposes that RNA transcripts, including circRNAs, messenger RNAs (mRNAs), and long non-coding RNAs (lncRNAs), contain miRNA response elements (MREs) which compete among themselves for miRNA binding to regulate the expression of each other. 12 , 28 , 29 Although there is still some debate about the clinical relevance of ceRNA, 30 many studies have adopted the ceRNA mechanism to explain the action of lncRNA 31 , 32 and their relationship with disease. 33 , 34 Meanwhile, in this study, we tested the difference among different co-expression thresholds, and the functional enrichment results of the network obtained under different thresholds were very similar (see Figures S8–S13 ). This suggests that the analysis pipeline of this study is robust to a certain extent. In this study, we try to identify key lncRNAs and explore their potential function in endometriosis patients by construction of a ceRNA network. The endometrial tissue samples are inner membrane in the secretory phase including two endometrial gland cells and stromal cells. A large number of existing studies 35 , 36 have shown that the occurrence and development of endometriosis is not only due to the role of glandular cells, but also the joint role of stromal cells, so we finally extracted glandular stromal cells, not specific to endometrial glandular cells to do further research. Firstly, lncRNA-associated ceRNA networks were constructed with combined DEmiRNAs, DEmRNAs and DElncRNAs with their interaction relationship of negative co-expression patterns and databases-sourced, which comprised two independent networks. Functional enrichment analysis found mRNAs in network 1 were enriched in immune-associated functions such as myeloid leukocyte activation, lymphocyte activation, and activation of immune response, and mRNAs in network 2 were enriched in development and differentiation associated functions such as gland development, cell-cell junction, cell division, and regulation of epithelial cell differentiation. Then, significant ternary relationships (DEmiRNA-DEmRNA-DElncRNA) were selected, and the top degrees of lncRNA in these were thought to be the key lncRNAs, and the mRNA interactions with these key lncRNAs were explored in depth. As a result, the top degrees of lncRNAs in the significant ternary relationship for network 1 were LINC01140, MSC-AS1, HAGLR, CKMT2-AS1, JAKMIP2-AS1 and AL365361.1 , and for network 2 were PAX8-AS1, MIR17HC, PART1, HOXA-AS3, PLAC4 and LINC00511. Functional enrichment analysis showed the mRNAs’ interaction with these top degree lncRNAs were enriched in positive regulation of phagocytosis, myeloid leukocyte activation, cell differentiation and tissue remodeling for module 1, and enriched in negative regulation of cell proliferation, blood vessel development and regulation of epithelial cell differentiation for module 2. Although lncRNA-associated research is limited, the key lncRNAs obtained above have been found to be related to some other diseases, for example, down-regulation of LINC01140 is associated with adverse features of breast cancer, 37 and significantly correlates with overall survival in gastric cancer patients. 38 High MSC-AS1 expression was correlated with poorer prognosis in patients with pancreatic ductal adenocarcinoma. 39 HAGLR acted as a microRNA-143-5p sponge to regulate epithelial-mesenchymal transition and metastatic potential in esophageal cancer. 40 Also, AL365361.1 was found to contribute to the modification of tumor immune microenvironment (TIME) in patients with head and neck squamous cell carcinoma. 41 Han et al. found expression quantitative trait loci in long non-coding RNA PAX8-AS1 are associated with decreased risk of cervical cancer. 42 Wang and Xu have found lncRNA PART1 promotes breast cancer cell progression by directly targeting miR-4516. 43 Also, HOXA-AS3 were found to integrate NF-κB signaling to regulate endothelium inflammation. 44 Especially, LINC00511 was found to play a major role in cancer in the reproductive system of women. LINC00511 accelerates proliferation and invasion in cervical cancer through targeting miR-324-5p/DRAM1 axis, 45 and inhibition of LINC00511 via transcription factor RXRA-regulated PLD1 can promote cell autophagy and apoptosis in cervical cancer. 46 LINC00511 can also mediate the effects of ESR1 on proliferation and invasion in ovarian cancer through miR-424-5p and miR-370-5p. 47 This means intensive study may focus on LINC00511 to reveal the pathway of reproductive system cancer and endometriosis in women. Previous studies have noticed some variation in endometriosis when researching the potential function of these key lncRNAs. Matthias and Menger have demonstrated the engraftment of endometriotic lesions is typically associated with angiogenesis, i.e. the formation of new blood vessels from pre-existing ones, and this angiogenic process underlies the complex regulation by angiogenic growth factors and hormones, 48 which supports our findings of key lncRNAs function in blood vessel development and regulation of hormone levels. Moreover, 63-kDa membrane protein (p63), which plays an important role in regulating epithelial proliferation and differentiation, were differentially expressed in different endometriotic lesions, 49 and Lin et al. explained integrin-mediated aberrant adhesion is also essential for establishment of endometriotic lesions, 50 both of which support our findings of lncRNAs function associated with negative regulation of cell proliferation and cell adhesion mediator activity. Also, some immune-associated genes were elucidated as related to the progression of endometriosis. Sun et al. demonstrated that M-MDSCs (the predominant type of myeloid-derived suppressor cells) recruited and activated by CCR9/CCL25 play a crucial role in the pathogenic progression of endometriosis, 51 which was consistent with our results about key lncRNA function of myeloid leukocyte activation. In our previous study, we identified global transcriptome abnormalities in eutopic endometria, 52 and further integrated microRNA and mRNA paired expression profiling to attain deregulated microRNA-transcription factor gene regulatory networks to find the crucial miRNAs, TFs and genes in ovarian endometriosis. 14 In the present study, key lncRNAs in endometriosis were surveyed by selecting top degrees of lncRNAs in a significant ternary relationship of the lncRNA-associated ceRNA network, and functional enrichment analysis was performed on the mRNAs interaction with these lncRNA to explore their functions. There is still some debate about the clinical relevance of ceRNA, which could be a potential limitation about the work. However, the present work still may provide novel and crucial insight for further study on molecular mechanisms involved in the development and proliferation of endometriosis.