Integrated Analysis of Hub Genes and MicroRNAs in Human Placental Tissues from In Vitro Fertilization-Embryo Transfer

Objective Supraphysiological hormone exposure, in vitro culture and embryo transfer throughout the in vitro fertilization-embryo transfer (IVF-ET) procedures may affect placental development. The present study aimed to identify differences in genomic expression profiles between IVF-ET and naturally conceived placentals and to use this as a basis for understanding the underlying effects of IVF-ET on placental function. Methods Full-term human placental tissues were subjected to next-generation sequencing to determine differentially expressed miRNAs (DEmiRs) and genes (DEGs) between uncomplicated IVF-ET assisted and naturally conceived pregnancies. Gene ontology (GO) enrichment analysis and transcription factor enrichment analysis were used for DEmiRs. MiRNA-mRNA interaction and protein-protein interaction (PPI) networks were constructed. In addition, hub genes were obtained by using the STRING database and Cytoscape. DEGs were analyzed using GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Differentially expressed miRNAs were validated through qRT-PCR. Results Compared against natural pregnancies, 12 DEmiRs and 258 DEGs were identified in IVF-ET placental tissues. In a validation cohort, it was confirmed that hsa-miR-204-5p, hsa-miR-1269a, and hsa-miR-941 were downregulation, while hsa-miR-4286, hsa-miR-31-5p and hsa-miR-125b-5p were upregulation in IVF-ET placentas. Functional analysis suggested that these differentially expressed genes were significantly enriched in angiogenesis, pregnancy, PI3K-Akt and Ras signaling pathways. The miRNA-mRNA regulatory network revealed the contribution of 10 miRNAs and 109 mRNAs while EGFR was the most highly connected gene among ten hub genes in the PPI network. Conclusion Even in uncomplicated IVF-ET pregnancies, differences exist in the placental transcriptome relative to natural pregnancies. Many of the differentially expressed genes in IVF-ET are involved in essential placental functions, and moreover, they provide a ready resource of molecular markers to assess the association between placental function and safety in IVF-ET offspring.


INTRODUCTION
In vitro fertilization-embryo transfer (IVF-ET) is the main method of assisted reproductive technology (ART), and can improve the success rates of infertility treatment (1,2). The use of this procedure has been steadily increasing worldwide and for example, in China, the average number of ART treatments performed now exceeds 700,000 annually (3). The vast majority of children successfully delivered through ART are healthy (4). Nonetheless, there are increased risks of several pregnancyrelated complications, including gestational diabetes (5), preeclampsia (6), placenta previa (7), abnormal placental growth, preterm delivery (8), and low birth weight (9). Some reports have suggested that the adverse perinatal outcomes occur due to IVF-ET procedures such as supraphysiological estrogen level during stimulation, in vitro culture, and microscopic manipulation (10,11). However, the underlying causes of IVF-ET-associated complications are largely unknown, although many conceptually appear related to placental vascular complications and the resulting effects on fetal development.
MicroRNAs (miRNAs) represent a major class of non-coding RNAs consisting of single-stranded RNAs of approximately 18-25 nucleotides in length (12). They function as negative gene regulators by binding to the 3' UTR of messenger RNAs (mRNAs), to either prevent protein translation or to direct the mRNA towards degradation (13). The placenta expresses many ubiquitous as well as specific miRNAs and a growing body of evidence proposes that miRNAs function as important regulators of placental development (14). Here, various miRNAs have been shown to control the differentiation, replication, apoptosis, invasion/migration and angiogenesis of trophoblasts, indicating the widespread contribution of miRNAs to the placental growth (15,16). For example, miR-346 and miR-582-3p down-regulate the expression of endocrine gland-derived vascular endothelial growth factor (EG-VEGF) and inhibit trophoblast cell invasion and migration (17). miR-191 inhibits angiogenesis by activating the nuclear factor-kB (NF-kB) signaling pathway (18). It was also found that miR-29b inhibited trophoblast invasion and angiogenesis by suppressing vascular endothelial growth factor (VEGF) expression (19). In addition, animal experiments in mice revealed that miR-450a-3p played a role in inhibiting cell proliferation, promoting apoptosis and interfering with embryonic development through regulating the target gene Bub1 (20). However, the impact of specific microRNAs in IVF-ET placental tissue, and their potential impact on related gene regulatory networks has to date been poorly investigated.
This study aimed to understand the differences in placental genomic expression profiles comparing IVF-ET and natural pregnancy-derived placentas, and the link between IVF-ET manipulation and placental structure and function. We used next-generation sequencing techniques to analyze of miRNA and mRNA expression profiles in IVF-ET and natural gestational placental tissue. Based on co-expression analysis and online prediction. we established a miRNA-mRNA regulatory network comprising 10 miRNAs and 109 mRNAs together with a PPI network comprised of ten hub genes. Furthermore, we verified that miR-204-5p, miR-1269a and miR-941 were downregulated in IVF-ET placentals thereby proposing these as key regulators involved in the effects of IVF-ET on placental development and function. Moreover, our study constitutes a verified resource for enabling further investigation into the transcriptomic and mechanistic differences between IVF-ET and naturally conceived pregnancies with the goal of providing new targets to assess the safety of IVF-ET in the clinic.

