Thyroid tumors and normal thyroid tissues from flash-frozen samples or optimal cutting temperature (OCT) compound-embedded samples were obtained from different institutions. This study was performed in accordance with the Declaration of Helsinki and collection of patient tissues and associated data was done in agreement with the Ethics Committees of Jules Bordet Institute (CE1978, CE2970) with informed consent of patients, Sapienza University of Rome, Pitié-Salpêtrière Hospital and with the Mayo Clinic Institutional Review Board protocol. Samples and associated data from the Institute – Instituto Português de Oncologia de Lisboa Francisco Gentil (IPOLFG) – were obtained in compliance with all applicable laws, including a written consent, and the study was approved by the institute’s Ethics Committee. Written informed consent was obtained from all patients of the tumor tissue bank – Tumorothèque Centre de Ressources Biologiques des Hospices Civil de Lyon. Samples and associated data were obtained from the tumor tissue bank – Tumorothèque ALLIANCE CANCER de Lille – that operates under the authorization AC-2018-3110 granted by the French ministry of research. Prior to scientific use, patients were appropriately informed and asked to consent in compliance with the French regulations.

All ATC and PDTC diagnoses were histologically re-evaluated and confirmed by an expert thyroid pathologist (MD), by examination of corresponding FFPE slides or, whenever possible, of serial sections of OCT slides. PDTC were defined following the Turin proposal (39).

A total of 140 samples were collected for two-dimensional gel (2D-gel) electrophoresis analysis (42 ATC, 30 PDTC, 8 follicular thyroid carcinoma (FTC), 9 oncocytic FTC, 23 papillary thyroid carcinoma (PTC) (1 metastasis), 7 lymph node metastases (LNM) and 21 normal thyroid tissues). Twenty-one samples were discarded due to signs of wide necrotic tissue (more than 80% of tissue area; 2 ATC), misclassification (1 ATC and 1 oncocytic FTC), bad RNA quality (RNA integrity number lower than 3)/protein quality (3 PTC, 8 ATC, 1 normal and 3 oncocytic FTC) or insufficient protein quantity (1 PTC and 1 PDTC). Eighty-four of 119 samples were selected for RNA-seq, and following an exploratory analysis, 27 cases were further excluded from analysis due to: gross contamination/misclassification (3 normal with neoplastic cells, 1 PDTC with no neoplastic signs, 2 ATC with more than 70% of PTC cells, 1 ATC with 60% of stroma area); sample duplicate (1 normal); low RNA quality (RNA integrity number lower than 4.5; 1 ATC and 2 PDTC); high proportion of PCR duplicates (more than 0.5% of the library; 1 FTC and 2 oncocytic FTC); or very low/undetectable 2D-gel signal (4 ATC, 6 PDTC and 2 PTC). In summary, 98 samples were included in the 2D-gel proteomic analysis and 57 samples were included in the transcriptomic analysis.

The antibodies used in this work are listed in Supplementary Table S2 . Equal amounts of whole‐cell extract proteins or immunoprecipitates were separated by SDS-PAGE and immunodetected. For 2D‐gel electrophoresis, cells were lysed in a buffer containing 7 M urea and 2 M thiourea. Ground tissues obtained from cryogrinding of flash-frozen tumor tissues and frozen tumor slides embedded in OCT (at least, 10 sections of 10 μm per sample) were solubilized as described (26). Proteins were separated by isoelectric focusing on immobilized linear pH gradient strips (pH 5 to 8, Bio‐Rad, Hercules, CA, USA) before separation by SDS-PAGE. Chemiluminescence images of the samples were acquired on films or with a Vilber‐Lourmat Solo7S camera and quantified using the Bio1D software (Vilber-Lourmat, Marne-la-Vallée, France). The profile of CDK4 separated by 2D-gel electrophoresis has been characterized previously (13, 26). Quantification of the spot 2 and spot 3 volumes (corresponding to the two main modified forms of CDK4) were done from 16‐bit scans of the 2D‐gel immunoblots, with unsaturated signals. After linear correction of the background, the volume ratio (spot3/spot2) was used to define three types of CDK4 modification profiles: a profile A (for absent) was attributed when the ratio was below 0.02; a profile L (for low) was attributed when the ratio lied between 0.02 and 0.50; and a profile H (for high) was given when the ratio was equal or above 0.50.

Co-immunoprecipitations were performed as described (13, 21). RB-kinase activity of immunoprecipitated CDK complexes was measured by in vitro incubation with ATP and a fragment of RB, as described (21, 36).

