The lncRNA-miRNA-mRNA ceRNA network in mural granulosa cells of patients with polycystic ovary syndrome: an analysis of Gene Expression Omnibus data.

BACKGROUND: Polycystic ovary syndrome (PCOS) is one of the most common endocrine abnormalities in women of reproductive age. In this study, we set out to construct a molecular long non-coding RNA (lncRNA)-microRNA (miRNA)-messenger RNA (mRNA) network according to the competitive endogenous RNA (ceRNA) theory and obtain insights into the related biological characteristics and pathways.
METHODS: We downloaded two gene expression profile datasets of mural granulosa cells (MGCs) of women with PCOS and healthy women without PCOS (GSE84376 and GSE106724) from Gene Expression Omnibus (GEO) DataSets. Using GEO2R, we identified the mRNAs and non-coding RNAs with differential expression. The DIANA-microT-CDS algorithm was applied to predict the genes targeted by the differentially expressed miRNAs. The lncRNA-miRNA interactions were predicted using DIANA-LncBase v2. Then, we constructed the lncRNA-miRNA-mRNA network. The Database for Annotation, Visualization and Integrated Discovery (DAVID) was employed to identify the functions and enriched pathways of the genes. Subsequently, STRING was used to construct the protein-protein interaction (PPI) network. cytoHubba in Cytoscape was used to rank the hub genes, and finally, PPI modules were screened with Cytoscape MCODE.
RESULTS: There were 462 mRNAs, 2,464 lncRNAs, and 55 miRNAs which showed differential expression between the MGCs of patients with PCOS and those of healthy controls. Based on the PPI analysis, differentially expressed genes (DEGs) were significantly enriched in retinol metabolism, drug metabolism-cytochrome P450, malaria, the Hippo signaling pathway, and glycine, serine, and threonine metabolism. The ceRNA network contained 71 lncRNA nodes, 14 miRNA nodes, and 69 mRNA nodes, as well as 167 edges. We identified some novel genes and non-coding RNAs that might be involved in PCOS, including CD163, MRC1, VSIG4, CCL2, CCR2, SPP1, hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-1231, hsa-miR-3609, and hsa-miR-4433b-3p.
CONCLUSIONS: This study identified a novel lncRNA-miRNA-mRNA network based on the ceRNA mechanism in PCOS. Some novel genes and non-coding RNAs that may be involved in the occurrence and development of PCOS were excavated, including CD163, MRC1, VSIG4, CCL2, CCR2, SPP1, hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-1231, hsa-miR-3609, and hsa-miR-4433b-3p. However, our findings need to be validated by in vivo and in vitro experiments. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Polycystic ovary syndrome (PCOS) is a highly prevalent hormone disorder in women of childbearing age as well as a prominent risk factor for female infertility. It is characterized by hyperandrogenism, ovulatory dysfunction, polycystic ovarian morphology, and metabolic disorders (including obesity, insulin resistance, and diabetes) (1,2). Despite the lack of clarity about PCOS at the etiological level, the majority of scientists agree that genetic factors play a pivotal role in the development and maintenance of the disorder.

Protein-coding genes comprise <2% of the human genome, while non-coding RNAs make up >90% of genomic transcription products (3). Non-coding RNAs include microRNAs (miRNAs) and long non-coding RNAs (lncRNAs), with lengths measuring 18–25 nucleotides and >200 nucleotides, respectively. At one point, scientists believed that lncRNAs were transcriptional “noise”, but increasing evidence has revealed that they are involved in various human diseases, such as cancers and reproductive system diseases (4,5). LncRNAs carry out an important molecular role by interacting with miRNAs, messenger RNAs (mRNAs), and even proteins. Among these mechanisms, lncRNAs function as competitive endogenous RNAs (ceRNAs) by communicating with mRNAs via competition for shared miRNAs (6,7).

In recent years, many lncRNAs and miRNAs were found to be associated with PCOS, such as lncRNA HOTAIR, miR-23a, and miR-23b (8,9). Fu et constructed a lncRNA-miRNA-mRNA ceRNA network in letrozole-induced PCOS rat model and verified it by quantitative real-time polymerase chain reaction (qRT-PCR) experiment (10). Microarray analysis of tissues from women with and without PCOS allows the differentially expressed mRNA and non-coding RNA expression profiles to be determined.

