Developing Models to Predict BRAFV600E and RAS Mutational Status in Papillary Thyroid Carcinoma

1. Introduction

Thyroid carcinomas are among the most prevalent malignancies of the endocrine system, with a notable rise in incidence over the last few decades [1]. The emergence of well-differentiated thyroid tumors such as papillary thyroid carcinoma (PTC) and follicular thyroid carcinoma (FTC) has been linked to alterations in various genes, including BRAF, RAS, RET, and recently discovered gene fusions such as EIFIAX, RET, NTRK1/3, ALK, PAX8-PPARG, RGADA, FGR2, and LTK genes [2]. PTC, which constitutes 80–85% of overall thyroid carcinoma cases, has particularly received an advanced exploration of its genomic landscape [3]. BRAFV600E and RAS mutations are the two most prevalent gene mutations detected in PTC, with prevalence rates of 25–82.3% [4,5] and 11.5–20% [3,6,7], respectively. Being mutually exclusive, these driver gene mutations ultimately led to the same incongruous activation of the mitogen-activated protein kinase (MAPK) pathway [8]. This pathway involves the sequential activation and phosphorylation of RAS, RAF, MEK, and ERK, which are important in the regulation of cellular growth and apoptosis. Dysregulation of the gene associated with this pathway can affect cellular function and promote tumorigenesis. Elevated pERK1/2 expression, as detected through immunohistochemistry staining or western blot analysis, has been considered a proxy indicator for heightened MAPK pathway activity in various malignancies [9]. Interestingly, the expression of pERK1/2 has reportedly differed between BRAF-mutated and RAS-mutated tumors, with the former being more elevated [2]. In addition, RAS-mutated PTC demonstrates concurrent activation of P13K/AKT signaling as well as MAPK activation [2]. These signaling differences result in distinct phenotypes of PTCs, which are characterized by varied clinical and histopathological findings.

Among several BRAF mutations that have been identified in PTC, BRAFV600E constituted most cases. This T1799A point mutation results in the replacement of valine with glutamate and has emerged as a significant clinical determinant due to its association with heightened disease aggressiveness [10]. PTC tumors typically display an inert behavior, with a considerable proportion of patients attaining a survival rate of ten years [11]. However, tumors harboring the BRAFV600E were associated with increased mortality rates, higher rates of recurrence, and resistance to radioiodine treatment [12,13]. Several studies have also linked this mutation to various aggressive pathology features, such as perithyroidal extension, node metastases, and advanced clinical stage [14,15,16]. In contrast to the BRAFV600E mutation, tumors possessing the RAS mutation were associated with less aggressive pathology characteristics, involving follicular histology variant [12], encapsulated tumors [17], minimal disease invasion [15,18], and lower risk of recurrences [12]. Among the three isoforms of the RAS gene, NRAS has been known as the most prevalent gene mutation that is closely related to PTC [19]. While it is less common in the American and European populations, RAS mutation was reported more frequent in the Asian population [20].

A recent discovery by The Cancer Genome Atlas (TCGA) has been able to classify PTC based on its two major driver gene mutations, which are BRAF-like and RAS-like tumors [2]. Identifying PTCs into BRAF-like and RAS-like tumors during diagnostic work-up is essential, not only for determining the most precise and targeted treatment but also to comprehend the biological behavior between the two. Hence, this study aimed to explore the differences between BRAF-like and RAS-like tumors concerning clinicopathology and pERK1/2 expression and further establish predictive models for BRAFV600E and RAS mutations in PTC.

2. Materials and Methods

2.1. Study Design and Population

This study retrospectively enrolled PTC patients who had undergone total thyroidectomy at Cipto Mangunkusumo National Hospital and MRCCC Siloam Hospital between January 2019 and September 2022. We excluded cases with aggressive features such as a high mitosis index and necrosis. The clinical information, including age, gender, and clinical stage, was procured from medical records. Three licensed pathologists gathered the histopathological data blindly, involving tumor size, PTC-nuclear score, capsule, histology variant, multifocality, lymphovascular invasion (LVI), extrathyroidal extension (ETE), and node metastases. Histology variants of PTCs were further classified into non-aggressive (classic and follicular) and aggressive (tall cell, oncocytic, and solid) [13].

2.2. pERK1/2 Immunohistochemistry Examination

