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

Insufficient endometrial receptivity is a major factor leading to implantation failure (IF), and the traditional way of morphological observation of endometrium cannot determine the condition of receptivity sufficiently. Considering that long-noncoding RNAs (lncRNAs) regulate endometrial receptivity and competing endogenous RNA (ceRNA) mechanism works in plenty of biological processes, ceRNA is likely to function in the pathology of IF. In the present study, we aim to construct an implantation failure related lncRNA-mRNA network (IFLMN), and to identify the key lncRNAs as the candidates for predicting endometrial receptivity. The global background network was constructed based on the presumed lncRNA-miRNA and miRNA-mRNA pairs obtained from lncRNASNP and miRTarBase. Differentially expressed genes (DEGs) of IF were calculated using the data of GSE26787, and then re-annotated as differentially expressed mRNAs (DEMs) and lncRNAs (DELs). IFLMN was constructed by hypergeometric test, including 255 lncRNA-mRNA pairs, 10 lncRNAs, and 212 mRNAs. Topological analysis determined the key lncRNAs with the highest centroid. Functional enrichment analyses were performed by unsupervised clustering, GO classification, KEGG pathway, and co-expression module analyses, achieving six key lncRNAs and their ceRNA sub-networks, which were involved in immunological activity, growth factor binding, vascular proliferation, apoptosis, and steroid biosynthesis in uterus and prepared endometrium for embryo implantation. Sixteen endometrial samples were collected during mid-luteal phase, including 8 recurrent implantation failure (RIF) or recurrent miscarriage (RM) women and 8 controls who conceived successfully. Quantitative real-time PCR was performed to compare the expression of the above six lncRNAs, which validated that the expression of all these lncRNAs was significantly elevated in endometrium of RIF/RM patients. Further studies are needed to investigate the underlying mechanism, and the lncRNAs may be developed into predictive biomarkers for endometrial receptivity.

Introduction

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.

Materials and Methods

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.

Figure 1

Workflow of the present study. First, the lncRNA-miRNA interact pairs were obtained from lncRNASNP, and the miRNA-mRNA interact pairs were obtained from miRTarBase. The pairs shared the same miRNA were combined to make up lncRNA-miRNA-mRNA triple pairs and the global triple network. Second, gene expression data were downloaded from Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified by Bioconductor limma package with fold change > 2 and adjusted P < 0.05. The DEGs were re-annotated and identified as differentially expressed mRNAs (DEMs) and lncRNAs (DELs). Third, IFLMN were extracted from the global network by hypergeometric test with P < 0.01. Fourth, the functional enrichment was mined by unsupervised clustering, GO, and KEGG annotation. Fifth, the topological features of IFLMN were analyzed and four key lncRNAs were selected. Sixth, construction of the sub-network with key lncRNAs and enrichment analysis. Seventh, extracting the co-expression modules with bi-directional clustering and enrichment analysis. Eighth, validation of the expression of key lncRNAs in endometrium with real-time PCR.

Construction of the global lncRNA-miRNA-mRNA triple network

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.

Screening of differentially expressed genes (DEGs) and re-annotation

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.

Construction of the IFLMN

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.

Hierarchical clustering

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.

Enrichment analysis, including GO and KEGG pathway analysis

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.

Topological analysis and selection of key genes

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.

Prediction of transcription factors (TFs)

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.

Construction of ceRNA sub-networks

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.

Bidirectional clustering

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.

Sample collection

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.

RNA Extraction and Reverse Transcription

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℃.

Quantitative Real-Time PCR

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.

Table 1

Primer sequences for real-time RT-PCR.