Tissue Collection and Ethics
Placental tissue samples were collected from women who underwent caesarean deliveries after IVF-ET assisted (n=3) or natural conceived (n=3) pregnancies. Inclusion criteria included full-term singleton delivery after IVF-ET, with age between 20 and 35 years, 37-42 gestational weeks, infant birth weight between 2500 g and 4000 g, and uncomplicated pregnancies. Three strictly matched natural pregnancies were selected as controls with matching parameters: delivery, maternal age, parity, and gestational duration ( Table 1). A validation cohort of 8 uncomplicated IVF-ET and 8 normal conception patients were similarly collected ( Table 2). We used an equal number of male and female placentals for both discovery and validation cohorts and selected tissue from the middle placenta throughout to minimize sampling bias. All placental tissues were rinsed extensively with ice cold PBS and stored at -80°C until later RNA extraction. Tissue collection was approved by the Ethics Committee of the Third Affiliated Hospital of Zhengzhou University with written informed consent provided by all patients prior to sample collection. the manufacturer's instructions. RNA concentration and integrity were then verified using an Agilent Bioanalyzer 2100 (Agilent technologies Santa Clara, US) before subjecting the samples to library preparation and sequencing. The concentration and size of the constructed libraries were measured using a Qubit ® 2.0 Fluorometer (Life Technologies, USA) and Agilent 2100 Bioanalyzer, respectively. The samples were prepared according to the HiSeq 2500 User Guide, and the flow cell with the cluster was loaded on the Illumina HiSeq 2500 (50-bp single-end FASTQ reads). The number of sequencing reads per sample was at least 10M and the proportion of bases with mass greater than 20 was greater than 95%, and the quality control meets the requirements of data analysis.

