Athymic BALB/c-nu/nu (nude mice) were acquired from Jackson Laboratory. Mice were provided with autoclaved food and water ad libitum. The study was approved by the Animal Care and Use Committee of the institution and was conducted in compliance with the Public Health Service Guidelines for the Care and Use of Animals in Research.

Four murine thyroid cancer cell lines derived from genetically engineered mouse models of PTC, FTC, PDTC, and ATC were established from primary tumors: PTC with Braf ^V600E mutation (PTC), FTC with Kras ^G12D mutation (FTC), PDTC with both Kras^G12D and Cdkn2a^null mutations (PDTC), and ATC with both Braf ^V600E and Trp53 ^null mutations (ATC). The establishment of PTC, FTC, and ATC strains was described previously (15–18). PDTC strain was established by cross-breeding among Kras^G12D , TPO-Cre, and Cdkn2a^null (strain 01XE4 obtained from The NCI Mouse Repository, https://frederick.cancer.gov/Resources/Repositories/nci-mouse-repository/MouseModels/AvailableStrains#/StrainDetails/79). PDTC was developed from a 13-month-old mouse with both Kras^G12D and Cdkn2a^null mutations. The PDTC cell line was established from the tumor. Thyroid origin was confirmed by genotyping. The cell lines were maintained in DMEM/F12 growth medium containing 10% fetal bovine serum, 100 units/mL penicillin, and 100 μg/mL streptomycin.

To establish pulmonary metastatic thyroid cancer cell lines, 1 × 10⁶ PTC, FTC, PDTC, or ATC cells were injected into the tail vein of five nude mice for each group. Six weeks after injection, pulmonary metastatic tumors were collected aseptically from the mice using blunt dissection, then mechanically dissociated by mincing and passing through a 40-μM mesh sterile screen, and suspended in DMEM/F12 growth medium for 3 months with a total of six passages to eliminate contaminated stromal fibroblasts, lymphocytes, and microphages present in the tumor cell culture. The primary cells were considered permanent cell lines after six passages. The established metastatic cell lines were named PTC-Met1, FTC-Met1, PDTC-Met1, and ATC-Met1. They were re-injected (1 × 10⁶ cells) to a new group of nude mice (n = 5 for each group) via tail vein for enrichment of cells with high metastatic potential. Three weeks following injection, lung metastatic tumors were harvested and propagated in DMEM/Ham’s F12 growth medium for 3 months with at least six passages. The metastatic cell lines were named PTC-Met2, FTC-Met2, PDTC-Met2, and ATC-Met2. The experimental procedures are summarized in Figure 1 .

RNA sequencing (RNA-Seq) was used for quantification of differentially expressed genes (DEGs) between primary (PTC, FTC, PDTC, or ATC) and metastatic thyroid cancer cell lines: PTC-Met1, FTC-Met1, PDTC-Met1, and ATC-Met1, or PTC-Met2, FTC-Met2, PDTC-Met2, and ATC-Met2 cell lines. Total RNA from cell lines was isolated, and libraries were constructed using an Illumina (San Diego, CA, USA) TruSeq RNA Library Prep kit according to the manufacturer’s procedure. Sequencing was performed on Illumina HiSeq 4000 with at least 20 million clean reads. The significant DEGs were selected based on the following criteria: Log2 fold change >2, false discovery rate (FDR) <0.001, and p-value from difference test <0.01. Gene list annotation and enrichment of biological pathways were performed using Metascape (https://metascape.org/gp/index.html#/main/step1).

The common DEGs were analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING; ver.12.0, https://string-db.org/) database for the construction of a protein–protein interaction (PPI) network with a confidence score of ≥0.4. The PPI network was then analyzed and visualized using Cytoscape software (ver.3.10.2, https://cytoscape.org). For the hub gene identification, CytoHubba, a plug-in of Cytoscape, was used. Genes with the degree of a node >10 (with more than 10 interacting genes) were considered hub genes.

The association of metastasis-related genes with OS was performed by the Kaplan–Meier analysis using TCGA-THCA mRNA expression dataset (n = 498) and cBioPortal For Cancer Genomics (https://www.cbioportal.org/).

Prism was used in statistical analysis. Chi-squared test was used when the data were not normally distributed, and the t-test was used when the data were normally distributed. p-Value <0.05 was considered significant. Time-to-event endpoints were visualized using the Kaplan–Meier curves, and differences in OS were compared using Cox’s proportional hazards model.

Significant DEGs were identified from RNA-Seq results. A total of 4,269 DEGs were found in the PTC (PTC vs. PTC-Met2), 4,219 DEGs in FTC (FTC vs. FTC-Met2), 1,897 in PDTC (PDTC vs. PDTC-Met2), and 3,301 in ATC (ATC vs. ATC-Met2) datasets ( Figure 2A ). A total of 130 common DEGs were found after the integration of four datasets, including 105 upregulated and 25 downregulated genes (Log2 fold change >2) present in all four Met2 cell lines ( Figure 2B ; Supplementary Table 1 ). The expression levels of most DEGs were increased in Met2 as compared to Met1 cells. Many DEGs not detected in Met1 cells also appeared in Met2 cells ( Supplementary Table 1 ). These data indicate the enrichment of DEGs in Met2 cells and suggest that these DEGs may be involved in DM. For example, “don’t eat me” signal expression was either no change or mildly increased in Met1 cells but significantly increased in Met2 cells such as Cd274 (PD-L1), Cd52, and Tbxas1 ( Table 1 ), indicating that these genes may play a significant role in immune evasion of metastatic cancer cells. Strictly speaking, Tbxas1 is not a “don’t eat me’ signal”, but it indirectly helps cancer cells evade immune elimination via platelet activation and aggregation (19). Met2 cells also took 3 weeks less than Met1 cells (3 vs. 6 weeks) to form lung metastasis, indicating enrichment of tumor cells with high metastatic potential.

The PPI network for the 130 DEGs was constructed after they were imported to STRING ( Figure 3A ). According to the ranking generated by the Degree algorithm, 32 hub genes were identified ( Figure 3B ). Their full names and related functions are shown in Supplementary Table 2 . Many of them have been reported to be involved in metastasis in different cancer types. Gene ontology and pathway enrichment analysis revealed the top 20 gene clusters in the regulation of cytokine production, inflammatory and negative regulation of immune response, neutrophil migration and phagocytosis, and MAPK cascade ( Figure 4A ). Visualization of these gene clusters by the PPI network demonstrated significant cross interactions among different pathways ( Figure 4B ). The gene list of enriched clusters and their involvement in multiple signaling pathways are shown in Supplementary Table 3 .

The association of common DEGs with OS of PTC patients was analyzed using TCGA-THCA dataset (n = 498). Seven genes were identified, and their overexpression was associated with poor OS (ARG1, TF, CD52, ITGBL1, PTGIR, TENM3, and AS3MT) ( Figure 5A ). The distribution of seven-gene expression signature among PTC samples is shown in Figure 5B . As shown in Table 2 , the expression of this group of genes was higher in Met2 cells than in Met1 cells, indicating that they may participate in DM. Since most of these genes were overexpressed in less than 2% of samples except for ITGBL1 ( Figure 5B ), we analyzed them as a group (13% samples, 63/498) for prognosis prediction and clinical attributes. The seven-gene expression signature was significantly associated with poor OS (p < 0.0001, Figure 5C ). The patients with seven-gene expression signature were associated with old age at diagnosis [56 vs. 45 years of median age (p < 0.0001)], late disease stage [38% (24/63) vs. 20% (85/435) in stage III and 24% (15/63) vs. 9% (39/435) in stage IV tumors (p < 0.0001)], tall cell variant [18% (11/63) vs. 6% (24/435) (p < 0.01)], and higher hypoxia score [−34 vs.− 38 of medium (p < 0.0001)] ( Figure 5D ). The overexpression of ITGBL1 was found in 8% of samples (40/498). Its overexpression was also associated with old age at diagnosis [7 vs. 46 years of median age (p < 0.0001)], late disease stage [43% (17/40) vs. 20% (92/458) in stage III and 30% (12/40) vs. 9% (42/458) in stage IV tumors (p < 0.0001)], tall cell variant [23% (9/40) vs. 6% (26/458) (p < 0.0001)], and higher hypoxia score [−35 vs. −42 of median value (p < 0.01)] ( Figure 5E ). Thus, ITGBL1 overexpression may be the most useful single gene marker for poor prognosis prediction.

We next compared mutational load between samples with or without seven-gene expression signature using TCGA-THCA dataset. We chose 38 genes whose mutations have been reported to be involved in carcinogenesis including genes known to be driver mutations in thyroid carcinogenesis. Among samples with seven-gene expression signature, 26 samples had more than one driver mutation (47%, 26/55) ( Figure 6A ). Deep homozygous deletions were frequently found in samples with BRAF ^V600E mutation ( Figure 6A ). Among patients without seven-gene expression signature, only 24 patients had more than one driver mutations (7%, 24/343) ( Figure 6B ). The difference is highly statistically significant (p < 0.00001). These data suggest that the seven-gene expression signature may reflect underlying increased chromosomal instability. Indeed, as shown in Figure 6C , patients with seven-gene expression signature had a higher rate of fraction of genome altered (FGA; the percentage of genome that has been affected by copy number gains or losses) than those without: 0.0399 ± 0.012 (mean ± SEM, n = 62) vs. 0.0153 ± 0.003 (n = 434, p < 0.01). Aneuploidy score was also increased in patients with seven-gene expression signature than those without: 1.60 ± 0.57 (mean ± SEM, n = 60) vs. 0.77 ± 0.14 (n = 405). These genetic changes may drive the metastatic transformation of tumor cells, resulting in distant metastasis.