Hematoxylin/eosin (HE), KI67 and p16 immunohistochemistry (IHC) was performed with FFPE tissue sections on a fully automated BenchMark Ultra IHC system (Ventana, Roche Diagnostics, Basel, Switzerland) with the UltraView Universal DAB Detection Kit (Roche Diagnostics), using standard routine protocols. Otherwise, FFPE sections from Mayo Clinic center were processed as described (49) and probed with a different set of antibodies ( Supplementary Table S2 ). Pictures were acquired with a NanoZoomer digital scanner (Hamamatsu Photonics, Hamamatsu, Japan) at 40× magnification or in an Aperio AT2 scanner (Leica Biosystems, Wetzlar, Germany) at 40x magnification (samples from IPOLFG institute).

The QuPath software (50) (version 0.4.0) was used for immunohistochemistry scoring analysis. For each section, the most appropriate tumor areas (without necrosis or fibrosis) were delimited and the positive cell detection command was applied using the sum of optical density, and the nuclear minimum area for detection was set between 10-25 µm². A 0.2 intensity threshold for the nuclear mean DAB optical density was applied for KI67 positive detection. Three intensity thresholds were set for the cellular mean DAB optical density of p16 staining, categorizing cells into weak (0.2 ≤ threshold < 0.4), moderate (0.4 ≤ threshold < 0.6) and strong (threshold ≥ 0.6) positive detection. The H-score was calculated by the software, which is a weighted sum (1 for weak, 2 for moderate and 3 for strong) of the percentage of cells in each intensity class, ranging from 0 to 300.

Response curves and statistical analyses were performed using GraphPad Prism version 6.0 (GraphPad Software, La Jolla, CA, USA). Half-maximal inhibitory concentrations of the cell proliferation (GI50) and 95% confidence intervals were estimated by fitting the data to the logarithm of concentration vs. normalized response model, with standard slope (Hill slope of -1.0) and using a least-squares fitting method. The two-sided unpaired Student’s t-test was used for comparison between two groups. Multiple group comparisons were done using the Kruskal–Wallis test followed by Dunn’s multiple comparison tests. Correlations were evaluated with the Pearson correlation (two-sided). Comparison of Kaplan–Meier survival curves was performed using the Mantel–Cox log-rank test. The Fisher’s exact test (two-sided) was used for the analysis of contingency tables. A p-value ≤ 0.05 was considered statistically significant.

To evaluate the proportion of patient thyroid tumors that would be potentially responsive or intrinsically resistant to CDK4/6 inhibition, we determined the CDK4 modification profile and the presence of activated T172-phosphorylated CDK4 from a cohort of fresh-frozen samples (n=98) representative of the different subtypes of thyroid tumors ( Figure 1A , Supplementary Figure S1 and Table S3 ). CDK4 was detected by immunoblotting from whole protein extracts separated by 2D-gel electrophoresis. CDK4 was thus resolved by its charge into three main forms as first characterized in thyroid primary cultures (13, 21, 51) and in other tumor types (26, 27). The most basic form (spot 1) is the native CDK4. The most acidic form (spot 3) increases in response to in vitro proliferation stimulus of various cells including thyrocytes (13, 21, 24, 35, 51–53) and has been identified as the highly regulated T172-phosphorylated CDK4 form, using several approaches including T172-phosphospecific antibodies (11, 22, 26, 54). To evaluate the relative level of CDK4 activation, we compared the presence of this T172-phosphorylated form to another modified CDK4 form (spot 2). This latter form does not incorporate [32P] phosphate and binds to p16 but only weakly to cyclins D (13). The abundance ratio of the T172-phosphorylated form (spot 3) over the form separated in spot 2 allowed us to define three types of profiles (L for low, H for high and A for absent phosphorylated CDK4), as previously reported (26, 27). All WDTC exhibited the T172-phosphorylated form of CDK4 (n=29). This form was also present in all LNM (n=7) and in one available PTC distant metastasis, with higher proportion of profiles H in these groups of samples (50%), than in primary WDTC (25%). In advanced dedifferentiated tumors, phosphorylated CDK4 was present in 95% of PDTC (n=20) but only in 70% of ATC (n=23). On the other hand, as expected with quiescent tissue, normal thyroid samples (n=18) predominantly lacked CDK4 phosphorylation.