LncRNAs	Genes	Primer forward	Primer reverse	Product size
NONHSAT083203.2	CECR3	5'-CCTCCTGGCAGAAAAGTTCCC-3'	5'-AAGGCAGAAGGACCAGTAGG-3'	150bp
NONHSAT212577.1	ST7-OT3	5'-GGCCTCTAAGAAACTTCGTGT-3'	5'-TGTCCTTCTTCAGGCTAGGG-3'	110bp
NONHSAT035952.2	DHRS4-AS1	5'-ATCGATAGGGCTGCACAGAG-3'	5'-TGTCCTCAGACCTTACACGC-3'	75bp
NONHSAT193031.1	C22orf34	5'-ATGCTCACAGGGATGTCACG-3'	5'-AACCCATTCCAAGCGGATCA-3'	165bp
NONHSAT053761.2	RAMP2-AS1	5'-GATCTATTTCGGGGGCTCGG-3'	5'-CGCTTCTCCTCCCTGGTTTT-3'	144bp
NONHSAT025064.2	PNCT_HSA157732	5'-TAGCCAGCTGAATTGCCACA-3'	5'-CGCTGGAGATCGTTCAGGAAA-3'	75bp
	ACTB	5'-CGCGAGAAGATGACCCAGAT-3'	5'-ACAGCCTGGATAGCAACGTA-3'	71bp

For PCR, the primers were designed using Primer-blast (https://www.ncbi.nlm.nih.gov/tools/primer-blast). All PCR reactions were repeated in triplicate for both target genes and internal controls.

Results and Discussion

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.

Construction of the implantation failure related lncRNA-mRNA network

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 2A.

Figure 2

Implantation failure related lncRNA-mRNA network (IFLMN) and enrichment analysis. A. The blue concave quadrilateral represents lncRNA and the green diamond represents mRNA. Grey edges indicate lncRNA-mRNA interactions. Size of the node represents the degrees of the nodes in IFLMN. The lncRNAs marked with purple labels are the four most central lncRNAs. B. GO analysis of IFLMN. C. KEGG pathway analysis of IFLMN. The dotted lines emphasize the annotations that have been proved to work in endometrial receptivity.

The functional characteristics of IFLMN

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 2B. The result of KEGG pathway revealed that these genes might be involved in the nine pathways listed in Table 3 and Figure 2C (P ≤ 0.05).

Table 2

Enriched GO terms in IFLMN.

GO Term	P-Value	Enrichment score
C-4 methylsterol oxidase activity	2.14E-02	100.00
polymeric immunoglobulin receptor activity	2.14E-02	100.00
low-affinity phosphate transmembrane transporter activity	2.14E-02	100.00
adenosine kinase activity	2.14E-02	100.00
S-adenosylmethionine-dependent tRNA (m5U54) methyltransferase activity	2.14E-02	100.00
oligopeptidase activity	2.14E-02	100.00
25-hydroxycholecalciferol-24-hydroxylase activity	2.14E-02	100.00
1-alpha,25-dihydroxyvitamin D3 24-hydroxylase activity	2.14E-02	100.00
triplex DNA binding	2.14E-02	100.00
cytosolic dipeptidase activity	2.14E-02	100.00
alanylglutamate dipeptidase activity	2.14E-02	100.00
dynorphin receptor activity	2.14E-02	100.00
stearoyl-CoA 9-desaturase activity	6.78E-04	66.67
tryptophan-tRNA ligase activity	4.26E-02	50.00
hydroxymethylglutaryl-CoA synthase activity	4.26E-02	50.00
prostaglandin-endoperoxide synthase activity	4.26E-02	50.00
interferon-gamma receptor binding	4.26E-02	50.00
glycine hydroxymethyltransferase activity	4.26E-02	50.00
L-allo-threonine aldolase activity	4.26E-02	50.00
mannose-6-phosphate isomerase activity	4.26E-02	50.00
N-acetyl-beta-glucosaminyl-glycoprotein 4-beta-N-acetylgalactosaminyltransferase activity	4.26E-02	50.00
adenylylsulfate kinase activity	4.26E-02	50.00
sulfate adenylyltransferase (ATP) activity	4.26E-02	50.00
adipokinetic hormone receptor activity	4.26E-02	50.00
phosphatidylinositol-3,5-bisphosphate 5-phosphatase activity	2.23E-03	40.00
neurotrophin binding	7.89E-06	36.36
neuroligin family protein binding	3.32E-03	33.33
TAP1 binding	3.32E-03	33.33
MHC class Ib protein binding	3.32E-03	33.33
peptide-transporting ATPase activity	3.32E-03	33.33
nerve growth factor binding	3.32E-03	33.33
cadmium ion binding	1.93E-04	33.33
NADPH-hemoprotein reductase activity	9.68E-03	20.00
calcium:cation antiporter activity	1.40E-02	16.67
co-SMAD binding	1.64E-02	15.38
ephrin receptor activity	3.45E-02	10.53
FMN binding	3.81E-02	10.00
acetylgalactosaminyltransferase activity	1.23E-02	8.57
4 iron, 4 sulfur cluster binding	2.31E-02	6.82
calcium-dependent phospholipid binding	3.09E-02	6.12
Rab guanyl-nucleotide exchange factor activity	4.01E-02	5.56
intramolecular oxidoreductase activity	4.21E-02	5.45
oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen	2.07E-04	5.20
growth factor binding	6.89E-03	4.41
serine-type endopeptidase inhibitor activity	4.49E-02	4.00
iron ion binding	5.78E-03	3.98
divalent inorganic cation transmembrane transporter activity	7.17E-03	3.83
glycosaminoglycan binding	4.47E-03	3.72
organic acid binding	4.06E-02	3.02
substrate-specific transmembrane transporter activity	3.06E-02	1.87
transition metal ion binding	7.56E-03	1.84
cation binding	3.21E-05	1.66
metal ion binding	5.65E-05	1.66
anion binding	8.11E-03	1.57
G-protein coupled receptor activity	2.19E-02	0.33

Table 3

Enriched KEGG Pathways in IFLMN.

GO Term	P-Value	Enrichment score
Steroid biosynthesis	2.39E-02	10.00
Biosynthesis of unsaturated fatty acids	3.10E-02	8.70
Glyoxylate and dicarboxylate metabolism	4.46E-02	7.14
Mineral absorption	2.47E-02	5.77
Basal cell carcinoma	2.85E-02	5.45
VEGF signaling pathway	3.42E-02	5.08
Axon guidance	1.89E-02	3.43
Phagosome	3.83E-02	3.25
Proteoglycans in cancer	3.58E-02	2.96

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 2B). 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 2B). 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 2C). 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 2C). 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.

