Integrative analysis of the immune-related ceRNA network in fetal growth restriction based on weighted gene co-expression network analysis

Immune disorders lead to placental dysfunction and fetal growth restriction (FGR), but current research on the immune regulation mechanisms of FGR and the involvement of competing for endogenous RNA (ceRNA) is insufficient. Therefore, this study aimed to construct an immune-related ceRNA network to predict FGR onset risk. Microarray data from the GEO database was used for gene expression and immune infiltration analyses. Weighted gene co-expression network analysis (WGCNA) was used to screen immune-related module genes with differential expression (DE) in FGR. A ceRNA network was constructed by integrating long non-coding RNA (lncRNA)-mRNA, lncRNA-microRNA (miRNA), and miRNA-mRNA relations. The diagnostic values of key genes in the network and their relationships with immune cell were confirmed in a validation cohort. By comparing FGR and normal samples, DE mRNAs, miRNAs, lncRNAs, and four types of immune cells with different infiltration levels were obtained. WGCNA then revealed 236 immune-related DE mRNAs that were involved in hormone secretion and immune cell differentiation. Based on co-expression analysis and miRNA prediction, we initially constructed a ceRNA network to screen several immune-related genes as potential diagnostic biomarkers of FGR, whose superior predictive performances were further confirmed by receiver operating characteristic curves. Among them, NEURL1 and ODF3B were found to positively correlate with M1 macrophages and may participate in the immunoregulation of FGR. From the perspective of ceRNA mechanism, we constructed an immune-related regulatory network for the first time wherein key genes are initially proposed as potential diagnostic biomarkers of FGR to involve in M1 macrophage-mediated immunoregulation.


Introduction
Fetal growth restriction (FGR), known as a common complication of pregnancy, is characterized by a fetal growth rate below the normal growth potential for a given gestational age [1,2]. It is a risk factor for adverse perinatal outcomes and may increase the risk of intrauterine fetal death three to seven times [3]. FGR can also cause a series of perinatal problems and postnatal complications, including neurodevelopmental disorders, metabolic syndrome, and cardiovascular and cerebrovascular diseases [4,5]. The etiology of FGR is complex and controlled by maternal, placental, or fetal factors. Furthermore, genetic factors play important roles, among which fetal-and placentalrelated genes have been identified in the etiology of FGR [6]. The complex etiology dictates the need for preventive and diagnostic testing for FGR to develop interventions and viable postpartum management.
FGR is associated with impaired cellular and humoral immunity, whereas placental dysfunction may be caused by an altered balance of immune cells [7]. Natural immune cells, including macrophages, dendritic cells, and natural killer cells, can influence placental development through immune regulation, growth factor supply, and maintenance of the uterine vascular system [8,9]. A high-throughput RNA sequencing analysis based on the FGR sheep model showed that the transcriptome responsible for B cell dysfunction was involved in immune response, WNT signal transduction, and adaptive stress response pathways [10]. Another transcriptome analysis screened candidates that were differentially expressed (DE) in the placenta of FGR patients, most of which participated in inflammatory and immune disorders [11]. Therefore, it is clinically significant to evaluate the potential diagnostic markers of FGR from the perspective of immune regulation.
Non-coding RNA is an important component of the human transcriptional genome and plays a regulatory role via various biological processes without encoding proteins [12]. Long non-coding RNAs (lncRNAs) serve as micro-RNA (miRNA) sponges to participate in competitive binding, thereby weakening the inhibitory impact of miRNAs on targeted mRNAs. They are known as competing endogenous RNAs (ceRNAs) to involve in biological functions and disease progression [13]. However, the mechanism of lncRNA-mediated ceRNA in FGR has not been studied and then attracted our interest. Therefore, this study sought to construct an immune-related ceRNA network to identify biomarkers that can be used for the diagnosis and risk prediction of FGR. We also explored the potential regulatory relationships between key genes and immune cells to elucidate the immune response mechanism of FGR.

Data collection
The RNA-seq of GSE114691 and miRNA-seq of GSE114349 datasets were captured from the gene expression omnibus (GEO) database [14]. These two datasets originated from the same batch of samples and contained 21 control and 18 FGR placenta samples (Homo sapiens), which were detected using a GPL11154 Illumina HiSeq 2000 platform (Homo sapiens). Furthermore, the GSE24129 dataset containing eight normal and eight FGR placenta samples was also collected from the GEO database as an external validation cohort.

