Integrated Transcriptomic Analysis of the miRNA–mRNA Interaction Network in Thin Endometrium

Although the thin endometrium (TE) has been widely recognized as a critical factor in implantation failure, the contribution of miRNA–mRNA regulatory network to the development of disease etiology remains to be further elucidated. This study performed an integrative analysis of the miRNA–mRNA expression profiles in the thin and adjacent normal endometrium of eight patients with intrauterine adhesion to construct the transcriptomic regulatory networks. A total of 1,093 differentially expressed genes (DEGs) and 72 differentially expressed miRNAs (DEMs) were identified in the thin adhesive endometrium of the TE group compared with the control adjacent normal endometrial cells. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses showed that the DEGs and the target genes of DEM were significantly enriched in angiogenesis, cell growth regulation, and Wnt signaling pathway. Multiple hub genes (CAV1, MET, MAL2, has-mir-138, ARHGAP6, CLIC4, RRAS, AGFG1, has-mir-200, and has-mir-429) were identified by constructing the miRNA–mRNA regulatory networks. Furthermore, a miRNA–mRNA pathway function analysis was conducted, and the hub genes were enriched in the FoxO signaling pathway, cell growth regulation, inflammatory response regulation, and regulation of autophagy pathways. Our study is the first to perform integrated mRNA-seq and miRNA-seq analyses in the thin adhesive endometrium and the control adjacent normal endometrial cells. This study provides new insights into the molecular mechanisms underlying the formation of thin endometrium.


INTRODUCTION
The endometrium is an indispensable factor for implantation and pregnancy, and an increase in endometrial thickness promotes an increased pregnancy rate. An endometrial thickness of <7 mm is usually regarded as sub-optimal for embryo transfer and results in a decreased probability of pregnancy (Shufaro et al., 2008). For patients with Asherman's syndrome (AS), repetitive curettage or invasive endometritis disrupts endometrial regeneration, thus resulting in a fibrotic and thin endometrium (TE) (Azizi et al., 2018). Patients suffering from thin or fibrotic endometrium are more susceptible to abnormal menstruation and, particularly, fertility impairments, such as a decreased pregnancy rate, unfavorable pregnancy outcomes, or recurrent pregnancy loss (Du et al., 2020). The current AS treatments aim to increase endometrial regeneration with low-dose aspirin, exogenous estrogen, vitamin E, vaginal sildenafil citrate, cytokines, and colony-stimulating factors (CSFs). Nonetheless, these treatments are unable to attain a satisfactory clinical response in many patients with TE (Azizi et al., 2018). The definite etiology and physiopathology of thin endometrium remain largely unclear at present. Therefore, studies aiming to explore the related molecular mechanism of TE are urgently needed to guide disease therapy in the future.
A transcriptomic analysis is essential for understanding the occurrence and pathogenic mechanism of thin endometrium. Only one existing study has reported the global transcriptomic abnormalities in thin endometrium at the mid-luteal phase (Maekawa et al., 2017). The study compared the transcriptomic profiles between three patients and three normal subjects using the Gene Chip Human Genome U133 Plus 2.0 Array platform. Finally, 318 genes were upregulated in the thin endometrium, while 322 genes were downregulated. According to that study, implantation failure induced by thin endometrium might be related to the abnormal activation of the inflammatory environment, together with an abnormally reduced oxidative stress (OS) response. Nonetheless, researchers have not clearly elucidated the underlying mechanism of endometrial regeneration dysfunction in patients with thin endometrium. Additionally, more studies are needed to comprehensively characterize how thin endometrium affects the transcriptomic profiles.
MicroRNAs (miRNAs) are non-protein-coding RNA molecules with short (20-25) nucleotides. miRNAs bind to target mRNAs for transcription and translation regulation, including mRNA degradation, cleavage, or translational repression (Shukla et al., 2011;Li et al., 2019). miRNAs have been deemed to participate in the regulation of various cellular processes, including cellular proliferation, differentiation, apoptosis, and angiogenesis (Laurent, 2008;Nicoli et al., 2010;Hong et al., 2019). Recently, more and more miRNAs were found to be associated with endometrial receptivity (Altmae et al., 2013), endometrial stromal cell differentiation (Qian et al., 2009), embryo development (Laurent, 2008), and implantation (Paul et al., 2019). The expression of miR-27a-3p and miR-124-3p was downregulated in the endometrium of chronic endometritis (Di Pietro et al., 2018). The expression of hsa-miR-449a, hsa-miR-3135b, and hsa-miR-345-3p could promote endometrium receptivity in preparation for in vitro fertilization and embryo transplantation (Mu et al., 2020). miR-30 and miR-200 family members have been repeatedly recognized as important miRNAs in the regulation of endometrial receptivity (Rekker et al., 2018). Aberrant miR-200 expression may negatively regulate endometrial development and decidualization (Jimenez et al., 2016) and plays an important role in regulating normal endometrial development and disorders such as endometriosis and endometrial cancer (Panda et al., 2012). However, few studies have investigated the effect of miRNA on thin endometrium. The dysfunction of endometrium cells in TE and how miRNAs regulate the pathogenesis of TE remain to be elucidated.
Our article aimed to identify the miRNA-mRNA networks and molecular pathways in women experiencing intrauterine adhesion (IUA) and to provide additional insights into the underlying transcriptomic mechanisms by performing RNA-Seq. The differentially expressed miRNA-mRNA regulatory axis along with the gene pathway-function network interactions in thin endometrium was constructed. Our findings supply a basis to better investigate the biological mechanisms of thin endometrium and facilitate the formulation of molecular targeted treatments for thin endometrium.