A co-expression module in IFLMN and the putative transfer factors

The 222 differentially expressed genes included in IFLMN were clustered. As shown in Figure 3A, 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.

Figure 3

Custer analysis of implantation failure related lncRNA-mRNA network (IFLMN) and its co-expression network. A. The clustering of the IFLMN genes. Unsupervised clustering of IFLMN genes could discriminate between fertile and IF groups. B. The Venn of predicted transcription factors (TFs) of the four co-expression lncRNAs. There are 25 common TFs and TFs involved in the sexual hormone functions were marked in brown. C. GO enrichment analysis of co-expression mRNAs. D. KEGG pathway analysis of co-expression mRNAs. The dotted lines emphasize the annotations that have been proved to work in endometrial receptivity.

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 3B. Among them, five TFs were involved in the sex steroid functions, and they were marked in brown in Figure 3B, 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 3C, 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 3D, 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 3C). 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 3D). 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.

Topological characteristics of IFLMN

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 4A. 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 4B).

Table 4

The 10 top genes with highest degrees, betweenness, and closeness.

N	Genes	Degree	Genes	Betweenness	Genes	Closeness
1	NONHSAT212577.1	57	NONHSAT212577.1	20696.65	NONHSAT083203.2	0.001773
2	NONHSAT083203.2	52	NONHSAT083203.2	19838.70	NONHSAT212577.1	0.001650
3	NONHSAT193031.1	37	NONHSAT193031.1	14863.68	NONHSAT035952.2	0.001634
4	NONHSAT035952.2	30	NONHSAT035952.2	12343.36	NONHSAT193031.1	0.001524
5	NONHSAT192631.1	24	NONHSAT192631.1	8979.926	PAPOLA	0.001511
6	NONHSAT025064.2	22	NONHSAT025064.2	8438.824	CD44	0.001493
7	NONHSAT053761.2	14	SCRG1	5681.178	PIGR	0.001462
8	NONHSAT128926.2	8	PIGR	4668.291	SCRG1	0.001445
9	NONHSAT129152.2	8	CD44	4418.550	EPHB1	0.001425
10	MSMO1	3	NONHSAT053761.2	4236.610	SLC25A24	0.001425

