Construction of implantation failure related lncRNA-mRNA network and identification of lncRNA biomarkers for predicting endometrial receptivity.

Embryo viability, endometrial receptivity, and embryo-maternal cross-talk are three determinants of successful embryo implantation. Preimplantation genetic screening (PGS) is used to exclude the embryos carrying aneuploidy to facilitate pregnancy, but the clinical pregnancy rate is not as satisfying as supposed. A retrospective cohort study in 2016 analyzed 5471 cycles with PGS and found that transfer of euploid embryos achieved a clinical pregnancy rate around 50% 1 . Another retrospective cohort study in US investigated 998 PGS cycles from 2011 to 2015 and found that the implantation rate was still lower than 60% after the genetic screening 2 . A study compared 88 recurrent implantation failure (IF) patients with PGS vs. without PGS, and the results showed that the clinical pregnancy rates were similar between these two groups 3 . All of these studies remind us that insufficient endometrial receptivity and failure in synchronization between embryo and endometrium is a leading cause of recurrent IF (RIF). Improving receptivity of endometrium can improve the outcome of in vitro fertilization and embryo transfer (IVF-ET). Since the traditional way to determine the condition of receptivity is simply by morphological observation of endometrium, more effective way for endometrial receptivity assessment is required. Long non-coding RNAs (lncRNAs) are a kind of RNAs that transcript longer than 200 nt nucleotides without translated proteins. It has been widely recognized that lncRNAs have plenty of functions in various endometrial biological processes. A series of lncRNAs, such as CASC2, LINC00672, TUG1, and TUSC7, was discovered to act in endometrial cancer as tumor suppressor and promote sensitivity to chemotherapy 4 - 7 . Overexpression of lncRNA H19 can stimulate the steroidogenic acute regulatory protein (StAR) so as to control the ovarian steroid production 8 . Comparison of gene expression profiles of eutopic and ectopic endometrium of adenomyosis and endometriosis showed that a bundle of lncRNAs were dysregulated in ectopic endometrium 9 , 10 . Furthermore, it has been reported that lncRNAs play an important role in the occurrence of IF. Animal study showed that the expression of a series of lncRNAs differed in pregnancy and non-pregnancy endometrium, which might play a vital role in facilitating embryo implantation 11 . The expression of H19 was down-regulated in the endometrium of RIF 12 . The expression of AK124742 was down-regulated together with PSMD6 in the cumulus cells of nonpregnant group compared with the pregnant group 13 . The above studies indicate that dysregulation of lncRNAs may affect the endometrial receptivity. The competing endogenous RNA (ceRNA) theory hypothesizes that lncRNAs can regulate other RNA transcripts by competing for shared miRNAs 14 . Since the introduce of this mechanism, it has been validated that ceRNA works in the development of various tumors, and ceRNA networks has been constructed in breast cancer 15 , pancreatic adenocarcinoma 16 , lung cancer 17 , hepatocellular carcinoma 18 , glioblastoma 19 , gastric cancer 20 and so on. Considering the role of lncRNAs in the regulation of endometrial receptivity and the function of ceRNA mechanism, we try to build the implantation failure related lncRNA-mRNA network (IFLMN) and to explore the key lncRNAs that can serve as biomarkers indicating endometrial receptivity.

