A total of 428 patients aged from 4 to 49 years at the time of surgery who were operated on for PTC at the State Institution “V.P. Komisarenko Institute of Endocrinology and Metabolism of the National Academy of Medical Sciences of Ukraine” (IEM, Kyiv) from 1990 to 2017 were enrolled. All patients were ≤18 years of age at the time of the Chernobyl accident, lived in the six most ¹³¹I contaminated northern regions of Ukraine, and thus belonged to the high-risk group for the development of radiogenic thyroid cancer.

The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the IEM Bioethics Committee (protocols N 22-KE of 26 April 2018 and N 31-KE of 27 February 2020), the Chernobyl Tissue Bank (CTB, project N001-2020), and the Ethics Committee of Nagasaki University (protocol 20130401–7 of 1 July 2021, the latest update). Informed consent was obtained from all subjects involved in the study or their guardians (for minors).

The primary analysis of histological specimens stained with hematoxylin/eosin, and PTC diagnosis was made in IEM by two experienced pathologists (TB and LZ). The pathological diagnosis was based on the 4th edition of the WHO histological classification (27). Most cases were reviewed by the international pathology panel of the Chernobyl tissue bank (CTB) project (28, 29). The diagnosis of PTC was confirmed in all cases. TNM categories were determined according to the 8th edition of the TNM classification (30).

PTCs were characterized by tumor size, dominant histological growth pattern, oncocytic changes, and major characteristics of tumor invasiveness (11, 31, 32). We also used an integrative clinicopathological variable, the “invasiveness score,” which is the arithmetic sum of every instance of multifocality, lymphatic/vascular invasion, any extrathyroidal extension (i.e., minimal or gross), N1 and M1, either isolated or in combination with other(s) for each tumor (32, 33). Thus defined, the invasiveness score ranged from 0 (no invasive feature presents) to 5 (all features present).

Clinical information was retrieved from the database of IEM. During the follow-up, patients received neck ultrasounds, and those thyroidectomized also received serum thyroglobulin tests on a 6-month to 1-year basis. Recurrence was defined as a tumor focus detected and treated not earlier than 6 months after the primary surgery. Radioiodine-refractory (RAI-R) recurrences were determined according to the guidelines of the American Thyroid Association (34).

Immunohistochemical (IHC) staining for BRAF^V600E (LZ, TIR) was performed as described before (33). In brief, we used a mouse monoclonal anti-BRAF (mutated V600E) antibody (VE1) ab228461(Abcam) at a 1:100 dilution and the Novolink Polymer Detection System (250T) (Leica RE7140-K) to detect IHC reaction product. IHC staining was evaluated by three qualified observers (L.Z., T.B., T.I.R.), and full agreement was achieved; there were no specimens interpreted by any observer as potentially false-negative or false-positive. A close correlation between the results of the VE1-based IHC for BRAF^V600E and molecular methods of the detection of the BRAF^V600E mutation at the DNA level has been reported in a meta-analysis (35) and confirmed in our previous study using formalin-fixed paraffin-embedded material (36). Therefore, we assumed the BRAF^V600E positivity on IHC was indicative of the BRAF^V600E mutation.

The proliferative activity of tumors was evaluated by IHC using Ki67 antibody (clone MIB-1; DAKO, Glostrup, Denmark, 1:100 dilution) in a Ventana BenchMark ULTRA instrument. The Ki67 labeling index (Ki67 LI) was determined with the image-analyzing software (CountσCell, Ki67 antigen Semi-Auto Counter, Seiko Tec LTD, Fukuoka, Japan) in a total of ~1,000 PTC cells (LZ). Image analysis was performed in a blind for the BRAF^V600E status manner.

