Bioinformatics Analysis Identifies Potential Related Genes in the Pathogenesis of Intrauterine Fetal Growth Retardation.

Background: Intrauterine growth retardation (IUGR) affects approximately 10% to 15% of all pregnancies worldwide. IUGR is not only associated with stillbirth and newborn death, but also the delay of cognition in childhood and the promotion of metabolic and vascular disorders in adulthood. Figuring out the mechanism of IUGR is rather meaningful and valuable.
Methods: Datasets related to IUGR were searched in the Gene Expression Omnibus website. Principal component analysis (PCA) was used for normalization. Differential expressed genes (DEGs) were screened out using the ggpot2 tool. DEGs were used to conduct Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses, and protein-protein interaction (PPI) analysis. IUGR related genes were searched in the OMIM website to look for the intersection with the DEGs. The DEGs were analyzed for tissue-specific expression by the online resource BioGPS. The results were displayed through volcano map, Venn map, box plot, heat map, and GSEA enrichment plots drawn by R language packages.
Results: Eleven DEGs were screened out of 2 datasets. One hundred ninety-five genes related to IUGR in OMIM were retrieved. EGR2 was the only intersection gene that was found in both groups. Genes associated with placental tissue expression include COL17A1, HSD11B1, and LGALS14. Molecular functions of the DEGs are related to the oxidoreductase activity. The following 4 signaling pathways, reactome signaling by interleukins, reactome collagen degradation, Naba secreted factors, and PID NFAT tfpathway, were enriched by GSEA. Two critical modules comprising 5 up-regulated genes (LEP, PRL, TAC3, MMP14, and ADAMTS4) and 4 down-regulated genes (TIMP4, FOS, CCK, and KISS1) were identified by PPI analysis. Finally, we identified 6 genes (PRL, LGALS14, EGR2, TAC3, LEP, and KISS1) that are potentially relevant to the pathophysiology of IUGR.
Conclusion: The candidate down-regulated genes LGALS14 and KISS1, as well as the up-regulated genes PRL, EGR2, TAC3, and LEP, were found to be closely related to IUGR by bioinformatics analysis. These hub genes are related to hypoxia and oxidoreductase activities in placental development. We provide useful and novel information to explore the potential mechanism of IUGR and make efforts to the prevention of IUGR.

Introduction

Intrauterine growth retardation(IUGR) is defined as an estimated fetal weight or abdominal circumference that is less than the 10th percentile for gestational age. IUGR is strongly associated with stillbirth and newborn death,[2-4] as well as the delay of cognition in childhood and the promotion of metabolic and vascular disorders in adulthood.[5,6] The etiology of IUGR can be divided into 3 categories: maternal, fetal, and placental, with the latter causing suboptimal uterine–placental perfusion. The placenta is the fetus’s principal contact with the mother, and it plays a vital role in fetal development and growth by allowing substrate transfer and moderating the maternal immune response to prevent immunological rejection of the conceptus. Most data shows that placental insufficiency is the most common pathology associated with IUGR. Placental restriction and insufficiency have been linked to several placental changes in IUGR pregnancies, including altered placental growth and substrate transport capacity, increased apoptosis and autophagy, and increased glucocorticoid action. Unfortunately, the mechanism of placental insufficiency is still unknown. With the development of bioinformatics, there are lots of open databases established. We performed a genome-wide gene expression analysis to identify the differentially expressed genes (DEGs) between IUGR and normal pregnancy using the placental microarray datasets GSE12216 and GSE147776 retrieved from Gene Expression Omnibus (GEO). A series of bioinformatic analyses were utilized for exploring the potential mechanism of IUGR. We attempted to identify genes that are potentially related to the pathogenesis of IUGR and to explore the diagnostic value of these screened-out genes.

Materials and Methods

Microarray data

