Induction of the BRAFV600E Mutation in Thyroid Cells Leads to Frequent Hypermethylation

Article information

Clin Exp Otorhinolaryngol. 2022;15(3):273-282
Publication date (electronic) : 2022 May 4
doi : https://doi.org/10.21053/ceo.2022.00206
1Department of Surgery, Inha University Hospital, Inha University College of Medicine, Incheon, Korea
2Department of Surgery, Seoul National University Hospital, Seoul National University College of Medicine, Seoul, Korea
3Cancer Research Institute, Seoul National University College of Medicine, Seoul, Korea
4Transdisciplinary Department of Medicine and Advanced Technology, Seoul National University Hospital, Seoul, Korea
5Department of Surgery, Seoul National University Boramae Medical Center, Seoul, Korea
6Department of Surgery, Seoul National University Bundang Hospital, Seongnam, Korea
Corresponding author: Kyu Eun Lee Department of Surgery, Seoul National University Hospital, Seoul National University College of Medicine, 101 Daehak-ro, Jongno-gu, Seoul 03080, Korea Tel: +82-2-2072-0393, Fax: +82-2-766-3975 E-mail: kyueunlee@snu.ac.kr
Received 2022 February 10; Revised 2022 March 14; Accepted 2022 March 21.

Abstract

Objectives.

The BRAFV600E mutation is a major driver mutation in papillary thyroid cancer. The aim of this study was to elucidate the correlation between DNA methylation and gene expression changes induced by the BRAFV600E mutation in thyroid cells.

Methods.

We used Nthy/BRAF cell lines generated by transfection of Nthy/ori cells with the wild-type BRAF gene (Nthy/WT cells) and the V600E mutant-type BRAF gene (Nthy/V600E cells). We performed gene expression microarray and DNA methylation array analyses for Nthy/WT and Nthy/V600E cells. Two types of array data were integrated to identify inverse correlations between methylation and gene expression. The results were verified in silico using data from The Cancer Genome Atlas (TCGA) and in vivo through pyrosequencing and quantitative real-time polymerase chain reaction (qRT-PCR).

Results.

In the Nthy/V600E cells, 199,821 probes were significantly hypermethylated, and 697 genes showed a “hypermethylation-downregulation” pattern in Nthy/V600E. Tumor suppressor genes and apoptosis-related genes were included. In total, 66,446 probes were significantly hypomethylated, and 227 genes showed a “hypomethylation-upregulation” pattern in Nthy/V600E cells. Protooncogenes and developmental protein-coding genes were included. In the TCGA analysis, 491/697 (70.44%) genes showed a hypermethylation-downregulation pattern, and 153/227 (67.40%) genes showed a hypomethylation-upregulation pattern. Ten selected genes showed a similar methylation-gene expression pattern in pyrosequencing and qRT-PCR.

Conclusion.

Induction of the BRAFV600E mutation in thyroid cells led to frequent hypermethylation. Anticancer genes, such as those involved in tumor suppression or apoptosis, were downregulated by upstream hypermethylation, whereas carcinogenic genes, such as protooncogenes, were upregulated by hypomethylation. Our results suggest that the BRAFV600E mutation in thyroid cells modulates DNA methylation and results in cancer-related gene expression.

INTRODUCTION

Thyroid cancer is one of the most common solid organ cancers, and its incidence has increased dramatically in the last decade [1,2]. Papillary thyroid carcinoma (PTC) is the most common histologic subtype and is associated with the BRAFV600E somatic mutation, which is known to be associated with poor prognostic factors such as lymph node metastasis, extrathyroidal extension, and advanced stage [3,4]. The BRAFV600E mutation is the highest signal activator of the mitogen-associated protein kinase (MAPK)/ extracellular signal-regulated kinase (ERK) pathway. It activates a carcinogenic signal cascade in thyroid cells and makes them more proliferative [3]. Recent whole-genome sequencing performed by The Cancer Genome Atlas (TCGA) project further revealed that BRAFV600E mutant-type PTC showed significant upregulation of the pERK-DUSP (Erk transcriptional program) pathway, possibly due to the insensitivity of BRAFV600E to ERK inhibition feedback [5].

The molecular mechanism of gene expression alterations resulting from the BRAFV600E driver mutation is not yet fully understood. Beyond genomic transcription, epigenomic changes, such as DNA methylation, or microRNAs can act as important factors that modulate gene expression [6]. DNA hypermethylation or hypomethylation in CpG islands can downregulate or upregulate downstream gene expression [7]. In thyroid cancer, some laboratory studies have suggested that specific aberrant DNA methylation is found in cancer-related genes such as PTEN and RASSF1A [6,8]. Previous studies have reported genomewide integrated analyses of methylation and expression arrays and the sequencing of thyroid cancer tissue or thyroid cancer cell lines [9-15]. However, previous studies could not determine the actual consequences of the primitive BRAFV600E mutation because the cell lines used in previous studies closely resembled dedifferentiated thyroid cancers, such as anaplastic thyroid cancer, rather than the precancerous state [16].

Our aim in this study was to investigate changes in DNA methylation and gene expression in thyroid cells initiated by the primitive BRAFV600E mutation. We used the Nthy/BRAF cell line, which was generated by transfection of either wild-type BRAF (Nthy/WT) or mutant BRAF (Nthy/V600E) into Nthy/ori cells at the authors’ institution [17]. We performed gene expression and methylation microarray analyses and analyzed two types of data to find correlations between DNA methylation and gene expression induced by the BRAFV600E mutation in normal thyroid cells.

MATERIALS AND METHODS

Nthy/BRAF cell line and gene expression microarray