Tissue Sample Collection
Eight females aged 20-40 years old, with a history of severe IUAs (Grade III-V) as diagnosed by hysteroscopy at the Reproductive Medicine Center of The First Affiliated Hospital of the University of Science and Technology of China, were enrolled in the study. The severity of IUAs was determined according to the American Fertility Society classification system (1988 version) (1988). Scores of 9-12 represented severe adhesions. The thickness of the endometrium was determined through vaginal ultrasound (at mid-luteal phase) as the maximum distance between endometrial interfaces, and the endometrial thickness in all patients was <7 mm. The sample information is described in Table 1. The endometrial tissue from the IUA (TE group) and adjacent normal endometrium tissues (AJ-CN group) from eight patients with severe IUAs (Grade III-V) were analyzed in the present study. This study was approved and monitored by the Human Research Ethics Committee of the First Affiliated Hospital of the University of Science and Technology of China. Each patient was required to provide a written informed consent prior to participation in this study. Endometrial tissues were sampled at mid-luteal phase during the menstrual cycle. Afterwards, the collected endometrial tissue samples were rinsed with saline to remove blood and then stored in liquid nitrogen at −80 • C until subsequent RNA isolation.

RNA Isolation and Library Construction
Total tissue RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), pooled equally, and reverse-transcribed into cDNAs using the QuantiTect Reverse Transcription Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's Frontiers in Genetics | www.frontiersin.org specific instructions. The quantity and quality of the extracted RNA were measured using Nanodrop (Thermo Scientific). The cDNA library was constructed using KAPA Stranded RNA-Seq Library Preparation Kit (Illumina) following the manufacturer's protocol. The synthetic cDNAs were end-repaired by polymerase and ligated with "A-tailing" base adaptors. Suitable fragments were selected for polymerase chain reaction (PCR) amplification to construct the final cDNA library. The final double-stranded cDNA samples were verified with Agilent 2100 Bioanalyzer (Agilent Technologies). Sequencing was performed on an Illumina HiSeq 4000 sequencing platform with 150-bp paired-end sequencing. Then, the combined RNA samples were separated using 15% (w/v) denaturing polyacrylamide gel electrophoresis. Subsequently, miRNA fragments with a size of ∼18-28 nt were separated by gel extraction, followed by RNA purification. The total RNA of each sample was used to prepare the miRNA sequencing library, which included the following steps: (1) 3 ′adaptor ligation, (2) 5 ′ -adaptor ligation, (3) cDNA synthesis, (4) PCR amplification, and (5) size selection of ∼135-155-bp PCR-amplified fragments (corresponding to ∼15-35 nt small RNAs). Libraries were quantified and validated with Agilent 2100 Bioanalyzer (Agilent Technologies). Thereafter, the small RNA library was sequenced using Illumina Hiseq 4000 (Illumina, San Diego, CA, USA), with a configuration of 50 cycles single reads according to the manufacturer's recommendations. All sequencing procedures were performed by Kang Chen Bio-tech (Shanghai, China).