In order to evaluate the molecular features associated with the different CDK4 profiles, we analyzed a sub-cohort of 84 cases by RNA-seq (complete gene expression data provided in Data Set 1 ). An exploratory classification of the gene expression profiles by PCA revealed misclassified outliers ( Supplementary Figure S2A ). After excluding the samples with major technical issues and/or clear misclassification (as detailed in section 2.1), further analysis was done for 57 selected cases ( Supplementary Table S4 ). As expected (37, 55–57), we observed that samples were mainly separated according to their tumor sub-types ( Figure 1B ). No clear segregation of samples according to the CDK4 modification profile was observed; nevertheless, most PTC samples with CDK4 profile H were closer to ATC samples ( Figure 1B ). Evaluation of the proliferation with the MKI67 gene expression ( Figure 1C ), of the thyroid differentiation (TDS, Supplementary Figure S2B ) and of the epithelial-mesenchymal transition scores (EMT, Supplementary Figure S2C ) validated the identity of ATC cases. As expected, these were the most proliferative samples and had lost thyroid identity with acquisition of mesenchymal features. PDTC and PTC samples although being very heterogeneous, showed higher MKI67 expressions and lower differentiation than normal tissues. Paradoxically, PDTC tended to have a higher differentiation score than PTC, as reported by others (37). Also, worth mentioning is that the only PDTC case with CDK4 profile A (sample JPI25) had the highest MKI67 expression ( Figure 1C ) and the second lowest thyroid differentiation score of PDTC group ( Supplementary Figure S2B ). Still, its EMT score ( Supplementary Figure S2C ) was similar to other PDTC [in contrast to sample JPI73 that is probably a misclassified ATC ( Supplementary Figure S2A )].

Search for mutations (synonymous mutations excluded) in the main genes involved in thyroid cancer and in cell cycle, revealed TP53 mutations as the major oncogenic driver for ATC, being present in 65% of cases ( Supplementary Table S5 ). For PDTC, NRAS and TP53 mutations were present in 21% and 37% of the tumors, respectively. As previously observed (32), these mutations were mutually exclusive and together affected more than half of the PDTC. PTEN mutations were also present in 21% of PDTC. Interestingly, these appeared to co-occur with TP53 mutations (in 3 of 4 PTEN mutated cases). Mutations in other genes were involved in less than 20% of the ATC and PDTC. BRAFV600E was the only driver mutation identified in PTC, lymph node or PTC metastases, being present in 88% of samples. Among the five tumor samples with no phosphorylated CDK4 (profile A), we identified two samples with mutated RB1 (a nonsense and a frameshift mutation). Three of these profile A samples had TP53 mutations.

As in our previous studies in breast tumors (26) and pleural mesotheliomas (27), we observed the lack of CDK4 T172-phosphorylation in some thyroid tumors to be associated with high expressions of CDKN2A (p16), CCNE1 (cyclin E1) and E2F1. Transcript levels of CDKN2A, E2F1 and CCNE1 were indeed elevated in ATC and PDTC samples with profile A of CDK4 ( Figures 1D-F ). High p16 cellular levels are expected to prevent the activating phosphorylation of CDK4 by impairing its binding to cyclins D. By contrast, in normal quiescent thyroid tissue, lack of CDK4 phosphorylation was likely due to the absence of mitogenic stimulation (13, 21) because it was associated to a particularly low p16/CDKN2A expression, an intriguing feature already reported by others (35, 58–61). Overexpression of p16/CDKN2A in cancer cells has been generally related to RB1 loss of function (28, 62, 63). However, three of these profile A samples (with high CDKN2A expression) were ATC without any detectable RB1 mutation (samples JPI21, JPI27 and JPI74). Lack of RB1 mutation was confirmed with targeted DNA-sequencing ( Supplementary Figure S3 ) in two samples (JPI21 and JPI27, although the latter showed RB1 hemizygous loss). Two other cases of profile A ATC, which could not be RNA-sequenced, also revealed an absence of RB1 mutations or loss. The copy number profile of the profile A PDTC (JPI25) was very complex ( Supplementary Figure S3 ) and we could not evaluate if the very low RB1 variant allelic fraction of 7% (compared to the TP53 variant allelic fraction of 25%) suffice to explain the high expression of CDKN2A and lack of CDK4 phosphorylation. Indeed, the RB alteration might be a subclonal event in this sample. Nevertheless, samples with profile A had a reduced RB1 expression, comparable to the one observed in normal tissue samples ( Figure 1G ).