The microarray expression profiling datasets, GSE12216 and GSE147776, were downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The dataset GSE12216 evaluated the global placental gene expression profiles in 8 IUGR and 8 healthy pregnancies based on the GPL2986 ABI Human Genome Survey Microarray Version 2 platform. The GSE147776 dataset was created by the GPL20844 Agilent-072363 SurePrint G3 Human GE v38x60K Microarray 039494 [Feature Number Version] platform. Consider changing the experiment to the GSE12216 dataset, the experiment contained 36 placental samples consisting of IUGR (n = 7), pre-eclampsia (n = 7), pre-eclampsia and IUGR (n = 6), and normal pregnancy (n = 8). We only included the subgroups of IUGR and normal pregnancy in our analysis. In combination, there were 15 IUGR and 16 healthy placentas were included in our study. The gestational age at birth of datasets GSE147776 was 36.8 ± 0.58 week, while that of GSE12216 was 236 ± 24 days, the data from both sets had comparability. The correlated annotation files for the platform were also downloaded from the GEO. Figure 1 displayed the overall research design.

Figure 1.

The overall research designs. The data were downloaded from GEO and OMIM databases.

Differential expression analysis

For normalization and principal component analysis (PCA), as well as to display the findings as heatmaps and volcano plots, the R package “ggplot2, version 3.3.3” was used. Quantitative performance was evaluated with normalization to allow data from different sets to be compared. The differences between groups became visible when the samples of each group were screened. Differential expressed genes (DEGs) were screened out using the ggpot2 tool. Each sample’s genes were kept if they met the following criteria: (1) a |log2 (fold-change)| > 1 and (2) a P-value < .05. The top 20 up-regulated genes were presented as a heatmap.

Screen for the IUGR related hub genes