Differential Expression Analysis
The DESeq2 package of R version 3.5.2 (http://www.r-project. org) was used to define differentially expressed miRNAs (DEmiRs) and genes (DEGs) with the Bayesian method used to correct batch effects. MiRNAs and mRNAs with statistical significance between the IVF-ET and control groups were selected according to the threshold criterion of fold change (FC) >1.5 and P < 0.05. Volcano maps were created in the R studio using the plot packages to illustrate the differential expression of DEmiRs and DEGs.

GO Enrichment Analysis for the Targets of Transcription Factors
DEmiRs were uploaded to FunRich software to screen for upstream transcription factors, which is primarily used for functional enrichment and interaction network analysis of genes and proteins, as well as enrichment targets for transcription factor pathways (21). For interaction network analysis between miRNAs, gene/mRNA, and transcription factors, gene ontogeny (GO) enrichment analysis was also used (22,23).

Construction of the miRNA-mRNA Regulatory Networks
The miRWalk V2.0, StarBase and TargetScan databases were used to predict the target mRNAs of the miRNAs identified as DEmiRs. Subsequently, the predicted DEGs were matched with the experimentally determined DEGs to develop miRNA-mRNA regulatory networks with visualization using Cytoscape software (http://www.cytoscape.org/) (24). All node degrees, proximity and presence of the regulatory network were simultaneously computed.

GO and KEGG Enrichment Analyses
For gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs, we used the Database for Annotation, Visualization and Integrated Discovery (DAVID) and used ggplot2 packages in the R studio to identify significantly altered biological processes (BPs), cellular components (CCs), molecular functions (MFs) and pathways associated with DEGs (P < 0.05).

Protein-Protein Interaction (PPI) Network Analysis and Hub Gene Identification
Differentially expressed gene data was uploaded to the STRING database (http://www.string-db.org/) (25). Interactions with a composite score of >0.4 were considered significant. The target genes in the PPI network act as nodes and the line from two nodes indicates relevant interactions. Cytoscape software was used to visualize the PPI network. We screened the top 10 genes with the highest degree of correlation to the others as hub genes with the CytoHubba plugin of Cytoscape (26).  One mg total RNA was reverse transcribed into cDNA using the ReverTra Ace qPCR RT Kit (Toyobo, Japan) according to the manufacturers' instructions. Specific primers were used to synthesize the cDNA of miRNAs. qRT-PCR reactions were performed with the indicated primers (

Statistical Analysis
Data were presented as means ± standard deviations (SD) for quantitative variables and the Student's t-test used to assess differences between groups. Otherwise, for discrete variables the Mann-Whitney U-test was used. A value of P < 0.05 was regarded as statistically significant. SPSS 21.0 and GraphPad Prism 6.0 were used for analysis.

Verification of the DEmiRs by qRT-PCR in Placental Tissues
To ensure the veracity of the high throughput sequencing analysis, it was necessary to validate our findings using alternative methodology and samples. On the basis of subsequent bioinformatics (see below), we selected 6 of the 12 DEmiRs (3 downregulated and 3 upregulated, respectively) and analyzed their expression using qRT-PCR in a validation cohort of 8 IVF-ET and 8 normal conception placentas. Instructively, we found that the relative expression levels of hsa-miR-204-5p, hsa-miR-1269a, and hsa-miR-941 were significantly downregulated whereas hsa-miR-4286, hsa-miR-31-5p and hsa-miR-125b-5p were all upregulated, respectively, in IVF-ET compared to the control placentas ( Figures 3A-F). These results suggest that the dysregulation of these miRNAs commonly occurs in IVF-ET pregnancies but further verification will be required to support the general veracity of the global sequencing data and analysis of the study.

Transcription Factor Enrichment and GO Enrichment Analysis
Transcription factors often represent the critical final step in signal transduction pathways. To investigate the enrichment of transcription factor targets likely associated with the genetic landscape of the IVF-ET placenta, we filtered out the top 10 transcription factors most closely linked to the DEmiRs. In deceasing probability order, this analysis identified SP1, EGR1, SP4, KLF7, ONECUT1, MYF5, TCF3, NFIC, SRF, and HOXB4 ( Figure 4A). It indicates a regulatory interaction between these transcription factors and DEmiRs. Intriguingly, SP1 was apparently able to regulate majority of the DEmiRs.
In concert with these findings, GO enrichment analysis indicated that the top 5 biological progress (BP) terms with the most enriched targets of the DEmiRs involved signal transduction; regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism, transport, apoptosis, and regulation of immune response ( Figure 4B). GO enrichment terms associated with molecular function (MF) indicated most of the genes were involved in transcription factor activity, transporter activity, RNA binding, Receptor signaling complex scaffold activity, and GTPase activity ( Figure 4C). The top 5 enriched cellular component (CC) terms were nucleus, cytoplasm, Golgi apparatus, lysosome and membrane raft ( Figure 4D).

Construction of miRNA-mRNA Regulatory Networks
Reliable identification of miRNA targets is still an imprecise process, but it is widely appreciated that predictions can be improved using the outputs of multiple algorithms (28).
Consequently, we employed three databases (miRWalk V2.0, StarBase and TargetScan) to analyze the potential impact of the DEmiRs on the placental transcriptome. Screening miRNA targets based on the overlapping results of the three databases and the intersection with DEGs. This analysis yielded paired interactions between 109 DEGs with 10 of the 12 identified DEmiRs. The network of miRNA-mRNA interactions was visualized in Cytoscape ( Figure 5) and the target genes of the DEmiRs are listed in Table 5. Notably, among these, the four downregulated DEmiRs, particularly hsa-miR-204-5p, hsa-miR-1269a and hsa-miR-941, formed the most extensive interactive network with multiple gene targets while the upregulated DEmiRs aligned with a more discrete set of target genes.

Functional Enrichment Analysis of the DEGs
Independent of the preceding analysis, we utilized ggplot2 and enrichment analysis to profile GO annotations and KEGG pathways associated with the differentially expressed genes. These predictions would allow an improved understanding of the functional impact of the genes dysregulated in IVF-ET. As illustrated, the top 10 enriched GO and KEGG terms are presented in categories of biological processes (BPs), cellular components (CCs), molecular functions (MFs) and defined KEGG pathway identifiers ( Figure 6).
Notably, this analysis revealed significant enrichment for BP entries aligned with angiogenesis, pregnancy, cell adhesion, positive regulation of transcription from the RNA polymerase II promoter and positive regulation of angiogenesis ( Figure 6A). Furthermore, the protein binding, poly (A) RNA binding and transcriptional activator activity accounted for the majority of MF terms ( Figure 6B) while the most enriched CCs were extracellular exosome, membrane, and extracellular region ( Figure 6C). The top 10 most highly enriched KEGG classifications included the PI3K-Akt and Ras signaling pathways along with focal adhesion ( Figure 6D).

Construction of Protein-Protein Interaction (PPI) Networks
Using the DEGs we next built a PPI network using the online STRING database and tools in Cytoscape. The network was mapped to a limit of 148 DEGs (Figure 7). Node size is proportional to the degree of the node itself. Edge width is proportional to the combined degree between genes. We assessed their degree of connectivity and identified 10 hub genes, and the genes are listed in Table 6.

Biological Analysis of the Hub Genes
Highly-connected genes with a network are considered master regulatory elements otherwise known as hub genes (29). We used the cytoHubba plugin of Cytoscape to reveal the ten most strongly related interactions amongst the DEGs in IVF-ET. This approach generated 10 nodes with 37 edges with the most likely hub genes consisting of EGFR, FOS, SERPINE1, LEP, HGF, EGR1, SPP1, HNRNPA2B1, IGF2 and ENG ( Figure 8A). In addition, KEGG analysis of the top 10 enrichment pathways were identified (Figures 8B, C).

DISCUSSION
As one of the main tools of ART, IVF-ET can improve the success rate of infertility treatments for many affected couples.   Nonetheless, despite its use for several decades and widespread acceptance, there are still uncertainties associated with IVF-ET in terms of perinatal risks and offspring health. On this basis we hypothesized that complications associated with IVF-ET may be associated with the altered expression of genes that regulate placental development and function. Consequently, we performed genome-wide miRNA and mRNA analyses comparing placentas from IVF-ET assisted and naturally conceived pregnancies. The choice of samples from uncomplicated births was deliberate, as large gene differences from pathological births would be expected while more fundamental changes could be revealed by the experimental design used. Foremost we considered the epigenetic regulatory mechanisms involving miRNAs as these are recognized effectors of placental-mediated complications during pregnancy (15). Indeed, there has been strong interest in miRNAs as predictive biomarkers for the detection of pathologies in pregnancy with the villous trophoblast being a major source of miRNAs found in maternal circulation (30). Such changes in circulating miRNAs would naturally reflect the changes occurring within the placenta. Our analysis produced a  -941  Down  28  RBM25, MPHOSPH8, SOCS1, SEC62, HP1BP3, JUP, NFE2L1, THY1, FLT1, GOLGA4, TRIQK, C6orf89, NFYC, PPP1R13L,  ARHGAP26, IGF2, BCL6, LIFR, ANGPT2, FOSB, SCIN, PHACTR2, RRBP1, ARL5A, MME, ERBB3, UBE2V1, ENG  hsa-miR-4286  Up  2  RBP4, HMSD  hsa-miR-31-5p  Up  5  GCSH, UCK2, HMSD, MXRA7, RPS4Y1  hsa-miR-224-5p  Down  5  SYNCRIP, ARHGAP26, ARL5A, SERPINE1, SLC7A8  hsa-miR-204-5p  Down  30  PAPPA2, WAC, SUB1, PROCR, SLCO2A1, HAPLN1, HP1BP3, LYNX1, NFE2L1, THY1, ACP2, CBR1, COBLL1, SLC7A8,  CREB3L2, ARID5B, SRSF11, TRIQK, SYNCRIP, ARHGAP26, BCL6, XPO7, SLC38A1, ATP6V1C2, ZFAT, NDEL1, SH3BP5,  RAP1B, SERPINE1, CYP19A1  hsa-miR-193b-5p  Up  5  UCK2, DNASE1L3 shortlist of 12 miRNAs that were differentially expressed in IVF-ET compared to placental tissue from natural pregnancies. Three of the 12 differentially expressed miRNAs were confirmed in an independent cohort to be specifically downregulated in IVF-ET placentas, i.e., hsa-miR-204-5p, hsa-miR-1269a and hsa-miR-941. Notably, all three microRNAs have been previously implicated in placental dysfunction and potential fetal growth complications although the underlying mechanisms generally remain unclear. The downregulation of miR-204-5p is associated with fetal growth abnormalities and is enriched in related biological pathways and has also been found to be regulated in the presence of adverse pregnancy-related outcomes (31). For instance, miR-204-5p can inhibit angiogenesis by regulating pro-angiogenic genes such as ANGPT1 and members of the VEGF family (32). miR-1269a was found to be a risk factor for ectopic pregnancy, and currently known risk factors include assisted reproductive technologies such as in vitro artificial insemination and hormonal stimulation (33). It has also been suggested that miR-941 is expressed in trophoblast cells and involved in insulin-related intracellular signaling pathways such as Wnt signaling, phosphoinositide-3kinase, TGF-b signaling, and PPAR-gamma (34). In addition, miR-941 has been shown to target Keap1 to activate the Nrf2 signaling pathway, which in turn protects human endometrial cells from oxygen and glucose deprivation-re-oxygenation induced oxidative stress and programmed necrosis (35). These correlates provide a compelling rationale for the functional significance of these miRNAs in IVF-ET fetal growth.
Given the central importance of miRNA expression in the execution of the transcriptional programs, we predicted transcription factors that might regulate these DEmiRs. The top ranked transcription factor was specificity protein 1 (SP1), a zinc finger transcription factor that binds to a variety of GCrich motifs and regulates the expression and function of miRNAs as well as the expression of genes associated with embryonic development and differentiation (36,37). SP1 has been shown to regulate the placental glucocorticoid barrier by repressing the expression of 11b-hydroxysteroid dehydrogenase type 2, leading to fetal growth restriction (FGR) (38). One study identified the interaction of miR-331-3p with SP1 and the interaction of miR-331-3p and miR-1908-5p with glycosyltransferases as a novel mechanism for ABO blood group regulation (39). The second ranked hit in the TF analysis was EGR1 which has been previously implicated in follicular development, ovulation, corpus luteum formation and placental angiogenesis (40), and plays a key role in placental implantation (41). Notably, EGR1 was also identified in the construction of the PPI network.
Here we screened key genes altered in IVF-ET placental tissue based on mRNA next-generation sequencing data and online tools to identify ten hub genes. Among these, EGFR signals through Src-and ERK-mediated pathways activated by VEGFR2. Placental trophoblast cells are enriched in EGFR (42) and activation of EGFR regulates the proliferation, migration and invasive capacity of extravillous trophoblast cells (43). FOS is another molecule involved in angiogenesis. FOS belongs to the transcription factor-activated protein 1 (AP-1) superfamily, which is responsible for a variety of cellular processes, including proliferation, differentiation, apoptosis, hypoxia, angiogenesis and steroidogenesis (44,45), as demonstrated in trophoblast cells (46). LEP, an important metabolic hormone, is highly expressed in the placenta and regulates placental, fetal growth and angiogenesis (47). SPP1, also called osteopontin (OPN), is located in the cytoplasm of placental syncytiotrophoblast and capillary endothelial cells and is considered a marker of placental bed remodeling. Placental development occurs in a hypoxic environment and can stimulate angiogenesis through upregulation of the vascular endothelial growth factor inhibitor of fibrinolytic plasminogen activator 1 (SERPINE1) (48).
Considering the important role of these key genes in placental development and angiogenesis, it was instructive to consider how these gene networks were impacted by altered miRNA expression. We used the intersection of miRNA and mRNA sequencing data to construct a microRNA-mRNA network. This analysis revealed more complexity in predicted interactions for the downregulated DeMiRs compared to their upregulated counterparts. Indeed, the downregulated DeMiRs display a far more extensive repertoire of target genes, particularly, hsa-miR-204-5p, hsa-miR-1269a and hsa-miR-941, for which, as discussed above, have tangible functional links to different aspects of placental regulation. Nonetheless, more work is needed to determine which gene targets are most functionally important but further understanding of the underlying mechanisms may help to reduce pregnancy complications and improve offspring safety.
Our current study has its limitations. The discovery cohort was small and we only validated selected miRNAs in the placentas of 8 IVF-ET assisted conception patients and 8 normal controls, which may reduce the reliability of our findings. Moreover, while we were diligent to sex match the placentals and minimize sampling bias by location, placental tissue is known to be divided into maternal, intermediate and fetal zones with different gene expression profiles (49) and interplacental differences are inevitable. Moreover, the molecular interactions and proposed functional relationships proposed by the bioinformatic analyses need formal verification in both experimental and clinical settings. Thus, more work is needed to explore the specific functions of differentially expressed genes in IVF-ET placentas and their mechanisms.

CONCLUSIONS
In summary, we constructed a miRNA-mRNA regulatory network to regulate the expression of genes essential for IVF-ET placental development and function using bioinformatic analysis. Our data reveal both previously identified miRNAs and mRNAs associated with placental dysfunction and pregnancy complications along with novel candidates. As such these data represent a ready resource for subsequent investigations into the effects of IVF-ET on placental development and function and may provide a basis for future prevention and treatment of adverse perinatal outcomes.

DATA AVAILABILITY STATEMENT
The data presented in the study is deposited in the NCBI repository, accession number GSE186149.

ETHICS STATEMENT
Tissue collection was approved by the Ethics Committee of the Third Affiliated Hospital of Zhengzhou University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
SY conducted the majority of the experiments and prepared the manuscript. WZ, MM, and SS were involved in the analysis and interpretation of the data. CY, RZ, NZ, SR, and HW collected samples of placenta and clinical data. RT contributed to the revision and critical discussion of the article. YG was engaged in research design, article drafting and critical discussion. All authors contributed to the article and approved the submitted version.