Data preprocessing and gene annotation
For the RNA-seq data in GSE114349, mRNAs and lncR-NAs were annotated according to the Ensembl ID in the count matrix. Genes annotated with "protein_coding" were retained as mRNAs, whereas genes annotated with "3prime_ overlapping_ncrna", "antisense", "sense_intronic", "lin-cRNA", "sense_overlapping", and "processed_transcript" were retained as lncRNAs. For the GSE24129 dataset, the normalized and log2-transformed expression matrix was obtained from the GEO database. The annotation file was also downloaded from the detection platform. After matching with the probes, gene expression levels were calculated.

Analysis of DE mRNAs, DE miRNAs, and DE lncRNAs
For the original expression matrix of RNA-seq and miRNAseq, the voom function from the Limma package (version 3.10.3) [15] was used to normalize the count value, followed by the difference analysis of gene expression between the FGR and control groups using the classical Bayesian method. Then, DE mRNAs, DE miRNAs, and DE lncRNAs that exhibited statistical significance were selected using the thresholds of |log fold change (FC)|> 1 and p value < 0.05, as adjusted by the Benjamini & Hochberg (BH) method.

Immune infiltration analysis
CIBERSORT [16,17] was used to calculate the infiltration abundances of 22 types of immune cells in each sample. By comparing FGR and normal samples, the differences in the infiltration abundance of immune cells between the two groups were analyzed using the Wilcox test and the statistical significance was set at p < 0.05.

Weighted gene co-expression network analysis (WGCNA) of immune-related genes
Based on the mRNA expression matrix, the median absolute deviation (MAD) of each mRNA was calculated. After removing the 25% mRNA ranked lowest in MAD, R.WGCNA (version 1.61) [18] was used to construct a scale-free co-expression network. Based on clustering and dynamic tree cut algorithms, highly correlated mRNAs were aggregated into a module with the parameters set as minModuleSize = 30 and MEDissThres = 0.3. Considering the infiltrate score of differential immune cells as traits, we calculated the correlation between modules and traits and obtained key module genes related to immune cells. The intersection of immune-related genes and DE mRNAs was then determined to obtain immune-related DE mRNAs.

Functional and pathway enrichment analyses
ClusterProfiler (version 3.8.1) [19] was used for gene ontology biological process (GO-BP) [20] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [21] pathway enrichment analyses of immune-related DE mRNAs. Enriched terms with BH adjusted p < 0.05 were considered as having statistical significance.

Co-expression analysis of DE mRNAs and DE lncRNAs related to immune cells
By one-to-one correspondence of samples, the Pearson correlation coefficients of immune-related DE mRNAs and DE lncRNAs were calculated, and significant relations were set at BH adjusted p < 0.05 and |r|> 0.8. To explore potential functions and pathways involving lncRNAs, R.clusterProfiler was used for enrichment analyses of lncRNA-targeted mRNAs in these relation pairs.

Prediction of miRNAs
The online tool. mirwalk3.0 [22] was used to predict the upstream miRNAs of immune-related DE mRNAs. With thresholds of binding probability > 95% and binding site position = 3UTR, the miRNA-gene relation pairs were obtained. The Pearson correlation analysis was further performed to select negative regulatory miRNA-mRNA relations at standards of BH adjusted p < 0.05 and r < -0.3.
To predict the downstream miRNAs of DE lncR-NAs, DIANA tools-LncBase v.2 [23] was applied with a score > 0.6. The obtained miRNAs were intersected with DE miRNAs to generate lncRNA-miRNA pairs.

Construction of a ceRNA network
By integrating the negative regulatory miRNA-mRNA relation pairs and lncRNA-miRNA relation pairs obtained above, the lncRNA-miRNA-mRNA regulatory axis was obtained. Combined with co-expression relations of immune-related DE mRNA-DE lncRNAs, a ceRNA network was finally generated using Cytoscape (version 3.6.0) [24].

Validation of genes in the ceRNA network
The gene expression levels in the ceRNA network were extracted from normal and FGR samples from the GSE24129 dataset. Then, receiver operating characteristic (ROC) curves were created to validate their diagnostic values. The CIBERSORT algorithm was used to screen immune cells with differences in infiltration abundance in the validation cohort. After selecting immune cells that matched the trend of differences in the training set, Spearman correlation coefficients between genes in the ceRNA network and key immune cells were calculated to observe their potential relationships.