We previously reported Nthy/ori cell lines expressing the BRAF gene with or without V600E mutation [17]. Two types of Nthy/BRAF cells were presented: Nthy/WT (Nthy cells with wild-type BRAF gene) and Nthy/V600E (Nthy cells with V600E mutanttype BRAF gene). Cells were grown in RPMI-1640 (Biowest, Riverside, MO, USA) supplemented with 10% FBS (Biowest), 2 mM GlutaMAX (Gibco; Thermo Fisher Scientific, Waltham, MA, USA), and 100 U/mL penicillin-streptomycin (Gibco) in a humidified atmosphere of 5% CO2 at 37°C. They were cultured in 100-mm dishes until confluent monolayers were reached, and RNA was extracted using an easy-spin total RNA extraction kit (iNtRON Biotechnology, Seoul, Korea). Total RNA was quantified using a Nanodrop ND-1000 spectrophotometer. The Illumina HumanHT-12 v4 Expression Bead Chip (Illumina Inc., San Diego, CA, USA) microarray services were provided by Macrogen (Macrogen Inc., Seoul, Korea). Two arrays were conducted for each cell type (2xNthy/WT and 2xNthy/V600E).

Methylation microarray

For DNA extraction, Nthy/WT and Nthy/V600E cells were cultured until a confluent monolayer was formed. Cells were harvested by trypsin-EDTA treatment, and DNA was extracted using QIAamp DNA Mini Kit (Qiagen, Hilden, Germany). DNA concentration was quantified using a spectrophotometer (NanoDrop, Wilmington, DE, USA). After DNA extraction, methylation microarray services were provided by Macrogen (Macrogen Inc., Seoul, Korea) using Illumina Infinium Methylation EPIC Bead Chip Kits (850K; Illumina Inc., San Diego, CA, USA). Three arrays were conducted for each cell type (3xNthy/WT and 3xNthy/V600E).

Bioinformatic analysis

Microarray data were analyzed according to two groups: Nthy/WT and Nthy/V600E. Detailed methods for gene expression microarray analysis are described elsewhere [17]. To identify differentially expressed genes (DEGs), a moderated t-test using “limma” package was applied [18]. Benjamini-Hochberg method was applied to correct the false positive rate from multiple comparisons. A log fold change of 2.0 was used as the cutoff value for significant DEGs.

Methylation array data were analyzed by “ChAMP” package provided by Bioconductor [19]. To identify differentially methylated probes (DMPs), “limma” method was also applied [18]. False positive rate correction was performed with the Benjamini-Hochberg test. Adjusted P-value under 0.05 was considered statistically significant for both gene expression and methylation analysis. Genes with inverse methylation-expression patterns were identified according to the following criteria: differentially hypermethylated probes in the promoter region with downregulated gene expression or promoter region hypomethylation with upregulated gene expression (Fig. 1). The promoter region was defined as CpG islands in TSS1500, TSS200, the 5´-untranslated region (UTR) and the 1st exon area.

Fig. 1.

Schematic diagram of the correlation analysis between differentially methylated probes (DMPs) in promoters and differentially expressed genes in exons. UTR, untranslated region.