Due to an alternative open reading frame (ORF), CDKN2A translates into two different proteins (p16 and p14). In contrast to breast tumors (26), high CDKN2A expression measured by RNA-seq analysis in thyroid tumors, did not entirely correlate with elevated p16 expression and lack of CDK4 phosphorylation. Among PDTC and ATC, eleven samples had elevated CDKN2A mRNA expression levels, but displayed CDK4 profiles H or L. In most of these cases, over-expressed CDKN2A mRNA seemed to mainly encode the p14 protein (only reads covering exon 1β, specific for p14 isoform, were observed in IGV) or would encode a truncated p16 protein (p.Q50* mutation). If we account for this by correcting CDKN2A transcription (as detailed in Materials and Methods), elevated expression of p16/CDKN2A was exclusively found in samples lacking phosphorylated CDK4 ( Figure 1H ). Immunohistochemical staining corroborated the absence of p16 protein in samples with p14-driven CDKN2A transcription ( Figure 1I ). By contrast, profile A samples with similar CDKN2A expression levels had strong p16 staining. Importantly, in all such cases, KI67 staining was similarly elevated in the p16-high tumor area, further demonstrating that p16 accumulation did not preclude tumor cell proliferation in these tumors.

To validate the high levels of p16 protein as a useful biomarker to identify tumors with no phosphorylated CDK4 (profile A), we performed immunohistochemistry staining in FFPE samples from 13 ATC, 17 PDTC, 12 PTC including 2 LNM, and 3 normal thyroid tissues ( Supplementary Table S6 and Figure S4 ). No p16 staining was observed in the normal thyroid samples. All six profile A tumors (5 ATC and 1 PDTC) were p16-positive, 3 of which with very strong staining. However, 9 out of 19 profile L and 7 out of 15 profile H, including most PTC, were also p16-positive in a significant proportion (more than 15%) of cells. The remaining cases (mostly profile L and H ATC and PDTC) had no or very faint staining, and positivity appeared rather in the stroma. Staining of serial cuts for the proliferation marker KI67 revealed that profile A tumors had a higher proportion of KI67-stained cells and, more importantly, that p16-positive areas were co-stained for KI67. By contrast, no such co-expression of p16 and KI67 was observed in most profile H and L tumors, including PTC samples in which the p16-positivity was also particularly heterogeneous ( Supplementary Table S6 and Figure S4 ). Only one profile L ATC case (JPI84, L* in Supplementary Table S6 ) exhibited areas with co-expression of both proteins, indicating that it could be a heterogeneous sample with profile A areas, which could also explain its very low ratio of CDK4 spot 3/spot 2.

We previously reported that the CDK4 modification profile of breast tumors and cell lines can be predicted using the expression values of 11 genes (26). The tool was initially developed with Affymetrix microarrays. Our ongoing work (manuscript in preparation) indicates that the prediction tool can be run with gene expression profiles quantified with other technologies (RNA-seq, AmpliSeq) and adapted to RNA extracted from FFPE tumor samples. Therefore, we first tested whether the tool based on the genes ( Supplementary Table S7 ) selected for the CDK4 profile prediction in breast tumors and the corresponding RNA-seq breast references could correctly predict the profiles of thyroid tumors, using raw CP20M gene expression values. The accuracy (proportion of predictions matching the observations) reached only 57.4% when applied to the prediction of the three CDK4 modification profile in thyroid tumors ( Supplementary Table S8 ). Five samples were erroneously predicted as A profiles. Four of them were samples with exclusive contribution of exon1β (coding for p14) to the expression of the locus. In addition, profiles H or L were confounded in 19 tumor samples (7 L profiles predicted as H and 12 H profiles predicted as L). Upon use of the CDKN2A CP20M gene expression values corrected by the contribution of the exon1α to the expression of the locus, the overall accuracy raised to 68.5% for the prediction of the three profiles. As risk evaluation by distinguishing H from L profiles may not be practically useful in advanced thyroid cancer, we next checked whether a binary prediction of tumors as A or non-A (L or H) profiles might be more accurate. When a tumor was predicted to have A profile by a highest correlation with the A profile reference compared to the H or L profile references, the accuracy of the prediction, based on the exon1α CDKN2A contribution, raised to 98.2% ( Supplementary Table S8 ). All observed A profile tumors were correctly predicted and only one sample (the PTC lymph node metastasis JPI38) was falsely predicted as A profile.