As shown in Figure 1 , we constructed the IFLMN and then analyzed the functional characteristics of IFLMN. The key lncRNAs involved in IF were identified, their related modules were investigated, the functional enrichments were studied, and at last the expression of key lncRNAs was validated in endometrium. The predicted lncRNA-miRNA interact pairs were obtained from lncRNASNP2 database 21 , and the predicted miRNA-mRNA interact pairs were obtained from miRTarBase 22 . Then we constructed the global triple network as the background network to identify gene interactions. The gene expression data of GSE26787 was downloaded from Gene Expression Omnibus (GEO) database, which was based on the platform of GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array). It compared the endometrial expression in five patients who underwent RIF after IVF/ICSI and five fertile women 23 . Data were imported into the R-studio and normalized with RMA algorithm 24 . Bioconductor limma package was applied to identify the DEGs between fertile and IF groups with fold change (FC) > 2 and adjusted P < 0.05 25 . The DEGs were re-annotated by HG-U133_Plus_2 Annotations (CSV format, Release 36, 7/12/16) and identified as differentially expressed lncRNAs (DELs) and mRNAs (DEMs). Only the lncRNAs defined in NONCODE database were annotated as lncRNAs. All these DELs and DEMs were mapped into the global triple network, and then lncRNA-miRNA-mRNA triplets were extracted by hypergeometric test with P < 0.01 26 , and were selected as the lncRNA-mRNA ceRNA pairs and constructed the IFLMN. Unsupervised hierarchical clustering was performed by Multiple Experiment Viewer (MeV V4.9) software. The data were normalized and processed using hierarchical clustering with Pearson correlation as distance metric selection and average linkage method. The genes with similar expression profiles were grouped together and selected as the candidates for further co-expression sub-network. Hierarchical clustering was performed on IFLMN and on the lncRNA related mRNAs. Gene Ontology (GO) terms of biological process and molecular function were used to determine the function of target genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was applied to explore the significant pathways of DEGs. We performed GO classification and KEGG pathway analyses using Cytoscape V3.5.1 with the ClueGo V2.3.5 plug-in. The p-value was calculated by two-sided hypergeometric test and Benjamini-Hochberg adjustment. GO terms and KEGG pathways were identified as enriched in target genes when P ≤ 0.05. We analyzed all of the target genes in IFLMN and ceRNA sub-networks using GO and KEGG pathway analyses, and the results were visualized with Cytoscape. The topological analysis of DELs and DEMs was conducted to identify the central nodes of IFLMN. The topological parameters, including degree, betweenness, and closeness, were calculated by Cytoscape with CentiScaPe 2.2 plug-in. The top-ranked genes in degree, betweenness, and closeness were compared and the overlapped genes in these three topological parameters were selected as the key genes for further ceRNA analysis. To characterize the potential association between the key lncRNAs, the TFs that might regulate the lncRNAs were identified. The sequences of presumed promoters were downloaded from UCSC genome browser database. PROMO software was used to detect the predicted TFs with Maximum matrix dissimilarity rate < 3. The Venn diagram was used to obtain the overlapping TFs targeting all the key lncRNAs. The key lncRNAs and its first mRNA neighbors clustered together were extracted from the IFLMN. The miRNAs connecting the lncRNA and mRNA pairs were extracted from the global triple network and used to construct the lncRNA-miRNA-mRNA triple pairs. The ceRNA sub-networks were constructed by all these competing triplets and were visualized with Cytoscape software. The gene co-expression modules were extracted from IFLMN by bidirectional clustering. The Pearson correlation coefficient (PCC) between lncRNA and mRNA was calculated and lncRNA-mRNA pairs with strong correlations (PCC > 0.6, P < 0.05) were selected for further bidirectional clustering. There were 10 lncRNAs, 163 mRNAs, and 192 lncRNA-mRNA pairs involved for further bidirectional clustering. These genes were hierarchically clustered on the basis of Pearson correlation, and co-expression modules were explored by computing distances with an average linkage matrix. This study was approved by the Ethics Committee of the Second Affiliated Hospital, School of Medicine, Zhejiang University. All participants signed informed written consents. Participants were recruited from the Second Affiliated Hospital of Zhejiang University School of Medicine and the Women's Hospital of Zhejiang University School of Medicine. Eight patients that failed to get pregnant after ≥ 3 times of embryo transfers or who experienced ≥ 3 times of pregnancy losses during 6-12 weeks were recruited in RIF group. Eight patients that conceived after their first embryo transfer were included in control group. Exclusion criteria were as follows: hysteromyoma, adenomyosis, endometrial polyps, intrauterine adhesion, endometriosis, or hydrosalpinx. Endometrial specimen was collected at mid-luteal phase with a SAP-1 device (Saipujiuzhou, Beijing, China). The specimens were put in 300 μl RNAiso Plus and stored at -80℃ until use. Total RNA was extracted from endometrial tissue using RNAiso Plus (TaKaRa, Dalian, China, 9109) according to manufacturer's instructions. The concentration and purity of RNA were determined with spectrophotometer (Nanodrop 2000, Thermo Scientific). RNA samples were stored at -80℃. Total RNA with good purity and integrity was reverse transcribed using PrimeScript™ RT Master Mix (TaKaRa, Dalian, China, RR036A) on the Mastercycler (Eppendorf AG 22331 Hamburg, Germany) under the following conditions: 37℃ for 15 minutes, 85℃ for 5 seconds, and 4℃. cDNA samples were stored at -20℃. Real-time RT-PCR analysis was run on 7500 Fast Real-Time PCR System (Applied Biosystems) using a TB Green™ Premix Ex Taq™II kit (TaKaRa, Dalian, China, RR820A), according to the manufacturer's protocol. The primer sequences for real-time RT-PCR were listed in Table 1 . The following parameters were used: 95℃ for 30 seconds; 40 cycles of 95℃ for 5 seconds and 60℃ for 34 seconds. All the experiments were repeated in triplicate. The housekeeping gene ACTB was used as an endogenous reference. The primer sequences used in this study were listed in Table 1 . The relative mRNA expression levels were calculated using the 2 -△△Ct method 27 .