The individual ¹³¹I thyroid radiation doses (the absorbed dose estimates in mGy) were calculated for each patient in the Dosimetry and Radiation Protection Department of the State Institution “National Research Center for Radiation Medicine of the National Academy of Medical Sciences of Ukraine” (NRCRM, Kyiv) using dosimetric models which include the system of ecological iodine transport and iodine biokinetics (“TD-CTB”) depending on the availability of direct measurements of thyroid activity in May-June 1986 in the person or in other persons in the same or adjacent settlements (37).

The probability of causation of a tumor by radiation exposure in a person of a given sex and age after a certain period of latency was determined using the US NIH/NCI Division of Cancer Epidemiology and Genetics' Interactive RadioEpidemiological Program—Probability of Cancer Causation from Radiation Version 5.7.1 software (38, 39). This software uses “Personal Information” such as gender, birth year, diagnosis year and cancer model (here, the “Thyroid (193)”), and “Dose Exposure Information” such as exposure year (here, 1986), exposure rate (here, the acute), radiation type (here, the electrons E > 15 keV as 90% of ¹³¹I beta-decay has the energy of 606 keV), organ dose (here, Constant) and parameter 1 (the thyroid dose in cSv; since radiation weighting factor for the beta-particles is 1, the equivalent doses were considered to be equal to the absorbed doses) as input variables. The output is the values of the “Assigned Share (Probability of Causation)” that range from the 1st to the 99^th percentile based on 10,000 random-seeded simulations. The higher POC value reflects the higher likelihood of cancer development due to radiation exposure. We used the median value (the 50th percentile) for calculations.

The Fisher's exact test, Fisher-Freeman-Halton exact test, and Cochran-Armitage test were used for univariate analysis of categorical data; the Mann-Whitney test was used to compare continuous data between two independent groups. Logistic regression models were adjusted for sex. Models with very small numbers of outcomes (<5 per cell) were conducted using Firth's approach to bias-reducing penalized maximum likelihood fit. Categorial variables with more than two response levels were assessed using multinomial logistic regression. Multivariate linear regression models were applied to continuous dependent variables. The BRAF^V600E effect in relation to the period of latency was estimated using survival analysis methods. The Kaplan-Meier method, the proportional hazard (Cox), and extended proportional hazard models were used. Computation and plotting of the results of the model with time-varying coefficients were performed with a SAS macro “coxtvc” (40). Multivariate models of the development of RAI-R recurrent PTC metastases in time were developed by non-automatic variable selection in the Cox proportional hazard model using minimization of the Akaike information criterion method. The integrated time-dependent area under curve was calculated for each model.