ATC with CDK4 profile H were significantly more proliferative than profile L cases (p = 0.0063), as judged by the higher expression of MKI67 ( Figure 2A , Supplementary Figure S5A ). A similar tendency could be seen for PDTC (p = 0.0854). However, we noticed that only a subset of profile H PDTC samples had higher MKI67 expression levels than profile L samples ( Figure 2A , Supplementary Figure S5B ). By contrast in PTC, several tumors displayed a H profile, but none was observed to be more proliferative than L profile ones ( Figure 2A , Supplementary Figure S5C ). Expression of MKI67 was indeed significantly correlated with phosphorylated CDK4 levels in ATC (p = 0.0462) and PDTC (p = 0.0490) but not in PTC (p = 0.4329) ( Supplementary Figures S5A-C ). Nevertheless, PTC with profile H had significantly lower thyroid differentiation than profile L (p = 0.0266), a tendency that was not clearly associated to BRAF mutations ( Figure 2B ). In each of the ATC and PDTC groups, samples with profile H showed significantly higher expression of E2F1 than samples with profile L (p = 0.0288 and p = 0.001, respectively; Figure 2C ). CCNE1 and RB1 expressions did not differ significantly between the two profiles ( Supplementary Figures S5D, E ).

PDTC patients, for which tumors exhibited a CDK4 profile H, had significantly shorter overall survival than those with profile L tumors (p = 0.0407; Figure 2D ). In ATC, the median survival tended to be lower for profile H (2 months vs 7.5 months for profile L), but more cases are needed to prove a true association (p = 0.1590; Figure 2E ). Similarly, PDTC and WDTC patients associated to profile H were more likely to have distant metastases (89% vs 50%; Figure 2F ) and nodal metastases (83% vs 33%; Figure 2G ), respectively. Therefore, even if it was not associated with a higher proliferation in WDTC, a higher CDK4 activation appeared associated to features of higher aggressiveness.

Unfortunately, complete clinical data could not be retrieved for several profile A tumor patients and several ATC patients were lost to follow-up. Remarkably, one of these profile A-ATC patients (JPI21) is still alive more than 8 years after diagnosis and chemotherapy.

Both commercial (n=17) and patient-derived cell lines (64) (n=4) from different subtypes of tumors were analyzed in order to characterize their response to two CDK4/6i - palbociclib and abemaciclib. Only three cell lines were intrinsically resistant to the drugs ( Figure 3A ). For the remaining cells, including all ATC-derived lines, both inhibitors reduced S-phase entry (as measured by BrdU incorporation into DNA), in a dose-dependent way. Palbociclib’s half-maximal inhibitory concentrations of the cell proliferation (GI50) were similar to the ones reported in other types of cancer cells and had no apparent correlation with the cell line subtype or driver mutations ( Figure 3B ). Nevertheless, six cell lines had a significant residual proliferation (defined by more than 20% of cells in S-phase, relatively to the control) at the maximum dose of 1µM of palbociclib ( Figure 3A ). Assessment of cell viability by other methods, dependent on mitochondrial dehydrogenases activity (MTT) or on the cellular protein content (SRB), showed less pronounced effect ( Supplementary Figure S6 ), stressing the mainly cytostatic properties of CDK4/6i. CDK4 immunoblotting from 2D-gel electrophoresis revealed that T172 phosphorylation of CDK4 was present in all the sensitive cell lines ( Figures 3C, D ). By contrast, the three completely resistant cells lacked the (active) phosphorylated CDK4 form ( Figure 3C ). Further analysis by western blotting showed that these resistant cell lines had no RB or no phosphorylated RB expression, in accordance with the mutations previously reported (65) (a frameshift mutation in T243, a nonsense mutation in FTC-238 and two amino acids deletion in FTC-133) and confirmed by our RNA-seq data ( Figure 3D ; Supplementary Table S2 ). The absence of functional RB in these cells probably accounted for the elevated expressions of p16 and Cyclin E1 observed in Figure 3D . In almost all sensitive cell lines, p16 was affected either by mutation, deletion or promoter methylation of the p16-coding transcript and thus was not expressed ( Figures 3D, E ). Nevertheless, CDKN2A transcripts were detected in most sensitive cell lines without reported CDKN2A gene deletion. Analogously to what we observed in profile L and H tumors with high CDKN2A transcript levels, CDKN2A expression was exclusively due to the alternative ORF p14 expression in these cell lines ( Figure 3E ). Exon 1β, which is unique to p14 transcript, was expressed in more than half of the cell lines and, consistently, p14 was detected by western blotting ( Figure 3E ). On the other hand, CDKN2A exon 1α, which is exclusively present in p16 transcript, was only expressed in 5 cell lines. In two of them (the sensitive HTC-C3 and THJ-29T cells), CDKN2A is mutated (H66R+D108N and G55Afs*91, respectively) and p16 protein was not detected. Cyclin E1 and CDK4 expressions were elevated in all the resistant cell lines and were more variable in sensitive ones ( Figure 3D ). High residual proliferation was seen in cell lines with higher expression of Cyclin E1 (WRO, CAL-62 and IHH-4) or CDK4 (THJ-11T). At the transcript level, CDK4 expression showed a significant negative correlation with the sensitivity of the cell lines and the expression of both CDK4 and CCNE1 were significantly correlated with the residual proliferation of the cell lines ( Supplementary Figure S7 ). The basal proliferation of the cell lines was also significantly correlated with the residual proliferation. By contrast, the transcript levels of CDK6, RB1 or CCND1 (cyclin D1) were not correlated with the residual proliferation ( Supplementary Figure S7 ).

