Identification and bioinformatics analysis of overlapping differentially expressed genes in depression, papillary thyroid cancer and uterine fibroids.

It is hypothesized that there may be common characteristics between the genetic regulatory networks of different diseases. To identify these potential similarities, analysis of overlapping differentially expressed genes (DEGs) in several diseases, which are believed to be associated in traditional Chinese medicine (TCM) was performed in the present study. The gene expression profiles associated with depression, papillary thyroid carcinoma (PTC) and uterine fibroids (UF) were preliminarily analyzed using Gene Expression Omnibus 2R tools. Gene Ontology enrichment analysis, Kyoto Encyclopedia of Genes and Genomes pathway analysis and protein-protein interaction network analysis of the overlapping DEGs in depression, PTC and UF was performed. The results indicated that multiple genes, including activating transcription factor 3 and WSC domain containing 2 and the phosphoinositide 3 kinase/protein kinase b signaling pathway and its downstream effectors may be common factors associated with depression, PTC and/or UF. The neuroendocrine functions of the hypothalamic-pituitary-ovarian axis and hypothalamic-pituitary-thyroid axis were also identified as being mutually associated with depression, PTC and/or UF. However, due to the limitations of DNA microassays, it is recommended that future studies take epigenetics into consideration. Further transcriptomic, methylomic and metabolomic analyses of depression, PTC and UF are also required to identify and elucidate the key associated biomarkers. In conclusion, the results of the current study shed light on the potential genetic interconnections between depression, PTC and UF, which may be beneficial for understanding their underlying coregulatory mechanisms and contributing to the development of homeotherapy based on bioinformatics prediction.

Introduction

Molecular biology methods, including DNA microarrays, are being increasingly used for high-throughput gene expression analysis and to assess the impact of genetic variation on phenotypes. As well as for the diagnosis and prognostic prediction of a single disease, DNA microarray analysis may be used to determine similarities among multiple diseases via comparative analysis, particularly when combined with Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and protein-protein interaction (PPI) analyses (1,2).

Recently, the incidence rates of thyroid cancer (TC) in women have increased in several countries, including the United States and Australia (3). Although this may be associated with over-diagnosis, it is leading to an increase in the number of studies focusing on TC (4). Uterine fibroids (UF) are common benign neoplasms that occur in the uterus (5). Thyroid dysfunction may develop via the hypothalamic-pituitary-ovarian (HPO) axis due to heavy menstrual bleeding in females diagnosed with UF (6). Traditional Chinese Medicine (TCM) suggests that there is an association between TC and UF, and that these conditions are also associated with depression. However, to the best of our knowledge, there has been no systematic attempt to identify potential associations between differentially expressed genes (DEGs) associated with depression, TC and UF. Therefore, the present study analyzed the potential overlapping DEGs between these conditions using DNA microarray data, in order to biologically interpret their potential genetic associations.

In the present study, original datasets from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/), including gene expression data series GSE12654 (depression), GSE3678 [papillary TC (PTC)] and GSE593 (UF), were evaluated. Each dataset was analyzed using GEO2R (www.nci.nlm.nih.gov/geo/geo2r/) (7), which is a GEO online tool for analyzing microarray data. The data of patients with depression, PTC or UF were initially compared with corresponding control groups to identify DEGs. Subsequently, the overlapping DEGs between two or three of the pathologies were screened and then subjected to GO, KEGG, PPI network analyses using various databases. By combining bioinformatics analysis methods, the present study aimed to gain insight into the potential underlying genetic interconnections between depression, PTC and UF.

Materials and methods

Microarray data

Gene expression profiles obtained from DNA microarrays were downloaded from the NCBI GEO database. The following datasets were obtained: GSE12654 (depression, n=11; control, n=15), GSE3678 (7 PTC samples, 7 control samples) and GSE593 (5 UF samples, 5 control samples). The relevant datasets were retrieved under the disease Medical Subject Headings for Homo sapiens.

Identification of DEGs

