Identification of differentially expressed genes in Cushing's disease by integrated bioinformatics analysis

doi:10.21203/rs.3.rs-29947/v1

Background

Cushing's disease is a rare and little-known disease, and the individualization of drug treatment varies greatly. Studies have shown that the gene expression profile of Cushing's disease is related to its clinical characteristics. Therefore, the study aims to identify key differential genes between the age and size of tumors through bioinformatics technology, thus providing a theoretical basis for personalized targeted therapy of Cushing's disease.

Methods

Downloading the gene expression microarray (GSE93825) data from the Gene Expression Omnibus (GEO) database and obtaining differentially expressed genes (DEGs) of different tumor sizes and ages through GEO2R. The DAVID database, Cytoscape and String platforms were utilized for functional enrichment analysis and protein-protein interaction (PPI) network analysis on selected differential genes.

Results

First, 96 DEGs were identified between macroadenoma (MAC) and microadenoma (MIC), which initially proved the different gene expression characteristics between them. Second, a total of 2128 DEGs were identified in MAC age group. The top five hub genes of the PPI network were GNGT2, LPAR3, PDYN, GRM3, and HTR1D. A total of 16 DEGs were identified in MIC age group. In addition, 88 DEGs were identified in younger MAC and MIC groups. The top five hub genes included LEP, PTGS2, STAT6, CXCL12, and ITPKB. 299 DEGs were identified in senior MAC and MIC groups. The first five hub genes were CCR7, LPAR2, CXCR5, ADCY3, and TAS2R14. By virtue of DAVID and Cytoscape software, the function enrichment analysis and core module analysis were performed successfully.

Conclusions

In summary, our research shows through bioinformatics analysis that different gene expression profiles of Cushing's disease are related to the size and age of the tumor, which may provide new insights into the molecular pathogenesis of Cushing's disease. These hub genes may be used for accurate diagnosis and treatment of Cushing's disease.

Translational Medicine

ACTH-secreting adenomas

Cushing's disease

Microarray

Differentially expressed genes

Bioinformatics

Cushing’s disease (CD) is a rare disease caused by pituitary adrenocorticotropic hormone (ACTH) adenoma, characterized by severe chronic hypercortisolemia, accounting for about 14% of all pituitary adenoma and about 70% of Cushing’s syndrome ^[1], most of which are microadenoma (MIC) less than 1 cm in diameter, and macroadenoma (MAC) accounts for only 10%–20% ^[2]. Patients with CD often have multiple serious complications due to hypercortisolemia, such as hypertension, diabetes, hyperlipidemia, osteoporosis, and mental depression. Untreated patients with CD often die from severe cardio-cerebrovascular diseases and severe infections ^[2]. The treatment of CD is mainly surgical treatment of trans nasal butterfly and ACTH pituitary adenoma ^[3], drug treatment ^[4], pituitary radiotherapy and bilateral adrenalectomy as adjuvant treatment methods. However, due to the complex etiology of CD, and the presence of patients who cannot be operated for some reason or need to use drugs to control blood cortisol concentration during radiation therapy, drug treatment is still a very important treatment ^[4]. Current drug treatments include pituitary tumor targeting drugs, glucocorticoid receptor antagonists, etc.

However, there is increasing evidence that ACTH-secreting adenoma has significant differences in response to the main and adjunct regulators of ATCH secretion, which obviously translates into different responses to these drugs ^[5]. In addition, studies have reported that somatic mutations and individual polymorphisms are related to differences in tumor size and secretion patterns ^[6,7]. The different gene expression profiles of ACTH-secreting adenoma are related to their clinical characteristics. Therefore, this study mainly evaluated the expression profiles of different genes in ACTH-secreting adenoma to analyze its correlation with clinical characteristics, such as tumor size and age, paving the way for personalized clinical treatment strategies.

In recent years, gene expression profile data has increased rapidly, and the use of bioinformatics methods to explore gene expression profile data has become a new research hotspot. As an effective technique for obtaining gene data on a large scale, gene expression microarrays have been used to collect gene chip expression profile data and study gene expression profiles in many human cancers. These microarrays provide a new method for studying tumor-related genes, and provide promising prospects for molecular prediction, drug-based molecular targeting and molecular therapy ^[8].

Therefore, in this study, we integrated and analyzed the data of gene expression profiling microarray (GSE93825) through a series of bioinformatics tools. Different gene expression profiles in ACTH-secreting adenoma were identified DEGs related to tumor size and age. Constructed a protein-protein interaction network and revealed the hub gene. In this way, we expect to find differentially expressed genes related to the clinical characteristics of tumors in CD, and thus provide a theoretical basis for personalized targeted therapy strategies.

Acquisition of Microarray data

DNA microarray is a new technology that can analyze the genome and gene expression characteristic map, including oligonucleotide chip, cDNA chip and genome chip, often divided into the following two modes: one is to fix the target DNA on the support, suitable for analyzing a large number of different target DNA. Another method is to fix a large number of probes on the support, which is suitable for analyzing different probe sequences of the same target DNA ^[9].