We next explored the impact of palbociclib treatment for 4 d (tested at 1 µM), on selected cell cycle related protein expressions ( Figure 4A ). At this drug concentration, no alteration was seen in the resistant cells upon palbociclib treatment. This suggests that any effect observed in responsive cells was specific. In sensitive cells, we observed diminished RB expression/hyper-phosphorylation and decreased levels of Cyclin A2 and E2F1 ( Figure 4A ). Palbociclib-treated cells exhibited elevated levels of CDK4, Cyclins D and E1. These observations are consistent with a cell cycle arrest at G1 phase. Drug-exposed cells observed at the microscope were flat, enlarged and seemed to have increased cellular complexity (not shown), suggestive of cells undergoing senescence. Palbociclib treatment induced higher senescence-associated β-galactosidase (SA-β-gal) activity in all the sensitive cell lines ( Figure 4B ). However, this increase was detected in more than 30% of cells only in three cell lines (THJ-29T, HTC-C3 and BHT-101). Aside from decreased levels of PARP1 [E2F target already described as down regulated by CDK4/6 inhibition (16)], treated cells showed no appreciable evidence for increased apoptosis (no cleaved PARP1 nor cleaved caspase 3) or autophagy (minimally changed p62 and LC3B levels) ( Supplementary Figure S8 ). This is in agreement with previous studies, where no or minimal induction of apoptosis was seen in thyroid cancer cell lines (66–68). Permanent arrest in the three cell lines with elevated SA-β-gal activity was further evidenced by reduced S-phase entry during prolonged palbociclib treatment (up to 10 d), which was maintained even after drug washout ( Figure 4C ). However, in most of the other sensitive cell lines, palbociclib did not induce a durable cell cycle arrest. Indeed, the proportion of cells entering into S-phase raised when the treatment was prolonged, or rapidly increased once the compound was withdrawn ( Figure 4C ).

The 11-gene expression signature used to predict the CDK4 profile in tumors was further evaluated in these 21 thyroid cell lines ( Supplementary Table S9 ), as performed for the tumors in section 3.2. Upon use of the CDKN2A gene expression values corrected by the contribution of the exon1α to the expression of the locus, the accuracy of the prediction of A or non-A CDK4 profiles raised to 90.5% ( Supplementary Table S10 ). The three profile A cell lines were correctly predicted and only two cell lines were falsely predicted to have a profile A. These two cell lines (HTC-C3 and THJ-29T) have mutated CDKN2A with high expression of the locus. These data further validate the use of the gene expression signature to predict the presence of CDK4 phosphorylation and hence the potential sensitivity to CDK4/6 inhibition. Nevertheless, in case of high expression of CDKN2A, potential mutations of the locus will have to be identified to exclude this confounding factor.