Following analysis of the different gene expression profiles in GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/), analyzed data were downloaded in SOFT format, converted into XLS files and screened using Microsoft Office Excel 2017 (Microsoft Corporation, Redmond, WA, USA). Genes with a log2 (fold-change) ≥1.0 and a P-value <0.05 were identified as DEGs. DEGs for depression, PTC and UF were listed and the overlapping DEGs for the following combinations were subsequently screened: Depression and PTC; depression and UF; PTC and UF; and depression, PTC and UF (Fig. 1A).

Figure 1.

Venn diagram, Volcano plot and heat map. (A) Venn diagram based on the overlapping DEGs between depression, PTC and UF compared with controls. (B) Volcano plot representing the DEGs of PTC (GSE3678), satisfying the criteria of log2 (fold-change) value >2 or <-2 and P<0.05. Significantly expressed genes are represented as red dots. (C) Heat map of the 25 most upregulated and downregulated genes in PTC represented by the red dots in (B). Higher values represent upregulation. GSM85215-85221, control group; GSM85222-85228, PTC group. PTC, papillary thyroid carcinoma; UF, uterine fibroids; DEGs, differentially expressed genes.

Volcano plot and heat map analyses

Data normalization in microarray dataset analysis is essential to remove systematic variations (8). As a common statistical transformation, a Z-score was calculated in Microsoft Office Excel 2017 for each individual data point in a population using the following formula: Raw intensity-mean intensity/standard deviation (9). In the present study, the raw microarray data for PTC was transformed to obtain corrected expression intensities as aforementioned. Subsequently, each probe (gene) in PTC was represented by an individual dot in a volcano plot and a heat map of highly significantly expressed genes was drawn using HemI 1.0 software (http://hemi.biocuckoo.org/; Fig. 1B and C).

GO and KEGG analyses of overlapping DEGs

To identify enriched functionally associated gene groups, GO analysis was performed via the Database for Annotation, Visualization and Integrated Discovery, version 6.8 (david.ncifcrf.gov; Fig. 2) (10). As a widely applicable method, GO-based analysis aids in interpreting the biological functions and cellular components of genes and gene products (11). In addition, KEGG (www.genome.jp) provides a comprehensive database for assembling large-scale data on biosystems and for interpreting and defining molecular-level functions (12). Therefore, KEGG pathway analysis was selected via the Enrichr database (amp.pharm.mssm.edu/Enrichr/) to perform functional annotation of the overlapping DEGs of PTC and UF (Table I) (13). As few DEGs were identified between depression and the other two conditions (Fig. 1A), these data were unsuitable for GO, KEGG and PPI analyses.

Figure 2.

Gene ontology analysis of the overlapping differentially expressed genes between papillary thyroid carcinoma and uterine fibroids. -LN (FDR);=loge (napier constant)FDR.

Table I.

Kyoto Encyclopedia of Genes and Genomes pathway analysis of overlapping differentially expressed genes of papillary thyroid carcinoma and uterine fibroids.

Term	P-values	Genes
AGE-RAGE signaling pathway (hsa04933)	<0.0001	COL1A1; EGR1; THBD; JUN; COL3A1; STAT1; FN1; BCL2; COL4A5; SELE
Malaria (hsa05144)	<0.0001	HBB; SDC1; HBA2; CD36; SELE; MET
ECM-receptor interaction (hsa04512)	<0.0001	COL1A1; FN1; SPP1; COL4A5; SDC1; CD36
Axon guidance (hsa04360)	<0.0001	SEMA3D; SLIT1; PAK3; MET; EPHB1; EPHA3; GNAI1
Focal adhesion (hsa04510)	0.0002	COL1A1; JUN; FN1; SPP1; BCL2; COL4A5; PAK3; MET
Cytokine-cytokine receptor interaction (hsa04060)	0.0002	CCL14; BMP2; EDA; KIT; TNFRSF11B; LIFR; IL1RAP; INHBA; MET
Pathways in cancer (hsa05200)	0.0003	JUN; BMP2; STAT1; KIT; FN1; BCL2; COL4A5; FOS; MET; GNAI1; FGFR2
Complement and coagulation cascades (hsa04610)	0.0004	CFD; THBD; PROS1; CFI; TFPI
Protein digestion and absorption (hsa04974)	0.0007	COL1A1; COL3A1; COL13A1; COL4A5; PRSS2
Morphine addiction (hsa05032)	0.0007	GABRB3; GABBR2; PDE10A; PDE8B; GNAI1
Cell adhesion molecules (hsa04514)	0.0009	CLDN10; VCAN; CDH2; SDC1; NCAM1; SELE
Retinol metabolism (hsa00830)	0.0017	ADH1B; ALDH1A1; AOX1; HSD17B6
African trypanosomiasis (hsa05143)	0.0025	HBB; HBA2; SELE
Tryptophan metabolism (hsa00380)	0.0037	AOX1; CYP1B1; OGDHL
Osteoclast differentiation (hsa04380)	0.0038	JUN; STAT1; FOSB; TNFRSF11B; FOS

AGE-RAGE, advanced glycation end products-receptor for advanced glycation endproducts; ECM, extracellular matrix.

PPI network analysis

With large-scale data, the STRING database 10.0 (string.embl.de/cgi/) is typically used to consolidate known and predicted PPIs (14). In the present study, to screen the significantly interactive proteins, the DEGs with a minimum required interaction score ≥0.4 in STRING were selected. Molecular Complex Detection (MCODE) is an automated algorithm for identifying highly interconnected clusters in large PPI networks (15). The PPI networks were mapped via the MCODE App in Cytoscape software 3.4.0 (http://www.cytoscape.org/) (16).

Results

Fibrinogen-like protein 2 (FGL2); activating transcription factor 3 (ATF3); ADAM metallopeptidase domain 12 (ADAM12); WSC domain containing 2 (WSCD2); and mucin 1 (MUC1) were identified as overlapping DEGs associated with depression, PTC and UF (Fig. 1A). In addition to these five genes, dual specificity tyrosine-phosphorylation-regulated kinase 2 (DYRK2); RAB9B, member RAS oncogene family pseudogene 1; thyrotropin-releasing hormone receptor (TRHR); integrin α-6; and retinoschisin 1 were listed as the overlapping DEGs between depression and UF (Fig. 1A). Phosphodiesterase 5A (PDE5A); MYC induced nuclear antigen (MINA); angiotensin II receptor type 1 (AGTR1); and solute carrier family 14 member 1 (SLC14A1) were overlapping DEGs between depression and PTC (Fig. 1A). Notably, 157 DEGs for PTC and UF were identified. Of these, 152 directly overlapped between PTC and UF and 5 overlapped with PTC, UF and depression. This was a markedly higher number than the total overlapping DEGs between the depression and UF/PTC subsets (Fig. 1A). The volcano plot and heat map of DEGs associated with PTC are presented respectively in Fig. 1B and C. However, the heat maps of UF from the GSE593 dataset and depression from the GSE12654 dataset were not assessed as they have been previously studied (17,18).

GO and KEGG analyses

The DAVID database 6.8 was used to perform GO pathway analysis of the overlapping DEGs of PTC and UF. The GO term enrichment results all met the inclusion criteria of false discovery rate <0.01 and P<0.01. The top 3 GO results with the highest enrichment scores in each category (automatically generated by DAVID) were selected and presented in Fig. 2. In addition, the top 15 results of the KEGG pathway analysis with the criteria P<0.01 are presented in Table I, which also indicated the top enriched pathways associated with the overlapping DEGs of PTC and UF.

To screen the high confidence PPIs, a minimum required interaction score ≥0.4 in the STRING database was used as a basic criterion. The high confidence PPI network was then mapped in Cytoscape (Fig. 3). This PPI network indicated that JUN proto-oncogene (JUN); AP-1 transcription factor subunit; FOS proto-oncogene; AP-1 transcription factor subunit (FOS); BCL2, apoptosis regulator (BCL2); early growth response 1 (EGR1); bone morphogenetic protein 2 (BMP2); collagen type I alpha 1 chain; fibronectin 1; signal transducer and activator of transcription 1 (STAT1); KIT proto-oncogene receptor tyrosine kinase (KIT); and cadherin 2 were the top 10 significant hub nodes.

Figure 3.

Constructed PPI network for the overlapping differentially expressed genes between papillary thyroid carcinoma and uterine fibroids with PPI scores ≥0.4 in the Search Tool for the Retrieval of Interacting Genes/Proteins database. The top 10 hubs are represented as pink dots. PPI, protein-protein interaction.

Discussion

PTC is the predominant histological subtype of TC and accounts for >80% of all thyroid malignancies (19,20). In TCM it is believed that PTC is associated with ‘liver qi stagnation Zheng’ (21,22), which is thought to be primarily induced by psychological conditions, including depression and anxiety (23). In TCM, UF are thought to be induced by ‘blood stasis’ and are closely associated with ‘qi stagnation’ (24). However, it remains unknown whether there are common genetic factors between depression, PTC and UF. In the present study, the GSEs of microarray data obtained from patients with depression, PTC and UF were evaluated using GO, KEGG pathway and PPI network analyses.

DEGs of depression, PTC and UF were screened using the GEO2R tool in the GEO database. The results demonstrated a markedly higher number of overlapping DEGs between PTC and UF (157 genes) compared with that between PTC and depression (9 genes) and UF and depression (10 genes). Functional annotation implied that the DEGs between PTC and UF were primarily associated with cell adhesion/migration, extracellular structure organization, cell differentiation and cell proliferation/death. KEGG analysis manifested similar results and verified the mutual association of cell adhesion/migration-related genes. Subsequently, the 157 DEGs were screened in the STRING database according to interaction scores and the majority of proteins (100 genes) were identified to have mutual connections. A PPI network was constructed to reveal the intricate associations between these proteins. The 10 hub proteins were identified to serve a role in the following pathways: i) Immune response-relevant pathways, including toll-like receptor, B/T cell receptor and tumor necrosis factor pathways; ii) cancer-relevant pathways, including extracellular matrix (ECM)-receptor, janus kinase/signal transducers and activators of transcription (JAK-STAT) pathways; and iii) other pathways, including the phosphoinositide 3 kinase/protein kinase B (PI3K/AKT), cyclic adenosine monophosphate, advanced glycation endproducts-receptor for advanced glycation endproducts (AGE-RAGE) and mitogen-activated protein kinase 1 (MAPK) pathways. These data indicated the associations between the functions of the top 10 hub proteins.