mRNA Sequencing and Data Analysis
Raw data were pre-processed using Solexa CHASTITY and Cutadapt to remove adaptor sequences, ribosomal RNA, lowquality reads, and other contaminants that may interfere with assembly. The criteria for this filtering procedure were set as follows: (1) RNA 5 ′ and 3 ′ adapters were removed, respectively, (2) bases with a phred quality score below 20 were clipped from both ends of reads, (3) after low-quality bases were trimmed, reads containing over two "N" were discarded, (4) reads with a length shorter than 75 nt were discarded; and (5) the parameters for BWA v0.5.724 were set as recommended according to Fastq_clean instructions.
Then, the sequence quality was examined using FastQC v0.11.7. Afterwards, Hisat2 was utilized to align those trimmed reads to the reference genome. StringTie (version 1.2.3) was used to reconstruct the transcriptome. Fragments per kilobase per million (FPKM) values of genes were normalized with Ballgown using the default parameters. FPKM ≥0.5 (Cuffquant) was considered as statistically significant for the next DEG analysis. RNA sequencing data were deposited into the Gene Expression Omnibus (GEO accession number GSE160635).

miRNA Sequencing and Data Analysis
The miRNA sequencing data from TE group and AJ-CN group endometrium cells were analyzed by our previously published tool, DeAnnIso (Zhang et al., 2016). Briefly, after sequencing, Bowtie was used for mapping reads into the reference genome. The aligned reads had no more than "N" mismatches (0-3, default is 2) in the first "L" bases (≥5, default is 10) of the left end. Thereafter, those precursor sequence-matched reads were aligned to the pooled pre-miRNA databases (known pre-miRNAs in miRBase v21) using the BLAST. The default E-value was set to 0.01 for BLAST. All the detected isomiRs were aligned with their canonical miRNAs, the numbers of mapped reads that were defined as the raw expression levels of that miRNA. To correct for the difference in read counts between samples, the read counts were scaled to reads per million (RPM). Small RNA sequencing data were deposited into the Gene Expression Omnibus (GEO accession number GSE108966).

Differential Expression Analysis
After excluding the transcripts with a low count, genes with an FPKM or RPM ≥5 in at least one sample were used for the analysis. Fold change (FC) and P-value for Fisher's exact test was calculated and used when comparing the differentially expressed mRNAs (DEGs) and miRNAs (DEMs) between the two groups. The log 2 FC derived from the comparisons of the FPKM or RPM values of the TE group with the AJ-CN group is depicted (|Log 2 FC| ≥ 2) and P < 0.05 were selected as the cutoff criteria to identify significant DEMs and DEGs. Additionally, TargetScan (Garcia et al., 2011) and miRDB (Wang and El Naqa, 2008) were used to predict mRNAs targeted by DEMs.

Functional Annotations
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the online analysis tool of Annotation Visualization and Integrated Discovery (https://david.ncifcrf.gov/). The P-value for Fisher's exact test was calculated as a result of enrichment degree. GO term enrichment of biological processes or KEGG pathway annotations with a P-value cutoff of 0.05 were identified as an important term in this study.

Construction of the Protein-Protein Interaction Network
The Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/) was used to construct the PPI network. The obtained interactions included both the known and the estimated interactions. A requisite confidence value (pooled score >0.4) was used as the threshold. In addition, Cytoscape v3.7.1 was utilized to visualize the PPI network, and CytoHubba functions were employed to identify the hub genes. Genes with Gene significance >0.2, module membership >0.8, and P ≤ 0.05 were defined as hub genes.