Figure 4

Topological characteristics of implantation failure related lncRNA-mRNA network (IFLMN) and ceRNA sub-network of NONHSAT212577.1. A. Topological characteristics of IFLMN. There are four genes existed simultaneously in top 10 lists of degree, betweenness, and closeness. B. Distribution of degree. The top four genes with highest degree were four lncRNAs. C. Molecular function enrichment of NONHSAT212577.1. D. Biological process enrichment of NONHSAT212577.1. E. KEGG enrichment of NONHSAT212577.1. F. The co-expression subnet of NONHSAT212577.1. The dotted lines emphasize the co-expressed genes. G. The ceRNA network of NONHSAT212577.1. The dotted lines emphasize the annotations that have been proved to work in endometrial receptivity.

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 4C and 4D, and Table S1. KEGG analysis with P ≤ 0.05 found six enriched KEGG pathways, as shown in Figure 4E 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 4C and 4D), 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 4E), 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 4F and 4G.

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 5A and Table S3. KEGG analysis with P ≤ 0.05 found 8 enriched KEGG pathways, as shown in Figure 5B and Table S4.

Figure 5

The ceRNA sub-network of NONHSAT083203.2. A. Biological process enrichment of NONHSAT083203.2. B. KEGG enrichment of NONHSAT083203.2. C. The co-expression subnet of NONHSAT083203.2. The dotted lines emphasize the co-expressed genes. D. The ceRNA network of NONHSAT083203.2.

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 5C and 5D.

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 6A and Table S5. KEGG analysis with P ≤ 0.05 found 8 enriched KEGG pathways, as shown in Figure 6B and Table S6.

Figure 6

The ceRNA sub-network of NONHSAT193031.1. A. Biological process enrichment of NONHSAT193031.1. B. KEGG enrichment of NONHSAT193031.1. C. The co-expression subnet of NONHSAT193031.1. The dotted lines emphasize the co-expressed genes. D. The ceRNA network of NONHSAT193031.1.

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 6C and 6D.

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 7A and Table S7. KEGG analysis with P ≤ 0.05 found 4 enriched KEGG pathways, as shown in Figure 7B and Table S8.

Figure 7

The ceRNA sub-network of NONHSAT035952.2. A. Biological process enrichment of NONHSAT035952.2. B. KEGG enrichment of NONHSAT035952.2. C. The co-expression subnet of NONHSAT035952.2. The dotted lines emphasize the co-expressed genes. D. The ceRNA network of NONHSAT035952.2.

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 7C and 7D.

Co-expression modules

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 8A, 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 8B, 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 8C, Table S9).

Figure 8

The co-expression networks in implantation failure related lncRNA-mRNA network (IFLMN). A. The bidirectional clustering of IFLMN. The genes with green bar were included in Module 1 and the genes with blue bar were included in Module 2. B. Module 1 of ceRNA networks. Module 1 included 3 lncRNAs and 29 mRNAs. C. Module 2 of ceRNA networks. Module 2 included 2 lncRNAs and 13 mRNAs. D. GO classification of Module 1. E. KEGG pathway of Module 1. F. GO classification of Module 2. G. KEGG pathway of Module 2.

For Module 1, GO classification (Molecular function) with P ≤ 0.05 obtained 35 enriched GO terms, as shown in Figure 8D. KEGG analysis with P ≤ 0.05 found three enriched KEGG pathways, as shown in Figure 8E.

Phosphatidylinositol 3 kinase (PI3K)-Akt signaling pathway was enriched in Module 1 (Figure 8E). 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 8B, 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 8E). 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 8B, 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 8F. KEGG analysis with P ≤ 0.05 found six enriched KEGG pathways, as shown in Figure 8G.

Steroid biosynthesis pathway was enriched in Module 2 (Figure 8G). 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 8C, 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.

Validation of lncRNAs expression of endometrium

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.

Figure 9

Expression of lncRNAs in endometrium. The six key lncRNAs were assayed by real-time PCR on a sample set of 8 RIF patients and 8 controls. The expression levels of lncRNAs in RIF increased significantly compared to fertile controls. * P < 0.05.

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.

Conclusion

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.