In all, there were 375647 lncRNA-miRNA interaction pairs, 19422 lncRNAs, and 2588 miRNAs obtained from the lncRNASNP2 database. 380633 miRNA-mRNA interaction pairs, 2599 miRNAs, and 15050 mRNAs were downloaded from miRTarBase. As shown in Figure S1 , lncRNA-miRNA pairs and miRNA-mRNA pairs were combined and constructed the global triple network, which included 756280 edges, 19422 lncRNAs, 2599 miRNAs, and 15050 mRNAs. This global network was used as a background network to construct the IFLMN. The gene expression data of GSE26787 was obtained from GEO database, comparing the endometrial expression between five RIF patients and five fertile women. 1251 DEGs were identified with FC > 2 and adjusted P < 0.05, including 990 genes up-regulated and 261 genes down-regulated in IF group. Only the lncRNAs defined in NONCODE database were annotated as lncRNA, and the DEGs were re-annotated and identified as 77 DELs and 1029 DEMs. All these DELs and DEMs were mapped into the global triple network, and then 1750 lncRNA-miRNA-mRNA triplets were extracted by hypergeometric test with P < 0.01. Finally, 255 lncRNA-mRNA pairs were selected as the lncRNA-mRNA ceRNA pairs and constructed the IFLMN, including 10 lncRNAs and 212 mRNAs, as shown in Figure 2 A. Enrichment analyses of GO classification and KEGG Pathway were performed on the DELs and DEMs in IFLMN. GO molecular function analysis found 55 GO terms enriched ( P ≤ 0.05 adjusted with Benjamini-Hochberg), as shown in Table 2 and Figure 2 B. The result of KEGG pathway revealed that these genes might be involved in the nine pathways listed in Table 3 and Figure 2 C ( P ≤ 0.05). Polymeric immunoglobulin receptor activity and interferon-gamma receptor binding were found to be enriched in IFLMN, with polymeric immunoglobulin receptor (PIGR) and interferon gamma (IFNG) (Figure 2 B). The recent transcriptional profile study on endometrial receptivity proved that immunological activity and vascular proliferation belong to the most over-expressed biological functions during window of implantation (WOI) 28 . The expression of molecules with immune function changes during WOI 29 . Animal study revealed that the decidual immune cell profile is consistently maintained during successive pregnancies and the decidua immune environment in pregnancy is precisely regulated 30 . IFNG is a proinflammatory cytokine secreted by uterine natural killer cells in maternal endometrium and by trophoblasts during early pregnancy. IFNG plays critical roles in initiating endometrial vasculature remodeling and angiogenesis at implantation sites 31 , 32 . Meanwhile, peptide-transporting ATPase activity was enriched in IFLMN, including TAP1 and TAP2. IFNG induces transcription of antigen peptide transporter (TAP1 and TAP2) and involve in the preparation for implantation 32 , 33 . Growth factor binding was enriched in IFLMN, including KAZALD1, NGFR, HAP1, ITGB3, NTRK3, and PCSK6 (Figure 2 B). Endometrial receptivity is affected by the interaction between embryo and endometrium via cytokines, growth factors, hormones and so on 34 . A series of growth factors, such as insulin like growth factor (IGF) 35 , vascular endothelial growth factor (VEGF) 36 , 37 , transforming growth factor (TGF) 38 , are involved in the endometrial receptivity. IGF1R expression decreases after ormeloxifene treatment and inhibits embryo implantation 35 . Steroid biosynthesis pathway was enriched in IFLMN (Figure 2 C). It is reported that decidualization of endometrial stromal cells is characterized by an increase in aromatase activity and activates the biosynthesis and metabolism of estrogen 39 . Estrogen plays a critical role in regulating endometrial natural killer and endothelial cells within the endometrium during early pregnancy 39 . The variation of steroid biosynthesis pathway may affect the regional concentration of steroid in endometrium and disturb the environment of embryo implantation 40 . These results indicate that endometrial steroid biosynthesis plays a pivotal role in the establishment of a receptive endometrium environment. VEGF signaling pathway was enriched in IFLMN (Figure 2 C). VEGF is expressed in the endometrium and changes cyclically 41 . It is involved in angiogenesis, an important phenomenon during embryo implantation 42 , 43 . During IVF/ICSI, VEGF expression is significantly higher in the women that became pregnant compared with those nonpregnant 44 . The increased progesterone level may affect the VEGF angiogenic factors of endometrium and decrease endometrial receptivity 45 . The 222 differentially expressed genes included in IFLMN were clustered. As shown in Figure 3 A, unsupervised hierarchical clustering of fertile and IF samples indicated a unique gene expression profile between fertile and IF groups. Meanwhile, 27 genes revealed a co-expression module, including 4 DELs and 23 DEMs, with a notably up-regulated expression in IF. These four DELs were NONHSAT083203.2, NONHSAT025064.2, NONHSAT129152.2, and NONHSAT193031.1. According to the sequences of promoters, the TFs of the four DELs were predicted on PROMO with Maximum matrix dissimilarity rate < 3. There were 25 common TFs appeared in these four lncRNAs, as shown in Figure 3 B. Among them, five TFs were involved in the sex steroid functions, and they were marked in brown in Figure 3 B, including C/EBPbeta, C/EBPalpha, ERalpha, PRA, and PRB. During pregnancy progesterone exerts a profound influence on uterus via PRA and PRB 46 . The expression of CCAAT/enhancer binding protein beta (C/EBPbeta) is regulated by progesterone in stroma and by estrogen in uterine epithelium 47 . PRA and PRB collaborate with C/EBPbeta to regulate the expression of many target genes to control uterine epithelial proliferation, stromal differentiation, angiogenesis, and local immune response to render the endometrial receptivity and allow embryo implantation 46 , 49 . Absence of C/EBPbeta leads to a block of stromal cell proliferation after decidual stimulation, causing a complete failure of the implantation 50 . Loss of appropriate ERalpha down-regulation in mid-secretory phase may lead to some types of uterine receptivity defects 51 . In all, sex steroids, TFs, lncRNAs and down-stream effectors constitute a network regulating endometrial receptivity and prepare endometrium for implantation. GO biological process analysis with P ≤ 0.01 adjusted with Benjamini-Hochberg found 13 GO terms enriched, as shown in Figure 3 C, including positive regulation of interleukin-1 alpha biosynthetic process, elastic fiber assembly, mRNA pseudouridine synthesis, regulation of microtubule motor activity, positive regulation of axon regeneration, odontogenesis of dentin-containing tooth, odontoblast differentiation, thymocyte migration, diphosphoinositol polyphosphate catabolic process, diadenosine pentaphosphate catabolic process, diadenosine hexaphosphate catabolic process, and adenosine 5'-(hexahydrogen pentaphosphate) catabolic process. The result of KEGG pathway with P ≤ 0.07 revealed that these genes might be involved in the five pathways listed in Figure 3 D, including ECM-receptor interaction, Ribosome biogenesis in eukaryotes, IL-17 signaling pathway, TNF signaling pathway, and Rheumatoid arthritis. Positive regulation of interleukin-1 alpha biosynthetic process was enriched in this module (Figure 3 C). During endometrial stromal cell decidualization there is a significant modulation of endometrial cell receptivity to IL1 family members and IL1R, which suggests the involvement of IL1 family members in physiological processes of endometrial receptivity 54 , 55 . Extracellular matrix (ECM)-receptor interaction pathway was enriched in this module (Figure 3 D). Both human and animal studies indicate that not only DEGs between RIF and conception after IVF/ICSI, but also DEGs between receptive endometrium and pre-receptive phase endometrium, are mainly participated in ECM-receptor interaction pathway 56 , 57 . The topological characteristics of IFLMN were analyzed with degree, betweenness, and closeness. The genes in IFLMN with the highest degree, betweenness, or closeness were listed in Table 4 . There are four genes existed simultaneously on all the three top 10 lists, that is, four lncRNAs, as shown in Figure 4 A. The top four genes with the highest degree coincide with these four lncRNAs, and they are NONHSAT212577.1, NONHSAT083203.2, NONHSAT193031.1, and NONHSAT035952.2 (Figure 4 B). For NONHSAT212577.1, 57 mRNAs were extracted from the IFLMN as the first mRNA neighbors. GO classification with P ≤0.05 adjusted by Benjamini-Hochberg found three enriched GO terms for Molecular function and one GO term for Biological process, as shown in Figure 4 C and 4 D, and Table S1 . KEGG analysis with P ≤ 0.05 found six enriched KEGG pathways, as shown in Figure 4 E and Table S2 . Molecular function was enriched in TAP1 binding, MHC class Ib protein binding, and peptide-transporting ATPase activity. Biological process was enriched in cytosol to ER transport. Both of these two pointed to TAP1 and TAP2 (Figure 4 C and 4 D), indicating transcription of antigen peptide transporters TAP1 and TAP2 were crucial in endometrial receptivity. What's more, KEGG pathway analysis also pointed to TAP1 and TAP2 (Figure 4 E), confirming their importance and emphasized the role of immunity in endometrial receptivity. It is reported that TAP1 and TAP2 are up-regulated in pregnant animals 58 . TAP proteins may up-regulate VEGF-C and stimulate the capacity of natural killer cells to be non-cytotoxic, so TAP1 and TAP2 are involved in local immune suppression during pregnancy 59 . Silencing of TAP1 abolishes the VEGF-C participated immune protection 60 . IFNG significantly regulates TAP1, chemokines, enzymes, and TFs, induces uNK cell production of the chemokines, direct migration and invasion, and promote physiological angiogenesis in the placental bed 61 . The above indicates that TAP1 and TAP2 play a role in implantation. Clustering of NONHSAT212577.1 and the 57 first neighbor mRNAs, the genes clustered together with the lncRNA were extracted. Then the linked miRNAs were added, and a new sub-network was constructed. In all, 21 mRNAs, 31 miRNAs, and 82 edges were included in the ceRNA sub-network of NONHSAT212577.1, as shown in Figure 4 F and 4 G. For NONHSAT083203.2, 52 mRNAs were extracted from the IFLMN as the first mRNA neighbors. GO classification (Biological process) with P ≤ 0.05 adjusted by Benjamini-Hochberg found 11 enriched GO terms, as shown in Figure 5 A and Table S3 . KEGG analysis with P ≤ 0.05 found 8 enriched KEGG pathways, as shown in Figure 5 B and Table S4 . Clustering of NONHSAT083203.2 and the 52 first neighbor mRNAs, the bundle of genes clustered together with the lncRNA were extract. Then the linked miRNAs were added, and a new sub-network was constructed. In all, 9 mRNAs, 20 miRNAs, and 105 edges were included in the ceRNA sub-network of NONHSAT083203.2, as shown in Figure 5 C and 5 D. For NONHSAT193031.1, 37 mRNAs were extracted from the IFLMN as the first mRNA neighbors. GO classification (Biological process) with P ≤ 0.05 obtained 128 enriched GO terms, as shown in Figure 6 A and Table S5 . KEGG analysis with P ≤ 0.05 found 8 enriched KEGG pathways, as shown in Figure 6 B and Table S6 . Clustering of NONHSAT193031.1 and the 37 first neighbor mRNAs, the bundle of genes clustered together with the lncRNA were extract. Then the linked miRNAs were added, and a new sub-network was constructed. In all, 8 mRNAs, 11 miRNAs, and 27 edges were included in the ceRNA sub-network of NONHSAT193031.1, as shown in Figure 6 C and 6 D. For NONHSAT035952.2, 30 mRNAs were extracted from the IFLMN as the first mRNA neighbors. GO classification (Biological process) with P ≤ 0.05 obtained 43 enriched GO terms, as shown in Figure 7 A and Table S7 . KEGG analysis with P ≤ 0.05 found 4 enriched KEGG pathways, as shown in Figure 7 B and Table S8 . Clustering of NONHSAT035952.2 and the 30 first neighbor mRNAs, the bundle of genes clustered together with the lncRNA were extract, then the linked miRNAs were added, and a new sub-network was constructed. In all, 17 mRNAs, 5 miRNAs, and 39 edges were included in the ceRNA sub-network of NONHSAT035952.2, as shown in Figure 7 C and 7 D. The lncRNAs and mRNAs included in the IFLMN was further selected with Pearson correlation analysis. The lncRNA-mRNA pairs with correlation coefficient > 0.632 and P < 0.05 was used for bidirectional clustering. There were 10 lncRNAs, 163 mRNAs, and 192 lncRNA-mRNA pairs selected for further bidirectional clustering. As shown in Figure 8 A, two co-expression modules were found. Module 1 included 3 lncRNAs and 29 mRNAs, and then 50 miRNAs were extracted from the global triple network. At last, 246 edges were included in the module 1 ceRNA sub-network (Figure 8 B, Table S9 ). Module 2 included 2 lncRNAs and 13 mRNAs. Six miRNAs were extracted from the global triple network. At last, 28 edges were included in the module 2 ceRNA sub-network (Figure 8 C, Table S9 ). For Module 1, GO classification (Molecular function) with P ≤ 0.05 obtained 35 enriched GO terms, as shown in Figure 8 D. KEGG analysis with P ≤ 0.05 found three enriched KEGG pathways, as shown in Figure 8 E. Phosphatidylinositol 3 kinase (PI3K)-Akt signaling pathway was enriched in Module 1 (Figure 8 E). It is reported that embryo implantation is disturbed when PI3K-Akt signaling pathway is blocked 62 . Cell motility and invasion are crucial events for endometrial remodeling during embryo implantation. Steroid hormones and growth factors can activate PI3K-Akt signaling pathway and stimulate migration of endometrial stromal cells 63 . TNXB was involved in PI3K-Akt pathway. As shown in Figure 8 B, hsa-miR-1301-3p, as link between NONHSAT083203.2 and TNXB, may work as sponge of this interaction to involve in endometrial receptivity. Apoptosis signaling pathway was also enriched in Module 1 (Figure 8 E). The proliferation and apoptosis of the endometrium is the foundation of the preparation for the WOI 64 . High estradiol concentration can induce the apoptosis of endometrial glandular cells and impairs endometrial receptivity, leading to adverse reproductive outcome in IVF cycles 65 , which explains the reason why pregnancy rate of frozen cycle is better than that of fresh cycle. NGFR was involved in apoptosis pathway. As shown in Figure 8 B, hsa-miR-6785-5p, hsa-miR-6883-5p, hsa-miR-4728-5p, hsa-miR-149-3p, hsa-miR-4747-5p, hsa-miR-5196-5p, hsa-miR-7106-5p, and hsa-miR-3202 linked NONHSAT083203.2/NONHSAT025064.2 and NGFR, and they may work as sponges of these interactions to involve in endometrial receptivity. For Module 2, GO classification (Molecular function) with P ≤ 0.05 obtained 12 enriched GO terms, as shown in Figure 8 F. KEGG analysis with P ≤ 0.05 found six enriched KEGG pathways, as shown in Figure 8 G. Steroid biosynthesis pathway was enriched in Module 2 (Figure 8 G). As we emphasized above, steroid biosynthesis pathway may play an important role in endometrial receptivity. MSMO1 was involved in steroid biosynthesis pathway. As shown in Figure 8 C, hsa-miR-6880-5p, hsa-miR-650, and hsa-miR-3612, as links between NONHSAT053761.2 /NONHSAT035952.2 and MSMO1, may work as sponges of these interactions to involve in endometrial receptivity. There were no significant differences between RIF/RM group and control group in age, BMI (body mass index), endometrial pattern, and endometrial thickness ( P > 0.05). The identified candidate lncRNAs were assayed by real-time PCR on a sample set of 8 RIF/RM patients and 8 controls. As shown in Figure 9 , the expression levels of NONHSAT083203.2, NONHSAT212577.1, NONHSAT035952.2, NONHSAT193031.1, NONHSAT053761.2, and NONHSAT025064.2 in RIF/RM patients were significantly higher than those in controls ( P < 0.01). The increased amplitudes range from 4.17 to 38.68 folds, which is obviously a substantial elevation. Our previous bioinformatics analysis discovered that these six lncRNAs were upregulated in IF patients and might predict endometrial receptivity. Furthermore, we investigated the expression of these lncRNAs in endometrium and validated their upregulation in RIF patients. It is reported that in patients with better endometrial receptivity, the expression of ST7-OT3 was down-regulated, which supported out results 66 . Considering their role in IFLMN and their functional enrichment, it will be meaningful to explore their upstream regulation mechanism and their downstream effector mechanism. A gene expression microarray study of endometrium in 2016 identified a signature containing 303 genes predictive of RIF. Further validation in independent samples showed that the gene signature predicted RIF with 100% positive predictive value (PPV). Among them, there were 11 lncRNAs identified in NONCODE database, including NONHSAT007100.2, NONMMUT115606.1, NONHSAT076715.2, NONHSAT030758.2, NONHSAT169502.1, NONHSAT133392.2, NONHSAT084833.2, NONHSAT198142.1, NONHSAT131185.2, NONHSAT073810.2, NONHSAT039454.2 68 . A lncRNA and mRNA microarray study of endometrium in 2017 detected 197 dysregulated lncRNAs in RIF endometrium compared with the control group. Real-time PCR validated the expression of ENSG00000267194.1 (p7482), XLOC_005748 (p23230), XLOC_004399 (p22979), ENSG00000224259.1 (p1238), ENSG00000261645.1 (p2900), ENSG00000267549.1 p8859), and ENSG00000229385.1 (p9212) were elevated and the expression of ENSG00000259716.1 (p10433) was decreased 67 . These results indicate that lncRNA expression of RIF endometrium was different from that of normal receptive endometrium, suggesting that lncRNAs may regulate endometrial receptivity. However, the key lncRNAs in the present study were not in accordance with the lncRNAs identified in these two studies. The reason might be as follows: ethnic disparities, difference microarray platforms, and different bioinformatics methodology.

In the present study, we constructed the IFLMN, identified six key lncRNAs in IFLMN through topological analysis, clustering and co-expression analyses, and validated their increased expression in endometrium of RIF/RM patients. NONHSAT083203.2 (CECR3), NONHSAT212577.1 (ST7-OT3), NONHSAT035952.2 (DHRS4-AS1), NONHSAT193031.1 (C22orf34), NONHSAT053761.2 (RAMP2-AS1), and NONHSAT025064.2 (PNCT_HSA157732), together with their ceRNA sub-networks, may be involved in immunological activity, growth factor binding, vascular proliferation, apoptosis, and steroid biosynthesis in uterus and prepare endometrium for embryo implantation. These lncRNAs might be regulated by transcription factors that involved in sex steroid functions, such as PRA and PRB. Since some lncRNAs maintain stable in peripheral blood and are protected from endogenous RNases, lncRNAs can function as biomarkers for detection and diagnosis. These six lncRNA emerged as promising candidates to predict endometrial receptivity since their role in IFLMN and their functional enrichment. Further studies are needed to validate the above lncRNAs in uterine luminal fluid and peripheral blood, and they may be developed into predictive biomarkers for endometrial receptivity.

Supplementary figure and tables. Click here for additional data file.