Adenosine monophosphate-activated protein kinase signaling is upregulated in PTC (25). Additionally, the PI3K/AKT/mammalian target of rapamycin (mTOR), PI3K/AKT/forkhead box O (FOXO), PI3K/AKT/phosphatase and tensin homolog deleted on chromosome ten (PTEN) and MAPK pathways have been demonstrated to be associated with TC (26–28). The PI3K-AKT/mTOR pathway has also been identified as one of the most upregulated signaling pathways in UF (29). In addition, it has previously been demonstrated that sirtuin 1 (SIRT1) is significantly associated with major depressive disorder (MDD) (30) and SIRT1, 2 and 3 bind and induce FOXO protein deacetylation (31). However, whether PI3K/AKT and downstream factors, including mTOR, FOXO and PTEN, participate in interregulatory mechanisms between depression, PTC and UF by influencing transcription remains unknown. Furthermore, clarification on the role served by AGE-RAGE signaling in PTC and UF is required, although it has recently been suggested that it is associated with the augmentation of collagen production (32).

In the present study, cGMP-specific phosphodiesterase type 5 (PDE5A) was identified as an overlapping DEG between depression and PTC. The differential expression of PDE5A has been previously reported in PTC (33) and cGMP signaling inhibits MDD development in cases of neuroplasticity deficit (34). Genetic variation in AGTR1, which was also identified as a common DEG between depression and PTC, has been associated with depression and differences in frontotemporal morphology (35) and is a molecular marker for distinguishing between malignant and benign TC (36).