By searching for the keyword "ACTH-secreting adenoma" in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) belonging to the National Center for Biotechnology Information (NCBI), the present study obtained the gene expression profile data set with sequence number GSE93825; the platform of GSE93825 was GPL18281, Illumina HumanHT–12 WG-DASL V4.0 R2 expression beadchip, which included 11 MAC and 29 MIC specimens. The MAC age group included specimens of 6 younger and 3 senior patients; the MIC age group included specimens of 11 younger and 16 senior patients. The younger MAC and MIC group included 6 MAC and 13 MIC specimens; the older MAC and MIC group included 4 MAC and 14 MIC specimens. The platform and series matrix files had been downloaded as TXT files, and the data set information was shown in Table 1.

Table 1 Details of ACTH-secreting adenoma data in GEO

Sequence number of chip	GSE93825
Platform	GPL18281
Sample type	Pituitary human tissue
Sample	MACs and MICs
Younger group	20-40 years old
Senior group	40-60 years old
Reference	Cassarino MF et al (2018)

The gene expression profile data set with sequence number GSE93825 includes 11 MAC and 29 MIC specimens. MAC specimens include specimens of 6 younger and 3 senior patients; MIC specimens include specimens of 11 younger and 16 senior patients.

Screening of DEGs and drawing of heatmaps