Analysis of immune infiltration differences between FGR and normal samples
Immune cell proportions in each FGR and normal sample were estimated using CIBERSORT ( Fig. 2A). The results suggested that four types of immune cells significantly differed in their infiltration abundance (Fig. 2B). The proportions of memory B cells and resting memory CD4 T cells were significantly elevated in the FGR samples, whereas M0 and M1 macrophages showed significantly higher infiltration levels in the FGR group than in the control group.

Identification of immune-related DE mRNAs using WGCNA
Based on mRNA expression levels in each sample, a WGCNA was used to determine immune-related modules. In this study, a power of β = 4 (scale-free R 2 = 0.85) was selected as the soft threshold to ensure a scale-free network (Fig. 3A). Through hierarchical clustering and dynamic tree cut algorithms, modules with a dissimilarity coefficient < 0.3 were merged, and a total of ten modules were finally integrated (Fig. 3B). By calculating the relationships between modules and traits (grouping and immune cell infiltration level), we found that the purple, dark green, and light cyan modules were significantly correlated with traits (p < 0.05 and |r|> 0.7), as shown in Fig. 3C. Therefore, the genes in these three modules were identified as immune-related. Considering the overlapping of immune-related genes and DE mRNAs, 236 immune-related DE mRNAs were obtained (Fig. 3D).

GO and KEGG enrichment analyses on immune-related DE mRNAs
Enrichment analyses were carried out to observe the potential biological functions and pathways in which these were screened according to the following thresholds: |logFC|> 1 and BH-adjusted p < 0.05. Red triangles indicate upregulated genes, and blue squares indicate downregulated genes candidate genes may be involved. The results indicated that 236 immune-related DE mRNAs were mainly enriched in GO-BPs for response to lipopolysaccharide, hormone transport, and molecules of bacterial origin (Fig. 4A), as well as KEGG pathways for cell adhesion molecules (Fig. 4B).

Co-expression analysis of immune-related mRNAs and lncRNAs
To obtain immune-related lncRNA-mRNA co-expression pairs, the Pearson correlation coefficients of immunerelated DE mRNAs and DE lncRNAs were computed.
With thresholds of BH adjusted p < 0.05 and |r|> 0.8, a total of 648 lncRNA-mRNA relation pairs comprising 54 lncRNAs and 83 mRNAs were constructed (Fig. 5A). To further explore the functions of the lncRNAs in this network, we conducted enrichment analyses on their targeted mRNAs. Results suggested that these immune-related lncRNAs may participate in GO-BPs, including gonadotropin secretion, regulation of the endocrine process, and regulation of transmembrane receptor protein-related signals (Fig. 5B), as well as KEGG pathways of Th17 cell differentiation and cytokine-cytokine receptor interaction (Fig. 5C).

Generation of an immune-related ceRNA network
The upstream miRNAs of the immune-related DE mRNAs were predicted, and then 312 miRNA-mRNA relation pairs were obtained based on their negative correlations. The potential correlations between DE lncRNAs and miR-NAs were also predicted, and 711 lncRNA-miRNA pairs were acquired. The lncRNA-miRNA-mRNA regulatory axis was constructed by integrating miRNA-mRNA and lncRNA-miRNA relation pairs. Combined with immunerelated DE mRNA-lncRNA co-expression relation pairs, a ceRNA regulatory network comprising 21 mRNAs, 16 miRNAs, and 16 lncRNAs was established (Fig. 6). In this network, the RP11-488P3.1-hsa-miR-1271-5p-ARNT2 and CTC-550B14.7-hsa-miR-3175-MAP3K9 regulatory axes may play important roles because their nodes have larger degrees of connections.

Verification of key genes in the ceRNA network based on an external validation cohort
The expression values of genes from the above ceRNA network were extracted from each sample of the validation set GSE24129, and a total of 21 mRNAs and one lncRNA (UCA1) were matched. ROC curves were created for these 22 genes to validate their diagnostic values. As a result, the area under the curve (AUCs) of 17 of these genes was greater than 0.6, indicating their high specificity and sensitivity for diagnosing FGR (Fig. 7A). The difference in immune cell infiltration between the FGR and control groups was also analyzed in the validation cohort. The box plot indicates two types of immune cells with statistical significance, but only the M1 macrophage shows consistent results with the training set, which means their infiltration level was significantly higher in FGR samples than in the control group (Fig. 7B). By calculating the correlation coefficients between M1 macrophages and the 22 candidate genes, significant positive correlations were found between NEURL1 and ODF3B expressions and M1 macrophage infiltrations (Fig. 7C-D). Furthermore, we found that UCA1 was significantly positively correlated with the majority of mRNAs, which further confirms the co-expression relationships between candidate mRNAs and lncRNAs.

