frontiers in Endocrinology
Check for updates
Identification of Seven Aberrantly Methylated and Expressed Genes in Adrenocortical Carcinoma
He Xiao 11, Weixiang He11, Ping Chen11, Deqiang Xu1, Guang Zeng1, Zhuo Li1,2,3,4,5, Mingliu Huang1, Xinghuan Wang1, Michael E. DiSanto6 and Xinhua Zhang1*
1 Department of Urology, Zhongnan Hospital of Wuhan University, Wuhan, China, 2 Department of Urology, Shenzhen Sixth People’s Hospital, Shenzhen, China, 3 The Sixth Affiliated Hospital of Shenzhen University Health Science Center, Shenzhen, China, 4 Affiliated Shenzhen Sixth Hospital of Guangdong Medical University, Shenzhen, China, 5 Shenzhen Key Laboratory for Endogenous Infection, Shenzhen, China, 6 Departments of Biomedical Sciences & Surgery of Cooper Medical School of Rowan University, Camden, NJ, United States
Background: Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with an unfavorable prognosis and limited treatment options. Nevertheless, no clinically applicable molecular markers have been identified for the progression of ACCs. DNA methylation alterations were found to contribute to the development of ACC in recent decades.
Material and Methods: The aims of the current study was to identify the abnormally methylated differentially expressed genes (DEGs) in ACCs, and to elucidate the mechanistic basis for these changes. Analyses were conducted on gene expression and gene methylation profile datasets to identify the aberrantly methylated DEGs. The DAVID software was used to conduct the analyses of functional enrichment on screened genes. Finally, expression was validated, and the relationship between abnormally methylated DEGs and clinical features was determined via the Oncomine database and The Cancer Genome Atlas (TCGA). To further verify the altered expression and methylation status of our identified genes we also validated these changes at the tissue and cellular levels.
Results: We screened and identified 92 differentially expressed genes and 802 abnormally methylated genes. Furthermore, seven aberrantly methylated and dysregulated genes were identified and validated, along with a number of functional enriched pathways. Among these seven genes, the expression or methylation status is significantly correlated with different pathological stages and overall rates of survival. In validation, the expression of seven genes were significantly altered and five genes were hypermethylated in ACC.
Conclusions: Our study identified abnormally methylated DEGs and potentially affected pathways in ACCs, from which we could begin to understand the basic molecular mechanisms of these alterations. Moreover, these abnormally methylated genes might serve as therapeutic targets and biomarkers to allow ACC patients to be more precisely diagnosed and effectively treated.
Keywords: adrenocortical carcinoma, bio-informatics, methylation, differentially expressed genes, biomarkers
OPEN ACCESS
Edited by:
Louise Metherell, Queen Mary University of London, United Kingdom
Reviewed by:
Sonir Roberto Rauber Antonini, University of São Paulo, Brazil Amit Tirosh, Sheba Medical Center, Israel
*Correspondence: Xinhua Zhang zhangxinhuad@163.com
+ These authors have contributed equally to this work
Specialty section:
This article was submitted to Systems Endocrinology, a section of the journal Frontiers in Endocrinology Received: 24 October 2018 Accepted: 28 June 2019 Published: 12 July 2019
Citation:
Xiao H, He W, Chen P, Xu D, Zeng G, Li Z, Huang M, Wang X, DiSanto ME and Zhang X (2019) Identification of Seven Aberrantly Methylated and Expressed Genes in Adrenocortical Carcinoma. Front. Endocrinol. 10:472. doi: 10.3389/fendo.2019.00472
INTRODUCTION
Adrenocortical tumors have a relatively high incidence of about 1-10%, most of which are benign tumors discovered by accident (1). In contrast, Adrenocortical carcinoma (ACC), is the most aggressive among all types of adrenal tumors and has a quite low incidence with only 0.7-2 diagnosed patients in one million people every year (2, 3). Nevertheless, this cancer is very deadly, yielding an average 5 year survival rate under 35% (4-7). ACCs are generally discovered in grown-ups aged 40-50 years, while it may also happen to children, especially those with a typical germline mutation in tumor protein p53 (8). Additionally, the metastasis of ACCs usually happens simultaneously when the cancer is diagnosed, resulting in poor prognosis (9). For the majority of the patients, mitotane treatment and systemic adjuvant chemotherapy are necessarily required to reduce the excessive secretion of hormones and the progression of the tumor (10-12). However, the effects of these agents are generally poor. At present, surgery has been taken as the only truly curative treatment for this disease. If the ACC tumor cannot be resected completely, the remaining therapeutic choices show little effect on the improvement of patient survival (13). In recent years, many studies identified genetic and epigenetic alterations in ACCs and suggested them to be key factors for the occurrence and development of ACCs. Hence, it is of great significance to unravel the alterations of genes in ACCs that can lead to the discovery of useful therapeutic targets along with an increased ability to predict individual responses.
Ever since the beginning of the post-genomic era, studies on tumorigenesis mechanisms have no longer been limited to gene deletion and mutation. In recent years, epigenetics has begun to play important roles in the study of tumorigenesis. DNA methylation is an epigenetic modification that is not accompanied by any alteration of DNA sequence, and usually occurs in the context of CpG dinucleotides. There are two types of DNA methylation patterns. One is a CpG-island hyper- methylation occurring in particular gene promoter areas, which is always related to tumor suppressor gene transcriptional silencing. Another is a global hypo-methylation related to increased instability of the chromosomes, reactivation of transposable elements, and parental imprinting loss (14). The abnormal DNA methylation can impact tumorigenesis and tumor progression by influencing the regulatory patterns of oncogenes and anti-oncogenes, as well as the structure of chromatin at the level of transcription (15, 16). Meanwhile, an increasing number of studies have been suggesting a critical role of histone modification and DNA methylation in the occurrence of tumors. Additionally, there is evidence showing the intimate relationship between tumorigenesis and GpG-island hyper-methylation which causes the silencing of downstream genes (17-19). With regard to adrenocortical tumors, evidence has shown the involvement of H19 promoter DNA methylation alteration in ACCs (20). Therefore, the abnormal methylation genes might be rediscovered as potential therapeutic targets or bio-markers for ACC.
In recent years, high-throughput sequencing has become a powerful tool which could search for epigenetic alternations in
genes and provide more useful information for diagnosis and prognosis in oncology. Differentially expressed genes (DEGs) and differentially methylated genes (DMGs) between tumor samples and normal controls can both be obtained from microarray data via a high-throughput platform. This data could lead to the identification of important epigenetic targets in ACCs.
There are very few studies that have simultaneously examined gene-methylation-profiles and gene-expression-profiles and their link to ACC development. In the current study, bio-informatics was employed to conduct systematic analysis and integration of the data collected from a gene-methylation-profile microarray (GSE77871) and multiple gene-expression-profile microarrays (GSE75415, GSE90713, GSE12368, GSE14922, and GSE19750. The analysis process involved screening gene enrichment in pathways and gene ontology (GO), overlapping five datasets with the use of a Venn diagram, obtaining DMGs and DEGs with the use of R software and validating them in a TCGA database, Genotype Tissue Expression (GTEx) database, Oncomine database and a GSE49280 database. This current study thus, seeks to identify abnormally methylated DEGs, along with alterations in signaling pathways and GO among normal and tumor samples and using this information elucidate the molecular mechanisms of ACCs.
MATERIALS AND METHODS
Microarray Data Information
Raw gene expression profiles (GSE75415, GSE90713, GSE12368, GSE14922, and GSE19750) and a raw gene methylation profile (GSE77871) were obtained from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Datasets GSE12368 (21) and GSE19750 (22) were performed on the same platform Affymetrix Human Genome U133 Plus 2.0 Array (HG U133 Plus 2.0). Dataset GSE14922 (23), GSE75415 (24), and GSE90713 (25) were performed on the platform Agilent-014850 Whole Human Genome Microarray 4x44K, Affymetrix Human Genome U133A Array and Affymetrix Human Gene Expression Array, respectively. These five datasets were combined and analyzed to screen DEGs. Dataset GSE77871 (26) based on the platform of Illumina HumanMethylation450 was used to screen DMGs. In addition, RNA sequencing data of 79 ACCs samples was downloaded from The Cancer Genome Atlas (TCGA) database (https://genome-cancer.ucsc.edu/) for further validation. Gene expression data were based on Illumina Hiseq’s RNA sequencing technology. Six normal samples and 18 ACC samples were included in the data of the gene-methylation- profile microarray, while 26 normal samples and 132 ACC samples were included in the data of the gene-expression- profile microarray. In order to verify the association between gene methylation and expression, we also downloaded the GSE49280 data.
Data Processing for the Identification of DEGs and DMGs
Raw expression data was calculated using preprocessing procedures: RMA background correction, log2 transformation,
quantile normalization, and median polish algorithm summarization using the “affy” (27) package of R software (version 3.5.0; Bell Laboratories, formerly AT&T, now Lucent Technologies, Murray Hill, NJ, USA). Additionally, “sva” (28) R package was used to remove batch effects between five datasets. Probes were annotated by the Affymetrix annotation files. R package “limma” (29) was applied to select the DEGs between 132 ACC samples and 26 normal samples. The cut-off standards for DMGs and DEGs were B > 0.2, p < 0.05 and |logFC| > 1, p < 0.05, respectively. Subsequently, the wANNOVAR tool (http:// wannovar.wglab.org/) was used to uncover the identities of DMGs. Then, the Venn diagram tool (http://bioinfogp.cnb.csic. es/tools/venny/) was applied to overlap the DEG identification processes for five gene-expression profile datasets. Finally, the upregulated hypo-methylated genes were acquired from the overlapped upregulated and hypo-methylated genes. The downregulated hyper-methylated genes were similarly acquired.
Gene Ontology and Pathway Enrichment Analysis of DMGs
For the analysis of gene ontology enrichment, the gene list, which contained 471 hyper-methylated and 336 hypo-methylated genes, was submitted to the DAVID website https://david-d. ncifcrf.gov/. The GO outcome p-value-ranking was performed with p < 0.01 as statistically significant. Those cut-off values of statistical significance were selected for molecular function (MF), cell component (CC), and biological processes (BP). For the analysis of KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment, the cut-off standard was set as p < 0.05.
Validation of the Identified Genes Expression and Methylation Levels
The multidimensional and comprehensive key genomic change maps for a variety of cancers were contained in the Oncomine, GTEx and TCGA databases. The identified genes were further validated with these three databases. Gene expression was compared between ACC and normal samples. In particular, TCGA data were analyzed with cBioPortal. As a publicly accessible platform, which provides the functions of downloading, analyzing, and visualizing of massive datasets of cancer genomes for different types of cancers, the tool of cBioPortal helped immensely in examining the correlation of DNA methylation with mRNA expression in the studies on ACCs, as well as for the altered genetic information associated with DMGs.
Correlation of the Identified Genes Expression and Methylation With Pathological Stages
We downloaded the TCGA data from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/) which is a robust data-driven platform which provides the functionality of data searching and downloading for bioinformatics and cancer researchers. The boxplot was applied to show the identified
genes’ relationship with the different pathological stages. One- way ANOVA was used to determine whether the findings are statistically significant.
ACC Cell Lines and Human ACC Samples
Human ACC cell lines (SW13, NCI-H296R) were purchased from the Procell Life Science & Technology Co., Ltd. in Wuhan, China. The cell line was identified by the China Centre for Type Culture Collection in Wuhan, China. SW13 was incubated in DMEM medium (Gibco, China), supplemented with 10% fetal bovine serum (FBS) (Gibco, Australia) according to the ATCC instructions. NCI-H296R was grown and maintained in DMEM media supplemented with 1% insulin transferrin selenium (Gibco, China) and 2.5% Nu-Serum I (Gibco, China) in a standard humidified incubator at 37°℃ in a 5% CO2 atmosphere. ACC and adjacent normal adrenal tissues (n = 8) were obtained from patients undergoing laparoscopic unilateral adrenalectomy at Zhongnan Hospital of Wuhan University. Two pathologists independently confirmed the histological diagnosis. All specimens were immediately fixed in 4% PFA (paraformaldehyde) and stored in liquid nitrogen. Prior to their operation, none of the patients enrolled in the study receive any treatment. The use of these ACC specimens was approved by the Ethics Committee at Zhongnan Hospital of Wuhan University, and informed consent was obtained from all patients.
Promoter Methylation Analysis and 5-Aza-2’-deoxycytidine Treatment
Bisulfite sequencing PCR (BSP) of the six methylated gene was accomplished by Sangon Biotech Co., Ltd. (Shanghai, China) and primers for the analysis were designed using the free on-line software MethPrimer (30). ACC cells were treated with 5-aza-2’- deoxycytidine (50 mM, Sigma-Aldrich, St. Louis, MO) for 48 h with a change of culture medium every day.
Total RNA Extraction and Real-Time RT-PCR
Total RNA was isolated from the frozen tissues using Takara RNAiso Plus (Takara Bio. Inc., Otsu, Shiga, Japan) according to the manufacturer’s protocol. Genomic DNA (gDNA) was removed and cDNA was reverse-transcribed using Takara PrimeScript™M RT reagent Kit with gDNA Eraser (Takara Bio. Inc., Otsu, Shiga, Japan) in a T100TM Thermal Cycler System (BioRad, USA). The experimental protocol utilized was first gDNA removal (42℃, 2 min), followed by reverse transcription (37℃ 15 min, 85℃ 5 s). Subsequently, all samples were amplified by a 25 ul reaction volume in a CFX96TM Real-time PCR Detection System (BioRad, USA), using SYBR® Premix Ex TaqTM II (Takara Bio. Inc., Otsu, Shiga, Japan). All samples were run in triplicate. The seven identified genes were investigated. The amplification program was repeated for 40 cycles. Primer sequences are shown in Table 1. For relative quantification, gene expression was normalized to expression of GAPDH housekeeping gene and compared by 2-AACT method
Immunofluorescence
Tissues were sectioned in 10 um thick slices and thawed, mounted onto glass slides using a cryostat (Leica CM 1850, Wetzlar, Germany), air-dried, and fixed for 10 min in ice cold acetone. Slides were washed in PBS and incubated for 2 h in a mixture of PBS supplemented with 0.2% Triton X-100 and 0.1% bovine serum albumin, followed by incubation overnight with the primary antibody (1:100). The secondary antibody employed to visualize the localization of HOXA5 is Cy3-conjugated goat anti-rabbit IgG (1:1,000). DAPI was used for staining the nucleus. Visualization was done with a laser scanning microscope (Olympus, Tokyo, Japan).
Survival Analysis
The survival rate analysis was conducted based on the TCGA database for the assessment of the identified genes’ effects on the prognosis of ACC patients. Firstly, patients with mRNA data were classified in two different categories in accordance with each gene’s median expressions (low vs. high). Patients with methylation data were similarly analyzed. Secondly, analysis was
| Target Gene | Primer sequence | |
|---|---|---|
| NR2F1 | Forward | 5'-CATTTTTGGGCGATCTCCAG-3' |
| Reverse | 5'-CTCGCTGCCTTCTTCTTTCG-3' | |
| KCNQ1 | Forward | 5'-GTCCATGCAACAAGGTGGTC-3' |
| Reverse | 5'-GTTGTTCGTAGGTGGGCAGG-3' | |
| HOXA5 | Forward | 5'-GCGTGGAAGTGTTCCTGTCT-3' |
| Reverse | 5'-ACCGCTTGGAGTCACAGTTT-3' | |
| CD14 | Forward | 5'-AAGCACTTCCAGAGCCTGTC-3' |
| Reverse | 5'-TCGTCCAGCTCACAAGGTTC-3' | |
| PTGER4 | Forward | 5'-GATCAAATGCCTCTTCTGCCG-3' |
| Reverse | 5'-ACTGCTGATCTCCTTCAGCTC-3' | |
| CYP11B1 | Forward | 5'-TGATGGCTTTGGTGCATGTG-3' |
| Reverse | 5'-CCCAACCGTGATCTGTCTGT-3' | |
| ESM1 | Forward | 5'-TTGCTACCGCACAGTCTCAG-3' |
| Reverse | 5'-GCAGGTCTCTCTGCAATCCA-3' | |
| GAPDH | Forward | 5'-ATCCCATCACCATCTTCCAGGAG-3' |
| Reverse | 5'-CCTGCTTCACCACCTTCTTGATG-3' |
conducted on the patients with both mRNA and methylation data. Finally, we performed Kaplan-Meier survival analysis and the log-rank test by adopting the “survival” R package. One-way analysis of variance (ANOVA) and paired 2-tailed Student’s t- tests were used to analyze the statistical significance of differences of data.
RESULTS
Identification of the DEGs and DMGs
The R software was used to preprocess data and assess quality. Then the matrices of expression were obtained from datasets of GSE75415, GSE90713, GSE12368, GSE14922, and GSE19750 (containing 132 ACC samples and 26 normal samples, Table 2). The DEGs of five datasets are shown as volcano plots in Figures 1A-E. We overlapped the DEGs of five datasets. In total we identified 92 DEGs with 57 downregulated and 35 upregulated genes in common between the 5 sets (Figure 1F). The DMGs are shown as a heatmap in Figure 2. In total we identified 802 DMGs including 336 hypo-methylated and 466 hyper-methylated genes in GSE77871.
Analysis of Functional Pathway Enrichment for 802 DMGs
Gene ontology (GO) and KEGG analyses of hypo-methylated and hyper-methylated genes were displayed in Figure 2. GO analysis indicated 42 terms (biological process, cellular component and molecular function) enriched from 466 hyper-methylated genes (p <0.01). Among those, 18 terms were found most significant (p < 0.001). KEGG analysis found seven pathways were significant enriched (p < 0.05), including melanogenesis, maturity onset diabetes of the young, signaling pathways regulating pluripotency of stem cells, cholinergic synapse, transcriptional misregulation in cancer, retrograde endocannabinoid signaling, and wnt signaling pathway. On the other hand, GO analysis indicated 27 terms enriched (p < 0.01) from 336 hypo-methylated genes and 13 terms were showed most significant (p < 0.001). KEGG analysis found six pathways were significantly (p < 0.05) enriched, including graft-vs .- host disease, autoimmune thyroid disease, staphylococcus aureus infection, asthma, nicotine addiction, and olfactory transduction.
| Dataset | No. of samples (N/T) | Array types | Experiment type | Origin |
|---|---|---|---|---|
| Soon et al. GSE12368 | 6/12 | Affymetrix Human Genome U133 Plus 2.0 Array | mRNA | Endocrine-related cancer 2009/6/1 |
| Tombol et al. GSE14922 | 4/4 | Agilent-014850 Whole Human Genome Microarray 4x44K G4112F | mRNA | Endocrine-related cancer 2009/9/16 |
| Demeure et al. GSE19750 | 4/44 | Affymetrix Human Genome U133 Plus 2.0 Array | mRNA | Surgery 2013/6/28 |
| West et al. GSE75415 | 7/15 | Affymetrix Human Genome U133A Array | mRNA | Cancer Research 2007/1/15 |
| Weiss et al. GSE90713 | 5/57 | Affymetrix Human Gene Expression Array | mRNA | Oncotarget 2017/9/26 |
| Legendre et al. GSE77871 | 6/18 | Illumina HumanMethylation450 | DNA | PLoS ONE 2016/3/10 |
A
GSE12368
B
GSE14922
C
GSE19750
00
6
+
+
៛
~
~
2
logFC
logFC
0
0
logFC
0
N
“
N
7
7
7
”
9
0
2
4
6
8
0
1
2
3
4
5
6
0
5
10
15
20
25
30
35
-log10(P.Value)
-log10(P.Value)
-log10(P.Value)
D
E
F
G5575415
GSE75415
GSE90713
+
530
+
140
9
2
23
2
71
2
GSEDOT13
87
22
31
GSE 18750
29
102
logFC
1451
0
logFC
0
8
27
44
92
33
32
“
“
11
4
31
68
11
77
194
77
Y
7
35
81
218
1111
614
0
2
4
6
8
10
12
0
5
10
15
20
GSE12308
-log10(P.Value)
-log10(P.Value)
GSE 14022
Identification of Abnormally Methylated Differentially Expressed Genes
To explore the different expressions of those abnormally DMGs, the overlapping of hyper-methylated genes with downregulated genes and hypo-methylated genes with upregulated genes was performed. As shown in Figure 3, this analysis yielded six downregulated hyper-methylated genes and one upregulated hypo-methylated gene.
Assessment of Seven Genes in TCGA and Oncomine Database
The seven genes identified above were further validated in the Oncomine database and TCGA database. As shown in Figures 4, 5, significant differences were identified and determined to be identical to those differences found between tumor and normal samples in both databases regarding the expression of the six downregulated hyper-methylated genes and the one upregulated hypo-methylated gene.
Genetic Information and Methylation Statuses of the Seven Genes in TCGA
Besides methylation, other genic alterations were explored with the software of cBioPortal. As shown in Figure 6A,
it was observed that 40 out of the total 79 cases having genic alterations, such as missense, amplification, etc., occurred most often in PTGER4, CYP11B1, and NR2F1 (13, 15, and 21%, respectively). The correlation of seven genes expression with their DNA methylation was shown in Figure 6B. It was noted that the negative correlation of six downregulated hyper-methylated genes were most significant, indicating that methylation might be highly important in gene expression.
Correlation of Seven Genes Expression and Methylation With Pathological Stages
The Cancer Genome Atlas (TCGA) dataset was used to analyze the association of the identified genes with pathological stage. It was found that the less expressed 4 genes (CD14, CYP11B1, KCNQ1, and PTGER4) were significantly negatively correlated with pathological stage (Figure 7A). However, the correlation of the rest three genes (ESM1, HOXA5, and NR2F1) were not significantly different. Comparisons were also made among all genes’ methylation status with various pathological stages and it was noted only the hyper-methylated two genes (PTGER4 and NR2F1) were significantly positively correlated with pathological stage (Figure 7B).
Z-Score 4
0
=4
Normal
ACC
Hypermethylated
Hypomethylated
plasma membrane
transcription, DNA-templated-
lipopeptide binding
DNA-templated
antigen processing and presentation, exogenous lipid antigen via MHC class Ib
positive regulation of transcription from RNA polymerase Il promoter
calcium ion binding
endogenous lipid antigen binding
Integral component of plasma membrane -
exogenous lipid antigen binding
sequence-specific DNA binding
Graft=versus-host disease
setgence-specific DNA binding
transcription from RNA polymerase Il promoter -
sodium-independent organic anion transport
homophilic cell adhesion via plasma membrane adhesion molecules
negative regulation of transcription from RNA polymerase Il promoter
MHC class Il protein complex
antigen processing and presentation of peptide or polysaccharide antigen via MIHC class II
RNA polymerase Il core promoter proximal region sequence-specific DNA binding -
RNA polymerase Il core promoter proximal region sequence-specific binding
MHC class Il receptor activity
cell adhesion
startle response
nervous system development
vocalization behavior
RNA polymerase Il regulatory region sequence-specific DNA binding -
Autoimmune thyroid disease
getin cytoskeleton
Asthma
transmembraneecon
regulation of ion transmembrane transport skeletal system devpinnesont
synapse organization
Staphylococcus aureus infection
jon transferhennes
potassium ion transmembraneque
Nicotine addiction
voltage-gated potassium channel complex
regulation of ion transmembrane transport
potassium ion transport -
Transcriptional misregulere
voltage-gated potassium channel complex
pluriparency basan das
Signaling pathways regulating pluripotency of stem cells -
dendrite
chemical synaptic transmission
anterior/posterior pattern specification -
inner ear morphogenesis-
detection of chemical stimulus involved in sensory perception
Wnt signaling pathway-
Chainemayans-
Count
Count
cell junction
Melanogenesis-
· 50
transmembrane signaling receptor activity
synapse assembly-
·
50
.
100
odorant binding
Retrograde endocannabinoid signaling -
phagocytosis-
100
150
sensory perception of smell
delayed rectifier potassium channel activity-
integral component of plasma membrane
embryonic forelimb morphogenesis-
regulation of heart contraction
p.value
p.value
Olfactory transduction
runnic skeletal system development
olfactory receptor activity
calcium-dependent cel-cell adhesion via plasma membrane cell adhesion molecules-
embryonic hindlimb morphogenesis-
0.03
0.04
detection of chemical stimulus involved in sensory perception of smell
Maturity onset diabetes of the young-
0.02
0.03
- G-protein coupled receptor activity
cell fate determination
enteric nervous system development-
0.02
G-protein coupled receptor signaling pathway
regulation of transcription in ventral spinal cord interneuron specification -
0.01
0.01
plasma membrane
subpallium development-
integral component of membrane
0.0
0.1
0.2
0.3
0.0
0.2
0.4
0.8
GeneRatio
GeneRatio
Survival Analysis
On the basis of the clinical information, and data of RNA sequencing and methylation in TCGA, patients were classified in two different categories in accordance to each gene’s median expressions. Then the Kaplan-Meier survival curve was plotted. It revealed that the overall rates of survival for the patients with four hyper-methylated genes (CD14, HOXA5, KCNQ1, and NR2F1), one more expressed gene (ESM1) or three less expressed genes (CD14, HOXA5, and KCNQ1) were greatly lower, which were shown in Figures 8A,B. When their expression and methylation were combined, the overall rates of survival significantly negatively correlated with hypermethylation and low expression of five genes (CD14, HOXA5, KCNQ1, and NR2F1) along with hypomethylation and high expression of ESM1 (Figure 8C).
Validation of the Methylation Status of Six Genes in GSE49280
The correlation of six genes expression with their DNA methylation was shown in Figure 9. It was noted that the negative correlation of six downregulated hyper-methylated genes were most significant, indicating that methylation might be highly important in gene expression.
Validation of the Methylation Status of Six Genes in ACC Tissues and Cells
In order to confirm the methylation of the six hyper-methylated gene promoters in ACC, bisulfite sequence analyses were performed on three ACC samples and three adjacent normal ones, demonstrating that methylation of all CpG sites of five
A
6
KCNQ1
B
1
ESM1
GSE77871(hypermethylated)
PTGER4
HOXA5
GSE77871(hypomethylated)
CD14
NR2F1
CYP11B1
Down-regulated
Up-regulated
460
51
335
34
FIGURE 3 | Identification of aberrantly methylated differentially expressed genes. (A). Six hypermethylated and downregulated genes (KCNQ1, PTGER4, HOXA5, CD14, NR2F1, and CYP11B1) were identified from hypermethylated genes in GSE77871 and common downregulated genes in five mRNA expression profiles. (B). One hypomethylated and upregulated gene (ESM1) was identified from hypomethylated genes in GSE77871 and common upregulated genes in five mRNA expression profiles.
A
B
C
D
1.5
1.5
2.5
2.0
2.0
1.5
log2 median-centered intensity
1.0
log2 median-centered intensity
1.0
log2 median-centered intensity
log2 median-centered intensity
1.5
1.0
0.5
KCNQ1
0.5
PTGER4
P-value: 5.85E-9
1.0
HOXA5
CD14
P-value: 6.20E-17 t.Test: 14.348
P-value: 2.69E-11
P-value: 2.14E-11
t.Test: 7.905
t-Test: +10.194
0.5
t-Test: 9.224
Fold Change: - 1.772
Fold Change: - 1.616
0.5
Fold Change: - 1.867
Fold Change: - 1.841
0.0
0.0
0.0
0.0
.
.
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
E
F
G
1.5
2.5
1.5
2.0
log2 median-centered intensity
1.0
log2 median-centered Incensky
1.5
log2 median centered intensity
1.0
ESM1
P-value: 1.90€-8
1.0
t-Test: 6.912
0.5
NR2F1
CYP11B1
P-value: 1.06E-12
P-value: 7.60E-6
0.5
Fold Change: 1,425
t-Test: - 10.343
0.5
t-Test: +5.145
Fold Change: +1.653
0.0
Fold Change: - 1.929
0.0
0.0
.
0.5
.
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
-1.0
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
-0.5
Adrenal Cortex (10)
Adrenal Cortex Carcinoma (29)
gene promoters were more frequent in ACC rather than in normal samples (Figures 10A-E). Among them, the degree of methylation of CYP11B1 is low and there is no significant difference compared with normal tissues (Figure 10F). Next, to further investigate the relationship between methylation and expression of six genes, SW13 and NCI-H295R cells were cultured with the demethylating agent 5-Aza-2’-deoxycytidine, and quantitative real time RT-PCR indicated the restoration of seven genes’ expression in ACC cells (Figures 11A-F, 12A-F).
Validation of the mRNA Expression of Seven Genes in ACC and Normal Tissues
Expression of CYP11B1, PTGER4, HOXA5, CD14, NR2F1, KCNQ1, and ESM1 mRNA were determined using quantitative real time RT-PCR between ACC samples and normal ones. HOXA5 expression was most significantly downregulated at the transcription level (p = 0.0019) in ACC. Real time RT-PCR also showed that for the other six genes that NR2F1 (p = 0.034),
A
B
C
D
12
KCNQ1
*
PTGER4
*
HOXA5
*
10
CD14
*
4
log2(TPM+1) transformed expression
10
log2(TPM+1) transformed expression
8
log2(TPM+1) transformed expression
log2(TPM+1) transformed expression
8
0
00
6
6
6
2
+
+
+
-
~
2
~
0
0
..
0
(num(T)=77; num(N)=128)
ACC
(num(T)=77; num(N)=128)
ACC
(num(T)=77; num(N)=128)
ACC
ACC (num(T)=77; num(N)=128)
E
F
15
G
1
NR2F1
CYP11B1
ESM1
log2(TPM+1) transformed expression
0
log2(TPM+1) transformed expression
log2(TPM+1) transformed expression
CO
5
0
+
+
0
5
~
2
-
0
0
0
(num(T)=77; num(N)=128)
ACC
(num(T)=77; num(N)=128)
ACC
(num(T)=77; num(N)=128)
ACC
KCNQ1 (p=0.04), CD14 (p=0.038), PTGER4 (p=0.03), CY11B1 (p = 0.0037), and ESM1 (p = 0.0013) the mRNA expression was significantly altered in the ACC samples (Figure 13A).
Immunohistochemistry and Immunofluorescence Analysis of HOXA5 in ACC and Normal Tissues
As shown in Figures 13B-E, HOXA5 is localized in adrenal cortical cytoplasm and nucleus. In normal adrenal gland, HOXA5 was mainly present in adrenocortical cells nuclei and partly present in the adrenocortical cytoplasm. In ACC, the localization is similarly with that in normal ones, but in ACC, the fluorescence intensity is significantly lower than that of normal. That is to say, immunohistochemistry and immunofluorescence also confirm that the expression of the
hypermethylated gene HOXA5 in ACC is lower than that in normal tissues.
DISCUSSION
The current study, for the first time, identified six downregulated hyper-methylated genes and one upregulated hypo-methylated gene related to ACC. Additionally, these genes are found to be correlated with different pathological stages and overall rate of survival of ACCs. Our study suggests these abnormally methylated and deregulated genes could be used as therapeutic targets or biomarkers for ACC patients to be precisely diagnosed and effectively treated.
The six downregulated hyper-methylated genes were CYP11B1, PTGER4, HOXA5, CD14, NR2F1, and KCNQ1 and
A
50%
NR2F1
21%
CYP11B1
15%
Alteration Frequency
40%
PTGER4
13%
30%
KCNQ1
6%
20%
CD14
5%
ESM1
5%
10%
HOXA5
4%
Genetic Alteration
Missense Mutation (unknown significance)
Amplification
Deep Deletion
mRNA Upregulation
No alterations
B
Pearson-Correlation :- 0.4208 P-value:1.125e-04 Sample Size:(N=79)
Pearson-Correlation :- 0.5884 P-value:1.182e-08 Sample Size:(N=79)
Pearson-Correlation :- 0.4856 P-value:5.727e-06 Sample Size:(N=79)
Pearson-Correlation :- 0.6582 P-value:4.354e-11 Sample Size:(N=79)
NR2F1 mRNA expression (log2)
SU
CYP11B1 mRNA expression (log2)
&
PTGER4 mRNA expression (log2)
KCNQ1 mRNA expression (log2)
10
₫
0
0
0
0
00
₪
00
0
2
10
0
*
៛
5
0
0
2
*
-0.4
-0.2
0.0
0.2
0.4
-0.4
-0.2
0.0
0.2
0.4
-0.3
-0.2
-0.1
0.0
0.1
0.2
0.3
-0.2
0.0
0.2
0.4
NR2F1 methylation (beta value)
CYP11B1 methylation (beta value)
PTGER4 methylation (beta value)
KCNQ1 methylation (beta value)
Pearson-Correlation :- 0.6935 P-value:1.406e-12 Sample Size:(N=79)
Pearson-Correlation :- 0.3024 P-value:6.764e-03 Sample Size:(N=79)
Pearson-Correlation :- 0.4101 P-value:1.742e-04 Sample Size:(N=79)
CD14 mRNA expression (log2)
0
ESM1 mRNA expression (log2)
$
HOXA5 mRNA expression (log2)
₪
¢
0
8
ç
0
0
0
2
10
1
0
0
9
0
0
8
00
0
+
0
0
*
₪
-0.3
-0.2
-0.1
0.0
0.1
0.2
0.3
-0.4
-0.3
-0.2
-0.1
0.0
0.1
-0.4
-0.2
0.0
0.2
0.4
CD14 methylation (beta value)
ESM1 methylation (beta value)
HOXA5 methylation (beta value)
the one upregulated hypo-methylated gene was ESM1. These abnormal expressions of these genes were further confirmed in other databases of TCGA and Oncomine. It is well-known that DNA methylation can turn off the activity of certain genes, while demethylation induces gene reactivation and expression. Therefore, hyper-methylation could lead to gene downregulation, functionally similar to a tumor suppressor gene, while hypo-methylation could result in gene upregulation, functionally similar to oncogenes. Thus, the hyper-methylated upregulated genes and hypo-methylated downregulated genes were not further analyzed in the present study.
It is widely acknowledged that the imbalance between oncogene and tumor suppressor genes plays a key role in tumorigenesis. Dysregulation of genes identified in our study may be involved in the development of ACC. Indeed, these seven genes have unique functions. As an enzyme coming from the big family of the cytochrome P450 genes, CYP11B1 can
serve as a catalyzer for lots of drug metabolism reactions, as well as the synthesis of lipids like steroids and cholesterol. Studies have found the association of downregulated members from cytochrome P450 with the incidence of various cancers, such as renal, breast, and colorectal cancers (31). PTGER4 is from the family of G protein coupled receptors, which resides on 5p13.1, encoding one out of the four identified receptors for prostaglandin E2 (32, 33). HOXA5 (homeobox A5) is a member of highly-conserved HOX gene family that is responsible for controlling the differentiation of cells and the development of the embryo (34). CD14 (CD14 molecule) is one of the critical elements that compose the immune system of humans and it has certain functions in protecting human bodies from the microbial products in the environment (35). NR2F1 is a nuclear hormone receptor and transcriptional regulator (36-38). KCNQ1 is a voltage-gated potassium channel required for the repolarization phase of
A
Relative mRNA expression (log2)
15
Stage
Stage I
Stage II
10
Stage III
Stage IV
5
p=0.02713
p=0.01443
p=0.9730
p=0.5650
p=0.04504
p=0.5905
p=0.01345
CD14
CYP11B1
ESM1
HOXA5
KCNQ1
NR2F1
PTGER4
B
1.00
Beta value of DNA methylation
0.75
Stage
Stage I
0.50
Stage II
Stage III
Stage IV
0.25
0.00
p=0.3633
p=0.5559
p=0.2783
p=0.2209
p=0.1917
p=0.01297
p=0.001726
CD14
CYP11B1
ESM1
HOXA5
KCNQ1
NR2F1
PTGER4
FIGURE 7 | Correlation of the expression levels and methylation levels of identified genes with the pathological stages of ACC. The boxplots show the medians and dispersions of the samples of different pathological stages for each of the seven genes. One-way ANOVA was used for different pathological stages and a p < 0.05 was regarded as significant. (A) Boxplots of the identified genes transcription level in different pathological stages. (B) Boxplots of the identified genes methylation level in different pathological stages.
the cardiac action potential (39). Additionally, the upregulated hypo-methylated gene ESM1, encodes one secreted protein, whose expression mainly occurs within the endothelial cells of human kidneys and lungs (40). Although the specific functions of these seven genes in ACC have never been demonstrated, their mutation and expression dysregulation may have impacts on tumorigenesis, staging, and prognosis of ACC.
The current study also demonstrated a negative correlation between gene expression and pathologic stages. However, only the four less expressed genes (CYP11B1, PTGER4, CD14, KCNQ1) identified in the study were significantly associated with a higher pathology stage. Consistently, our study showed hypermethylation of two genes (PTGER4 and NR2F1) led to a
higher pathology stage. Finally, five genes (CYP11B1, PTGER4, CD14, KCNQ1, and NR2F1) were significantly correlated with pathology stages. These genes functionally act as tumor suppressor genes when the genes are hypermethylated. In survival analysis, the overall rates of survival for the patients with 4 hyper-methylated genes (CD14, HOXA5, KCNQ1, and NR2F1), one more expressed gene (ESM1) or three less expressed genes (CD14, HOXA5, and KCNQ1) were greatly lower. When their expression and methylation were combined, the overall rates of survival significantly negatively correlated with hypermethylation and low expression of five genes (CD14, HOXA5, KCNQ1, and NR2F1) along with hypomethylation and high expression of ESM1. In addition, altered DNA methylation patterns have been demonstrated to be the first detectable
A
Survival curve (p=0.008)
Survival curve (p=0.091)
Survival curve (p=0.048)
Survival curve (p=0.005)
Survival curve (p=0.042)
Survival curve (p=0.21)
Survival curve [p=0.259)
2
CD 14 high expression
9
CYP-1 18-1 hich expression
1
TA
ESM1 high expression
9
HOKAS high expression - HOMAG love expression
9
-
9
CDtd low expression
… CYP1181 low espression
… ESM1 low expression
KONG1 high expression
9
KONQt low eagression
NRJF1 high expression
NR2F1 low expression
PTGER4 high expression
PTOERd kovi expression
8
a
3
2
8
3
3
3
Survival tata
3
Survival tata
3
Servinal tatn
0
Survindl sata
8
3
Survival tuta
3
3
3
2
3
2
3
=
=
=
=
~
=
=
9
2
9
9
8
2
0
0
2
4.
₿
8
10
12
0
2
4
.
8
10
12
e
2
4
·
M
50
12
@
2
4
.
3
10
12
0
2
4
6
8
15
12
₿
2
#
·
4
»
12
a
2
4
.
4
10
12
Tine (roör!
Time lycer)
Tina (year)
Tino (ymär)
Time (poor]
Tira lyear)
Tinso (ymar)
B
Survival curve (p=0.004)
Survival curve (p=0.064)
Survival curve (p=0.125)
Survival curve (p=0.018)
Survival curve (p=0.001)
Survival curve (p=0.013)
Survival curve [p=0.981)
2
3
CD 54 Hypermethylation
9
TE
CYP-1 181 Hypermethylation
1
Ya
ESMi Hyperrathıyladior
9
… CYP1181 hyporastylation
… ESMi hypomethylation
Hyparmethylation - HOKAS typomethylation
B
-
KONQ1 Hypernethyarsice KCNQt hypomethylation
9
—
NRJF1 Hypermachylation
9
NR2F1 hypomethylation
PTGERA Hyperuttrylation
.-.- PTOERi typoinetrelacion
3
3
3
8
3
3
g
Sarvinal soft
3
Survival tata
2
3
Barvinal toft
3
Savingi satın
3
Survival tata
3
3
3
3
3
3
3
3
=
=
=
=
7
=
=
3
a
2
0
2
0
9
0
2
6
8
10
12
0
2
4
·
8%
10
12
2
2
4
.
M
50
12
a
2
4
6
3
10
12
0
2
4
.
3
10
12
0
4
.
4
10
12
#
2
2
·
3
10
12
Time (rødr)
Time lycer)
Tina (year)
Tine (ynär)
Time (poor]
Tira byeağ
Tinso (ymar)
C
Survival curve (p=0.002)
Survival curve (p=).067)
Survival curve (p=0.038)
Survival curve (p=0.006)
Survival curve (p=0.004)
Survival curve (p=0.025)
Survival curve [p=0.441)
2
“A
CD14 hyper & low expression
=
Kg
CYP-4101 hyper & low expression
9
O
” ESMi hyper & low expression
9
1 HORAS hyper & low napression
:
+
. NONG1 hyper & low expression
9
NR2F5 typer I love expression NR2Fs trypo & high expression
9
?
CD14 hypo & high expression
CYP-1101 hypo & high expression
-
. ESMff typo & high expression
*** HORAS hypa & high expression
- HONG1 hypo & high expression
-
PTGER4 hyper & low expression
- PTGERli hypa & high expression
3
H
a
3
3
3
3
3
8
Survival tata
3
Survival tiền
3
8
3
Survival tain
3
Survival tiên
3
=
2
O
¥
3
=
=
=
=
2
=
=
2
8
g
8
8
8
3
0
2
4
.
9
10
12
0
2
.
.
3
10
13
D
2
4
6
3
50
12
a
3
4
0
a
10
12
0
2
4
6
8
10
12
D
2
4
.
@
10
4
2
4
.
a
10
12
Tine (roar!
Tive (year)
Time (year)
Tine (ymary
Tanve (poor]
Time tyear
Tine (ymary
A
B
C
1
1
0.9
Spearman: - 0.44
Spearman: - 0.65
1
Spearman: - 0.48
(p = 4.4890-5)
0.9
(p = 1.03e-10)
(p = 5.911e-6)
0.8
Pearson: - 0.36
Pearson: - 0.50
(p = 1.1910-3)
(p = 3.347e-6)
Pearson: - 0.40
0.8
0.8
(p = 2.9330-4)
0.7
NR2F1: Methylation (HM450)
0.7
HOXA5: Methylation (HM450)
0.6
0.6
0.6
0.5
0,5
0.4
0.4
y = - 0.03x + 0.6
R2 = 0.13
0.4
0.3
=- 0.14x + 0.59
0.2
y = - 0.17x + 0.54
0.3
R= = 0.24
R= = 0.16
0.2
0.2
0-
0.1
0.1
0
0
-0.2
2
0
2
4
6
8
10
12
14
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
-1
.0.5
0
0.5
1
1.5
2
25
3
3.5
4
NR2F1: mRNA Expression z-Scores (RNA Seq V2 RSEM)
KCNQ1: mRNA Expression z-Scores (RNA Seq V2 RSEM)
HOXA5: mRNA Expression z-Scores (RNA Seq V2 RSEM)
D
E
F
0.9-
0.9
1
Spearman: - 0.71
Spearman: - 0.63
0.8
(p = 2.05e-13)
0.8
(p = 4.61e-10)
Spearman: - 0.69
(p = 1.350-12)
Pearson: - 0.67
Pearson: - 0.39
0.7
(p = 2.040-11)
(p = 3.5580-4)
0.8
Pearson: - 0.55
(p = 1.50e-7)
0.7
CD14: Methylation (HM450)
0.6
PTGERA: Methylation (HM450)
CYP11B1: Methylation (HM450)
0.6
0.6
.
.
0.5
0.5
= - 0.05x + 0.69
R= = 0.16
0.4
0.4
y =- 0.14x + 0.61
0.4
R = 0.44
0.2
-0.08x + 0.61
0.3
.
R° = 0.31
0.3
0.2
0
0.2
0.1
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
0.1
-1
0
1
2
3
4
5
6
7
8
0
2
4
6
8
CD14: mRNA Expression z-Scores (RNA Seq V2 RSEM)
PTGER4: mRNA Expression z-Scores (RNA Seq V2 RSEM)
CYP11B1: mRNA Expression z-Scores (RNA Seq V2 RSEM)
A
TSS
B
TSS
C
T88
GC Percentage
20 40 60 80
GC Percentage
20 40 000
DC Percentage
20 40 60
0
0
500 bp
1000 bp
1500 bp
2000 bp
1
2500 bp
3000 bp
3500 Bp
4000 bp
4500
CpGl
500 lbp
1000 bp
1500 bp
2000 bp
2500 bp
3009 bp
3200 bp
I
I
4000 bp i
4500
0
o
CpGI
11
500 bp M
Y
I
1000 bp
1500 bp
2000 bp
2500 bp
3000 hp
3:00 hp
4000 bp
4500
NR2F1, Chromosome 5 (from 93580337ep to 53584837bp)
KCNQ1, Chromosome 11 (from 2441991bp to 2446491bp)
HOXAS, Chromosome 7 [from 27142168bp to 27146868bp)
Normal 1
Normal 1
Normal 1
ACC 1
ACC 1
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 1
Normal 2
Normal 2
-OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 2
ACC 2
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 2
00000000OOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 2
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 3
ACC 3
-000000000000000000000000000000000OOO
ACC 3
CpG
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 96 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
CpG
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 26 27 28 29 30 31 32 33 34 35 36 Unmethylated Methylabad
CpG
Bisulfite sequence (from 793 - 1256bp)
Unmethylated
Methylated
Bisulfite sequence (from 2645 - 3582bp)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 Bisulfite sequence [from 2493 - 2902bg) Methylated
Unmethylated ☐
D
E
F
TSS
TSS
TSS
GC Percentage
20,40,60,00,
GC Percentage
20, 40 00, 00
GC Percentage
2
0
a
0
CpGI
I
500 bp
1000 bop
1500 bp
2000 bp
2500 bp
3000 bp
3300 bp
4000 bp
1 11 1000 ER EN ITI
4500
CpG
o
I
500 bp
1000 bp
1500 bp
2000 bp
3000 bp
3500 bp
4000 bp 4500
CpCl
500
bp
I
1000
1500 bp
2000 bp
2500 bp
3500 bp
4000 bp
4500
CDI4, Chromosome 5 (from 140632201bp to 140636701bp)
PTGER4, Chromosome 5 (from 48676930bp to 40681430bp)
CYP11B1, Chromosome 8 (from 142878320bp to 142882820bp)
Normal 1
-000OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal1 .OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 1
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 1
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 1
-OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO -OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 1
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 2
Normal 2
Normal 2
000000000000000000OOOOOOOOOOOOOOOOOO
ACC 2
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 2
-OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 2
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
Normal 3
-000000000000000000OOOOOOOOOOOOOOOOOO
ACC 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
ACC 3
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
CpG
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 Bisulfite sequence (from 3735 - 41550p) Unmethylated Methylated
CpQ
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 20 27 28 29 30 31 32 33 34 36 38 Bisulfite sequence (from 2407 - 2964bp)
CpG
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 Methylated Bisulfite sequence |from 3200 - 3500bp) Unmethylated
Unmethylated
Wetirylated
A
B
C
p = 0.0006
p = 0.0012
p = 0.0008
Relative expression of NR2F1
10-
p = 0.0227
Relative expression of KCNQ1
15-
p = 0.0410
Relative expression of HOXA5
15-
p = 0.0182
8
p = 0.0090
10
p = 0.0027
6
10-
p = 0.0009
4
5.
5
2
0
0
0
0
0.5
1
0
0.5
1
0
0.5
1
5 Aza (um)
5 Aza (um)
5 Aza (um)
D
E
p = 0.0003
F
p = 0.0005
p = 0.0904
Relative expression of CD14
15-
p = 0.0238
Relative expression of PTGER4
10-
p = 0.0034
Relative expression of CYP11B1
2.5
p = 0.1584
8.
2.0
10.
p = 0.0087
p = 0.0064
p = 0.6006
6-
1.5.
5.
4.
1.0
2-
0.5
0
0
0.0
0
0.5
1
0
0.5
1
0
0.5
1
5 Aza (um)
5 Aza (um)
5 Aza (um)
neoplastic changes associated with tumorigenesis. Furthermore, these DNA changes may be easily detected in the plasma or serum of cancer patients, thereby highlighting the potential of DNA methylation as a novel non-invasive molecular marker for
cancer diagnosis and prognosis (41, 42). Kim et al. showed the rate of 3 year survival of the patients with acute myelocytic leukemia depends on their levels of HOXA5 methylation (43). In our study, survival analysis results demonstrated that HOXA5
A
p = 0.0004
B
p = 0.0012
C
p = 0.0056
Relative expression of NR2F1
10-
p = 0.0352
Relative expression of KCNQ1
10
p = 0.0261
Relative expression of HOXA5
15-
p = 0.0310
8
p = 0.0336
8
p = 0.0088
10-
6
6.
p = 0.0052
4.
4.
5.
2
2.
0
0
0
0
0.5
1
0
0.5
1
0
0.5
1
5 Aza (um)
5 Aza (um)
5 Aza (um)
D
E
F
p = 0.0731
p = 0.0008
Relative expression of CD14
10-
p = 0.0460
Relative expression of PTGER4
p = 0.0006
p = 0.1036
15.
Relative expression of CYP11B1
p = 0.0441
2.5.
p = 0.4871
8.
p = 0.0207
2.0
6.
10
p = 0.0027
1.5-
4.
5.
1.0-
2.
0.5-
0
0
0.0
0
0.5
1
0
0.5
1
0
0.5
1
5 Aza (um)
5 Aza (um)
5 Aza (um)
was closely related with overall survival, indicating that HOXA5 could be a predictor of poor prognosis in ACC patients. In addition, compared to the patients with lesser expression of ESM1, those with highly expressed ESM1 were reported to have an obviously worse relapse-free survival in triple- negative breast cancer (44). In our study, we found that ESM1 was hypomethylated and upregulated in ACCs, indicating that aberrant methylation in ACCs might lead to the upregulation of ESM1, resulting in ACC tumorigenesis. In addition, some other genes like KCNQ1, CD14, and NR2F1 were significantly related with overall survival, which might be able to be used as a prognosis indicator for ACC. Besides the methylation, other genetic alterations, such as missense mutation, amplification and deep deletion could also contribute to tumorogenesis. Therefore, not all hypermethylated genes identified in the present study are indeed associated with a higher pathological stage or a poorer survival.
Our quantitative real-time PCR and sodium bisulfite sequencing data strongly indicates that five genes (NR2F1, KCNQ1, HOXA5, CD14, and PTGER4) undergo methylation induced regulation of gene expression in ACC tissues. ACC tissues have a higher number of methylated sites of 36 CpG sites as compared to the normal samples. In contrast, expression of these five genes in normal adrenal tissues were higher than that in ACC tissues. We found that the degree of methylation of the CYP11B1 gene is inconsistent with
our prediction, but its low expression is consistent with our prediction. Because of the multiple ways to inhibit gene expression, this difference may be due to other pathways leading to low expression of CYP11B1 in ACC. However, we have noticed that hypermethylation of the CYP11B1 gene and the loss of its expression may be related to the occurrence of hypercortisolemia (45). Interestingly, among the seven genes we screened, CYP11B1 and HOXA5 are also among the adrenal specific expression of the gene signature defined by Zheng et al. (46). In addition, as Rechache et al. discovered, we also screened out KCNQ1 as one of the genes that are hypermethylated and downregulated in ACC (47). We included the dataset GSE49280 (48) in our research and found that the negative correlation of six downregulated hyper-methylated genes were significant, indicating that methylation might be highly important in gene expression. In order to further explore the expression and localization of HOXA5 gene, we performed novel immunohistochemistry and immunofluorescence experiments in ACC and adjacent normal tissues, and found that the gene is mainly expressed in the nucleus.
In conclusion, we identified seven genes hypo/hyper methylated and dysregulated, which could play important roles in the development and progression of ACCs. Additionally, some genes were found to negatively correlate with different pathological stages or overall rates of survival of ACCs. We
A
p= 0.034
6-
p= 0.04
p = 0.0019
B-
6
Relative expression of NR2F1
Relative expression of KCNQ1
Relative expression of HOXA5
7
C
5
4
ã
3
4-
2
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
3
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
P = 0.038
p = 0.03
p= 0.0037
p= 0.013
7
12
Relative expression of CD14
Relative expression of PTGER4
₪
Relative expression of CYP11B1
11
Relative expression of ESM1
M
6
10
9
3
-
8
4
7
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
2
Normal adrenal gland (n=4)
Adrenal cortical carcinoma (n=4)
B
C
Normal adrenal gland
Adrenal cortical carcinoma
D
E
Normal adrenal gland
Adrenal cortical carcinoma
DAPI + HOXA5
3
DAPI + HOXA5
have verified these seven genes by various experimentation and have demonstrated that they are differentially expressed in ACC. Among six hypermethylated genes we identified in screening, five of them (NR2F1, KCNQ1, HOXA5, CD14, PTGER4) were found hypermethylated through BSP experiments. Our novel
data suggests these abnormally methylated genes could be used as therapeutic targets and biomarkers for ACC patients to be precisely diagnosed and effectively treated. However, further studies are required to validate the functional activities and molecular mechanism of these genes on ACC.
ETHICS STATEMENT
The use of these ACC specimens was approved by the Ethics Committee at Zhongnan Hospital of Wuhan University, and informed consent was obtained from all patients.
AUTHOR CONTRIBUTIONS
HX, PC, and XZ conceived and designed the study. HX, XW, and WH performed the analysis procedures. ZL, GZ, and DX analyzed the results. HX, PC, MH, MD, and XZ contributed to the writing of the manuscript. All authors reviewed the manuscript.
REFERENCES
1. Grumbach MM, Biller BM, Braunstein GD, Campbell KK, Carney JA, Godley PA, et al. Management of the clinically inapparent adrenal mass (“incidentaloma”). Ann Intern Med. (2003) 138:424-9. doi: 10.7326/0003-4819-138-5-200303040-00013
2. Kebebew E, Reiff E, Duh QY, Clark OH, McMillan A. Extent of disease at presentation and outcome for adrenocortical carcinoma: have we made progress? World J Surg. (2006) 30:872-8. doi: 10.1007/s00268-005- 0329-x
3. Kerkhofs TM, Verhoeven RH, Van der Zwan JM, Dieleman J, Kerstens MN, Links TP, et al. Adrenocortical carcinoma: a population-based study on incidence and survival in the Netherlands since 1993. Eur J Cancer. (2013) 49:2579-86. doi: 10.1016/j.ejca.2013.02.034
4. Luton JP, Cerdas S, Billaud L, Thomas G, Guilhaume B, Bertagna X, et al. Clinical features of adrenocortical carcinoma, prognostic factors, and the effect of mitotane therapy. N Engl J Med. (1990) 322:1195-201. doi: 10.1056/NEJM199004263221705
5. Icard P, Goudet P, Charpenay C, Andreassian B, Carnaille B, Chapuis Y, et al. Adrenocortical carcinomas: surgical trends and results of a 253-patient series from the French association of endocrine surgeons study group. World J Surg. (2001) 25:891-7. doi: 10.1007/s00268-001-0047-y
6. Abiven G, Coste J, Groussin L, Anract P, Tissier F, Legmann P, et al. Clinical and biological features in the prognosis of adrenocortical cancer: poor outcome of cortisol-secreting tumors in a series of 202 consecutive patients. J Clin Endocrinol Metab. (2006) 91:2650-5. doi: 10.1210/jc. 2005-2730
7. Allolio B, Fassnacht M. Clinical review: adrenocortical carcinoma: clinical update. J Clin Endocrinol Metab. (2006) 91:2027-37. doi: 10.1210/jc.2005-2639
8. Allolio B, Hahner S, Weismann D, Fassnacht M. Management of adrenocortical carcinoma. Clin Endocrinol. (2004) 60:273-87. doi: 10.1046/j.1365-2265.2003.01881.x
9. Volante M, Buttigliero C, Greco E, Berruti A, Papotti M. Pathological and molecular features of adrenocortical carcinoma: an update. J Clin Pathol. (2008) 61:787-93. doi: 10.1136/jcp.2007.050625
10. Lubitz JA, Freeman L, Okun R. Mitotane use in inoperable adrenal cortical carcinoma. JAMA. (1973) 223:1109-12. doi: 10.1001/jama.1973.03220100011003
11. Pommier RF, Brennan MF. An eleven-year experience with adrenocortical carcinoma. Surgery. (1992) 112:963-70.
12. Bukowski RM, Wolfe M, Levine HS, Crawford DE, Stephens RL, Gaynor E, et al. Phase II trial of mitotane and cisplatin in patients with adrenal carcinoma: a southwest oncology group study. J Clin Oncol. (1993) 11:161-5. doi: 10.1200/JCO.1993.11.1.161
13. Fassnacht M, Allolio B. Clinical management of adrenocortical carcinoma. Best Pract Res Clin Endocrinol Metab. (2009) 23:273-89. doi: 10.1016/j.beem.2008.10.008
FUNDING
This work was supported by the National Natural Science Foundation of China under Grant Nos. 81160086, 81270843, and 81770757. Study funders have no responsibility for manuscript preparation, publishing decision, data collection, data analysis, or study design.
ACKNOWLEDGMENTS
We would like to express our appreciation for the pre-dominant technical support of Ms. Shanshan Zhang and Miss Danni Shan. Also, we would like to thank Kanehisa Laboratories who have developed the KEGG database which has been greatly helpful in our study.
14. Kulis M, Esteller M. DNA methylation and cancer. Adv Genet. (2010) 70:27- 56. doi: 10.1016/B978-0-12-380866-0.60002-2
15. Farkas SA, Milutin-Gasperov N, Grce M, Nilsson TK. Genome-wide DNA methylation assay reveals novel candidate biomarker genes in cervical cancer. Epigenetics. (2013) 8:1213-25. doi: 10.4161/epi.26346
16. Towle R, Truong D, Hogg K, Robinson WP, Poh CF, Garnis C. Global analysis of DNA methylation changes during progression of oral cancer. Oral Oncol. (2013) 49:1033-42. doi: 10.1016/j.oraloncology.2013.08.005
17. Das PM, Singal R. DNA methylation and cancer. J Clin Oncol. (2004) 22:4632- 42. doi: 10.1200/JCO.2004.07.151
18. Jones PA, Baylin SB. The epigenomics of cancer. Cell. (2007) 128:683-92. doi: 10.1016/j.cell.2007.01.029
19. Wright KD, Gilbertson RJ. To infinium, and beyond! Cancer Cell. (2010) 17:419-20. doi: 10.1016/j.ccr.2010.04.020
20. Gao ZH, Suppola S, Liu J, Heikkila P, Janne J, Voutilainen R. Association of H19 promoter methylation with the expression of H19 and IGF-II genes in adrenocortical tumors. J Clin Endocrinol Metab. (2002) 87:1170-6. doi: 10.1210/jcem.87.3.8331
21. Soon PS, Gill AJ, Benn DE, Clarkson A, Robinson BG, McDonald KL, et al. Microarray gene expression and immunohistochemistry analyses of adrenocortical tumors identify IGF2 and Ki-67 as useful in differentiating carcinomas from adenomas. Endocr Relat Cancer. (2009) 16:573-83. doi: 10.1677/ERC-08-0237
22. Demeure MJ, Coan KE, Grant CS, Komorowski RA, Stephan E, Sinari S, et al. PTTG1 overexpression in adrenocortical cancer is associated with poor survival and represents a potential therapeutic target. Surgery. (2013) 154:1405-16. doi: 10.1016/j.surg.2013.06.058
23. Tombol Z, Szabo PM, Molnar V, Wiener Z, Tolgyesi G, Horanyi J, et al. Integrative molecular bioinformatics study of human adrenocortical tumors: microRNA, tissue-specific target prediction, and pathway analysis. Endocr Relat Cancer. (2009) 16:895-906. doi: 10.1677/ERC-09-0096
24. West AN, Neale GA, Pounds S, Figueredo BC, Rodriguez Galindo C, Pianovski MA, et al. Gene expression profiling of childhood adrenocortical tumors. Cancer Res. (2007) 67:600-8. doi: 10.1158/0008-5472.CAN-06-3767
25. Weiss ID, Huff LM, Evbuomwan MO, Xu X, Dang HD, Velez DS, et al. Screening of cancer tissue arrays identifies CXCR4 on adrenocortical carcinoma: correlates with expression and quantification on metastases using (64)Cu-plerixafor PET. Oncotarget. (2017) 8:73387-406. doi: 10.18632/oncotarget.19945
26. Legendre CR, Demeure MJ, Whitsett TG, Gooden GC, Bussey KJ, Jung S, et al. Pathway implications of aberrant global methylation in adrenocortical cancer. PLoS ONE. (2016) 11:e0150629. doi: 10.1371/journal.pone.0150629
27. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy-analysis of affymetrix genechip data at the probe level. Bioinformatics. (2004) 20:307-15. doi: 10.1093/bioinformatics/btg405
28. Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLOS Genet. (2007) 3:1724-35. doi: 10.1371/journal.pgen.0030161
29. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007
30. Li LC, Dahiya R. MethPrimer: designing primers for methylation PCRs. Bioinformatics. (2002) 18:1427-31. doi: 10.1093/bioinformatics/18.11.1427
31. Halberg RB, Larsen MC, Elmergreen TL, Ko AY, Irving AA, Clipson L, et al. Cyp1b1 exerts opposing effects on intestinal tumorigenesis via exogenous and endogenous substrates. Cancer Res. (2008) 68:7394-402. doi: 10.1158/0008-5472.CAN-07-6750
32. Inoue H, Takamori M, Shimoyama Y, Ishibashi H, Yamamoto S, Koshihara Y. Regulation by PGE2 of the production of interleukin-6, macrophage colony stimulating factor, and vascular endothelial growth factor in human synovial fibroblasts. Br J Pharmacol. (2002) 136:287-95. doi: 10.1038/sj.bjp.0704705
33. Honda T, Segi-Nishida E, Miyachi Y, Narumiya S. Prostacyclin-IP signaling and prostaglandin E2-EP2/EP4 signaling both mediate joint inflammation in mouse collagen-induced arthritis. J Exp Med. (2006) 203:325-35. doi: 10.1084/jem.20051310
34. Liu WJ, Huang MX, Guo QL, Chen JH, Shi H. Effect of human cytomegalovirus infection on the expression of Hoxb2 and Hoxb4 genes in the developmental process of cord blood erythroid progenitors. Mol Med Rep. (2011) 4:1307-11. doi: 10.3892/mmr.2011.576
35. Martinez FD. CD14, endotoxin, and asthma risk: actions and interactions. Proc Am Thorac Soc. (2007) 4:221-5. doi: 10.1513/pats.200702-035AW
36. Richardson AL, Wang ZC, De Nicolo A, Lu X, Brown M, Miron A, et al. X chromosomal abnormalities in basal-like human breast cancer. Cancer Cell. (2006) 9:121-32. doi: 10.1016/j.ccr.2006.01.013
37. Hou J, Lambers M, den Hamer B, den Bakker MA, Hoogsteden HC, Grosveld F, et al. Expression profiling-based subtyping identifies novel non-small cell lung cancer subgroups and implicates putative resistance to pemetrexed therapy. J Thorac Oncol. (2012) 7:105-14. doi: 10.1097/JTO.0b013e3182352a45
38. Smits BM, Haag JD, Rissman AI, Sharma D, Tran A, Schoenborn AA, et al. The gene desert mammary carcinoma susceptibility locus Mcs1a regulates Nr2f1 modifying mammary epithelial cell differentiation and proliferation. PLoS Genet. (2013) 9:e1003549. doi: 10.1371/journal.pgen.10 03549
39. Gaston V, Le Bouc Y, Soupre V, Burglen L, Donadieu J, Oro H, et al. Analysis of the methylation status of the KCNQ1OT and H19 genes in leukocyte DNA for the diagnosis and prognosis of Beckwith-Wiedemann syndrome. Eur J Hum Genet. (2001) 9:409-18. doi: 10.1038/sj.ejhg.5200649
40. Zhang J, Liu L, Wang X, Huang Q, Tian M, Shen H. Low-level environmental phthalate exposure associates with urine metabolome
alteration in a Chinese male cohort. Environ Sci Technol. (2016) 50:5953-60. doi: 10.1021/acs.est.6b00034
41. Hao X, Luo H, Krawczyk M, Wei W, Wang W, Wang J, et al. DNA methylation markers for diagnosis and prognosis of common cancers. Proc Natl Acad Sci USA. (2017) 114:7414-9. doi: 10.1073/pnas.1703577114
42. Xu RH, Wei W, Krawczyk M, Wang W, Luo H, Flagg K, et al. Circulating tumour DNA methylation markers for diagnosis and prognosis of hepatocellular carcinoma. Nat Mater. (2017) 16:1155-61. doi: 10.1038/nmat4997
43. Kim SY, Hwang SH, Song EJ, Shin HJ, Jung JS, Lee EY. Level of HOXA5 hypermethylation in acute myeloid leukemia is associated with short-term outcome. Korean J Lab Med. (2010) 30:469-73. doi: 10.3343/kjlm.2010.30.5.469
44. Sagara A, Igarashi K, Otsuka M, Kodama A, Yamashita M, Sugiura R, et al. Endocan as a prognostic biomarker of triple-negative breast cancer. Breast Cancer Res Treat. (2017) 161:269-78. doi: 10.1007/s10549-016-4057-8
45. Kometani M, Yoneda T, Demura M, Koide H, Nishimoto K, Mukai K, et al. Cortisol overproduction results from DNA methylation of CYP11B1 in hypercortisolemia. Sci Rep. (2017) 7:11205. doi: 10.1038/s41598-017-11435-2
46. Zheng S, Cherniack AD, Dewal N, Moffitt RA, Danilova L, Murray BA, et al. Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer Cell. (2016) 30:363. doi: 10.1016/j.ccell.2016. 07.013
47. Rechache NS, Wang Y, Stevenson HS, Killian JK, Edelman DC, Merino M, et al. DNA methylation profiling identifies global methylation differences and markers of adrenocortical tumors. J Clin Endocrinol Metab. (2012) 97:E1004- 13. doi: 10.1210/jc.2011-3298
48. Assie G, Letouze E, Fassnacht M, Jouinot A, Luscap W, Barreau O, et al. Integrated genomic characterization of adrenocortical carcinoma. Nat Genet. (2014) 46:607-12. doi: 10.1038/ng.2953
Conflict of Interest Statement: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Copyright @ 2019 Xiao, He, Chen, Xu, Zeng, Li, Huang, Wang, DiSanto and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.