Of the overlapping DEGs associated with depression, PTC and UF in the present study, ATF3 directly interacts with several key cancer-associated proteins, including p53, E6 and androgen receptor (37) and similarly, MUC1 is aberrantly overexpressed in human breast cancer and other types of cancer (38,39). ATF3 is an endoplasmic reticulum stress transcription factor, which has also been demonstrated to be significantly associated with depression- and addiction-like behavior in rats (40). In addition, as one of the 25 most downregulated genes in PTC, WSCD2 was identified to be significantly associated with extraversion of personality traits in a meta-analysis of genome-wide association studies (41). The overlapping genes identified between depression and PTC/UF are a minority and may predominantly depend on the small number of DEGs for depression. However, analysis of these DEGs may help to elucidate the underlying influence of depression on PTC/UF.

In the PPI network analysis, among the most significant hub proteins identified in the network were JUN and FOS, which induce the transcription of proteins and are associated with cell activation, proliferation, differentiation and death (42). Further functions of AP-1 include inflammatory and apoptotic activities (43,44). The significant hub proteins BCL2 and STAT1 have been identified as closely associated with cancer (45,46) and the PPI network also revealed interactions with JUN, FOS, KIT, EGR1 and BMP2. The results of the present study indicated that the top 10 hub proteins may serve a role in the AGE-RAGE, PI3K-Akt, ECM-receptor, TNF and JAK-STAT signaling pathways, which corresponds with the results of the GO analysis (with regards to cell proliferation, migration, adhesion and death) and indicates potential associations between PTC and UF.

