Nthy/BRAF^WT (human thyroid cell lines which harbor wild type BRAF gene [Nthy/WT]) and Nthy/BRAF^V600E (Nthy/V600E) are newly established cell lines expressing either BRAF^WT or BRAF^V600E by lentiviral transduction of human thyroid cell line Nthy-ori 3-1 [28,29]. Nthy/vector was not included because its functional activities such as colony formation and invasion were not observed in the previous study [28]. Also, in vivo tumorigenesis is also not observed in any of Nthy-ori 3-1, Nthy/vector, and Nthy/WT (data not included). Thus, Nthy/WT and Nthy/V600E were used in this study and grown in Roswell Park Memorial Institute (RPMI) 1640 supplemented with 10% fetal bovine serum (FBS, #S1480, Biowest, Nuaillé, France), 2 mM GlutaMAX (#35050061, Gibco, Waltham, MA, USA), and 100 U/mL penicillin-streptomycin (#15140122, Gibco) in a humidified atmosphere of 5% CO₂ at 37°C. All cell lines were maintained for a month and replaced with another identically passaged vial. A Short Tandem Repeat (STR) analysis was performed by Korean Cell Line Bank (Seoul, Korea) using AmpFLSTR identifiler PCR Amplification Kit (#4322288, Applied Biosystems, Waltham, MA, USA) and validated that all the cell lines have the same STR profiles as that of ‘Nthy-ori 3-1’ (RRID: CVCL_2659). The absence of mycoplasma in the cell lines was demonstrated by Korean Cell Line Bank using e-Myco plus Mycoplasma PCR Detection Kit (#25237, iNtRON Biotechnology, Seongnam, Korea).

Before experiments, cells were incubated in phenol red-free RPMI 1640 with 10% FBS for 24 hours followed by of treatment of 100 nM E2 (#E2257, Sigma, St. Louis, MO, USA) for 24 hours. The concentration was determined based on other studies and the results of our invasion assay [10-12,14]. Cell lysates were prepared using RIPA buffer (#RC2002-050-00, Biosesang, Seongnam, Korea) containing protease inhibitor cocktail (#8340, Sigma) and phosphatase inhibitor cocktail 2 and 3 (#p5726 and #p0044, Sigma). The protein concentrations were determined using BCA method and resuspended in sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) loading buffer (#SF2002-110-00, Biosesang) at 1 μg/μL followed by incubation at 95°C for 5 minutes. For 10% SDS-PAGE, 15 μg (for ERα) or 25 μg (for ERβ) of the protein sample was loaded and separated in Tris-Glycine running buffer. Polyvinylidene difluoride membrane was used to transfer the proteins and then blocked in 5% nonfat dry milk at room temperature for 30 minutes. Subsequently, they were incubated in primary antibodies at 4°C overnight: rabbit monoclonal antibody ‘ab108398’ (1:4,000, clone: EPR4097, Abcam, Cambridge, MA, USA) for ERα, rabbit polyclonal antibody ‘ab3576’ (1:4,000, Abcam) for ERβ, rabbit monoclonal antibody ‘#2118’ (1:10,000, clone: 14C10, Cell Signaling Technology, Danvers, MA, USA) for glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as an endogenous control. Detection was performed using Luminata Forte Western HRP Substrate (#ELLUF0100, Millipore, Burlington, MA, USA) following 2 hours of incubation in the secondary antibody at room temperature: anti-rabbit HRP-conjugated antibody ‘#7074’ (1:4,000, clone: 14C10, Cell Signaling Technology). Densitometry analysis was performed using Image Lab (Bio-Rad Laboratories, Hercules, CA, USA). ERα and ERβ were normalized by GAPDH first, and fold change by E2 treatment was calculated in each experiment. The ratio of ERα and ERβ was calculated by using the GAPDH-normalized ERα and ERβ values, and the relative fold change in Fig. 1B is visualized based on the ratio of the normalized expression in E2-treated group over control group to show the effect of E2 on expression of ERs depending on BRAF status. Statistical tests were performed between Nthy/WT and Nthy/V600E with four biological replicates.

BD cell culture inserts with 8-μm pore for 24-well plates (#353097, BD Bioscience, Bedford, MA, USA) were used and coated with 40 μL of 0.5 mg/mL growth factor-reduced Matrigel (#356230, BD Bioscience) for invasion assay. Matrigel was equilibrated for 12 hours at 37°C in 5% CO₂ before cell loading. Cells were incubated in phenol red-free complete RPMI medium supplemented with 2% dextran charcoal-stripped FBS for 48 hours before the experiments and harvested using TrypLE Express without phenol red (#12604-021, Gibco). The harvested cells were resuspended in RPMI without FBS at 1E6 cells/mL, and 100 μL of the cell suspension was loaded onto the culture inserts. The bottom chambers were filled with RPMI with 2% FBS. Either dimethyl sulfoxide (DMSO) or E2 was added to both the upper and bottom chambers and incubated for 24 hours at 37°C in 5% CO₂. After the incubation, the cells on the upper side of the culture inserts were removed using cotton swabs and the media in the bottom chambers were carefully aspirated. To measure the amount of migrated or invaded cells under the culture inserts and on the bottom wells quantitatively, CellTiter-Glo Luminescent Cell Viability Assay (#G7572, Promega, Madison, WI, USA) was used based on the manufacturer’s instruction. The measured values (relative light unit) were displayed with fold change by dividing by negative controls (0 nM of E2 in each cell line), and statistical significance was tested after log2 transformation.

Soft-agar assay was performed using CytoSelect 96-Well Cell Transformation Assay kit (#CBA-130, Cell Biolabs, San Diego, CA, USA). Nthy/WT and Nthy/V600E cells were seeded at about 7,000 cells per well in 24-well plates in a top layer of 0.4% agarose on a base layer of 1.2% agarose. Culture medium containing DMSO or E2 was added to each well and cultured for 2 weeks changing the media every 2 to 3 days. After about 2 weeks, the number of cells was measured quantitatively by using the fluorescent dye included in the assay kit following the manufacturer’s instruction. The detected relative fluorescence units using SpectraMax M5 (Molecular Devices, San Jose, CA, USA) were normalized by the no-cell controls, and fold changes of E2-treated groups over the control group in each experiment were visualized.

Nthy/WT and Nthy/V600E were seeded at 5E4 cells per each well in 24-well culture plate and incubated until all wells reach 100% confluence. Each well was wounded by scratching with a SPLScar scratcher (#201925, SPL Life Science, Pocheon, Korea) and washed with Dulbecco’s phosphate-buffered saline (DPBS) with Ca²⁺ and Mg²⁺ to remove the debris followed by filling up 400 μL of complete media with 10% FBS and phenol red. After confirming that wounds are evenly formed without a chunk of cells, we added 100 μL of media containing either DMSO, 50 nM of diarylpropionitrile (DPN; #1494, ERβ-selective agonist, Tocris Bioscience, Bristol, UK), or 5 μM of methylpiperidino pyrazole (MPP; #13863, ERα-selective antagonist, Cayman Chemical, Ann Arbor, MI, USA). The concentrations were determined based on the previous studies using these molecules [30-32]. We marked the clearly scratched area under the plates before incubation and took pictures after 0 and 6 hours of incubation. In the two matched images, we made two representative borders for before and after incubation with black and red colors, respectively. The difference of widths between the borders was measured for each condition to compare the relative migrated area.

To analyze expression pattern of ESR1 and ESR2 in thyroid cancer patients depending on BRAF^V600E status and age group, we downloaded clinicopathologic, somatic mutation, and mRNA expression count data of 365 female PTC patients that both clinical and RNA-sequencing data are available from The Cancer Genome Atlas (TCGA) data portal. In the dataset, mRNA sequencing experiments were performed by University of North Carolina using a HiSeq 2000 (Illumina, San Diego, CA, USA). Gene expression counts were expressed as RNA-sequencing by expectation maximization (RSEM) with upper quartile normalization, while RSEM results pertaining to the ESR1 and ESR2 genes were extracted from RNA-sequencing data. The ‘ESR ratio’ was defined as ESR1 to ESR2 (ESR1/ESR2) RSEM count. Data pertaining to somatic mutations of BRAF genes were obtained from the mutation calling file provided by Broad Institute and Baylor College of Medicine using an Illumina Genome Analyzer. We have previously reviewed pathologic documents downloaded from TCGA to correct some inaccurate data. Expression of ESR1, ESR2, and ESR ratio were analyzed by dividing into four groups depending on the two factors: BRAF^V600E status (excluded two patients with different types of BRAF mutation) and age 50 at diagnosis. Gene signatures such as thyroid differentiation score (TDS), BRAF^V600E-Ras score (BRS), and M2 macrophage signature (M2) were calculated using single-sample gene set enrichment analysis (ssGSEA) based on the markers reported by the TCGA group and expressed in log2- transformed z-score [33].

Statistical analyses were performed using GraphPad PRISM version 8.0.0 (GraphPad Software, San Diego, CA, USA) for experimental data and R version 4.0.3 (R Foundation for Statistical Computing, Vienna, Austria) for TCGA data. Data are presented as mean±standard deviation. The independent t test or Mann-Whitney test was used for comparison of continuous variables. Statistical significance was defined as two-sided P<0.05.

At the basal level in the described culture condition, Nthy/V600E expressed lower ERα levels and higher ERβ levels than Nthy/WT (Fig. 1A). However, following E2 treatment, Nthy/WT cells decreased the ERα/ERβ ratio by downregulating and upregulating the expression of ERα and ERβ, respectively (Fig. 1). In contrast, Nthy/V600E manifested an increased ERα/ERβ ratio through the upregulation of ERα and a slight decrease in ERβ expression in response to E2 (Fig. 1). As these cell lines have an identical genomic background except for BRAF^V600E expression, BRAF^V600E is likely to be involved in the downstream estrogen signaling pathways in this human thyroid cancer cell line.

The E2 treatment increases the migration and invasion ability of thyroid cancer cell lines, which are especially important in the initial process of cancer metastasis. There was little change in the migratory ability of Nthy/WT based on the treatment concentration of E2, whereas Nthy/V600E was more sensitive to the E2 concentration (Fig. 2A). Moreover, invasion was significantly increased in Nthy/V600E, but not in Nthy/WT (Fig. 2A). Despite no significant difference in migration, it was evident that Nthy/V600E was more susceptible to E2-induced increase in migration and invasion than was Nthy/WT. It is consistent with the expression changes of ERs in response to E2 treatment: increased and decreased ERα/ERβ ratio in Nthy/V600E and Nthy/WT, respectively.

Anchorage-independent growth is considered to be a key phenotype of cancers [34]. In particular, anchorage-independent growth plays a critical role in the late stage of metastasis because this ability is closely associated with anoikis resistance and stemness of cancers [35]. E2 increases stem-like features and anchorage-independent growth of breast cancer cell lines and primary thyrocytes from nodular goiters [14,36]; therefore, we evaluated the role of E2 in the anchorage-independent growth of thyroid cell lines depending on their BRAF^V600E status.

Nthy/WT formed only a few colonies and grew rarely, and E2 treatment slightly inhibited anchorage-independent growth (Fig. 2B). In contrast, E2 enhanced the growth of Nthy/V600E in soft-agar, forming more and larger colonies in a dose-dependent manner (Fig. 2B). These results corresponded to the changes in ERα and ERβ expression following E2 treatment, which showed an increased and decreased ERα/ERβ ratio in Nthy/V600E and in Nthy/WT, respectively. This suggests that BRAF^V600E contributed to the enhancement of the tumor-promoting effects of E2 by modulating the expression of ERα and ERβ in response to E2 treatment.

In order to confirm that the differential responsiveness BRAF^V600E to E2 is mediated by the altered expression of ERα and ERβ, we tested the efficacy of ERα-selective antagonist “MPP” and ERβ-selective agonist “DPN” in reducing the migration ability via a wound healing assay. We cultured Nthy/WT and Nthy/V600E in complete media with phenol red and 10% FBS to ensure the proper functioning of both ERα antagonist and ERβ agonist. The representative pictures were selected based on the mean values in each condition.

As Nthy/WT showed relatively low ERα and high ERβ expression compared to Nthy/V600E in the presence of E2, DPN significantly inhibited the migration of Nthy/WT whereas MPP did not. Thus, ERβ agonist was more effective in Nthy/WT (Fig. 3). Unlike Nthy/WT, DPN did not significantly inhibit migration, whereas MPP significantly reduced migration in Nthy/V600E, which expresses higher ERα and lower ERβ concentrations in the presence of E2 (Fig. 3). As both cell lines were tested under identical conditions, it was clear that BRAF^V600E rendered the ERα signaling pathway more dominant over ERβ in Nthy under E2-treated conditions, with a resultant E2-induced enhancement of cancer-promoting responses.

In the TCGA data, the overall ESR1 expression is higher in BRAF^V600E group than in BRAF^WT group, and women who are younger than 51 years expresses higher ESR1 levels than those older than 50 years (Fig. 4A). Both BRAF^V600E and age less than 51 years are associated with a higher ESR1 expression, although the effect of age is stronger in BRAF^V600E group than in BRAF^WT group (744.5 vs. 426.0, P=0.0129 in BRAF^V600E; and 502.3 vs. 389.3, P=0.5057 in BRAF^WT) (Fig. 4A). The overall ESR2 expression is lower in BRAF^V600E group compared to BRAF^WT group, but significant differences by age were not observed in both groups (8.902 vs. 8.049, P=0.2054 in BRAF^V600E group; and 11.46 vs. 12.90, P=0.333 in BRAF^WT group) (Fig. 4A). Nevertheless, the ESR ratio of the young group is significantly higher in BRAF^V600E group, but not in BRAF^WT group (127.6 vs. 68.58, P=0.0257 in BRAF^V600E group; and 58.21 vs. 41.39, P=0.3214 in BRAF^WT group) (Fig. 4A). In correlation analysis between ESR genes and BRS, BRAF^WT group does not show strong correlations regardless of age groups although BRAF^V600E-like signature is positively correlated with higher ESR1 and lower ESR2. However, a different correlation pattern is observed in BRAF^V600E group. In the younger age group with BRAF^V600E, BRAF^V600E-like signature is positively correlated more strongly than in the BRAF^WT groups (Fig. 4B). On the other hand, the older age group with BRAF^V600E had the opposite pattern showing BRAF^V600E-like signature is associated with lower ESR1 and higher ESR2 (Fig. 4B).

As BRAF^V600E-like gene signature has been reported to be associated with poor differentiation (low TDS) and activated MAPK pathway signatures by TCGA group, high ESR1 and low ESR2 also showed correlation with the prognosis-related gene signatures including EMT-related genes and M2 macrophage signature (Fig. 4C, Supplemental Fig. S1). These results suggest that BRAF^V600E might be associated with the regulation of ESR1 and ESR2 expression, which would be contributed by the higher level of estrogen in the group of younger women. Taken together, BRAF^V600E might enhance the effect of estrogen on thyroid cancer cells by modifying the responsiveness.