Mural granulosa cells (MGCs), which surround follicles, provide nutrients and growth regulators to oocytes during their development (11). Besides, MGCs participate in the arrest of oocytes in the meiotic phase through their secretion of oocyte maturation inhibitor (OMI) (12). Thus, dysfunction of MGCs may contribute to the aberrant folliculogenesis observed in PCOS. For the current research, the lncRNA, miRNA, and mRNA expression profiles of MGCs in women with and without PCOS were obtained. Subsequently, we established the lncRNA-miRNA-mRNA network in PCOS based on the ceRNA theory. We aimed to provide novel insights into the lncRNA-miRNA-mRNA interaction network and its related pathways in PCOS, and to provide theoretical perspectives to inform future research on PCOS occurrence and development.

We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/atm-21-2696).

Methods

Microarray data

The National Center for Biotechnology Information (NCBI)’s Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/) was searched for gene expression profile datasets of MGCs from women with PCOS and healthy women without PCOS. The GEO repository is a public database that stores chip, second-generation sequencing, and other high-throughput sequencing data. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The keywords for the search included “polycystic ovary syndrome”, “PCOS”, “granulosa cells”, and “homo sapiens”. GSE84376 and GSE106724 were selected for further analyses (). The GSE84376 dataset (platform: GPL16384, Affymetrix Multispecies miRNA-3 Array) contains 13 PCOS samples and 15 non-PCOS samples. The GSE106724 dataset (platform: GPL21096, Affymetrix Human Genome U133 Plus 2.0 Array Affymetrix Human Clariom D Assay) consists of 8 PCOS samples and 4 non-PCOS samples.

Figure 1

Flowchart of data processing and analysis. LncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of differentially expressed mRNAs and non-coding RNAs

We screened the differentially expressed mRNAs and non-coding RNAs between the PCOS and non-PCOS samples using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r). We defined RNAs with |log FC| >1 and P<0.05 as being differentially expressed.

Prediction and construction of the miRNA-mRNA network

Target genes of differentially expressed miRNAs were predicted using the DIANA-microT-CDS algorithm (http://www.microrna.gr/microT-CDS) (13). Target genes and miRNAs show contrasting alteration patterns in their differential expression. Cytoscape (http://www.cytoscape.org/) was employed to visualize the miRNA-mRNA regulatory network.

Constructing the lncRNA-miRNA-mRNA ceRNA network

LncRNA-miRNA interactions were predicted using DIANA-LncBase v2 (14). The lncRNAs and the miRNAs are differentially expressed with opposite expression alteration patterns. Then, the co-expressed competing triplets were assembled to obtain differentially expressed lncRNA-miRNA-mRNA networks Cytoscape was used for network visualization.

Functional and pathway enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation of the differentially expressed genes (DEGs) related to PCOS was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) (15). GO analysis examines genes with respect to their molecular functions (MFs), biological processes (BPs), and cellular components (CCs). We defined P<0.05 as indicating statistical significance for both GO terms and KEGG pathways.

Protein-protein interaction (PPI) network generation and module analysis