The expression of pERK1/2 in this present study were classified into high and low expression based on the cut-off point 10% established by Gomes et al [21]. The standard immunohistochemical evaluation procedures were used to assess the expression of pERK1/2. Positive and negative control was included in each specimen. Colon adenocarcinoma paraffin blocks as a positive control were taken from the routine control archives in our institution. Three-mm thickness of unstained slides were cut and rinsed under running water for 2 minutes following deparaffinization and rehydration. In a de-cloaking chamber at a temperature of 96 °C for 25 minutes, antigen retrieval was carried out using pH 9 Tris-EDTA buffer. After 3 minutes of washing in PBS pH 7.4, a blocking solution containing 10% normal horse serum (Thermo Fisher Scientific, Inc.) was administered for 10 minutes at room temperature to block non-specific protein. Each slide was incubated with rabbit monoclonal anti-Phospho-p44/42 MAPK (Erk1/2) (20G11; Cell Signaling Technology, Danvers, MA, USA) at a dilution ratio of 1:600. Subsequently, each slide was incubated with the PolyVue Plus Mouse/Rabbit Enhancer (Diagnostic Biosystems, USA) for 15 minutes followed by PolyVue Plus Mouse/Rabbit HRP Label for 15 minutes. The slides were repeatedly washed before being incubated to diluted diaminobenzidine chromogen buffer substrate for 1 minute at room temperature. Mayer's hematoxylin was used for a 10-second counterstaining procedure at room temperature. Each slide was examined under a light microscope (Leica Microsystems GmbH) and photographed in five representative fields at x400 magnification with a minimum of 500 tumor cells for each case. Tumor cells-stained brown in the nucleus counted as positive. The quantitative evaluation of pERK1/2 expression was performed by counting the proportion of cells stained positively using ImageJ software version 1.51 (National Institutes of Health). The Kappa inter-observer analysis indicated an agreement of 0.879 (p<0.001), which is near-perfect.

2.3. Mutational Analysis

The tumor specimens were subjected to mutational analysis using Sanger sequencing to detect BRAFV600E mutation as well as N/H/K-RAS codon 12, 13, and 61 mutations. The procedures were performed according to the methodology outlined in the previous study [22].

2.4. Statistical Analysis

The statistical software package SPSS version 20 was utilized for the purpose of data processing. Bivariate analyses were conducted utilizing the Chi-square and Mann-Whitney U test. Clinicopathological variables that showed a P value <0.05 during bivariate analysis were considered as significant variables and were subsequently added to a multivariate analysis. Binary logistic regression test was employed to conduct the multivariate analysis, utilizing a backward: conditional method. The model's goodness of fit was evaluated by conducting the Hosmer-Lemeshow test. A significance level of 0.05 or higher is indicative of a reliable predictive model. The fittest model results in selected clinicopathological variables which act as the predictor for each BRAFV600E and RAS mutational status. The development of a scoring value for each predictor involved the formulation of the coefficient B and S.E. as displayed in the regression test. Following the implementation of the scoring system for the study sample, an analysis of the receiver operating characteristic (ROC) curve was conducted. An area under the receiver operating characteristic curve (AUC) exceeding 0.7 indicates a satisfactory level of diagnostic precision. The scoring wizard tool was utilized to evaluate the probability of each total score in every predictive model. To evaluate the applicable combinations of BRAFV600E and RAS model prediction, a multinomial logistic regression test and internal validation was performed.

3. Results

3.1. Baseline characteristics

The study retrospectively collected PTC patients who had undergone total thyroidectomy from January 2019 to September 2022, with an initial recruitment of 527 patients. A total of 305 patients were excluded from the study for multiple reasons as illustrated in Figure 1. The study consisted of a total of 222 participants.

A predominant demographic profile of this study included females (162, 73%), diagnosed under the age of 55 years (163, 73.4%), at clinical stage 1 (166, 74.8%). Histopathological features are dominated by cases characterized by tumor size less than 4 cm (166, 74.8%), nuclear score of 3 (152, 68.5%), lack of capsules (139, 62.6%), non-aggressive variants (152, 68.5%). The histology variants were dominated by follicular variants (87, 39.2%), followed by classic (65, 29.3%), tall cells (52, 23.9%), oncocytic (10, 4.4%), and solid (7, 3.2%), respectively. The presence of multifocality (168, 75.7%), absence of LVI (135, 60.8%), lack of ETE (158, 71.2%), and lack of node metastases (134, 60.4%) constituted the majority of cases in this study. The expression of pERK 1/2 was quantified to range from 0.2% to 99%, with a median value of 5%. Based on the 10% cut off point [21], this study was dominated with low pERK1/2 expression (145, 65.3%).