We have observed that MAPK/ERK activity was increased upon palbociclib treatment, as judged by the elevated phosphorylation levels of ERK in most sensitive cell lines, and this mirrored the effect seen for Cyclin D1 ( Figure 5A ). In 12 of 18 sensitive cell lines, palbociclib also induced increased phosphorylation of p70S6 kinase and an electrophoretic mobility shift of this protein to a slower-migrating form, suggestive of hyperphosphorylation (69). By contrast, a palbociclib effect on 4EB-P1 and ribosomal protein S6 (RPS6) was less apparent ( Figure 5A ). We have thus investigated the possible synergism between CDK4/6 inhibition and BRAF and/or MEK inhibition using a combination of dabrafenib (BRAF inhibitor) and trametinib (MEK inhibitor). This dabrafenib/trametinib combination is the only targeted therapy approved by the FDA for the treatment of patients with locally advanced or metastatic ATC harboring BRAF V600E mutation (70). In a subset of seven ATC-derived cell lines carrying different oncogenic drivers and having distinct responses to palbociclib ( Supplementary Table S11 ), triple combined treatment potentiated drugs action, leading to complete suppression of proliferation ( Figure 5B ). A strong synergistic effect was obtained, as judged by the calculated combination indexes (CI below one) ( Figure 5C ). Double combination of trametinib and palbociclib was equally effective against cell lines possessing mutated oncogenes other than BRAF. The synergism was reached even in cells displaying partial response to palbociclib alone (8505C and CAL-62 lines). Identical drug cooperation was observed using either palbociclib or abemaciclib ( Supplementary Figure S9 ).

To evaluate the long-term effect of the combinations, we performed clonogenic assays at two different concentration schemes, in two RAS- and in two BRAF-mutated cell lines ( Figure 6A ). Tested concentrations were chosen in accordance to the concentration-response curves for each of the drugs and to known pharmacodynamics data [100-200 nM dabrafenib (71); 15-30 nM trametinib (72, 73); 200-260 nM palbociclib (74)]. Hence, higher palbociclib concentration was tested (1µM instead of 500 nM) for the most palbociclib-partially resistant line CAL-62, and trametinib concentration was unchanged between the two treatment schemes for cell lines with BRAF V600E mutation (highly effective at 5 nM). Cell lines treated for 10 d had their ability to form colonies reduced by each of the drugs alone (for cell lines with mutant RAS, this was seen only with higher concentrations). However, when the compounds were removed, there was a complete or almost complete regrowth (except for HTC-C3 line, which was more sensitive to palbociclib and trametinib than other cell lines). In combined treatments (either trametinib with dabrafenib, trametinib with palbociclib or the three drugs), clonogenic potentials of the cell lines were drastically reduced during the continuous presence of drugs ( Figure 6A ). Importantly, following drugs withdrawal, only regimens combining palbociclib to BRAF and/or MEK inhibitors showed the most efficient and durable responses, completely abolishing the regrowth of most cell lines. For RAS mutant lines, this effect was observed only with higher concentration of drugs, while the triple combination was already effective at lower concentrations in BRAF mutant lines. The same effects were also seen for cell lines possessing other oncogenes ( Supplementary Figure S10 ). We observed that SA-β-gal activity was higher following combined treatment than with each of the drugs alone ( Figure 6B ), whereas neither cleaved caspase 3 nor cleaved PARP1 were changed (data not shown). Thus, inhibition of ERK pathway synergizes with CDK4/6i, potentiating the proliferation-suppressive properties of each of the drugs and strongly impairing the long-term growth of thyroid cancer cells and the development of resistance.

We next explored the molecular mechanisms behind the synergistic inhibition of cell cycle progression and cell proliferation in response to combined targeting of CDK4/6 and MEK/ERK pathways. We observed that a 3 d treatment with the combination of the drugs at low concentrations, similar to the ones used for clonogenic assays ( Figure 6A ), already allowed to efficiently arrest the cell proliferation by repressing RB phosphorylation and Cyclin A2 expression ( Figure 7A ). In addition, the combination with trametinib and dabrafenib prevented the upregulation of both Cyclins D and Cyclin E1, as well as the activation of ERK and mTOR pathways (as inferred by the decreased phosphorylation levels of ERK, p70S6K1, RPS6 and 4EB-P1), that were observed upon palbociclib treatment ( Figure 7A ). These effects were also maintained or further enhanced when 250 nM of palbociclib was used instead of 62.5 nM ( Supplementary Figure S11 ).