Additional adjustment for age at operation was performed to evaluate the impact of this variable as a confounder in separate models. Matching 1:1 the BRAF^V600E-positive to the BRAF^V600E-negative PTCs by age at operation (±2 years) was performed with the SAS macro “match” (https://git hub.com/Jiangtang/Programming-SAS/blob/master/UserMacros/ mayo/match.sas) using the “optimal” method. For analyses of matched groups, the univariate related samples tests, i.e., the Wilcoxon signed-rank test for continuous data, McNemar test for 2x2 contingency tables, and Bowker test for categorical variables with several response levels (SAS PROC FREQ with the “agree” option), and clustered log-rank test (41) for the period of latency and recurrence-free survival were applied. Multivariate analysis included linear regression for continuous data and conditional logistic regression analysis for categorical variables with a dichotomous response (SAS PROC LOGISTIC with the STRATA statement). The models included age at operation as an explanatory variable to control for residual confounding.

Calculations were performed using version 9.4 of SAS (SAS Institute, Cary, NC, USA) or IBM SPSS Statistics Version 24 software (International Business Machines Corp., Armonk, NY, USA). All tests were two-sided; p < 0.05 was considered statistically significant.

All data collected for or generated during this study are presented in Figure 1. Cytoplasmic expression of BRAF^V600E (indicative of the BRAF^V600E mutation at the DNA level) was observed in 136/428 (31.8%) of cases. As shown in Table 1, the proportion of male patients was lower in the BRAF^V600E-positive PTCs (OR = 0.609, p = 0.037); therefore, the regression models were adjusted for sex. The distributions of the BRAF^V600E-positive and BRAF^V600E-negative PTCs by age at operation, age at exposure, latency, and ¹³¹I thyroid dose shown in Figure 2 strongly suggest substantial differences between the two PTC groups. Indeed, the BRAF^V600E-positivity was associated with significantly older age at operation (b = 14.346, p < 0.001), age at exposure (b = 4.858, p < 0.001), and longer period of latency (b = 9.681, p < 0.001). The onset of PTCs with different BRAF status in time (i.e., after a certain latency) assessed in a Cox model did not meet the assumption of proportional hazards (Supplementary Figures 1A–C). We, therefore, used an extended Cox model with time-dependent BRAF status, which performed adequately (Supplementary Figure 1D). Parameters of the model confirmed the delay of the BRAF^V600E-positive PTC development in time (HR = 0.007, p < 0.001) and their accelerated failure time behavior (HR = 1.181, p < 0.001 for the BRAF status^*latency variable, Table 1).

The BRAF^V600E-positivity was associated with lower ¹³¹I thyroid dose (b = −560.333, p = 0.006) and POC (b = −27.534, p < 0.001). With increasing patient age at exposure, POC was decreased in all dose ranges, both for the BRAF^V600E-positive and -negative PTCs (Supplementary Figure 2).

The BRAF^V600E-positive PTCs were characterized by a smaller size (b = −7.119, p < 0.001) and higher frequencies of microcarcinomas (OR = 4.788, p < 0.001). Among microcarcinomas, there were no BRAF^V600E-positive and only 2 BRAF^V600E-negative incidentalomas. The higher frequencies were observed for oncocytic changes in tumor cells (OR = 4.820, p < 0.001), multifocal growth (OR = 2.539, p < 0.001), higher Ki67 LI (b = 1.747, p < 0.001) due to more frequent Ki67 LI from 5 to 10% (OR = 2.553, p < 0.001), while the frequency of fully encapsulated PTCs was lower (OR = 0.297, p = 0.014) (Table 1). Of interest, the size of the BRAF^V600E-positive PTCs did not correlate with Ki67 LI (b = −0.008, p = 0.938) while there was a significant inverse correlation in the BRAF^V600E-negative group (b = −0.127, p = 0.027); the difference in these regression coefficients was statistically significant (p_het = < 0.001).

With regard to tumor architecture, the dominant papillary growth pattern was most common for the BRAF^V600E-positive PTCs (OR = 4.715, p < 0.001). Follicular and solid-trabecular dominant components were less frequent (OR = 0.302, p < 0.001; and OR = 0.523, p = 0.003, respectively). Note that the dominant structural component generally coincides with histological PTC subtype/variant or implies conventional PTC if a tumor has a mixed structure (Supplementary Table 1).

Only 17/136 (12.5%) of the BRAF^V600E-positive PTCs belonged to rare histological variants, of which 9 (6.6%) were the tall cell and 8 (5.9%) were the Warthin-like variants (Supplementary Table 1). These tumors were characterized by pronounced oncocytic changes not only in the cells of the primary tumor but also in the cells of lymph node metastases (Figure 3). Among the BRAF^V600E-negative PTCs, rare histological variants were even less common 13/292 (4.4%) of cases (p = 0.004 for frequency comparison with the BRAF^V600E-positive PTCs). The most frequent was the diffuse sclerosing variant, 11/292 (3.8%) of cases; the tall cell and hobnail variants were represented by one case (0.3%) each. The distribution of particular rare histological variants in the BRAF^V600E-positive and BRAF^V600E-negative PTCs was statistically significantly different (p < 0.001).

The BRAF^V600E-positive PTCs displayed less frequent invasive features: lymphatic/vascular invasion (OR = 0.303, p < 0.001), extrathyroidal extension (OR = 0.417, p < 0.001), regional (OR = 0.518, p = 0.003) and distant metastases (OR = 0.091, p = 0.001), and the integrative invasiveness score (OR = 0.485, p < 0.001; Table 1). To evaluate the impact of rare histological variants on tumor invasiveness, we compared PTCs of common variants only (i.e., without rare variants). The effects of BRAF^V600E were very similar to the results presented in Table 1: less frequent lymphatic/vascular invasion (OR = 0.308, p < 0.001), extrathyroidal extension (OR = 0.425, p < 0.001), regional (OR = 0.514, p = 0.004) and distant metastases (OR = 0.117, p = 0.004), and the integrative invasiveness score (OR = 0.495, p < 0.001). Hence, rare histological variants did not markedly contribute to the differences in tumors invasiveness observed between the BRAF^V600E-positive and BRAF^V600E-negative PTCs.

Within the BRAF^V600E-positive group, we did not find signs of higher aggressiveness of PTCs of rare variants as compared to those of PTCs of common variants (the strongest OR = 1.505, p = 0.442; all 95% CIs for ORs included the value of 1). In contrast, in the BRAF^V600E-negative group, PTCs of rare variants had distant metastases more frequently (OR = 4.517, p = 0.014) and higher invasiveness score (OR = 2.867, p = 0.039); associations for other invasive features were statistically non-significant, although all ORs were >1 (data not shown).

Total thyroidectomy was the principal extent of operation in more than 90% of cases in both BRAF^V600E-positive and BRAF^V600E-negative groups (Table 1), although statistically, it was more frequent in the former (OR = 2.690, p = 0.047). Lymph node dissection was performed in about half of the cases without an overall difference between the groups (OR = 0.777, p = 0.230), but the central dissection was more frequent (OR = 2.487, p < 0.001) and the lateral less frequent in the BRAF^V600E-positive group (OR = 0.402, p < 0.001). Most patients (>80%) in each group received postoperative RAI treatment (OR = 1.111, p = 0.704).

The median follow-up of 3.7 years in the BRAF^V600E-positive group was shorter than that of 9.7 years in the BRAF^V600E-negative group (b = −6.129, p < 0.001), likely due to the longer period of latency of the BRAF^V600E-positive tumors. In the course of follow-up, 6/136 (4.4%) of the BRAF^V600E-positive PTCs and 9/290 (3.1%) of the BRAF^V600E-negative PTCs (OR = 1.421, p = 0.516) developed recurrent lymph node metastases not earlier than 6 months after the primary surgery. Lymph node metastases were the only type of recurrences observed in this study, and all of them were reoperated. We confirmed the fully concordant BRAF^V600E status of the primary and recurrent tumors in all 15 recurrent cases (Cohen's κ = 1.000, p = 0.001). Analysis of disease-free survival in the Cox model showed that despite the relatively small absolute difference in the frequency of recurrences, the BRAF^V600E positivity was a risk factor for recurrence (HR = 5.334, p = 0.034 adjusted for sex, tumor size, N and M categories, extrathyroidal extension, multifocality, lymphatic/vascular invasion, extent of thyroid operation, lymph node dissection, and RAI treatment, which are the factors potentially affecting the chance of recurrence; Figure 4). The effect of BRAF^V600E was nearly independent of other variables as judged from the HR = 5.075 value after adjustment for sex only (see Table 1).

Comparative characteristics of the primary tumors and recurrent metastases with different BRAF statuses are presented in more detail in Supplementary Table 2. No difference in tumor size was observed, yet of note, among all recurrent PTCs there was only one BRAF^V600E-negative and zero BRAF^V600E-positive microcarcinomas, and all primary tumors were non-encapsulated. The recurrent BRAF^V600E-positive primary tumors displayed oncocytic changes and multifocality more frequently than the BRAF^V600E-negative PTCs.

It is unlikely that rare PTCs or those associated with higher tumor aggressiveness histological subtypes were overrepresented among the recurrent PTCs. There were only one (of a total of eight, 12.5%) Warthin-like and one (of a total of nine, 11.1%) tall cell variants among the BRAF^V600E-positive PTCs and two (of a total 11, 18.2%) diffuse-sclerosing variants among the BRAF^V600E-negative primary PTCs that recurred.

Of importance, all 6 (100%) recurrent BRAF^V600E-positive metastases were radioiodine-refractory (RAI-R) while there were only 2/9 (22.2%) of such among the BRAF^V600E-negative recurrences (Table 1 and Figure 3). The difference in frequencies of the RAI-R recurrent metastases was statistically significant (OR = 7.113, p = 0.019). In view of clinical significance of RAI-R recurrences, we created several multivariate models (Table 2). Histopathological characteristics associated with RAI-R in these models were greater tumor size, oncocytic changes, multifocality, and the BRAF^V600E-positivity. As seen in Table 1, most of these variables are closely related to BRAF^V600E.

No deaths were documented in the group of 136 patients with the BRAF^V600E-positive PTCs, and 4 deaths were registered among 292 patients (1.4%) with BRAF^V600E-negative PTCs (p = 0.312) from 9 to 25 years after surgery, but none was due to PTC progression.

During exploratory data analysis, we noticed that age at operation was strongly correlated with the BRAF status (which was the focus of this study, Table 1) and a number of baseline or tumor characteristics (which were considered the outcomes; Supplementary Table 3). Correlations of age at operation with the outcomes, however, could not be construed as causative, indicating that this variable may be a confounder. To account for its effect, we used two statistical solutions.

First, we adjusted the multivariate regression models for age at operation. The BRAF^V600E positivity retained significant associations with older age at exposure, a longer period of latency, lower POC, smaller tumor size, higher frequency of oncocytic changes, higher KI67 LI, higher frequency of dominant papillary growth pattern, and lower frequency of the follicular pattern (Supplementary Table 4). No statistically significant associations with aggressive features were seen. Among clinical characteristics, only a higher frequency of total thyroidectomy was observed. The associations were lost for radiation dose to the thyroid, large-size tumors, multifocality, tumor capsule, lymphatic/vascular invasion, extrathyroidal extension, regional and distant metastasis, and the integrative invasiveness score (see also Table 1 footnote “c”). No associations were longer seen with lymph node dissections, recurrence-free survival, and RAI treatment response.

Second, we performed fuzzy (± 2 years) 1:1 BRAF^V600E-positive-BRAF^V600E-negative PTCs matching by age at operation. On multivariate analysis, the association with smaller tumor size (and microcarcinomas), papillary growth pattern, and higher Ki67 LI remained statistically significant (Supplementary Table 5), implying these are inherent properties of the BRAF^V600E-positive radiogenic PTCs. The number of lost associations was even greater, including those with baseline characteristics, thyroid radiation dose and POC, aggressive tumor features, and clinical data.

The associations between BRAF^V600E and different characteristics that were lost after the varying-stringency adjustments for age at operation should be interpreted with caution since a potent confounder, when introduced in a multivariate statistical model, may mask/obscure the effects of causative explanatory variables. We, therefore, believe that the data presented in Table 1 (i.e., without adjustment for age at operation) describe the BRAF^V600E correlations with baseline and clinicopathological characteristics of radiogenic PTCs adequately.

The BRAF^V600E-positive PTCs and BRAF^V600E-negative PTCs displayed statistically significant differences in the probability of causation (POC) of a tumor due to radiation exposure (Table 1 and Figure 5). We therefore further examined whether comparative characteristics of tumors with different BRAF statuses were changing with regard to the POC level. For this purpose, we calculated the BRAF^V600E effect on PTC characteristics for different POC quartiles and then determined whether the linear trend for changes across POC quartiles was statistically significant.