3.2. Bivariate analysis: correlation between clinico-histopathology characteristics with BRAFV600E and RAS mutational status

One hundred and sixteen cases that did not exhibit the BRAFV600E and RAS mutation are further designated as control. Two distinct bivariate analyses were performed to compare the BRAFV600E-mutated and control, as well as the RAS-mutated and control. The results were outlined in Table 1 and Table 2. The bivariate analysis conducted on BRAFV600E-mutated group revealed a significant correlation between BRAFV600E mutation and nuclear score 3 (P=0.001; OR=4.1; 95% CI=1.7–9.4), the absence of tumor capsules (P=0.001; OR=3.2; 95% CI=1.5–6.7), tall cell variant (P=0.001; OR=8.9; 95% CI=3.5–22.6), aggressive histology variants (P=0.001; OR=2.9; 95% CI=1.5–5.4), the presence of ETE (P=0.01; OR=2.3; 95% CI=1.2–4.4), the presence of node metastases (P=0.008; OR=2.3; 95% CI=1.2–4.3), and high pERK1/2 expression (P=0.008; OR=2.4; 95% CI=1.2–4.8). There was a significant association between RAS mutation and the follicular histology variant (P=0.001; OR=2.6; 95% CI=1.1–6.2), the non-aggressive histology variants (P=0.001; OR=17; 95% CI=2.2–128.6), as well as high pERK1/2 expression (P=0.001; OR=7.6; 95% CI=3.5–16.7).

3.3. Multivariate analysis: establishing the BRAFV600E prediction model

A multivariate analysis using logistic regression test was conducted to examine the association between BRAFV600E mutation and multiple variables. Nuclear score of 3, the absence of tumor capsules, aggressive histology variants, and high pERK1/2 expression were identified as predictive factors contributing to the presence of BRAFV600E mutation. As indicated in Table 3, the predictor variables were assessed individually to determine their respective score for the development of a BRAFV600E prediction model. The nuclear score of 3, the lack of tumor capsules, and aggressive histology variants each contribute a score of 1. pERK1/2 expression level exceeding 10% corresponds to a score of 2.

The Hosmer-Lemeshow goodness-of-fit test yielded results indicating that the logistic regression model exhibited a favorable level of calibration (X²=0.128, P> 0.05). The receiver operating characteristic (ROC) curve's area under the curve (AUC) was determined to be 0.734 with P-value <0.001 and a 95% CI of 0.661-0.807 (Figure 2). This finding suggests that the logistic regression model exhibits a favorable level of discrimination. Probability, sensitivity, and specificity values for every outcome of the overall score was summarized in Table 4. Based on the results of probability analysis, it was determined that the highest probability, amounting to 82%, is associated with a total score of 5. This indicates that if PTC achieves a score of 5, there is an 82% likelihood of the occurrence of the BRAFV600E mutation.

3.4. Multivariate analysis: establishing the RAS mutation prediction model

Based on the multivariate analysis of RAS mutational status, the predictor variables which may be included for the development of the RAS prediction model were follicular histology variant and pERK1/2 expression exceeding 10%. Each corresponding variable gives a score of 1 based on the B/SE value as summarized in Table 5.

The result of the Hosmer-Lemeshow goodness-of-fit test suggests that the RAS mutation-logistic regression model demonstrated a satisfactory level of calibration (X2=0.174, P>0.05). The AUC of the ROC curve was found to be 0.8, with a P-value <0.001 and a 95% CI of 0.702-0.854 (Figure 3). Probability, sensitivity, and specificity values for every outcome of the overall score generated was summarized in Table 6. According to the findings of the probability analysis, it has been ascertained that the highest probability, comprising 70%, is linked to a cumulative score of 2. This suggests that in the case of PTC obtaining a score of 2, there is a probability of 70% for the presence of the RAS mutation.

3.5. Internal validation: applicating BRAFV600E and RAS mutation prediction model to study samples

Both prediction models were used to internally validate all study samples. The results showed that a sample capable of fulfilling two prediction models had varying probabilities. Consequently, we established four possible combination outcomes based on the score obtained from the combination of BRAFV600E and RAS mutation model (Figure 4).