We then evaluated the effect of the combined treatments on the formation and specific kinase activities of CDK4 complexes, by performing co-immunoprecipitation assays from cells treated for 6 d ( Figure 7B , upper panel), with the high concentration scheme of drugs tested in the clonogenic assays ( Figure 6A ). The in vitro RB-kinase activity of CDK4 complexes (assayed in the absence of palbociclib, hence reflecting the presence of activated CDK4) was increased by palbociclib mostly in CAL-62 cells (the more palbociclib-resistant cell line) and to a lesser extend in HTC-C3, 8505C and C-643 cell lines. This increase was prevented by the co-treatment with trametinib in CAL-62 cells or by the dabrafenib/trametinib combination in the BRAF-mutated cells ( Figure 7B , upper panel). In CAL-62 (and more weakly in 8505C cells), the dramatic increase of RB-kinase activity of CDK4 induced by palbociclib was associated with a similarly increased formation of cyclin D3-CDK4 complexes, which were devoid of p21 and thus did not require the p21 assembly activity. This observation exactly corresponds to our previous description of a paradoxical effect of palbociclib that specifically stabilized activated cyclin D3-CDK4/6 complexes, which could become hyperactive upon drug removal, because their activity was not restricted by binding to CIP/KIP inhibitory proteins (52). Similar to what we observed in Figure 4A , palbociclib-treated cells exhibited increased levels of either CDK4, Cyclin D3 and/or Cyclin D1, and reduced levels of p27 ( Figure 7B , lower panel). Conversely, trametinib and the dabrafenib/trametinib combination reduced the protein concentrations of CDK4, Cyclin D1 and Cyclin D3, while increasing p27 accumulation ( Figure 7B , lower panel). The p27 modulation is in agreement with the already described reduction of p27 levels upon expression of RAS, BRAF or RET/PTC mutants in thyroid cancer cells, which would be reversed by MEK inhibition (75). These effects, which were variably observed in the different cell lines, might also concur (respectively) to increase the formation of activated CDK4 complexes in response to palbociclib, and to explain the inhibition of this formation by BRAF/MEK inhibitors ( Figure 7B , upper panel). To further characterize the effect of the drugs on CDK4 activity, whole protein extracts from treated cell lines were resolved by 2D-gel electrophoresis ( Figure 7C ). Treatment with palbociclib increased the presence of all the CDK4 forms, and even more the proportion of the activated T172-phosphorylated form. Again, combination with MEK/BRAF inhibitors was able to counteract these effects, consistently with our previous observation that CDK4 phosphorylation depends on MEK activity (24).

In the same experiments, we similarly evaluated the effects of the combined treatments on the formation and RB-kinase activity of CDK2 complexes that were (co-) immunoprecipitated using either CDK2 or Cyclin E1 antibodies ( Figure 7B , middle panel). Intriguingly, the RB-kinase activity associated with Cyclin E1 complexes also increased in response to palbociclib treatment in CAL-62, HTC-C3 and 8505C cells, mirroring the situation observed in CDK4 complexes. Again, this paradoxical effect of palbociclib was prevented by BRAF/MEK inhibitors ( Figure 7B , middle panel). These effects were (respectively) associated with an increased abundance of Cyclin E1 and Cyclin E1-CDK2 complex in palbociclib-treated cells, and a reduction of the binding of CDK2 to Cyclin E1 together with an increased association of p27, in cells treated with MEK or MEK/BRAF inhibitors.

At variance with the situation observed in Cyclin E1 co-immunoprecipitates, the total RB-kinase activity associated with CDK2 complexes (immunoprecipitated using a CDK2 antibody, and thus comprising both Cyclin E1 and Cyclin A2 complexes) was not increased by palbociclib in CAL-62 cells and was reduced in the three other cell lines with the most complete inhibition in HTC-C3 (the most palbociclib-sensitive cells) ( Figure 7B , middle panel). In all the cell lines, accumulation of Cyclin A2 and its association with CDK2 were reduced in palbociclib-treated cells, compensating for the increased presence and activity of Cyclin E1-CDK2 complexes. Intriguingly, treatments with trametinib and with dabrafenib/trametinib used alone (respectively in CAL-62 and in 8505C), reduced Cyclin A2 expression ( Figure 7B , lower panel) but not its presence in CDK2 complexes ( Figure 7B , middle panel). Nevertheless, the most complete inhibitions of both CDK2 activity and CDK2-Cyclin A2 association were observed in cells treated with the combination of palbociclib and MEK/BRAF inhibitors ( Figure 7B , middle panel).

To summarize, the upregulation of Cyclins D with stabilization of activated (T172-phosphorylated) CDK4 complexes and the upregulation of Cyclin E1 with increase of its kinase activity are most probably the mechanisms that bypass CDK4/6 inhibition in thyroid cancer cells. Nonetheless, this can be prevented through combined MAPK/ERK inhibition, which reduces the accumulation of Cyclins D and Cyclin E1, the phosphorylation of CDK4, and increases p27 to further inhibit CDK2 activity.