To find functional terms of inversely correlated genes, Database for Annotation, Visualization, and Integrated Discovery (DAVID 6.7) and Swiss-Prot Protein Information Resource (SPPIR) were used. Adjusted P-value less than 0.1 was considered statistically significant. All statistical analyses were performed with the R programming language (https://www.R-project.org/).

Validation of inverse correlations using a dataset from TCGA

To validate genes with inverse correlation patterns using other public databases, we assessed somatic mutations, mRNA expression by RNA sequencing and DNA methylation by Illumina Infinium Methylation 450K from the TCGA data (https://gdc.cancer.gov/). Patients in the TCGA cohort were divided into two groups according to BRAF mutation status. For the inversely correlated genes derived from the microarray analysis, we checked mRNA expression and promoter site methylation patterns between BRAF wild-type PTC versus BRAF mutant-type PTC. If the same hypermethylation-downregulation or hypomethylation-upregulation pattern was observed, it was marked as “TRUE”; otherwise, it was marked as “FALSE.”

Validation of gene expression by quantitative real-time polymerase chain reaction

Cells were cultured in 100-mm dishes until confluent monolayers were formed, and total RNA was extracted using RNeasy Mini Kit (Qiagen) and quantified using Nanodrop One (Thermo Fisher). cDNA was transcribed from 1 µg of RNA using Legene Premium Express 1st strand DNA Synthesis System Kit (LeGene Biosciences, San Diego, CA, USA). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed using PowerTrack SYBR green master mix (Applied Biosystems, Waltham, MA, USA) and a QuantStudio 3 real-time PCR System (Applied Biosystems); all samples were run on the same 96-well plate. Relative mRNA expression levels were calculated by normalizing the value for each gene to that of the housekeeping gene ACTB (beta-actin). The primer sequences for qRT-PCR are described in Supplementary Table 1.

Validation of promoter methylation status by pyrosequencing

Nthy/WT and Nthy/V600E cells were cultured until confluent monolayers were formed. Genomic DNA was extracted using QIAamp DNA mini kit (Qiagen) according to manufacturer’s recommendations. The genomic DNA (gDNA) concentration was quantified using a Nanodrop One (Thermo Fisher), and >1 µg/20 µL of gDNA was provided for pyrosequencing analysis. Pyrosequencing was performed by Macrogen (Seoul, Korea) using PyroMark Q48 Autoptrep (Qiagen) system. Primers for pyrosequencing are described in Supplementary Table 2. EpiTect Fast DNA Bisulfite Kit (QIAGEN) was used for bisulfite treatment. The methylation rate was compared between the two types of pyrosequencing data before and after bisulfite treatment. For plotting the pyrosequencing and qRT-PCR validation results, GraphPad Prism ver. 9.1.0 for Windows (GraphPad, San Diego, CA, USA) was used.

Ethical statement

This is not the study that deals with human or animal subjects, so it doesn’t need to gain the approval from Institutional Review Board/Institutional Animal Care and Use Committee (IACUC).

RESULTS

In total, 5,165 genes were differentially expressed between Nthy/WT and Nthy/BRAF cells; 2,441 DEGs were upregulated and 2,724 DEGs were downregulated in Nthy/V600E cells (Fig. 1) [17]. In the methylation array analysis, 266,267 probes were differentially methylated between the Nthy/WT and Nthy/V600E cells. The detailed number of DMPs according to the genomic structure and geographic region of the CpG island is shown in Fig. 2. In Nthy/V600E cells, 199,821 probes were hypermethylated and 66,446 probes were hypomethylated (Fig. 1).

Fig. 2.

Distribution of differentially methylated probes according to genomic location and CpG site. (A) Number of hypermethylated probes in Nthy/V600E cells; genomic structure. (B) Number of hypermethylated probes in Nthy/V600E cells; geographic regions from the CpG area. (C) Number of hypomethylated probes in Nthy/V600E cells; genomic structure. (D) Number of hypomethylated probes in Nthy/V600E cells; geographic regions from the CpG area. UTR, untranslated region; IGR, intergenic region.

Nthy/V600E cells showed 14,516 hypermethylated DMPs in the promoter region with CpG islands. There were 2,724 significantly downregulated genes in Nthy/V600E cells. According to genomic location, 697 genes showed an inverse correlation pattern (promoter hypermethylation with exonic downregulation) in Nthy/V600E cells (Supplementary Table 3). Table 1 lists the top 20 genes that showed promoter hypermethylation and gene downregulation, sorted by the number of hypermethylated probes of CpG islands in the promoter region. Apoptosis-related genes (FASKT, PTEN, EPB41L3), structural constituent-related genes (FBLN2), DNA repair genes (MAD2L2), and scavenger receptor genes (SCARF2) were included.

Top 20 genes with DNA hypermethylation above the exonic area and downregulation of exons in Nthy/V600E cells

Conversely, the number of hypomethylated DMPs in the promoter region in Nthy/V600E cells was 3,022, and 2,441 genes were significantly upregulated. Furthermore, 227 genes had an inverse correlation pattern of promoter hypomethylation and exon upregulation in Nthy/V600E cells (Supplementary Table 4). Table 2 shows the top 20 genes that have DNA hypomethylation and exon upregulation in Nthy/V600E cells. Cell cycle-related genes (CCND1, BCL2, PSMD5), MAPK activation-related genes (GHR), and transmembrane transporter-related genes (SLC6A15, BCL2, HOXB7, NPAS2, KCNK1, ARHGAP22) were included.

Top 20 genes that had an inverse correlation with hypomethylation and upregulation in Nthy/V600E cells

SP-PIR keywords obtained from 697 hypermethylated-downregulated genes in Nthy/V600E are listed in Table 3. Tumor suppressor genes (RASSF4, PRR5, PLA2G16, CADM1, EFNA1, PYCARD, BIN1, MN1, PTEN, and MTUS1) were downregulated in response to upstream hypermethylation in Nthy/V600E cells. Table 4 lists the SP-PIR terms corresponding to the 227 hypomethylated-upregulated genes in Nthy/V600E cells. Protooncogenes (EGFR, FGF5, CCND1, FLI1, MAFB, BCL2, ETV1, FOXO1, MECOM, CBFB) were upregulated in response to upstream hypomethylation in Nthy/V600E cells. In silico validation using TCGA data showed that 491 of 697 genes (70.44%) presented the same promoter hypermethylation with a gene downregulation pattern in PTC with the BRAFV600E mutation tissue compared with PTC tissue without the BRAFV600E mutation. Furthermore, 153 of 227 genes (67.40%) showed hypomethylation with a gene upregulation pattern in PTC with the BRAFV600E mutation. Details of the TCGA validation results for 924 genes are described in Supplementary Table 5.

SP-PIR keywords obtained from 697 hypermethylated-downregulated genes in Nthy/V600E cells

SP-PIR keywords from 227 hypomethylated-upregulated genes in Nthy/V600E cells

Fig. 3 shows the in vitro validation results for some inversely correlated genes (RUNX3, MEST, TP53INP1, RASSF4, and GPX3 among the hypermethylated genes with downregulation in Nthy/V600E cells; CCND1, BCL2, DUSP6, EGFR, and ZEB1 among the hypomethylated genes with upregulation in Nthy/V600E cells). Pyrosequencing was used for the methylation analysis, and qRT-PCR was used for mRNA expression analysis. RUNX3, MEST, TP53INP1, RASSF4, and GPX3 showed higher promoter methylation rates in Nthy/V600E cells than in control cells (Fig. 3A). RUNX3, MEST, RASSF4, and GPX3 also showed lower mRNA expression in Nthy/V600E cells (Fig. 3B). CCND1, BCL2, DUSP6, EGFR, and ZEB1 showed lower promoter methylation rates in Nthy/V600E cells, and these genes showed higher mRNA expression in Nthy/V600E cells (Fig. 3C and D).

Fig. 3.

Pyrosequencing and quantitative real-time polymerase chain reaction results for selected genes. (A) Methylation status of selected genes that showed hypermethylation and downregulation in Nthy/V600E cells. (B) mRNA expression status of selected genes that showed hypermethylation and downregulation in Nthy/V600E cells. (C) Methylation status of selected genes that showed hypomethylation and upregulation in Nthy/V600E cells. (D) mRNA expression status of selected genes that showed hypomethylation and upregulation in Nthy/V600E cells. Nthy/WT, Nthy/ori cells with the wild-type BRAF gene; Nthy/V600E, V600E mutant-type BRAF gene.

DISCUSSION

The somatic BRAFV600E mutation is a major driver mutation in thyroid cancer, particularly in differentiated PTC. The prevalence of the BRAFV600E mutation in thyroid cancer has been reported to be up to 30%–80%, and it is known to be associated with aggressive characteristics of thyroid cancer [3,20]. The basic molecular mechanism of the BRAFV600E mutation in thyroid cancer is activation of the MAPK/ERK pathway. The activation of MAPK-related gene expression in cells is considered to be a major tumorigenic mechanism of various cancers associated with enhanced cell cycle progression, uncontrolled division, and overproliferation [21,22].

Concurrent genomic and epigenomic changes initiated by the BRAFV600E mutation in thyroid cancer have been specifically demonstrated in previous research. Gene expression and DNA methylation changes in MAPK pathway genes are associated with specific clinicopathologic characteristics of thyroid cancer [10,23]. The TCGA project and other next-generation sequencing studies have also revealed that specific sets of gene expression alteration patterns are present in thyroid cancer, followed by the BRAFV600E mutation [5,24]. According to TCGA, BRAFV600E-like thyroid cancer involves upregulated ERK-related genes, including DUSP, which initiates Erk transcription and indicates a low thyroid differentiation score. The methylation signature was also different between BRAF-type and RAS-type thyroid cancers [5]. However, the intracellular mechanism by which the BRAFV600E driver mutation alters the expression of various genes in thyroid cells has not yet been clearly elucidated. A TCGA analysis showed two DNA methylation clusters in BRAFV600E-like thyroid cancer, but the molecular consequences of CpG methylation and gene expression changes were not evaluated [5].

DNA methylation is an epigenomic phenomenon in which a methyl group is attached to the fifth carbon of the cytosine residue in CpG dinucleotides. Clusters of CpG dinucleotides are called CpG islands and are generally located in promoter regions. Hypermethylation of promoter CpG islands represses the transcription of downstream exons and results in suppressed gene expression [7,25]. Hou et al. initially proposed that the BRAFV600E mutation may alter MAPK gene expression through DNA methylation [14]. They reported a functional association of the BRAFV600E mutation with methylation changes in thyroid cancer cells and suggested that the BRAFV600E mutation could facilitate epigenetic changes. They used the BCPAP and OCUT1 cell lines, which harbor the BRAFV600E mutation [14]. In the control group, they used the shRNA method to knock down the BRAF gene in the thyroid cancer cell lines. Studies have demonstrated a correlation between DNA methylation and gene expression in thyroid cancer according to BRAFV600E mutation status using patient-derived PTC tissue and normal thyroid tissue as controls [10,12,26-28]. These studies found that hypermethylation of tumor suppressor genes (PTEN, RASSF1A, TIMP3, SLC5A8, DAPK, RARβ2, etc.) and hypomethylation of oncogenes (VEGF, etc.) were identified in thyroid cancer with the BRAFV600E mutation. However, previous studies used thyroid cancer cells that had already differentiated into cancer rather than cells that were in the precancerous state or continuously differentiated from normal cells. Thus, those studies could not determine the order of the consequences of the primitive BRAFV600E mutation at the genomic, methylation, and expression levels. In particular, those studies were not able to clarify whether mutation or methylation came first. Although Hou et al. used BRAF knockdown PTC cell lines, they were still not representative of the initial period of carcinogenesis induced by the BRAFV600E mutation. As such, previous studies have limitations in explaining the cellular epigenomic and genomic consequences of the BRAFV600E mutation.

To overcome the abovementioned limitation, we used our BRAF gene-transfected Nthy/ori cell lines (Nthy/WT and Nthy/V600E) [17]. In our initial research, Nthy/WT cells showed similar biological behavior to Nthy/ori cells. In contrast, Nthy/V600E cells showed increased anchorage-independent growth and invasion ability and upregulated genes with the ERK/MAPK cascade over time, gradually becoming thyroid cancer cells. We suggest that Nthy/V600E cells are a good model to study early-stage biological consequences induced by the BRAFV600E mutation in thyroid cells. Therefore, we can safely conclude that the methylation changes in our study are processes that follow the initial BRAF mutation.

Our research showed that induction of the BRAFV600E mutation in Nthy/ori cells resulted in frequent DNA hypermethylation; 199,821 probes were hypermethylated in Nthy/V600E cells, while 66,446 probes were hypermethylated in Nthy/WT cells. In the correlation analysis between methylation and gene expression, 697 genes were hypermethylated and downregulated in Nthy/V600E cells. These genes included FASTK, a strong inducer of apoptosis and the immune response; PTEN, a tumor suppressor gene in many cancers; SGCE, whose expression is significantly increased by CpG methylation in gastric cancer; PAOX, whose inhibition results in the rapid death of progenitor cells; MEST, whose loss of expression with methylation is related to carcinogenesis; and RUNX3, whose site-specific hypermethylation predicts PTC recurrence. In the protein information analysis, tumor suppressor genes that repress carcinogenesis were also downregulated in Nthy/V600E cells. In contrast, 227 genes were hypomethylated and overexpressed in Nthy/V600E cells. The functions of these genes were related to transporting glucose and other sugars and hypermethylation in cancer (SLC6A15), blocking apoptotic death and differential methylation in breast cancer (BCL2), and altering cell cycle progression and contributing to tumorigenesis and hypomethylated-high gene expression (CCND1). Protooncogenes that promote carcinogenesis were significantly enriched in the protein information analysis.

We selected 10 genes for validation analysis as follows: RUNX3, MEST, TP53INP1, RASSF4, GPX3, CCND1, BCL2, DUSP6, EGFR, and ZEB1. RUNX is a tumor suppressor gene, whose methylation in thyroid cancer can predict the recurrence of PTC. DNA methylation following MEST gene expression knockdown has been proposed to facilitate thyroid cancer cell survival. For TP53IP1, loss of p53 function is important for immortalization in thyroid cancer cell lines. Methylation of Ras association domain family (RASSF) genes was reported in thyroid cancer. GPX3 is frequently methylated in human PTC tissue, and its expression is associated with tumor size and lymph node metastasis. Enhanced CCND1 gene expression has been proposed as a prognostic factor and possible mechanism for recurrence in PTC. Higher expression of the BCL2 gene was associated with poor survival in thyroid cancer. Increased expression of the DUSP6 gene with decreased DNA methylation was reported in PTC. Overexpression of EGFR in thyroid cancer is well known and associated with tumor aggressiveness and progression. Overexpression of ZEB1 was associated with aggressive tumor progression of PTC. For those genes, the same methylation alteration patterns were obtained in promoter sites by pyrosequencing, and inverse correlations of mRNA expression were also identified, as shown in Fig. 3. Therefore, we propose that our results for inversely correlated genes from our microarrays are reliable, as the important genes were validated by pyrosequencing and qRT-PCR (Fig. 3).

We checked the methylation and gene expression changes of significant genes in the TGCA data, part of in silico validation. We identified that 491 of 697 genes (70.44%) showed the same hypermethylation-downregulation pattern, and 153 of 227 genes (67.40%) showed hypomethylation-upregulation in BRAFV600E mutated PTC tissue (Supplementary Table 5). The remaining 30% of genes did not show significant methylation-gene expression changes. This difference is probably because one experiment was performed on cell lines and another experiment was performed on surgically removed cancer tissue. Although the cancer cell line has very similar properties to those of the cancer tissue from which it originated, it is inevitably genetically different from cancer tissue growing in the human body as a result of being grown in an in vitro environment. Therefore, three-dimensional cell line models such as organoids are being used in recent studies. We suggest that future research should investigate these issues using three-dimensional cell lines.

Based on our analysis, we suggest that initiation of the BRAFV600E mutation in thyroid cells can modulate cancer-related gene expression by altering promoter site CpG island methylation. The possible carcinogenic mechanism is illustrated in Fig. 4. There have been many reports of methylation changes caused by BRAF mutations in various cancers. However, there has been relatively little research on the mechanism through which the BRAF mutation directly and indirectly affects DNA methylation. One study suggested that the BRAF mutation results in the deregulation of the TET genes (TET1, TET2, and TET3), which encode DNA demethylases, leading to demethylation in colon cancer [29]. Another report suggested that the BRAF oncoprotein activates the transcriptional repressor MAFG to mediate CpG island methylation [30].

Fig. 4.

Possible carcinogenic mechanisms of the BRAFV600E mutation.

The BRAFV600E mutation suppresses antitumor-related gene expression by promoting hypermethylation, but it enhances the expression of genes with tumorigenic effects by hypomethylation. In particular, our findings strongly support this argument, unlike previous research, as we experimented with both gene expression and methylation arrays in normal thyroid cells with cancer progression induced by the BRAFV600E driver mutation, which is the highest trigger in the MAPK pathway. This study is also the first report that used the latest methylation array chip, Infinium MethylationEPIC BeadChip (850k) in a thyroid cancer study. Using this higher-density chip, we were also able to perform a more detailed analysis of promoter site methylation than previous studies that used 27K or 450K chips [10,27,28].

The results of our study are suggested to have the following clinical implications. If a test method using a DNA methylation biomarker related to the BRAF mutation is developed and applied in the clinical field, the diagnostic accuracy for ambiguous thyroid nodules could be improved. Thus, unnecessary re-testing or diagnostic surgery could be reduced. In addition, the BRAF mutation is expected to be used as a biomarker to determine the prognosis and future treatment policies for patients undergoing surgery for thyroid cancer. Finally, since methylation changes are epigenetic, meaning that they are easier to modify than DNA mutations or gene expression alterations, there may be the potential to develop a methylation-targeted therapeutic agent to control CpG islands located on carcinogenic genes.

This study has certain limitations. Although our cell lines are good models for studying the genomic consequences of the BRAFV600E mutation, they are not yet standardized worldwide. Further studies including previously known BRAF mutationdominant or RAS mutation-dominant thyroid cancer cell lines are needed to support our results [16]. In addition, our study could not conclude that methylation changes due to BRAF gene mutations can be reversed by BRAF gene suppression (e.g., with BRAF inhibitor treatment). Finally, in our study, we did not identify the direct mechanisms of how the BRAF mutation affects methylation in thyroid cancer cell lines, because our study dealt with associations using two types of microarrays. We intend to address these limitations in our subsequent research.

In summary, induction of the BRAFV600E mutation in normal thyroid cells changed gene expression by affecting the frequency of DNA methylation of promoter site CpG islands. This result suggests that the BRAFV600E mutation modulates DNA methylation to change gene expression in the carcinogenic cascade. In future research, this result will be used to identify potential diagnostic and therapeutic targets in thyroid cancer.

HIGHLIGHTS

▪ Induction of the BRAFV600E mutation in thyroid cells led to frequent hypermethylation, which was associated with the downregulation of anticancer genes, such as those involved in tumor suppression and apoptosis.

▪ Carcinogenic genes, such as protooncogenes, were upregulated by hypomethylation.

▪ Our results suggest that the BRAFV600E mutation in thyroid cells modulates DNA methylation and results in cancer-related gene expression.

Notes

No potential conflict of interest relevant to this article was reported.

AUTHOR CONTRIBUTIONS

Conceptualization: SJK, KEL. Data curation: JWY, SYH, KK, SJK, YJC, JYC. Formal analysis: JWY, KK. Funding acquisition: JWY, SJK. Methodology: YJC, JYC, KEL. Project administration: KEL. Visualization: JWY. Writing–original draft: JWY, SYH, KEL. Writing–review & editing: all authors.

Acknowledgements

This study was funded by research grant from Seoul National University Hospital (no. 04-2013-0770), research grants from Inha University Hospital and College of Medicine, and National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. NRF-2021R1G1A1091596).

Supplementary Materials

Supplementary materials can be found via https://doi.org/10.21053/ceo.2022.00206.

Supplementary Table 1.

Sequences of primers used for quantitative real-time polymerase chain reaction (qRT-PCR)

ceo-2022-00206-suppl1.xlsx
Supplementary Table 2.

Primer design for pyrosequencing

ceo-2022-00206-suppl2.xlsx
Supplementary Table 3.

Genes that had inverse correlations with hypermethylated downregulated genes in Nthy/V600E cells

ceo-2022-00206-suppl3.xlsx
Supplementary Table 4.

Genes that had inverse correlations with hypomethylated upregulated genes in Nthy/V600E cells

ceo-2022-00206-suppl4.xlsx
Supplementary Table 5.

Validation of the inverse correlation patterns of 924 genes in Nthy/V600E cells and The Cancer Genome Atlas papillary thyroid carcinoma (TCGA PTC) tissue with the BRAF mutation

ceo-2022-00206-suppl5.xlsx

References

1. Kitahara CM, Sosa JA. The changing incidence of thyroid cancer. Nat Rev Endocrinol 2016;Nov. 12(11):646–53.
2. Jung KW, Won YJ, Kong HJ, Lee ES, ; Community of PopulationBased Regional Cancer Registries. Cancer statistics in Korea: incidence, mortality, survival, and prevalence in 2015. Cancer Res Treat 2018;Apr. 50(2):303–16.
3. Xing M. BRAF mutation in thyroid cancer. Endocr Relat Cancer 2005;Jun. 12(2):245–62.
4. Li C, Lee KC, Schneider EB, Zeiger MA. BRAF V600E mutation and its association with clinicopathological features of papillary thyroid cancer: a meta-analysis. J Clin Endocrinol Metab 2012;Dec. 97(12):4559–70.
5. Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell 2014;Oct. 159(3):676–90.
6. Faam B, Ghaffari MA, Ghadiri A, Azizi F. Epigenetic modifications in human thyroid cancer. Biomed Rep 2015;Jan. 3(1):3–8.
7. Bird A. DNA methylation patterns and epigenetic memory. Genes Dev 2002;Jan. 16(1):6–21.
8. Xing M. Gene methylation in thyroid tumorigenesis. Endocrinology 2007;Mar. 148(3):948–53.
9. Chen YC, Gotea V, Margolin G, Elnitski L. Significant associations between driver gene mutations and DNA methylation alterations across many cancer types. PLoS Comput Biol 2017;Nov. 13(11):e1005840.
10. Ellis RJ, Wang Y, Stevenson HS, Boufraqech M, Patel D, Nilubol N, et al. Genome-wide methylation patterns in papillary thyroid cancer are distinct based on histological subtype and tumor genotype. J Clin Endocrinol Metab 2014;Feb. 99(2):E329–37.
11. Beltrami CM, Dos Reis MB, Barros-Filho MC, Marchi FA, Kuasne H, Pinto CA, et al. Integrated data analysis reveals potential drivers and pathways disrupted by DNA methylation in papillary thyroid carcinomas. Clin Epigenetics 2017;May. 9:45.
12. Kikuchi Y, Tsuji E, Yagi K, Matsusaka K, Tsuji S, Kurebayashi J, et al. Aberrantly methylated genes in human papillary thyroid cancer and their association with BRAF/RAS mutation. Front Genet 2013;Dec. 4:271.
13. White MG, Nagar S, Aschebrook-Kilfoy B, Jasmine F, Kibriya MG, Ahsan H, et al. Epigenetic alterations and canonical pathway disruption in papillary thyroid cancer: a genome-wide methylation analysis. Ann Surg Oncol 2016;Jul. 23(7):2302–9.
14. Hou P, Liu D, Xing M. Genome-wide alterations in gene methylation by the BRAF V600E mutation in papillary thyroid cancer cells. Endocr Relat Cancer 2011;Nov. 18(6):687–97.
15. Ozer B, Sezerman OU. A novel analysis strategy for integrating methylation and expression data reveals core pathways for thyroid cancer aetiology. BMC Genomics 2015;16(Suppl 12):S7.
16. Saiselet M, Floor S, Tarabichi M, Dom G, Hebrant A, van Staveren WC, et al. Thyroid cancer cell lines: an overview. Front Endocrinol (Lausanne) 2012;Nov. 3:133.
17. Kim BA, Jee HG, Yi JW, Kim SJ, Chai YJ, Choi JY, et al. Expression profiling of a human thyroid cell line stably expressing the BRAFV600E mutation. Cancer Genomics Proteomics 2017;Jan. 14(1):53–67.
18. 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;Apr. 43(7):e47.
19. Morris TJ, Butcher LM, Feber A, Teschendorff AE, Chakravarthy AR, Wojdacz TK, et al. ChAMP: 450k chip analysis methylation pipeline. Bioinformatics 2014;Feb. 30(3):428–30.
20. Nikiforova MN, Kimura ET, Gandhi M, Biddinger PW, Knauf JA, Basolo F, et al. BRAF mutations in thyroid tumors are restricted to papillary carcinomas and anaplastic or poorly differentiated carcinomas arising from papillary carcinomas. J Clin Endocrinol Metab 2003;Nov. 88(11):5399–404.
21. Hilger RA, Scheulen ME, Strumberg D. The Ras-Raf-MEK-ERK pathway in the treatment of cancer. Onkologie 2002;Dec. 25(6):511–8.
22. Peyssonnaux C, Eychene A. The Raf/MEK/ERK pathway: new concepts of activation. Biol Cell 2001;Sep. 93(1-2):53–62.
23. Giordano TJ, Kuick R, Thomas DG, Misek DE, Vinco M, Sanders D, et al. Molecular classification of papillary thyroid carcinoma: distinct BRAF, RAS, and RET/PTC mutation-specific gene expression profiles discovered by DNA microarray analysis. Oncogene 2005;Oct. 24(44):6646–56.
24. Yoo SK, Lee S, Kim SJ, Jee HG, Kim BA, Cho H, et al. Comprehensive analysis of the transcriptional and mutational landscape of follicular and papillary thyroid cancers. PLoS Genet 2016;Aug. 12(8):e1006239.
25. Jones PA, Baylin SB. The epigenomics of cancer. Cell 2007;Feb. 128(4):683–92.
26. Lee EK, Chung KW, Yang SK, Park MJ, Min HS, Kim SW, et al. DNA methylation of MAPK signal-inhibiting genes in papillary thyroid carcinoma. Anticancer Res 2013;Nov. 33(11):4833–9.
27. Rodriguez-Rodero S, Fernandez AF, Fernandez-Morera JL, CastroSantos P, Bayon GF, Ferrero C, et al. DNA methylation signatures identify biologically distinct thyroid cancer subtypes. J Clin Endocrinol Metab 2013;Jul. 98(7):2811–21.
28. Cai LL, Liu GY, Tzeng CM. Genome-wide DNA methylation profiling and its involved molecular pathways from one individual with thyroid malignant/benign tumor and hyperplasia: a case report. Medicine (Baltimore) 2016;Aug. 95(35):e4695.
29. Noreen F, Kung T, Tornillo L, Parker H, Silva M, Weis S, et al. DNA methylation instability by BRAF-mediated TET silencing and lifestyle-exposure divides colon cancer pathways. Clin Epigenetics 2019;Dec. 11(1):196.
30. Fang M, Ou J, Hutchinson L, Green MR. The BRAF oncoprotein functions through the transcriptional repressor MAFG to mediate the CpG Island Methylator phenotype. Mol Cell 2014;Sep. 55(6):904–15.

Article information Continued

Fig. 1.

Schematic diagram of the correlation analysis between differentially methylated probes (DMPs) in promoters and differentially expressed genes in exons. UTR, untranslated region.

Fig. 2.

Distribution of differentially methylated probes according to genomic location and CpG site. (A) Number of hypermethylated probes in Nthy/V600E cells; genomic structure. (B) Number of hypermethylated probes in Nthy/V600E cells; geographic regions from the CpG area. (C) Number of hypomethylated probes in Nthy/V600E cells; genomic structure. (D) Number of hypomethylated probes in Nthy/V600E cells; geographic regions from the CpG area. UTR, untranslated region; IGR, intergenic region.

Fig. 3.

Pyrosequencing and quantitative real-time polymerase chain reaction results for selected genes. (A) Methylation status of selected genes that showed hypermethylation and downregulation in Nthy/V600E cells. (B) mRNA expression status of selected genes that showed hypermethylation and downregulation in Nthy/V600E cells. (C) Methylation status of selected genes that showed hypomethylation and upregulation in Nthy/V600E cells. (D) mRNA expression status of selected genes that showed hypomethylation and upregulation in Nthy/V600E cells. Nthy/WT, Nthy/ori cells with the wild-type BRAF gene; Nthy/V600E, V600E mutant-type BRAF gene.

Fig. 4.

Possible carcinogenic mechanisms of the BRAFV600E mutation.

Table 1.

Top 20 genes with DNA hypermethylation above the exonic area and downregulation of exons in Nthy/V600E cells

Gene name Number of upstream hypermethylated probes Log FC T-value Adjusted P-value Selected gene function from the Gene Ontology database
SGCE 22 –0.872 –11.014 0.002 Calcium ion binding, cytoskeleton, cytoplasm
EPB41L3 15 –1.364 –19.244 0.001 Cell-cell junction, apoptotic process, regulation of cell shape
FASTK 12 –0.658 –5.482 0.017 Serine/threonine kinase activity, apoptotic signaling pathway
DPYSL4 11 –1.232 –15.216 0.001 Pyrimidine nucleobase catabolic process, axon guidance
PTEN 11 –0.495 –6.981 0.009 Apoptotic process, epidermal growth factor signaling pathway
KIAA1217 10 –1.020 –10.900 0.003 Cytoplasm, embryonic skeletal system development
NTF3 10 –0.398 –6.606 0.010 Activation of MAPK activity
FAM163A 10 –0.344 –5.730 0.015 Integral to membrane
PAOX 9 –1.329 –14.875 0.001 Cellular nitrogen compound metabolic process
RUNX3 9 –0.902 –11.945 0.002 Negative regulation of transcription and cell cycle
ZFP3 9 –0.392 –6.522 0.010 DNA-templated, regulation of transcription
PRMT6 9 –0.335 –4.078 0.040 Histone methylation, negative regulation of transcription
UNC45A 9 –0.294 –4.076 0.040 Cell differentiation, perinuclear region of cytoplasm
MEST 9 –0.250 –4.026 0.041 Endoplasmic reticulum, regulation of lipid storage
CYP24A1 8 –1.604 –4.462 0.031 Vitamin D receptor signaling pathway
FAM110A 8 –0.539 –8.772 0.005 Cytoplasm, microtubule organizing center
FBLN2 8 –0.466 –4.293 0.035 Extracellular matrix structural constituent
SAMD11 8 –0.397 –4.283 0.035 Negative regulation of transcription from RNA polymerase II
MAD2L2 8 –0.317 –4.440 0.031 Negative regulation of transcription, DNA repair
SCARF2 7 –1.081 –15.262 0.001 Scavenger receptor activity

FC, fold change; MAPK, mitogen-associated protein kinase.

Table 2.

Top 20 genes that had an inverse correlation with hypomethylation and upregulation in Nthy/V600E cells

Gene name Number of upstream hypomethylated probes Log FC T-value Adjusted P-value Selected gene function from the Gene Ontology database
CCND1 8 1.184 11.630 0.002 Mitotic cell cycle, positive regulation of protein phosphorylation
THBD 7 0.570 9.215 0.004 Negative regulation of platelet activation, leukocyte migration
NTNG1 6 0.594 7.708 0.007 Axonogenesis, anchored to plasma membrane
NPAS2 5 2.031 24.169 0.000 DNA binding transcription factor, signal transducer
SLC6A15 5 0.961 8.960 0.004 Proline:sodium symporter activity, neurotransmitter transporter
HOXB7 5 0.667 10.938 0.003 Sequence-specific DNA binding transcription factor activity
BCL2 5 0.666 11.088 0.002 G1/S transition of mitotic cell cycle, negative regulation of apoptotic process and mitotic cell cycle
TMEM108 5 0.439 7.422 0.007 Integral to membrane
SPRY4 4 2.731 18.910 0.001 Negative regulation of MAPK activity
SATB2 4 1.690 20.881 0.001 Negative regulation of transcription from RNA polymerase II
GAD1 4 1.259 14.122 0.001 Glutamate decarboxylation to succinate
KCNK1 4 1.077 16.265 0.001 Potassium ion transport, synaptic transmission
GHR 4 0.939 15.545 0.001 Activation of MAPK activity, activation of JAK2 kinase activity
ARHGAP22 4 0.897 9.592 0.004 Regulation of small GTPase mediated signal transduction
PSMD5 4 0.705 11.726 0.002 Signal transduction by p53 class mediator resulting in cell cycle arrest
SLAIN1 4 0.681 7.139 0.008 NA
HAS3 4 0.525 4.934 0.023 Positive regulation of transcription, hyaluronan biosynthetic process
RIPPLY2 4 0.435 6.315 0.011 Ossification, somitogenesis
MRPS14 4 0.295 4.582 0.029 Structural constituent of ribosome, mitochondrion
CPEB1 4 0.264 4.027 0.041 Cytoplasmic mRNA processing

FC, fold change; MAPK, mitogen-associated protein kinase; NA, not applicable.

Table 3.

SP-PIR keywords obtained from 697 hypermethylated-downregulated genes in Nthy/V600E cells

Term P-value Fold enrichment Gene
Glyoxylate bypass 0.069 28.370 IDH2, IDH1
Stress-induced protein 0.030 10.639 STK25, HSPB1, SERPINH1
Ehlers-Danlos syndrome 0.055 7.737 COL1A1, B4GALT7, ADAMTS2
Cholesterol biosynthesis 0.004 7.466 CYB5R3, TM7SF2, HMGCR, CYP51A1, FDFT1
Sterol biosynthesis 0.011 5.674 CYB5R3, TM7SF2, HMGCR, CYP51A1, FDFT1
Tyrosine-specific phosphatase 0.004 5.491 PTPRM, PTPRE, PTPN2, PTPRA, PTPRN, CDC25B
Dwarfism 0.010 4.480 COL9A2, FGFR3, PTH1R, DYM, COL1A1, FLNB
Steroid biosynthesis 0.013 4.256 CYB5R3, TM7SF2, HMGCR, CYP51A1, ACAT2, FDFT1
Homotetramer 0.087 3.783 ASS1, GUSB, ALDH2, FUCA1
Triple helix 0.094 3.661 COL9A2, COL6A1, COL1A1, COL16A1
Hydroxylysine 0.094 3.661 COL9A2, COL6A1, COL1A1, COL16A1
Phospholipid biosynthesis 0.059 3.377 CHKA, CPT1B, ISYNA1, CRLS1, ETNK2
Fatty acid biosynthesis 0.064 3.299 PTGIS, SCD, ELOVL2, FADS3, FADS2
Phosphoric monoester hydrolase 0.011 3.289 PTPRM, PTPRE, PTPN2, PTPRA, PGAM1, PPP3CC, PTPRN, CDC25B
Lipid synthesis 0.002 3.184 CYB5R3, TM7SF2, A4GALT, PTGIS, HMGCR, CYP51A1, SCD, ELOVL2, FADS3, FADS2, FDFT1
Copper 0.050 2.986 ATOX1, COX17, MT1X, MOXD1, APLP1, MT1F
Tyrosine-specific protein kinase 0.088 2.955 DDR1, FGFR4, FGFR3, ERBB2, ROR2
LIM domain 0.039 2.797 LIMS2, LMO1, PRICKLE1, LIMK2, PDLIM1, ISL1, TES
Collagen 0.019 2.688 CTHRC1, COL9A2, C1QL1, COL6A1, SCARA3, COL1A1, C1QL4, COL16A1, WDR33
GTP binding 0.087 2.541 RAB31, TUBB2B, GNB2, GNAI1, TUBB2A, RHEB
Heterodimer 0.029 2.479 INHBB, HLA-H, TUBB2B, TUBB2A, ITGB4, TUBA4A, PPP3CC, HLA-B, ADD3
Blocked amino end 0.043 2.467 SRI, CYB5R3, PFN1, GSTM3, SERPINB6, GNAI1, CSTB, COL6A1
Growth factor 0.017 2.382 INHBB, VEGFB, FGF18, PDGFB, NTF3, GDF6, NRG1, GDF15, VGF, MDK, NGF
Tyrosine-protein kinase 0.042 2.300 DDR1, FGFR4, FGFR3, CLK3, ERBB3, ERBB2, TESK1, ROR2, EPHB4
Prenylation 0.022 2.182 RND2, PALM, RAB31, PRICKLE1, MRAS, RHEB, RAB15, RAB6B, RHOD, RAB20, RASD2, GNG7
Protein phosphatase 0.040 2.182 MTMR3, DUSP18, PTPRM, PPM1E, PTPRE, PTPN2, PTPRA, PPP3CC, PTEN, CDC25B
Tumor suppressor 0.053 2.071 RASSF4, PRR5, PLA2G16, CADM1, EFNA1, PYCARD, BIN1, MN1, PTEN, MTUS1
Actin-binding 0.007 2.067 MYO1C, HIP1R, SPIRE1, TMSB10, MTSS1 L, PALLD, TPM2, TPM1, FLNB, DSTN, PFN1, EPB41L3, CORO1A, AIF1 L, ADD3, FHOD1, ADD1, SNTA1
NADP 0.049 2.000 BLVRA, ME1, TM7SF2, CYP24A1, HMGCR, CYP51A1, IDH2, IDH1, NADK, HSD11B1 L, FDFT1

SP-PIR, Swiss-Prot Protein Information Resource.

Table 4.

SP-PIR keywords from 227 hypomethylated-upregulated genes in Nthy/V600E cells

Term P-value Fold enrichment Gene
Proto-oncogene 0.001 3.767 EGFR, FGF5, CCND1, FLI1, MAFB, BCL2, ETV1, FOXO1, MECOM, CBFB
Chromosomal rearrangement 0.000 3.727 SATB2, CCND1, FLI1, BRD3, BCL2, MSI2, ETV1, FOXO1, MECOM, CHCHD7, HMGA1, CBFB
Cytoplasmic vesicle 0.020 2.912 SLC2A8, TM9SF1, ICA1, NPTX1, CADPS2, PRDX6, AGTRAP, PHLDA1
Homeobox 0.022 2.864 TSHZ3, IRX3, SATB2, MSX1, HOXB7, ZHX2, MKX, ZEB1
Developmental protein 0.000 2.447 TSHZ3, FZD8, IRX3, SATB2, EGFL7, ANO1, FOXA1, NTNG1, ANPEP, MECOM, PPDPF, SPRY4, ZIC2, ARHGAP22, NOTCH1, MSX1, HOXB7, SFRP1, POGK, RIPPLY2, MKX, TWIST2
Activator 0.007 2.333 MAFB, FOXA1, PPARG, KLF15, CPEB1, ZEB1, CITED4, NOTCH1, FLI1, GATA3, DNER, ETV1, NFIL3, ETV4
Repressor 0.029 2.191 IRX3, SATB2, MSX1, MAFB, FOXA1, ZHX2, ZNF503, CPEB1, ZEB1, NFIL3, TWIST2
DNA binding 0.099 2.039 MSX1, FLI1, HOXB7, PPARG, NR3C2, ZEB1, NR2F2, HMGA1

SP-PIR, Swiss-Prot Protein Information Resource.