Construction of the DEM-DEG Regulatory Network
TargetScan (http://www.targetscan.org/) and miRDB (http:// www.mirdb.org/) were utilized to preliminarily predict DEM target genes. The co-predicted targets were used for further GO and KEGG pathway enrichment analyses. The genes shared between DEM targets and DEGs were used to analyze the miRNA-mRNA pairs, which were maintained to construct the DEM-DEG regulatory network with Cytoscape. Differentially expressed target genes were chosen for GO and KEGG pathway analyses to investigate the miRNA-mRNA regulatory networks in TE.

Genome-Wide Patterns of the mRNA Transcriptomic Landscape
Using Illumina Hiseq 4000, 18,354,811 original RNA reads were obtained from the thin endometrial cells of patients with IUA, and 21,755,164 reads were obtained from adjacent normal endometrial cells. After removing adaptor sequences and lowquality reads, 18,288,140 (thin endometrial cells from patients with IUA) and 21,745,564 (adjacent normal endometrial cells) clean reads remained. Then, the genes were normalized to FPKM, and 15,561 genes were expressed in endometrial tissues from those eight women.
In the thin adhesive endometrial tissue of the TE group, 374 genes were upregulated, while 719 genes were downregulated compared to the control adjacent normal endometrial cells (Supplementary Tables 1, 2). The GO analysis of 1,093 DEGs identified many genes that were significantly enriched in the cell adhesion process (GO: 0007155, P = 1.92E-10), negative regulation of growth (GO: 0045926, P = 4.46E-05), angiogenesis (GO: 0001525, P = 4.63E-05), cell junction assembly (GO: 0034329, P = 2.99E-04), negative regulation of cell migration (GO: 0030336, P = 7.39E-04), the Wnt signaling pathway (GO: 0016055, P = 0.0017), and negative regulation of the BMP signaling pathway (GO: 0030514, P = 0.003) ( Table 2 and Figure 1). A blockade angiogenesis was considered as the FIGURE 2 | Kyoto Encyclopedia of Genes and Genomes pathway analysis of 1,093 differentially expressed genes between thin endometrium and adjacent normal endometrium. main pathological change in the scarred thin endometrium (Jiang et al., 2019). Moreover, this study identified several DEGs-related signaling pathways by performing KEGG pathway enrichment analysis, including the vascular smooth muscle contraction pathway, extracellular matrix-receptor interaction, focal adhesion, tight junction, cell adhesion molecules, calcium signal transduction pathway, p53 signal transduction pathway, and adherens junction pathway (Table 3 and Figure 2). The 1,093 DEGs were also compared with the primary associated changes identified in the transcriptome of the thin endometrium (Maekawa et al., 2017), and nine commonly upregulated genes (PDLIM3, FABP3, HIF3A, FILIP1, DPP6, MYOCD, PRKCB, ALDH1B1, and TRNP1) and 65 commonly downregulated genes were identified (Supplementary Table 3). The expression of MYOCD (myocardin), a cardiac-specific co-activator of serum response factor, was upregulated in thin endometrium, while ADAM12 (a disintegrin and metalloproteinase 12) expression was decreased in thin endometrium, and these genes are associated with the fibrosis process (Li et al., 2018;Mittal et al., 2019;Nakamura et al., 2019).

Genome-Wide Patterns of the miRNA Transcriptomic Landscape
First, clean reads were mapped to the human genome, and then those mapped reads were further matched to miRbase (V22). Notably, 7,004,583 reads (TE sample) and 5,717,874 reads (AJ-CN sample) were aligned to human pre-miRNAs. A total of 1,244 known miRNAs were altogether identified in our endometrial samples. According to the results of the miRNA-seq analysis, 72 known miRNAs were deemed to be DEMs between the thin adhesive endometrium of the IUA group and the control adjacent normal endometrial cells. Among these DEMs, five miRNAs were upregulated and 67 were downregulated compared with the control adjacent normal endometrial cells (Supplementary Table 4). The five upregulated and top 10 downregulated DEMs are shown in Table 4.
TargetScan and miRDB were used to characterize the putative target mRNAs of the 72 candidate DEMs in thin endometrium and to better illustrate the functions of DEMs. TargetScan and miRDB were employed to identify 812 common candidate target genes for the 15 DEMs (Supplementary Table 5). Then, GO and KEGG analyses were performed for the 812 target genes. GO enrichment analyses suggested that the target genes of multiple DEMs were associated with the regulation of angiogenesis, MAPK activation, negative regulation of cell migration, negative regulation of stress fiber assembly, positive regulation of epithelial cell proliferation, regulation of the canonical Wnt signaling pathway, and positive regulation of cell proliferation (Table 5 and Supplementary Figure 1). The KEGG pathways in which the DEM targeted genes are involved were discovered, which included the Ras signal transduction pathway, Hippo signal transduction pathway, MAPK signal  transduction pathway, PI3K-Akt signal transduction pathway, gap junction, p53 signaling pathway, Wnt signal transduction pathway, and ErbB signal transduction pathway (Figure 3 and Table 6). The PI3K/Akt signal transduction pathway is suggested to participate in endometrial regeneration induced by granulocyte macrophage-CSF therapy (Liu et al., 2020).

DEG-DEM Regulatory Network and Functional Assessment
For the establishment of the DEG-DEM regulatory network, 53 (21 upregulated genes and 32 downregulated genes) overlapping genes were discovered by comparing the target genes of DEMs (five were upregulated and 10 were downregulated) with DEGs, and they were deemed as consistently expressed genes (CEGs) (Figure 4). The STRING database was used to construct the PPI network using the CEG list. As shown in Figure 5, CAV1, MET, MAL2, has-mir-138, ARHGAP6, CLIC4, RRAS, AGFG1, has-mir-200, and has-mir-429 were the top 10 hub genes that interacted with the maximum number of nodes. Additionally, the gene pathway-function interactions were analyzed, and the identified hub genes showed significant enrichment in negative regulation of cell growth and inflammatory response regulation. For a better assessment of how this miRNA-mRNA regulatory network affected thin endometrium, a KEGG pathway analysis of CEGs was performed. The miRNA-mediated gene regulatory network in thin endometrium plays important roles in the regulation of the FoxO signaling pathway and the regulation of autophagy (Table 7).

DISCUSSION
IUA, which is characterized by endometrial fibrosis and thin endometrium, was always regarded as a major cause of female infertility and a major challenge to clinical therapy. Even through a surgical operation combined with hormone treatment, TE with severe endometrial injuries is difficult to restore. The previous transcriptomic microarray analysis discovered 318 upregulated genes and 322 downregulated genes in thin endometrium and revealed the abnormal activation of the inflammatory environment and an abnormal decrease in the OS response in thin endometrium (Maekawa et al., 2017). Current knowledge about the pathogenesis and involvement of miRNA-mRNA networks in thin endometrium is limited. In this study, gene expression patterns of thin endometrium along with the matched control endometrial tissues from women were explored, and we revealed the abnormal activation of the inflammatory environment and an abnormal decrease in the OS response in thin endometrium. To our knowledge, this study is the first to employ self-controlled transcriptomic analysis to investigate the regulatory functions of miRNA-mRNA networks of cells from the mid-secretory thin endometrium and adjacent normal endometrial cells.
As revealed in our results, some genes were abnormally expressed at the time of disease onset, revealing that thin endometrium may have occurred as a type of endometrial disorder due to the abnormal expression of genes within endometrial tissues prior to lesion occurrence. Indeed 1,093 genes were significantly differentially expressed in thin endometrium. A total of 74 DEGs associated with TE in our study were consistent with a previous study performed in thin and control endometrial samples using a microarray (Maekawa et al., 2017), including those that were up-regulated. Furthermore, our DEG functional enrichment analysis also revealed the involvement of angiogenesis and negative regulation of growth and cell migration in thin endometrium. Typically, during each menstrual cycle, angiogenesis promotes new blood vessel formation and is crucial for endometrial regeneration by supplying a vascularized and receptive endometrium for embryo implantation. Previous studies also show that the vascular endothelial growth factor (VEGF) could be a regulator of endometrial angiogenesis. Thus, the differential expression of VEGF and the blockade of angiogenesis in our study could be considered as pathological changes of the scarred thin endometrium (Jiang et al., 2019).
Interestingly, consistent with previous miRNA expression profiles reported for the recurrent implantation failure endometrium (Vilella et al., 2015;Rekker et al., 2018), some upregulated DEMs in TE in our study also belonged  (Panda et al., 2012). Through analyzing the interactions between DEMs and their targets, some vital pathways, including MAPK, p53, PI3K-Akt, and Wnt signal transduction, were found to participate in TE. As endometrial thickness has been recognized as an important indicator of endometrial receptivity (Ledee-Bataille et al., 2002), we thus assume that the abnormalities of these pathways may compromise the development of the endometrium. For example, rapid activation of PI3K/Akt signaling cascades by growth factors and estrogen is involved in the migration of normal endometrial stromal cells (Gentilini et al., 2007). However, the expression of DEGs in the PI3K/AKT pathway, including EFNA1, FGF9, LPAR1, CCNE2, SGK1, MET, IL6R, and PDGFRA, was decreased in thin endometrium, which suggests that the repair ability of the thin endometrium was impaired during the proliferative phase (Le et al., 2016). Similarly, the abnormal Wnt/beta-catenin signal pathway would also impair the proliferation of estrogendependent endometrial cells (Tepekoy et al., 2015).
In the present study, our miRNA-mRNA regulatory networks provided a complete profile for the underlying mechanism of thin endometrium formation, and the hub genes identified in the networks may play certain roles in the development of thin endometrium. CAV1 expression is associated with cell survival and proliferation (Zhao et al., 2013). MET, the receptor for insulin-like growth factor, potentially affects the functions of the endometrium (Satterfield et al., 2008). Therefore, the present study may provide useful information for understanding FIGURE 4 | The intersecting mRNAs between the common predicted target mRNAs and differentially expressed genes.
of the miRNA-mediated changes in mRNA expression in thin endometrium, and a further understanding of the functions of miRNA-mRNA networks can provide a new perspective for future studies examining potentially novel biomarkers and therapeutic targets. FIGURE 5 | Differentially expressed miRNA-differentially expressed gene regulatory network. The red and green colors denote upregulation and downregulation, respectively.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Materials.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committees on Human Research of the First Hospital Affiliated with University of Science and Technology of China. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
DL, ZW, and WT contributed to sample collection. XT and YM designed the experiment. SZ, LZ, and BX performed the experiment, data analysis, and manuscript preparation. XT and BX revised the manuscript. All authors contributed to the article and approved the submitted version.