Regarding systemic disease, the impact of endocrine disorder (ED) on metabolic and physiological functions should not be understated, particularly when considering the characteristic ‘domino effect’ of ED. Hormones and their receptors and enzymes and transporter proteins serve a role in the intricate reactions of the endocrine system, including the hypothalamus-pituitary-thyroid (HPT) axis. In particular, hypothalamic thyrotropin-releasing hormone stimulates pituitary thyrotropin release and 80% of the triiodothyronine (T3) in the cerebral cortex is derived from the local conversion of thyroxine by deiodination (47). In the present study, TRHR was identified as an overlapping DEG in depression and UF. Active T3 exerts its function through nuclear thyroid hormone receptors (THR) and notably, deiodinase, THRA and THRB are expressed in the endometrium (48). In addition, there is a clear association between depression and thyroid function with regard to regulation by the HPT axis (49). Furthermore, interconnections between the hypothalamic-pituitary-ovarian (HPO) and HPT axes have been reported; notably hypothyroidism is closely associated with hyperprolactinemia, as well as menstrual problems, including oligomenorrhea, amenorrhoea and polymenorrhea (50). Therefore, an imbalance of hormones acting at the HPT and HPO axes may have adverse effects on mental, thyroid and/or uterine health (51).

In the present study, although only DNA microassay data analysis was used, it appeared that depression, PTC and UF may be co-associated with multiple signaling mechanisms and ED. However, due to the epigenetic effects of environmental exposure, differences in gene levels do not adequately explain the development of disease. Numerous factors, including genetic diathesis and environmental influences have been reported to influence the vulnerability of individuals with depression (52). UF has multiple symptoms, including abnormal uterine bleeding, infertility and urinary incontinence (53). Researchers have inferred that there is a lack of evidence to prove that total thyroidectomy and radioactive iodine therapy improves the survival rate of patients (54). Additionally, hysterectomy precludes future fertility and alternative treatments may also affect fertility (55). At least 50% of UF cases (56) and almost all microcarcinomas of PTC are considered clinically asymptomatic (57). Similarly, subclinical depression is highly prevalent and the therapeutic effects of drug and psychological treatment remain unclear (58). However, in cases of long-term depression, UF or PTC, apart from direct hormone supplement/regulation, it is unknown whether other alternative therapies aimed at regulating the neuroendocrine system or signaling networks are able to replace resective surgery such as thyroidectomy, myomectomy/hysterectomy. Recently, the identification of informative predictive biomarkers, which may be used to guide therapeutic programs for various types of cancer, has increased (59,60). In the future, it is likely that insight from population-based empirical study of TCM combined with bioinformatics analysis will elucidate the potential of alternative therapies for depression, UF and PTC.