Discussion
FGR is the second leading cause of perinatal mortality and morbidity worldwide and has become a public health concern [25]. The importance of genetics in the etiology of FGR has been emphasized. A transcriptome analysis based on FGR rats found that thousands of transcripts were modified more than 2.5 times in the experimental group, suggesting a significant contribution of genetic alterations to the pathogenesis of FGR [26]. However, there is still a lack of appropriate diagnostic methods in clinics. Therefore, it is necessary to explore the cellular and molecular mechanisms based on the genetic background to identify effective biomarkers and explore therapeutic interventions. Considering the important role of immune cells in placental development, we constructed an immune-related ceRNA network to identify potential diagnostic biomarkers of FGR and regulatory pathways in which they may be involved in the FGR progression.
Based on lncRNA-mRNA, lncRNA-miRNA, and miRNA-mRNA relations, a ceRNA network was constructed. In this network, the RP11-488P3.1-hsa-miR-1271-5p-ARNT2 and CTC-550B14.7-hsa-miR-3175-MAP3K9 regulatory axes may play more important roles. A comparison of the epigenetic modification sites of FGR and healthy placenta revealed that the transcription factor ARNT2, which corresponds to motifs enriched in the differential acetylation region, was upregulated in FGR [27]. Studies based on animal models have shown that ARNT2 is also a key gene associated with the innate immune phenotype, and is involved in immune response pathways and metabolism, such as lymphocyte differentiation and activation [28]. MiR-1271-5p has been found to participate in toll-like receptor signaling pathways and the activation of immune cells [29]. Although no studies have confirmed the relationship between RP11-488P3.1, miR-1271-5p, and ARNT2, we hypothesized that the regulatory mechanism of ARNT2 and miR-1271-5p may affect the innate immune phenotype of FGR based on the above findings. In addition, miR-3175 participates in the MAPK-related signaling pathway, whereas MAP3K9 is an upstream activator of MAPK signaling [30,31]. Choi et al. also found the immunoreactivity of MAPK signaling in FGR rats [32]. Therefore, miR-3175-MAP3K9 may trigger the pathogenesis of FGR by regulating the immune activity of the downstream MAPK pathway. However, further experiments are still needed to confirm our hypothesis. Fig. 6 The construction of a ceRNA network by integrating miRNA-mRNA negative correlations, lncRNA-miRNA relations, and mRNAs-lncRNAs co-expression relation pairs. Green and red circles indicate downregulated and upregulated mRNAs, respectively; purple and yellow-green triangles indicate downregulated and upregulated miRNAs, respectively; and blue and pink diamonds indicate downregulated and upregulated lncRNAs, respectively By analyzing immune cell infiltration levels, the infiltration of M1 macrophages was significantly increased in the FGR group in this study. CD68 + M1 macrophages have been reported to be associated with placental dysfunction under fetal growth restriction [33]. Abnormal activation of macrophages in utero affects trophoblast function and placental development, leading to a variety of adverse pregnancy outcomes [34]. These findings support our conclusion that increased levels of abnormal activation and infiltration of M1 macrophages affect placental development. In this study, we also found that M1 macrophage infiltration was positively associated with the expression of NEURL1 and ODF3B. Of these, ODF3B is related to methylation and is involved in pathways of immune cell migration, proliferation, activation, and inflammatory activity [35]. A previous study also reported the regulation of the neuralized family member NEURL1 and the downstream Notch signaling pathway, which is also closely related to macrophage activation [36,37]. Using an independent dataset, we further indicated that the potential of NEURL1 and ODF3B in predicting FGR risks may be related to the transfer of macrophage polarity to the M1 phenotype.

Study limitations
Although this study found the possible immune-related ceRNA mechanism in FGR by means of a series of bioinformatics analyses, the lack of experimental verification is one of the non-negligible limitations in this study. Furthermore, the number of samples we included was under restriction due to the lack of publicly available chip data of FGR. In the follow-up research, mechanistic studies will be carried out to further confirm our findings and explore the competitive regulatory relationships between key mRNAs, miRNAs, and lncRNAs in FGR. We will also explore the role of macrophages in FGR and their relationships with NEURL1 and ODF3B expression based on clinical FGR placental samples.

Conclusion
To summarize, we constructed an immune-related ceRNA network for the first time to screen key genes which may predict the risk of FGR onset. Most candidate genes have been verified to have high specificity and sensitivity for the diagnosis of FGR. Among them, NEURL1 and ODF3B were