R software was utilized to convert the downloaded platform and matrix files, and deleted the unqualified data. The data was calibrated, standardized and log2 converted. The ID corresponding to the probe name was converted to the corresponding gene symbol and saved in the TXT file. GEO2R online software (http://www.ncbi.nlm.nih.gov/geo/geo2r) was used for differential gene expression analysis. GEO2R allowed users to compare different sample groups in the GEO series in order to screen DEGs under the entire experimental conditions. Samples with a ANOVA P-value of < 0.05 and a logarithmic multiple change (FC) of 0.1 / 0.2 were considered DEGs, and TXT results were saved for subsequent analysis.

Finally, the pheatmap software package of the R software was used to construct a heatmap and highlighted the regions where the differential genes were mainly concentrated.

KEGG and GO enrichment analyses of DEGs

The database DAVID (https://david.ncifcrf.gov/) for annotation, visualization and integrated discovery provides a comprehensive set of gene and protein functional annotation information, which is an important foundation for the successful analysis of any high-throughput gene function for researchers to reveal the biological significance behind a large number of genes ^[10]. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding advanced functions and biological systems from large-scale molecular data sets generated through high-throughput experimental techniques ^[11]. Genome Ontology (GO) is a major bioinformatics tool used to annotate genes and analyze the biological processes of these genes^[12].

The DAVID online tool was utilized to analyze the function and pathway enrichment of the proteins encoded by DEGs and annotate these genes. In this study, we analyzed the significant up-regulation and down-regulation of DEGs determined from the comprehensive microarray “ATCH-secreting adenoma” data. P < 0.05 was considered statistically significant.

PPI network construction and module analysis

Functional protein-protein interaction (PPI) analysis is essential to explain the molecular mechanism of key cellular activities in the process of canceration. The STRING database (http://string-db.org/) is commonly used to identify interactions between known proteins and predicted proteins. The results come from experimental data, databases, text mining and predictive bioinformatics data ^[13]. The Cytoscape software’s plug-in molecular complex detection (MCODE) can be used to cluster a given network based on topology to find densely connected areas^[14].

In this study, the STRING online database was used to construct a PPI network of DEGs, and interactions with a combined score > 0.4 were considered statistically significant. Subsequently, MCODE was performed in Cytoscape software to screen important modules in the PPI network. The selection criteria were as follows: MCODE score > 5, degree cutoff = 2, node score cutoff = 0.2, maximum depth = 100, k score = 2. Hub genes were selected with a connectivity > 10. The gene enrichment analysis of genes in a single module was performed by DAVID, and the significance threshold was P < 0.05.

Identification of DEGs in ATCH-secreting adenoma

Identification of DEGs in MAC and MIC group. According to GSE93825, a total of 96 DEGs were selected, including 25 genes with up-regulated expression and 71 genes with down-regulated expression, which were statistically significant (ANOVA P-value < 0.05, log FC > 1). The top 45 DEGs were made into a cluster heatmap (Fig. 1a). The heatmap and hierarchical cluster of these 45 genes initially indicated different gene expression characteristics between MAC and MIC.

Identification of DEGs in MAC age group. According to GSE93825, a total of 2128 DEGs were selected, including 1643 up-regulated genes and 485 down-regulated genes, which were statistically significant (ANOVA P-value < 0.05, log FC > 2). The top 45 DEGs were made into a cluster heatmap (Fig. 1b).

Identification of DEGs in MIC age group. According to GSE93825, a total of 16 DEGs were selected, including 2 genes with up-regulated expression and 14 genes with down-regulated expression, which were statistically significant (ANOVA P-value < 0.05, log FC > 1). DEGs were shown in Table 2.

Table 2 DEGs in MIC age group

Gene symbol	Log FC	Gene symbol	Log FC
IL26	-5.583311	BRIP1	-2.80321
OR52N4	-4.856554	GSTTP2	-2.74332
GLB1L3	-4.436037	FAM119A	-2.3704
HFE	-3.604451	ZC3HAV1L	-2.02566
SPN	-3.213325	PIGW	-1.91017
SHROOM4	-3.202547	LILRA6	-1.67649
HIST1H2BH	-2.994321	C1ORF103	2.06169
MIR2116	-2.933099	MAN1B1	3.913303

A total of 16 DEGs were selected in MIC age group, including 2 genes with up-regulated expression and 14 genes with down-regulated expression; ANOVA P-value <0.05.

The heatmap and hierarchical cluster of these 45 genes and Table 2 indicated that different gene expression characteristics existed between MACs of different ages and MICs of different ages.

Identification of DEGs in younger MAC and MIC group. A total of 88 DEGs were selected according to GSE93825, including 86 genes with up-regulated expression and 2 genes with down-regulated expression, which had statistical significance (ANOVA P-value < 0.05, log FC > 1). The top 45 DEGs were made into a cluster heatmap (Fig. 1c).

Identification of DEGs in senior MAC and MIC group. According to GSE93825, a total of 299 DEGs were selected, including 216 up-regulated genes and 83 down-regulated genes, which were statistically significant (ANOVA P-value < 0.05, log FC > 1). The top 45 DEGs were made into a cluster heatmap (Fig. 1d).

The heatmaps and hierarchical clusters of these 90 genes clearly indicated that there were different gene expression characteristics between MAC and MIC of the same age group.

KEGG and GO enrichment analysis of DEGs

To analyze the biological significance of these genes, we used DAVID and R software to reveal functional description, classification and location of DEGs.

KEGG and GO enrichment analysis in MAC and MIC group. GO analysis results showed that DEGs in MAC and MIC group were mainly enriched in biological processes. In the biological process group, genes were mainly enriched in positive transcription regulation, cell proliferation, cell biosynthesis, and macromolecular metabolism. As for molecular function, genes were mainly abundant in binding, including growth factor binding and chromatin binding. In addition, cell components were mainly concentrated in extracellular regions (Fig. 2).

KEGG and GO enrichment analysis in MAC age group. GO analysis results showed that DEGs in MAC age group mainly enriched the biological processes of cell surface receptor-linked signal transduction and nucleic acid biosynthesis. In the cell group, genes were mainly enriched in the inherent components of the cell membrane. KEGG pathway analysis showed that DEGs were mainly enriched in neuroactive ligand-receptor interactions. These results were shown in Table 3.

Table 3 GO and KEGG pathway analysis of DEGs in MAC age group

Term	Description	Count					P-value
Upregulated
GO:0042981	Regulation of apoptosis	10					0.010548
GO:0043067	Regulation of programmed cell death	10					0.01121
GO:0010941	Regulation of cell death		10				0.011466
GO:0016021	Integral to membrane		38			0.020699
GO:0031224	Intrinsic to membrane		38			0.036694
GO:0007166	Cell surface receptor linked signal transduction		15			0.038692
GO:0046873	Metal ion transmembrane transporter activity			5		0.049871
hsa04080	Neuroactive ligand-receptor interaction			6	0.048516
Downregulated
GO:0009113	Purine base biosynthetic process			2		0.030846
GO:0006144	Purine base metabolic process			2		0.047579
GO:0046112	Nucleobase biosynthetic process			2		0.047579

Select the top two hundred DEGs with the most significant P-value for GO and KEGG pathway analysis.

KEGG and GO enrichment analysis in MIC age group. GO analysis results showed that in the cell group, DEGs in MIC age group were mainly enriched in the basic part of the cell. The results were shown in Table 4.

Table 4 GO enrichment analysis of DEGs in MIC age group

Term	Description	Count	P-value
Downregulated
GO:0045178	Basal part of cell	2	0.018628

Use the DAVID online tool to analyze DEGs; ANOVA P-value <0.05.

KEGG and GO enrichment analysis in younger MAC and MIC group. GO analysis results showed that the DEGs in younger MAC and MIC group were mainly enriched in biological processes, cell composition and molecular functions. Genes were mainly enriched in the biological processes of positive regulation of transcription, positive regulation of nucleic acid metabolism, positive regulation of gene expression, organic matter response, positive regulation of nitrogen compound metabolism, positive regulation of cell biosynthesis and transcription regulation of RNA polymerase II promoter. The changes in molecular functions were mainly concentrated on calcium ion binding. The cellular composition of DEGs changed mainly in the extracellular area (Fig. 3).

KEGG and GO enrichment analysis in senior MAC and MIC group. GO analysis results showed that the DEGs in senior MAC and MIC group were mainly enriched in biological processes and cell composition. The biological processes of gene enrichment included the regulation of phosphorylation, the regulation of phosphate metabolism, the regulation of enzyme catalytic activity and the positive regulation of molecular functions. As for molecular functions, genes were mainly abundant in helicase activity. In the cell group, genes were mainly enriched in the organelle membrane and membrane system (Fig. 4).

PPI Network Analysis and hub gene selection

PPI network analysis and hub gene selection in MAC and MIC group. Based on the STRING database, the PPI network of MAC and MIC group was constructed; and the module analysis was performed by MCODE in the Cytoscape software, and the top two modules were selected (Fig. 5a). Functional analysis showed that the important core modules were mainly enriched in cell growth, endocrine processes and coagulation cascade pathways [see Additional file 1]. The top five hub genes were PLG, DCN, FBLN1, RNF41 and ISL1.

PPI network analysis and hub gene selection in MAC age group. Based on the STRING database, the PPI network of MAC age group was constructed; and the module analysis was performed by MCODE in the Cytoscape software, and the top two modules were selected (Fig. 5b). Important key modules showed functions including protein-coupled receptor protein signaling pathways, cell-surface receptor-linked signal transduction and protein catabolism [see Additional file 1]. The top five hub genes were GNGT2, LPAR3, PDYN, GRM3 and HTR1D.

PPI network analysis and hub gene selection in younger MAC and MIC group. Based on the STRING database, the PPI network of younger MAC and MIC group was constructed; and the module analysis was performed by MCODE in the Cytoscape software, and the top two modules were selected (Fig. 6a). Important core modules demonstrated negative regulation of stimulus response, angiogenesis and vascular development functions [see Additional file 1]. The top five hub genes included LEP, PTGS2, STAT6, CXCL12 and ITPKB.

PPI network analysis and hub gene selection in senior MAC and MIC group. Based on the STRING database, the PPI network of senior MAC and MIC group was constructed; and the module analysis was performed by MCODE in the Cytoscape software, and the top two modules were selected (Fig. 6b). Important key modules displayed functions including protein-coupled receptor protein signaling pathways, cell surface receptor-linked signal transduction and mRNA metabolism processes [see Additional file 1]. The top five hub genes included CCR7, LPAR2, CXCR5, ADCY3, and TAS2R14.

ACTH-secreting adenoma is a rare and little-known disease. More and more studies have shown that there are large differences in the secretion parameters and the response to drug treatment between adenomas that secrete ACTH ^[5]. Therefore, studying the variability of these tumors is necessary to identify clinical features related to tumor size and age and to help develop targeted therapies. Microarray technology allows us to explore the genetic changes of adrenocorticotropic tumors, and has been proven to be a useful method to identify new biomarkers in other diseases ^[15].

Our research design was original. Previous studies have identified different gene expression profiles in pituitary adenomas (including adrenocorticotropic tumors), but have not considered the tumor size classification and age of them ^[5]. On the basis of this study, a comprehensive bioinformatics method was used to analyze DEGs between MACs of different ages, MICs of different ages, and MAC and MIC of the same age, respectively.

Heatmap analysis revealed DEGs between different groups. First, this study initially identified different gene expression characteristics of MACs and MICs. Among the genes with the largest differences, we highlighted PLG and ISL1, which were reduced by 7-fold and 5-fold respectively; CLTCL1 and RNF41 increased by 2-fold and 1.5-fold respectively.

Second, we identified genes that were selectively over-expressed and under-expressed in younger MACs that showed significantly different gene expression characteristics than senior MACs. Among genes with the most differences, we highlighted MCPH1 and TARS2, which were reduced by 14-fold and 7.5-fold, respectively, and SLC17A6 and SGTB increased by 12-fold and 9-fold, respectively. However, it may be due to insufficient sample size. In this study, only 16 differential genes were identified between younger MICs and senior MICs, and it was impossible to analyze difference in gene expression by heatmap. Among genes with the greatest differences, we highlighted IL26 and OR52N4, which were reduced by 5.5-fold and 5-fold, respectively, and C1ORF103 and MAN1B1 were increased by 4-fold and 2-fold, respectively. Based on the above results, this study initially proved that the age of tumor patients may lead to different gene expression profiles of ACTH.

Next, we separately identified that there were significantly different gene expression characteristics between younger MACs and MICs, and between senior MACs and MICs. Among genes with the greatest differences at younger age, we highlighted GNRHR and HIST1H2BE, which were increased 10-fold and 5-fold, respectively, and COL25A1 and TLE1 were decreased 4-fold and 2-fold, respectively. Among genes with the greatest differences at senior age, we highlighted AGMAT and FLJ35220, which were increased by 12-fold and 7-fold, respectively, and PRAM1 and CLCN2 were decreased by 6-fold and 4-fold, respectively. The above results preliminarily proved that tumor size may also be one of the factors that lead to different gene expression profiles of ACTH.

PPI network of MAC age group illustrated the overview of its functional connections. Hub genes were also selected: GNGT2, LPAR3, PDYN and GRM3. Perhaps due to insufficient sample size, only 16 differential genes were selected in MIC age group, and genes with the largest differences were selected: IL26 and MAN1B1.

G protein subunit γ transduction protein 2 (GNGT2) belongs to G protein γ family. GNGT2 may mediate β-inhibitor 1 to induce Akt phosphorylation and NF-κB activation^[16]. When activity of NF-κB increased, it may inhibit apoptosis and promote tumorigenesis, angiogenesis, invasion and metastasis. It has been reported in the literature that expression of GNGT2 in esophageal cancer patient samples and cell lines was significantly upregulated, thereby activating NF-κB pathway, and promoting proliferation of esophageal cancer cells ^[17].

In many tumor cells such as human colon cancer and ovarian cancer, LPA enhanced cell motility, metastasis and invasion ability through lysophosphatidic acid receptor 3 (LPAR3). In addition, LPAR3 may also participate in the regulation of protein phosphorylation, which was required for anti-apoptotic function ^[18].

Metabolic glutamate receptor 3 (GRM3) is an inhibitory molecule on the surface of B cells and a subtype of metabolic glutamate receptor group II. It is a G protein-coupled receptor that can inhibit adenylate cyclase system and reduce formation of cAMP after activation. The abnormal level of GRM3 is related to tumor cell apoptosis associated with B cells. Apoptosis plays a key role in occurrence and development of tumors and autoimmune diseases, and resistance and susceptibility of various therapeutic agents. It has been shown in the literature that GRM3 was involved in apoptosis of B-cell-related tumors such as multiple myeloma and B-cell leukemia; GRM3 may mediate apoptosis via Foxo1, effectively inhibiting growth of mouse myeloma cell line SP 2/0 ^[19]. The genetic abnormalities of GRM3 were also often described in renal cell carcinoma and esophageal cancer.

The MAN1B1 gene product is called ER α–1,2-mannosidase (ERManI), an enzyme located in the Golgi complex of mammalian cells. ERManI was target of miR–125b, and miR–125b was a microRNA that was frequently down-regulated in many types of cancer. Up-regulated ERManI may prevent inappropriate secretion of misfolded glycoproteins, maintain protein homeostasis in mammalian secretory pathways, and promote cancer cell proliferation, migration, and invasion ^[20].

Animal studies have shown that predynorphin (PDYN) affected cognition and memory, especially in elderly animals. PDYN gene polymorphism played a role in the memory function of the elderly ^[21]. Another study found that expression of PDYN increased with age; during normal aging, increased expression of PDYN reduced expression and signaling of metabolic glutamate receptor (mGluR), which impaired cognitive function and increased Anxiety ^[22]; In addition, blockade of mGluR may also affect growth and migration of various tumor cells. Therefore, PDYN may be a differential gene related to the age of ACTH.

The study found that in the pancreatic tumor cell line, interleukin (IL) 26 immune cells infiltrated, resulting in phosphorylation of key cancer signaling pathways ERK1 / 2 and STAT3 pathways, and increased colony formation of tumor cells ^[23]. In another study, serum IL–26 levels were positively correlated with gastric cancer and its clinicopathological stage. Serum IL–26 levels in patients with gastric cancer were significantly higher than those in benign stomach disease group, and gradually increased with increase of gastric cancer clinical stage. IL–26 was still an independent risk factor for gastric cancer after adjustment for age, gender, carcinoembryonic antigen (CEA), CA125 and other risk factors ^[24]. Therefore, difference in expression of IL–26 may also be related to different clinical stages of ACTH.

PPI network of younger MAC and MIC group illustrated the overview of its functional connections, and hub genes were selected: LEP, PTGS2, STAT6, CXCL12 and ITPKB. Similarly, senior MAC and MIC group selected hub genes: CCR7 and ADCY3.

Adenylate cyclase 3 (ADCY3) is a membrane-associated protein widely expressed in human tissues and can catalyze the formation of cyclic 3’, 5’-adenosine monophosphate (cAMP). The study found that there was a significant correlation between upregulation of ADCY3 and Lauren’s intestinal type gastric cancer. ADCY3 overexpression may exert its tumor-promoting effect through cAMP / PKA / CREB pathway ^[25].

Transcription activator–6 (STAT6) is a member of the family of signal transduction and transcription activators. STAT6 signaling pathway can promote proliferation of colon cancer cells, possibly by regulating expression of p21 and p27 genes. It has been reported in the literature that upregulation of long non-coding RNA RP11–468E2.5 may inhibit JAK / STAT signaling pathway by targeting STAT5 and STAT6 genes, thereby inhibiting cell proliferation and promoting apoptosis in colorectal cancer ^[26].

The presence of chemokines and their receptors indirectly or directly regulates tumor cell invasion and metastasis. At present, there are more researches on chemokine 21 (CCL21) and its receptor 7 (CCR7) ^[27], chemokine CXCL13 and its receptor CXCR5 ^[28] have been proven to regulate invasion and metastasis of prostate cancer, colon cancer, breast cancer, gastric cancer, renal cell carcinoma and other malignant tumors, tumor-related angiogenesis and tumor cell movement, etc. CXCR5 is a member of CXC receptor family of chemokines, with 7 transmembrane domains, coupled with the GTP2 protein to form a transmembrane receptor; CCR7 is also a seven-transmembrane structure coupled to a heterotrimeric G protein. CXCR5 and CCR7 are highly expressed in various malignant tumor cells.

Although difference in methylation curve between invasive and non-invasive non-functional pituitary adenoma (NFPA) was very small. However, the study found that abnormal epigenetic disorders related to the invasion of specific genes included up-regulation of ITPKB and down-regulation of CNKSR1 in aggressive tumors ^[29].

In Cushing’s syndrome, increase in body fat-related LEP levels was not related to source of hypercortisolism, and visceral fat may be its main source. Hypothalamic-pituitary-adrenal axis dysfunction did not directly affect diurnal variation of plasma LEP levels ^[30]. However, in this study, differential expression of LEP affected gene profiles between MACs and MICs. Therefore, it is not clear why our results differ from previously published results.

In addition, there have been reports in the literature that compared with male rectal cancer patients, female patients had higher levels of adipokine and LEP. LEP levels were higher in female patients after a significant weight loss. The distribution of fat factors in patients with rectal cancer had a gender difference ^[31]. It suggested that distribution of LEP in adrenocorticotropic tumors may also have gender differences. Other studies will expand and confirm the above results.

Prostaglandin epoxide synthase 2 (PTGS2), also known as COX–2, is a key enzyme in prostaglandin biosynthesis and is highly expressed in many benign precancerous lesions and malignant tumors in humans ^[32]. Under normal physiological conditions, most tissue cells did not express it, but in the course of pathological reactions such as inflammation or tumor, its expression was rapidly up-regulated by certain cytokines, growth factors, inflammatory mediators, cancer promoting factors, hypoxia, hormones and other stimulating factors ^[33]. Prostaglandin, the major product of PTGS2, may inhibit cell apoptosis, promote cell proliferation, suppress immune surveillance, promote angiogenesis and play a key role in occurrence and development of tumors ^[34].

Nonetheless, few studies have investigated its role in pituitary adenoma, and only one study observed that compared with normal pituitary tissue, the expression levels of cyclooxygenase isoforms (COX–1 and COX–2) in pituitary tumors including CD were increased. Compared with MICs and non-invasive pituitary tumors, in MACs and invasive tumors, the expression levels of COX–2 and its product PGE2 increased significantly ^[35]. This was consistent with our finding that PTGS2 may cause different gene expression profiles with different tumor sizes.

Pituitary adenomas produced chemokine stromal cell-derived factor (SDF–1α / CXCL12) and its receptor CXCR4. One study showed that under hypoxic conditions, CXCL12 and CXCR4 were up-regulated simultaneously. CXCR4 expression was significantly increased in pituitary adenoma and was positively correlated with Knosp grade, while CXCR4 expression was higher in MACs and adenomas produced by growth hormone (GH) ^[36]. Therefore, CXCL12 gene may be differentially expressed between MACs and MICs.

Some limitations should be recognized in this study. The chip release time and the number of samples may be the biggest defects. Over time, the differential genes of CD may change for some reason. Besides, due to the undisclosed data, this study has not yet analyzed its clinical parameters and prognosis. Further molecular biology experiments are needed to confirm that different gene expression profiles of CD are related to tumor size and age.

In summary, our research indicated that the differential gene expression profiles of CD were related to their clinical characteristics, such as tumor size and age, through bioinformatics analysis of gene expression microarrays, which may provide a theoretical basis for the future study of molecular pathophysiology of CD and the differential response of drug therapy between individuals, and pave the way for targeted precision therapy.

CD: Cushing’s disease; GEO: Gene Expression Omnibus; DEGs: Differentially expressed genes; PPI: Protein-protein interaction; MAC: Macroadenoma; MIC: Microadenoma; ACTH: Adrenocorticotropic hormone; NCBI: National Center for Biotechnology Information; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Genome Ontology; MCODE: Molecular complex detection; ERManI: ER α–1,2-mannosidase; mGluR: Metabolic glutamate receptor; CEA: Carcinoembryonic antigen; cAMP: Cyclic adenosine monophosphate; NFPA: Non-functional pituitary adenoma; GH: Growth hormone.

Ethics approval and consent to participate

According to the database strategy, this study obtained access to de-identified linked datasets from the GEO database. For analyses of de-identified data from the GEO database, institutional review board approval and informed consent were not required.

Consent for publication

Not applicable.

Availability of data and materials

The dataset supporting the conclusions of this article is available in the GEO repository, [unique persistent identifier and hyperlink to dataset in http://www.ncbi.nlm.nih.gov/geo/ format]. GSE93825 dataset was downloaded from GEO and expression profiling arrays were generated using GPL18281, Illumina HumanHT–12 WG-DASL V4.0 R2 expression beadchip. We applied GEO2R online software (http://www.ncbi.nlm.nih.gov/geo/geo2r), DAVID online tool (https://david.ncifcrf.gov/), STRING database (http: // string- db.org/) and Cytoscape software to analyze and identify different gene expression profiles of Cushing’s disease.

All data generated or analyzed during this study are included in this published article.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was supported by National Natural Science Foundation of China (Grant number: 81670456), Major science and technology project for ‘Significant New Drugs Creation’ (2018ZX09711001–003–011), Innovation Fund for Medical Sciences (Grant No: 2017-I2M–1–011), and Beijing Municipal Natural Science Foundation (Grant number: 7162132).

Authors’ contributions

LTT conceived and executed the project, analyzed the research data and wrote the manuscript. LTT, WJR, YHY and RSY participated in conceptual design, data collection and basic analysis of the project. SL and DGH were involved in conceptual design of the project, research supervision, manuscript revision, financial support and final approval of the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank all participants of the group in the National center for Pharmaceutical Screening, Institute of Materia Medica, Chinese Academy of Medical Science and Peking Union Medical College.

Asa SL, Ezzat S. The pathogenesis of pituitary tumors. Annu Rev Pathol. 2009;4:97–126.
Kreitschmann-Andermahr I, Psaras T, Tsiogka M, Starz D, Kleist B, Siegel S, et al. From first symptoms to final diagnosis of Cushing’s disease: experiences of 176 patients. Eur J Endocrinol. 2015;172:285–9.
Biller BM, Grossman AB, Stewart PM, Melmed S, Bertagna X, Bertherat J, et al. Treatment of adrenocorticotropin-dependent Cushing’s syndrome: a consensus statement. J Clin Endocrinol Metab. 2008;93(7):2454–62.
Lacroix A, Feelders RA, Stratakis CA, Nieman LK. Cushing’s syndrome. Lancet. 2015;386(9996):913–27.
Cassarino MF, Ambrogio AG, Cassarino A, Terreni MR, Gentilini D, Sesta A, et al. Gene expression profiling in human corticotroph tumours reveals distinct, neuroendocrine profiles. J Neuroendocrinol. 2018;30(9):e12628.
Nakano-Tateno T, Tateno T, Hlaing MM, Zheng L, Yoshimoto K, Yamada S, Asa SL, Ezzat S. FGFR4 polymorphic variants modulate phenotypic features of Cushing disease. Mol Endocrinol. 2014;28:525‐33.
Perez-Rivas LG, Theodoropoulou M, Ferraù F, Nusser C, Kawaguchi K, Stratakis CA, et al. The gene of the ubiquitin‐specific protease 8 is frequently mutated in adenomas causing Cushing’s disease. J Clin Endocrinol Metab. 2015;100:E997‐E1004.
Petryszak R, Burdett T, Fiorelli B, Fonseca NA, Gonzalez-Porta M, Hastings E, et al. Expression Atlas update-a database of gene and transcript expression from microarray- and sequencing-based functional genomics experiments. Nucleic Acids Res. 2014;42(Database issue):D926-D932.
Marzancola MG, Sedighi A, Li PC. DNA microarray-based diagnostics. Methods Mol Biol. 2016;1368:161–78.
Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, Stephens R, Baseler MW, Lane HC, Lempicki RA. The DAVID Gene Functional Classification Tool: A novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8:R183.
Kanehisa M. The KEGG database. Novartis Found Symp. 2002;247:91–252.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.
von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258–61.
Bandettini WP, Kellman P, Mancini C, Booker OJ, Vasu S, Leung SW, Wilson JR, Shanbhag SM, Chen MY, Arai AE. MultiContrast Delayed Enhancement (MCODE) improves detection of subendocardial myocardial infarction by late gadolinium enhancement cardiovascular magnetic resonance: A clinical validation study. J Cardiovasc Magn Reson. 2012;14:83.
Kulasingam V, Diamandis EP. Strategies for discovering novel cancer biomarkers through utilization of emerging technologies. Nat Clin Pract Oncol. 2008;5(10):588–99.
Tang H, Wei P, Duell EJ, Risch HA, Olson SH, Bueno-de-Mesquita HB, et al. Genes-environment interactions in obesityand diabetes-associated pancreatic cancer: A GWAS data analysis. Cancer Epidemiol Biomarkers Prev. 2014;23(1):98–106.
Liu GM, Ji X, Lu TC, Duan LW, Jia WY, Liu Y, Sun ML, Luo YG. Comprehensive multi-omics analysis identified core molecular processes in esophageal cancer and revealed GNGT2 as a potential prognostic marker. World J Gastroenterol. 2019;25(48):6890‐901.
Hua X, Chen J, Wu L. Identification of candidate biomarkers associated with apoptosis in melanosis coli: GNG5, LPAR3, MAPK8, and PSMC6. Biosci Rep. 2019;39(1):BSR20181369.
Liu X, Zhang Y, Wang Z, Wang X, Zhu G, Han G, et al. Metabotropic glutamate receptor 3 is involved in B-cell-related tumor apoptosis. International Journal of Oncology. 2016;49(4):1469–78.
Pan S, Cheng X, Chen H, Castro PD, Ittmann MM, Hutson AW, Zapata SK, Sifers RN. ERManI Is a Target of miR–125b and Promotes Transformation Phenotypes in Hepatocellular Carcinoma (HCC). Plos One. 2013;8(8):e72829.
Kölsch H, Wagner M, Bilkei-Gorzó A, Toliat MR, Pentzek M, Fuchs A. Gene polymorphisms in prodynorphin (PDYN)are associated with episodic memory in the elderly. J Neural Transm (Vienna). 2009;116(7):897–903.
Ménard C, Tse YC, Cavanagh C, Chabot JG, Herzog H, Schwarzer C, Wong TP, Quirion R. Knockdown of Prodynorphin Gene Prevents Cognitive Decline, Reduces Anxiety, and Rescues Loss of Group 1 Metabotropic Glutamate Receptor Function in Aging. Journal of Neuroscience. 2013;33(31):12792–804.
Mayer P, Linnebacher A, Glennemeier-Marke H, Marnet N, Bergmann F, Hackert T, et al. The Microarchitecture of Pancreatic Cancer as Measured by Diffusion-Weighted Magnetic Resonance Imaging Is Altered by T Cells with a Tumor Promoting Th17 Phenotype. International Journal of Molecular Sciences. 2020;21(1):346.
Xue TH, Yang J, Song P, Zhou G. Investigation on correlations of serum IL–26 with diagnosis and staging of gastric cancer. J BUON. 2019;24(1):215–20.
Hong SH, Goh SH, Lee SJ, Hwang JA, Lee J, Choi IJ, et al. Upregulation of adenylate cyclase 3 (ADCY3) increases the tumorigenic potential of cells by activating the CREB pathway. Oncotarget. 2013;4(10):1791–803.
Jiang L, Zhao XH, Mao YL, Wang JF, Zheng HJ, You QS. Long non-coding RNA RP11–468E2.5 curtails colorectal cancer cell proliferation and stimulates apoptosis via the JAK/STAT signaling pathway by targeting STAT5 and STAT6. J Exp Clin Cancer Res. 2019;38(1):465.
Maolake A, Izumi K, Natsagdorj A, Iwamoto H, Kadomoto S, Makino T, et al. Tumor necrosis factor-α induces prostate cancer cell migration in lymphatic metastasis through CCR7 upregulation. Cancer Sci. 2018;109(5):1524‐31.
Hussain M, Adah D, Tariq M, Lu Y, Zhang J, Liu J. CXCL13/CXCR5 signaling axis in cancer. Life Sci. 2019;227:175‐86.
Kober P, Zbijewska J, Rusetska N, Maksymowicz M, Goryca K, Kunicki J, Bonicki W, Siedlecki JA, Bujko M. DNA methylation profiling in nonfunctioning pituitary adenomas. Molecular and Cellular Endocrinology. 2018;473:194–204.
Robaczyk MG. [Evaluation of leptin levels in plasma and their reliance on other hormonal factors affecting tissue fat levels in people with various levels of endogenous cotisol]. Ann Acad Med Stetin. 2002;48:283‐300.
Florescu A, Branisteanu D, Bilha S, Scripcariu D, Florescu I, Scripcariu V, Dimofte G, Grigoras I. Leptin and adiponectin dynamics at patients with rectal neoplasm - Gender differences. PLoS One. 2019;14(8):e0212471.
Khorshidi F, Haghighi MM, Nazemalhosseini ME, Azimzadeh P, Damavand B, Vahedi M, et al. The prostaglandin synthase 2/cyclooxygenase 2 (PTGS2/COX2) rs5277 polymorphism does not influence risk of colorectal cancer in an Iranian population. Asian Pac J Cancer Prev. 2014;15(8):3507‐11.
Cao Y, Prescott SM. Many actions of cyclooxygenase–2 in cellular dynamics and in cancer. Journal of Cellular Physiology. 2002;190(3):279–86.
Koji K, Keisuke O, Takahisa M. Role of prostaglandins in tumor microenvironment. Cancer & Metastasis Reviews. 2018;37(2–3):347–54.
Akbari N, Ghorbani M, Salimi V, Alimohammadi A, Khamseh ME, Akbari H. Cyclooxygenase enzyme and PGE2 expression in patients with functional and non-functional pituitary adenomas. BMC Endocr Disord. 2020;20(1):39.
Yoshida D, Koketshu K, Nomura R, Teramoto A. The CXCR4 antagonist AMD3100 suppresses hypoxia-mediated growth hormone production in GH3 rat pituitary adenoma cells. Journal of Neuro Oncology. 2010;100(1):51–64.

Identification of differentially expressed genes in Cushing's disease by integrated bioinformatics analysis

Status:

Version 1

Abstract

Figures