The keyword “intrauterine growth retardation” was searched at the OMIM website (https://www.omim.org/), which may provide us IUGR related hub genes. These genes had been identified by other researches. A Venn diagram was utilized to explore the intersection of hub genes from OMIN and the DEGs from our analysis.

Tissue-specific gene expression analysis

The tissue-specific pattern of mRNA expression can reveal crucial information regarding gene function. BioGPS is a full resource for learning about genes and protein functions. The DEGs were analyzed the tissue-specific expression by the online resource BioGPS (http://biogps.org/). Highly tissue specific transcripts were identified as those that were mapped to a particular tissue and met the following criteria: (1) The tissue-specific expression level was more than 10 times higher than the median, and (2) the second highest expression level was less than one-third of the highest level.

Functional enrichment analysis of DEGs

The DEGs were then processed on the R package ggplot2 for functional enrichment. The R package of “clusterProfiler, version 3.14.3” and “org.Hs.eg.db, version 3.10.0” were used to conduct Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses for statistically significant DEGs. The cut-off threshold of the analysis was set as P.adj < .1 and q value < 0.2. The results were presented as a bar plot via the “ggplot2” package in R.

Gene set enrichment analysis

We performed Gene Set Enrichment Analysis (GSEA) of the DEGs of each dataset via the “clusterProfiler, version 3.14.3” package to determine the biological pathway. As reference gene sets in GSEA, C2: curated gene sets from MSigDB collections were chosen (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Pathways with a false discovery rate (FDR) <0.25 and an adjusted P-value < .05 were considered to be significant enrichment. Gene set permutation was performed 5000 times for GSE147776 analysis, and 1000 times for GSE12216 analysis.

Protein-protein interaction (PPI) network analysis

To figure out the potential protein interactions, PPI networks were investigated using the STRING database (http://string-db.org) and displayed using the Cytoscape software (version 3.7.2). We mapped the DEGs from the datasets onto the PPI network and set a threshold of medium confidence (0.40) for interaction scores. The app Molecular Complex Detection (MCODE) from the Cytoscape software suite was used to identify the key modules.

Results

Data process and differentially expressed genes

Normalization was used to remove global differences from the differential samples so that they could be compared. As we illustrated, the median of each sample is essentially the same, indicating that the degree of normalization between samples is good (Figure 2A and B). Additionally, the PCA of GSE12216 showed the ratio of PC1 and PC2 is 29.5% (Figure 3A), while the combined proportion of GSE147776 is 36% (Figure 3B), suggesting that there were abundant DEGs in these 2 datasets. GSE12216 had 154 up-regulated genes and 74 down-regulated genes that were enrolled for further investigation (Figure 3C), while GSE147776 had 901 up-regulated genes and 916 down-regulated genes that were screened out (Figure 3D). The top 20 relative up-regulated DEGs of each group were exhibited by heatmaps (Figure 4A and B). The details of the top 20 up-and down-regulated DEGs are provided in Supplemental Table 1.

Figure 2.

The box plot of normalization: (A) Box plot of GSE12216 gene chip and (B) Box plot of GSE147776 gene chip.

Figure 3.

(A) PCA analysis plot of GSE12216 gene chip. (B) PCA analysis plot of GSE147776 gene chip. Red and blue spots represent samples from IUGR group and Normal group, respectively. (C) Volcano map of GSE12216 gene chip. (D) Volcano map of GSE2147776 gene chip. Red dots represent up-regulated genes and green represent down-regulated genes.

Figure 4.

(A) Heatmap of up-regulated genes in GSE12216 gene chip, (B) heatmap of up-regulated genes in GSE147776 gene chip, and (C) Venn diagram of DEGs in GSE12216 and GSE147776 gene chips and OMIM.

In the OMIM online catalog of human genes and genetic illnesses, we discovered 195 genes related to intrauterine growth retardation that have been published in the literature. By using a Venn diagram, only one intersection hub gene from OMIN and the DEGs, EGR2, was discovered (Figure 4C). The DEGs from the 2 datasets are listed as DHRS2, COL17A1, S100A12, ADAMTS4, HSD11B1, KHDRBS3, PRL LGALS14, SLC2A5, EGR2, and DNM1 in Supplemental Table 2.

Tissue-specific expression of genes

The screened out 11 DEGs were explored in a specific tissue or organ system using BioGPS. The top 2 associated tissues were recorded for each differential gene (Table 1). Genes associated with placental tissue expression include COL17A1, HSD11B1, and LGALS14.

Table 1.

Tissue-specific expressed genes identified by BioGPS.

Gene	The tissue-specific pattern of mRNA expression
	The first expression tissue	The second expression tissue
DHRS2	Kidney	Bonemarrow
EGR2	CD33+ myeloid	Thyroid
COL17A1	Bronchial epithelial cell	Placenta
S100A12	CD33+ myeloid	Bonemarrow
ADAMTS4	Ovary	Spinal cord
HSD11B1	Liver	Placenta
KHDRBS3	Testis interstitial	Testis Germ cell
PRL	Pituitary
LGALS14	Placenta
SLC2A5	Lymphoma Burkitts
DNM1	Prefrontal cortex	Amygdala

Functional and pathway enrichment of DEGs

For GO and KEGG analysis, the 11 intersection DEGs were all changed to Entrez ID. Three molecular functions (MF) have been enhanced when P.adj is less than .05 and q value is less than 0.2. GO:0016616, GO:0016614, and GO:0017124 are the identifiers (Figure 5A). Supplemental Table 3 contained the details of the MF enrichment.

Figure 5.

(A) Results of GO Molecular Function analysis of 11 DEGs. Enrichment plots of GSEA involved in (B and C) reactome signaling by interleukins, (D and G) reactome collagen degradation, (E and H) Naba secreted factors, and (F and I) PID NFAT tf pathway.

The DEGs from GSE147776 were utilized for GSEA. There were 192 datasets satisfied the FDR < 0.25 and P.adjust < .05 requirements. At the same time, the DEGs of GSE12216 were assessed on an equal footing, and 164 datasets were enriched as a result. Placental tissue correlated genes (COL17A1, HSD11B1, LGALS14) and EGR2 from OMIM were individually explored in the 2 datasets. The following 4 signaling pathways: reactome signaling by interleukins, reactome collagen degradation, Naba secreted factors, and PID NFAT tfpathway, were found in both datasets (Figure 5B-I).

PPI network analysis of DEGs

According to the STRING online database, a PPI network with 33 nodes and 68 edges was constructed, with an interaction score >0.4 (Figure 6A). The nodes represent genes, and the edges represent gene connections. The color red denotes genes that have been up-regulated, while blue nodes stand for down-regulating. To identify the major PPI network modules, we utilized the MCODE program in Cytoscape to conduct network gene clustering. Two critical modules comprising 5 up-regulated genes (LEP, PRL, TAC3, MMP14, and ADAMTS4) and 4 down-regulated genes (TIMP4, FOS, CCK, and KISS1) were identified, as illustrated in Figure 6B and C. Additionally, functional enrichment analysis revealed that these 9 genes were mostly involved in protease binding, metallopeptidase activity, and metalloendopeptidase activity (Supplemental Table 4).

Figure 6.

(A) Cytoscape network visualization of the 33 nodes and 68 edges that was obtained with interaction scores >0.4 according to the STRING online database. (B and C) Two key modules were identified by MCODE, which was used to identify hub gene.

Identification of hub genes

Genes linked to IUGR were found by the following methods after a rigorous inspection. Firstly, PPI network analysis revealed 2 critical modules containing 9 genes (LEP, PRL, TAC3, MMP14, ADAMTS4, TIMP4, FOS, CCK, and KISS1). Secondly, the tissue-specific gene expression analysis revealed that the placenta expressed COL17A1, HSD11B1, and LGALS14. Finally, according to the OMIN website, EGR2 is a hub gene. In addition, we manually identified 6 genes (PRL, LGALS14, EGR2, TAC3, LEP, and KISS1) that are potentially relevant to the pathophysiology of IUGR using the GeneCards database and literature review. In the HPA database (https://www.proteinatlas.org/), the expression of these 6 DEGs were checked. KISS1 was stained high in heathy placenta, while the genes: PRL, TAC3, and LEP were negatively stained. The immunohistochemical pictures of EGR2 and LGALS14 could not been searched (Figure 7).

Figure 7.

Protein expression in normal placenta: (A) KISS1, (B) PRL, (C) TAC3, and (D) LEP.

Discussion

The placenta serves as the primary link between the fetal and the maternal circulation, so placental insufficiency could cause the placenta to fail to provide an appropriate amount of substrates to the fetus, resulting in IUGR. A lot of placental genes have been reported to associate with IUGR to a certain degree, such as placental insulin-like growth factor 1 (IGF1), epidermal growth factor (EGF), endoglin, and vascular endothelial growth factor (VEGF-A).

In our study, 11 DEGs have been screened out which are strongly related to IUGR. GO function enrichment resulted in the oxidoreductase activity (GO:0016616, GO:0016614) and SH3 domain binding which is found in a great variety of intracellular or membrane-associated proteins. According to the annotation of the quickGO web (https://www.ebi.ac.uk). The go:0016616 could activate the go:0016614, which activates oxidoreductase activity, cascades, into molecular function. The antioxidant activity, oxygen carrier activity, ATP-dependent activity, binding and translation regulator activity were all included in the molecular function. There are lots of studies that have reported the relationship between IUGR and oxidoreductase activity. Reactive oxygen species and oxidative stress appear to be important factors in the physiological and pathological states of the IUGR placenta. Heme oxygenases (HO) have become essential regulators of cardiovascular function in recent decades, owing to their synthesis of physiologically active metabolites such as carbon monoxide, bilirubin, and elemental iron.[18,19] In the field of physiological and pathological placental function, notably, the protective role of HO-1 against IUGR has been shown.[20,21] Indoleamine 2,3-dioxygenase(IDO) activity was significantly lower in placenta with IUGR, suggesting that the importance of placental IDO during fetal development.[22,23]

Four GSEA pathways were enriched in our research, the first one is Naba secreted factors pathway. Three kinds of matrisome-associated proteins were characterized in this pathway by using bioinformatic pipelines identical to those used to describe the core matrisome: ECM (Extracellular Matrix)-affiliated proteins, ECM regulators, and secreted factors. ECM might regulate the angiogenesis of placenta to influence fetal development.[24-26] The second one is PID NFAT tfpathway, which is calcineurin-regulated NFAT-dependent transcription in lymphocytes. Five proteins (NFAT1-5) make up the NFAT transcription factor family. NFAT5 is a transcription factor, which has broader implications for development, immune function, and cellular stress responses. Dobierzewska et al had shown that NFAT5 is up-regulated in IUGR placental hypoxia and ischemia.

Reactome signaling by interleukins (IL) is the third enriched pathway. Interleukins are low-molecular-weight proteins that bind to cell surface receptors and act in an autocrine and/or paracrine way, influencing processes such as tissue growth and repair, hematopoiesis, and several levels of the host defense against pathogens. IL-6 was studied most. Placental IL-6 concentration was confirmed to relate to fetal growth.[29,30] Studies stated that higher concentrations of IL-6 were observed in the IUGR placenta,[31,32] in cord blood, While, down-regulation of IL-6 in IUGR was reported by Cecati et al. IUGR is linked to a reduction in IL-10 levels, and an increment in IL-1α levels. All of the above findings point to IL being an inflammatory factor in IUGR.

The last pathway that we enriched is reactome collagen degradation. There is evidence that changes in vascular and uteroplacental matrix metalloproteinases (MMPs) and collagen content could be corrected by angiogenic agents and MMP modulators, alleviating IUGR. According to these 4 pathways, inflammatory activities might incite ECM to alter placental maturation and affect fetal growth.

The genes that had been screened out in our study were divided into 2 groups. The down-regulated group contained the genes LGALS14 and KISS1. The remaining 4 genes, PRL, EGR2, TAC3, and LEP, were included in the up-regulated group. LGAS14 expression is strongly associated with placenta. LGALS14 is a strong inducer of T-cell apoptosis and could bind beta-galactoside and lactose. LGALS14 is a highly differentially expressed autosomal gene that regulates inflammation and the immune system, mediates cellular apoptosis and tissue development, facilitates metabolic processes, and regulates the cell cycle.[25,38] The expression of LGALS14 is down-regulated in adverse pregnancy.[39,40] This is consistent with our bioinformatical analysis of LGALS14. The other down-regulated gene is KISSS-1(Kisspeptin), which is a potent positive regulator of gonadotrophin-releasing hormone and leptin. KISS levels are lower in maternal serum in pregnancies associated with IUGR. Metastin, as a protein encoded by the KISS-1 gene, is significantly lower in maternal plasma with fetal growth impairment in the first trimester. This provides further evidence for KISS’s down-regulation. Further, we identified that PRL was up-regulated in IUGR placentas. A higher amount of PRL was found in the cord blood of neonates born to malnourished and anemic mothers. indicating an adaptive reaction on the part of the fetus to offset an in-utero growth disadvantage. The early growth response 2 (EGR2) is another up-regulated gene in our study. IL-6 is known as an inducer of EGR2, its up-regulation could promote the expression of EGR2 in IUGR. Furthermore, Rt-PCR confirmed that TAC3 was significantly increased in both maternal blood and placenta in severe FGR compared to normal pregnancy, and correlated with the severity of IUGR.[45,46] LEP encodes leptin, which is a hormone released by the placenta and is important in fetal growth throughout the entirety of pregnancy. LEP up-regulation in IUGR placentas was reported in another 3 studies.[48-50] Leptin is crucial in placental development and function. Abnormal trophoblast proliferation or invasion could be related to excess placental leptin release in IUGR pregnancy. Leptin could stimulate lipolysis in the placenta and impact free fatty acid availability to the fetus.

In conclusion, the candidate down-regulated genes LGALS14 and KISS1, as well as the up-regulated genes PRL, EGR2, TAC3, and LEP, were found to be closely related to IUGR by bioinformatics analysis. What’s more, immunohistochemical expression of the DEGs is opposite in healthy placenta to the genes expression in IUGR’s placenta. We believe that these 6 hub genes may have an impact on placental development and function based on our findings. Because of their hypoxia and oxidoreductase activities, LEP and PRL may influence the trophoblast invasion. Micro-array was used to screen these hub genes out of the placenta, which could lead to future investigations focusing on circulating placental RNA. If the hub gene’s expression matches that of circulating placental RNA, this could be a novel non-invasive technique to explore more about the mechanism of IUGR. These hub genes have potential value in developing biomarkers for the early diagnosis or detection of IUGR.

Click here for additional data file.

Supplemental material, sj-rar-1-evb-10.1177_11769343221112780 for Bioinformatics Analysis Identifies Potential Related Genes in the Pathogenesis of Intrauterine Fetal Growth Retardation by Chao Xiao, Yao Wang and Yuchao Fan in Evolutionary Bioinformatics