The PPI network of the DEGs was analyzed by STRING (https://string-db.org/). The cut-off standard was defined as >0.4. We visualized the PPI network using Cytoscape. Hub genes were ranked by Cytoscape cytoHubba (16). Cytoscape MCODE was applied to screen modules within the PPI network and identify the central genes among the DEGs (degree cutoff =2, max. depth =100, k-core =2, and node score cutoff =0.2) (16). The MCODE application in Cytoscape was utilized to examine the modules of the PPI network and identify the central genes among the DEGs (standard: degree cutoff =2, max. depth =100, k-core =2, and node score cutoff =0.2). Modules with anMCODE score >3 and nodes >5 were selected.

Results

Differentially expressed lncRNAs, miRNAs, and mRNAs

A total of 462 mRNAs, 2,464 lncRNAs, and 55 miRNAs were differentially expressed between the MGCs of patients with PCOS and those of healthy controls. Of them, 213 mRNAs, 392 lncRNAs, and 52 miRNAs were over-expressed, while 249 mRNAs, 2,072 lncRNAs, and 3 miRNAs were under-expressed.

Functional enrichment and KEGG pathway analysis of the DEGs

lists the GO terms and KEGG pathways of DEGs in GSE106724. The most significantly enriched BPs were cell adhesion, response to wounding, ascending aorta morphogenesis, angiogenesis, and negative regulation of insulin receptor signaling pathway. The KEGG pathways showing the most significant enrichment were retinol metabolism, drug metabolism—cytochrome P450, malaria, the Hippo signaling pathway, and glycine, serine and threonine metabolism.

Table 1

Top five GO and KEGG enrichment of DEGs in GSE106724

Category	Term	Gene function/pathway	Count	P value
BP	GO:0007155	Cell adhesion	21	2.0E–3
	GO:0009611	Response to wounding	7	2.1E–3
	GO:0035910	Ascending aorta morphogenesis	3	2.6E–3
	GO:0001525	Angiogenesis	13	3.0E–3
	GO:0046627	Negative regulation of insulin receptor signaling pathway	5	3.1E–3
CC	GO:0016021	Integral component of membrane	139	1.6E–3
	GO:0017053	Transcriptional repressor complex	6	6.5E–3
	GO:0005886	Plasma membrane	110	8.7E–3
	GO:0005667	Transcription factor complex	11	9.1E–3
	GO:0005884	Actin filament	6	1.3E–2
MF	GO:0031727	CCR2 chemokine receptor binding	3	1.3E–3
	GO:0042803	Protein homodimerization activity	28	3.6E–3
	GO:0005089	Rho guanyl-nucleotide exchange factor activity	7	5.8E–3
	GO:0008373	Protein kinase A binding	4	7.1E–3
	GO:0051018	Sialyltransferase activity	4	7.1E–3
KEGG pathway	hsa00830	Retinol metabolism	5	3.7E–2
	hsa00260	Glycine, serine and threonine metabolism	4	4.1E–2
	hsa00982	Drug metabolism—cytochrome P450	5	4.5E–2
	hsa05144	Malaria	4	7.2E–2
	hsa04390	Hippo signaling pathway	7	7.8E–2

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function.

PPI network of the DEGs and module analysis

The PPI network of the 462 DEGs obtained from the GSE106724 dataset was established using STRING. It was composed of 420 nodes and 556 edges (). By overlapping the top 20 genes according to the degree and maximal clique centrality (MCC) methods in cytoHubba, we identified CD163, MRC1, VSIG4, MSR1, UBR4, CCL2, CYBB, CCR2, SPP1, ESR1, LPL, and MGAM as hub genes.

Figure 2

PPI analysis. (A) The PPI network of DEGs from the GSE106724 dataset. Disconnected nodes in the network are hidden; (B) module 1 (score: 5); (C) module 2 (score: 4.8); (D) module 3 (score: 3.333); and (E) module 4 (score: 3.143). PPI, protein-protein interaction; DEGs, differentially expressed genes.

Four clusters from the PPI network were identified (). Module 1 comprised 23 nodes and 55 edges; module 2 comprised 6 nodes and 12 edges; module 3 consisted of 4 nodes and 5 edges; and module 4 contained 8 nodes and 11 edges. Enrichment analyses of module 1 demonstrated that the most significantly enriched BPs were protein ubiquitination involved in ubiquitin-dependent protein catabolic process, inflammatory response, and positive regulation of monocyte chemotaxis. For module 1, no KEGG pathway met the criterion for statistically significant enrichment (P<0.05). The most significantly enriched BPs in module 2 included oxygen transport, bicarbonate transport, blood coagulation, hydrogen peroxide catabolic process, and positive regulation of cell death. The KEGG pathways in module 2 were mainly enriched in African trypanosomiasis and malaria.

Analysis of the miRNA-mRNA network

We screened out 55 differentially expressed miRNAs from the GSE84376 dataset. Among them, 52 miRNAs were upregulated and 3 miRNAs were downregulated. There were 91 overlapping genes between the predicted target genes of these differentially expressed miRNAs and the DEGs from the GSE106724 dataset. These genes were related to 14 miRNAs. A miRNA-mRNA regulatory network comprising 124 miRNA-mRNA pairs was established (). The most enriched BPs included negative regulation of DNA binding, negative regulation of apoptotic process, fat pad development, protein ubiquitination involved in ubiquitin-dependent protein catabolic process, and the glucocorticoid receptor signaling pathway. No KEGG pathway showed significant enrichment.

Figure 3

miRNA-mRNA network. Rectangular nodes represent miRNAs, and oval-shaped nodes represent mRNA. Green represents downregulation and red shows co-upregulation. miRNA, microRNA; mRNA, messenger RNA.

Analysis of the lncRNA-miRNA-mRNA ceRNA network

We created a lncRNA-miRNA-mRNA ceRNA network on the basis of the relationships among differentially expressed lncRNAs, miRNAs, and mRNAs. As can be seen in , the ceRNA network contains 71 lncRNA nodes, 14 miRNA nodes, 69 mRNA nodes, and 167 edges.

Figure 4

The lncRNA-miRNA-mRNA ceRNA network. Diamond-shaped nodes represent lncRNAs; rectangular nodes represent miRNAs; and oval-shaped nodes represent mRNAs. Green represents downregulation and red shows co-upregulation. (A) The whole lncRNA-miRNA-mRNA ceRNA network; (B) the subnetworks of five hub miRNAs with their first neighbors. LncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA.

DEGs in the ceRNA network were subjected to functional analysis with DAVID. The most significantly enriched BPs included negative regulation of apoptotic process, protein ubiquitination involved in ubiquitin-dependent protein catabolic process, DNA repair, osteoblast development, and cellular sodium ion homeostasis. lists the top five GO terms. No KEGG pathway met the criterion for statistically significant enrichment (P<0.05).

Table 2

Top five GO enrichment of DEGs in lncRNA-miRNA-mRNA ceRNA network

Category	Term	Gene function/pathway	Count	P value
BP	GO:0043066	Negative regulation of apoptotic process	7	4.0E–3
	GO:0042787	Protein ubiquitination involved in ubiquitin-dependent protein catabolic process	4	1.4E–2
	GO:0006281	DNA repair	4	4.4E–2
	GO:0002076	Osteoblast development	2	5.5E–2
	GO:0006883	Cellular sodium ion homeostasis	2	6.2E–2
CC	GO:1903561	Extracellular vesicle	3	1.2E–2
	GO:0016324	Apical plasma membrane	5	1.6E–2
	GO:0005829	Cytosol	19	2.1E–2
	GO:0015629	Actin cytoskeleton	4	3.7E–2
	GO:0016020	Membrane	13	5.9E–2
MF	GO:0017080	Sodium channel regulator activity	3	5.9E–3
	GO:0005546	Phosphatidylinositol-4,5-bisphosphate binding	3	1.9E–2
	GO:0010314	Phosphatidylinositol-5-phosphate binding	2	4.3E–2
	GO:0070300	Phosphatidic acid binding	2	5.6E–2
	GO:0005545	1-phosphatidylinositol binding	2	7.0E–2

GO, Gene Ontology; DEGs, differentially expressed genes; lncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA; BP, biological process; CC, cellular component; MF, molecular function.

By applying the degree and MCC methods in cytoHubba, we discovered the top five hub miRNAs to be hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-1231, hsa-miR-3609, and hsa-miR-4433b-3p. The subnetworks of these five hub miRNAs and their first neighbors are presented in . We picked the subnetworks of four hub miRNAs (hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-3609, and hsa-miR-4433b-3p) for further analysis. These four miRNAs interacted with 34 lncRNAs and 42 mRNAs. The most enriched BPs included negative regulation of apoptotic process, osteoblast development, protein ubiquitination involved in ubiquitin-dependent protein catabolic process, pituitary gland development, and chondrocyte differentiation. The most enriched CCs included cytosol and extracellular vesicles. The most enriched MFs were sodium channel regulator activity and protein binding. No KEGG pathway met the criterion for statistically significant enrichment (P<0.05).

Discussion

PCOS constitutes one of the most prevalent endocrine disorders affecting women of childbearing age. Although PCOS is etiologically multifactorial, genetic abnormality is thought to be a key cause. Recent microarray studies have identified some differentially expressed lncRNAs, miRNAs, and mRNAs as being involved in PCOS. Based on the ceRNA hypothesis, lncRNAs are referred to as miRNA sponges, owing to their ability to bind and sequester miRNAs in order to exert a regulatory effect (17,18). Thus, in-depth insight into lncRNA-miRNA-mRNA interactions will improve the understanding of the occurrence and development of PCOS.

In this study, we identified a total of 462 mRNAs, 2,464 lncRNAs, and 55 miRNAs as being differentially expressed between the MGCs of women with PCOS and those of healthy women. Among the DEGs, CD163, MRC1, VSIG4, MSR1, UBR4, CCL2, CYBB, CCR2, SPP1, ESR1, LPL, and MGAM were identified as hub genes. CD163 is an endocytic scavenger receptor that serves as a marker of macrophage activation in adipose tissue. Circulating soluble CD163 is frequently associated with obesity and metabolic disorders (19). Møller et al. found an elevated serum level of soluble CD163 to be a predictive factor for type 2 diabetes mellitus in an age- and body mass index-independent manner (20). MRC1 is a conserved replication fork factor that stabilizes stalled forks and plays an essential role in checkpoint signal transmission (21). When the S-phase checkpoint becomes activated following DNA damage, MRC1 is critical for cell-cycle reentry (22). VSIG4 is a complement receptor belong to the immunoglobulin superfamily and is only expressed by a subset of macrophages which reside in tissues. A protective effect of the VSIG4-Fc fusion protein has been reported in experiments on the development of autoimmune arthritis, uveoretinitis, and hepatitis (23). MSR1 is a macrophage scavenger receptor that has been reported to participate in the maintenance of immunological tolerance (24). UBR4 belongs to the mammalian N-recognin family. Tasaki et al. found that mice with UBR4 deficiency died during embryogenesis and exhibited pleiotropic abnormalities, such as vascular development impairment in the yolk sac (25). CCL2 is a chemokine which promotes cell migration by inducing epithelial-mesenchymal transition in oral squamous cell carcinoma (26), and its receptor is CCR2 (27). SPP1 is a secreted non-collagenous protein that acts as a cytokine. It is likely important to cell-matrix interaction, and can regulate several cell behaviors and modulate leukocyte migration (28). ESR1 gene polymorphisms have been found to be involved in the phenotype of complications of PCOS (29). LPL is a key lipolysis-associated factor which is associated with obesity. Finally, MGAM appears to have great diagnostic value for type 2 diabetes mellitus (30).

Our lncRNA-miRNA-mRNA ceRNA network comprises 71 under-expressed lncRNA nodes, 14 over-expressed miRNA nodes, and 69 under-expressed mRNA nodes. The top five hub miRNAs included hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-1231, hsa-miR-3609, and hsa-miR-4433b-3p. Wang et al. found that hsa-miR-3135b in granulosa cells was significantly negatively correlated with follicle-stimulating hormone. Receiver operating characteristic (ROC) curve analysis revealed that the accuracy of hsa-miR-3135b in diagnosing PCOS was 0.760, with a sensitivity of 62.5% and a specificity of 85.7% (31). Zhang et al. found that miR-1231 was expressed at reduced levels in human glioma tissues and had a negative correlation epidermal growth factor receptor (EGFR) expression. Further, luciferase assay verified that miR-1231 directly targets EGFR, and when overexpressed, EGFR inhibited the suppression of the PI3K/AKT pathway and G1 arrest induced by miR-1231 (32). At present, little research has been done on hsa-miR-4649-3p, hsa-miR-3609, and hsa-miR-4433b-3p. This study also identified several 41 differentially expressed lncRNAs which interact with hub miRNAs in the ceRNA network, including NONHSAT018285.2, NONHSAT142592.2, NONHSAT112972.2, and NONHSAT012542.2; however, little research has been done on these lncRNAs.

Although our study provides important information for elucidating the development of PCOS and provides a theoretical basis for future research, these mechanisms have not been experimentally confirmed; further in vivo and in vitro experiments are needed to validate our findings.

Conclusions

This study identified a novel lncRNA-miRNA-mRNA network based on the ceRNA mechanism in PCOS. Some novel genes and non-coding RNAs that may be involved in the occurrence and development of PCOS were excavated, including CD163, MRC1, VSIG4, CCL2, CCR2, SPP1, hsa-miR-3135b, hsa-miR-4649-3p, hsa-miR-1231, hsa-miR-3609, and hsa-miR-4433b-3p. We have unearthed important information to elucidate the development of PCOS, which may provide a theoretical basis for future studies. However, in vivo and in vitro experiments are needed to validate our findings.

The article’s supplementary files as