The BRAFV600E prediction model results were classified into low BRAFV600E-scoring group (total score 0- 2) and high BRAFV600E-scoring group (total score 3-5) based on the specificity value of 65% as the middle threshold for identifying the BRAFV600E mutational status. The RAS prediction model results were classified into low RAS-scoring group (total score of 0-1) and high RAS-scoring group (total score of 2) using a specificity value of 91%.

Table 7 provides a summary of a multinomial analysis on four combination outcomes. Samples with low BRAFV600E-scoring group and low RAS-scoring group, as further known as combination 1, acted as the reference group since they had the greatest proportion of non-BRAFV600E and non-RAS patients. The combination 2 (adjOR=4.857, p=0.01, 95% CI=1.470-16.049), which included samples with low BRAFV600E-scoring group and high RAS-scoring group, was substantially linked to more occurrences of RAS mutation, considered as a RAS-like combination. A strong correlation exists between BRAFV600E mutation and combination 3, which includes samples with high BRAFV600E-scoring group and low RAS-scoring group (adjOR=3.091, p=0.001, 95% CI=1.594-5.995), further considered as BRAF-like combination. The combination 4 group, which contained samples with high BRAFV600E-scoring group and high RAS-scoring group, was shown to have significantly more RAS-mutated patients (adjOR=14.571, p=0.001, 95% CI=4.095-51.855).

4. Discussion

Various genetic alterations have been identified as contributing factors to the development of PTC. The most observed genetic changes in PTC are BRAFV600E and RAS mutations [23]. These driver gene mutations are mutually exclusive and contribute to the aberrant activation of the MAPK pathway [24]. The different signaling cascades associated with BRAFV600E and RAS mutations give rise to the specific phenotypic and behavioral characteristics observed in PTC. Tumors harboring BRAFV600E mutations exhibit a greater propensity for aggressiveness, characterized by an increased likelihood of disease recurrences [25], mortality [12], and resistance to radio-ablation [13]. Conversely, tumors carrying RAS mutations tend to display more indolent behavior [18]. TCGA has emphasized the importance of categorizing PTCs into two distinct subtypes, namely BRAF-like and RAS-like, according to their distinctive biological behaviors [2,26]. This present study aimed to develop a predictive model for BRAFV600E and RAS mutations in PTC using various histopathological features, including the novel PTC-nuclear score and pERK1/2 expression.

Histopathological factors that are known to contribute to the disease's aggressiveness are the presence of tumor multifocality, vascular invasion, perithyroidal soft tissue invasion, and node metastases [27]. The present study provides more evidence for prior research [28,29] that has established an association between these parameters and the BRAFV600E mutation status. A phospho-specific antibody that detects pERK1/2 is also examined in this study to assess the activation of the MAPK pathway on a cellular level. It was documented that pERK1/2 expression exceeding 10% was associated with a higher risk of BRAFV600E mutation. Jung et al have discovered a correlation between BRAF-like tumors and high nuclear scores [16]. The findings of this current investigation align with those of a prior study, which demonstrated an association between PTC-nuclear score of 3 and the BRAFV600E mutation. We identified three features that emerged as significant predictors of the presence of the BRAFV600E mutation. These variables include a nuclear score of 3, the aggressive histology variants, the lack of tumor capsule, and an expression level of pERK1/2 greater than 10%. The scoring system was established utilizing the characteristics, as indicated in Table 1. The nuclear score of 3, the lack of a tumor capsule, and the aggressive histology variants each contribute a score of 1. If the expression of pERK1/2 exceeds 10%, it is assigned a score of 2. The current research demonstrates a positive correlation between higher total scores and an increased probability of observing a BRAFV600E mutation in PTC. The probability of the BRAFV600E mutation is highest at 82% when a set of total five scores is taken into account.

In comparison to the BRAFV600E mutation, tumors harboring the RAS mutation have been linked to a less aggressive pathological phenotype, characterized by a follicular patterned tumor, encapsulated tumors [17], less disease invasion [15,18] and a reduced likelihood of recurrence [18]. Our finding is in line with previous literature, in which RAS mutations were significantly more common in the follicular variant of PTCs. Follicular variant of PTCs is considered the non-aggressive histology variants, which comprises most of low-grade tumors. A significant difference between the signaling pathways of BRAFV600E-mutated and RAS-mutated tumors resides in the lower level of MAPK activity found in RAS-mutated tumors [2]. This present study, however, was able to display a significant association between the enhanced MAPK activity in both BRAFV600E and RAS-mutated PTCs, as displayed in pERK1/2 immunohistochemistry expression. During a multivariate analysis, it was determined that two features, namely the presence of the follicular variant and an expression level of pERK1/2 greater than 10%, were identified as significant predictors for RAS mutation. The scoring system was constructed based on the criteria listed in Table 2. Each variable is assigned a value of 1, resulting in a total score of 2, which contributes to a probability of 70% for RAS mutation.

The provided study sample exhibits the capacity to satisfy two distinct prediction models with differing probabilities for BRAFV600E and RAS mutations, posing challenges in determining mutational status. Hence, this work proposes the concurrent utilization of the BRAFV600E and RAS prediction models in routine clinical applications. All samples were applied to both the BRAFV600E and RAS prediction models for internal validation. The initial utilization of the BRAFV600E prediction model involved its categorization into two distinct groups: a low BRAFV600E-scoring group (score 0–2) and a high BRAFV600E-scoring group (score 3–5). The RAS prediction model consists of two distinct sample groups based on their total scores. One group is characterized by RAS scores ranging from 0 to 1, indicating low RAS-scoring group, while the other group has a uniform RAS score of 2, indicating high RAS-scoring group. Ultimately, four possible outcomes were established, denoted as combinations within the context of this investigation (Figure 4). There are four possible outcome groups: combination 1 involves cases with low scores in both the BRAFV600E and RAS groups, combination 2 involves cases with a low BRAFV600E score and a high RAS score, combination 3 involves cases with a high BRAFV600E score and a low RAS score, and combination 4 involves cases with high scores in both the BRAFV600E and RAS prediction models. The prevalence of combination 1 was seen to be highest among individuals with non-BRAFV600E and non-RAS mutations. Combination 2 exhibited the highest prevalence in samples where RAS mutations were detected, with a statistically significant positive correlation. In combination 3, most of the samples exhibited BRAFV600E mutations, which were shown to be statistically significant. Based on the obtained results, it was determined that the dominant parts of combination 2 and combination 3 were RAS-like and BRAFV600E-like PTCs, respectively. This discovery provides evidence in favor of the proposed hypothesis. In combination 4, a significant association is shown between RAS mutations and a 14-fold increased likelihood compared to non-BRAFV600E non-RAS mutations. However, no association is found with BRAFV600E mutations. On the one hand, this combination exhibits a proclivity toward RAS mutations. However, given that PTC with the BRAFV600E mutation tends to display a more aggressive nature, it is advisable to approach its interpretation with caution to avoid potential undertreatment. Undertreatment refers to a situation in which PTC has aggressive characteristics, although the management approach is based on low-risk criteria, thereby elevating the likelihood of disease recurrence or metastasis.

Hitherto, this study is the sole research that has constructed a predictive model pertaining to gene mutations in PTC, owing to the routine implementation of molecular examination in developed countries. Our findings can map the histopathology characteristics of PTC into BRAF-like and RAS-like tumor as a foundation of biological behavior of the tumor. This present study is limited as it does not include other variables such as mortality, recurrence, distant metastases, and therapy response. Additional external validation studies are required to further assess the predictive model, utilizing larger and more diverse samples as well as incorporating additional variables as previously mentioned.

5. Conclusions

Using clinico-histopathology features and pERK1/2 expression, two distinct predictive models for BRAFV600E and RAS mutational status in PTC were developed. The BRAFV600E prediction model consists of a PTC-nuclear score of 3 (score 1), a lack of capsules (score 1), the aggressive histology variants (score 1), and pERK1/2 expression >10% (score 2). The probability of the BRAFV600E mutation is highest at 82% when a set of total five scores was reached. The RAS prediction model consists of follicular variant (score 1) and pERK1/2 expression > 10% (score 1). BRAF-like tumors are those included in combination 3 (high BRAFV600E-scoring group and low RAS-scoring group), which exhibits a significant 3-fold increase in the BRAFV600E mutation. RAS-like tumors are those belonging to combination 2 (low BRAFV600E-scoring group and high RAS-scoring group) which possess a significant 4.8-fold increase in RAS mutation. These prediction models may serve as a fundamental basis for comprehending the distinct phenotypic and molecular characteristics of BRAF-like and RAS-like PTCs.