FAM-related genes were sourced from the MSigDB available at https://www.gsea-msigdb.org/gsea/msigdb, a comprehensive resource that curates gene sets to facilitate functional enrichment analysis and the biological interpretation of genomic data (12–15). The Cancer Genome Atlas (TCGA) is a comprehensive database that integrates genomic and clinical data across various cancer types, serving as a valuable resource for cancer research (16). The TCGA cohort includes 59 normal thyroid tissues and 513 PTC tissues. We used survival, clinical features, and expression data in the TCGA cohort. Gene Expression Omnibus (GEO) is a public repository designed for the storage and sharing of gene expression and functional genomics data from high-throughput experiments. We used the GSE29265 and GSE33630 cohorts. The GSE29265 cohort includes 20 normal and 20 PTC tissues (17). The GSE33630 cohort includes 45 normal and 49 PTC tissues (18, 19).

The scRNA-seq data of GSE232237 were obtained from the GEO, which included three normal thyroid, seven PTC, and five anaplastic thyroid cancer (ATC) cases (20). In this study, we selected seven PTC samples to perform our analysis. The R package “Seruat” was utilized to perform the single-cell analysis. According to cell-specific markers, we isolated eight types of cells, namely, thyrocytes (“TG,” “EPCAM,” “KRT18”), T cells (“CD3D,” “CD3E,” “CD8A”), B cells (“CD79A,” “MS4A1,” “IGHG1”), myeloid cells (“CD68,” “CD163,” “LYZ”), endothelial cells (“PECAM1,” “VWF,” “CLDN5”), plasma cells (“MZB1,” “SDC1,” “CD79A”), mast cells (“CST3,” “KIT,” “CPA3”), and fibroblasts (“FGF7,” “COL1A1,” “COL1A2”).

The RStudio (Version 4.4.1) software was used to construct the FAM-related gene prognostic signature. In this study, we employed several R packages, including “survival,” “survminer,” “heatmap,” “survivalROC,” “ComplexHeatmap,” “ggplot2,” “ggpubr,” “ggExtra,” and “limma.” To build the prognostic signature, the “survival” R package was employed for multivariate Cox regression analysis. Survival analysis utilized the “survival” and “survminer” R packages. Additionally, the “survival” R package was also used for both univariate and multivariate Cox regression analyses. A heatmap of three FAM-related genes was generated utilizing the “heatmap” R package. The ROC curve was plotted employing the “survivalROC” R package. The “ComplexHeatmap” R package was used to depict the association between clinical features and risk scores.

Gene set enrichment analysis (GSEA) is a computational method that evaluates whether predefined gene sets exhibit significant changes in expression across distinct biological states or conditions. In this study, the GSEA software was utilized to carry out function analysis. The plots for multiple GSEA were generated using the “plyr,” “ggplot2,” “grid,” and “gridExtra” R packages.

Mutation data were downloaded from the TCGA database. Mutation analysis mainly employed the “maftools” R package to generate waterfall plots, which visually represented the mutation frequency and types across samples, displaying gene-wise and sample-wise mutation distributions to highlight common mutations and patterns in genomic data.

The CIBERSORT algorithm, based on linear support vector regression (SVR), is a widely recognized and reliable machine learning technique for deconvolving expression matrices to evaluate the ratio of 22 distinct human immune cell subtypes (21). In this study, based on the CIBERSORT algorithm, we evaluated the immune cell abundance in each PTC sample from the TCGA database. In addition, we used two R packages: “preprocessCore” and “limma”.

Next, we utilized seven algorithms, namely, CIBERSORT-ABS (21), CIBERSORT (21), Tumor Immune Estimation Resource (TIMER) (22), QUANTISEQ (23), MCPcounter (24), xCell (25), and EPIC (26), to comprehensively analyze the immune landscape of two groups, which enabled us to assess various aspects of immune cell infiltration and characterize the immunological features in greater detail. We used the following eight R packages: “limma,” “scales,” “ggplot2,” “ggtext,” “reshape2,” “tidyverse,” “ggpubr,” and “pheatmap.” Additionally, immune function analysis was carried out utilizing five R packages: “limma,” “GSVA,” “GSEABase,” “ggpubr,” and “reshape2”.

In this study, all cell lines were purchased from the China Type Culture Collection (CTCC). Nthy-ori 3-1 is a normal thyroid cell line, while TPC-1, K1, KTC-1, and BCPAP are PTC cell lines. Nthy-ori 3-1 cells were cultured in Dulbecco’s modified eagle medium (DMEM). TPC-1, KTC-1, and BCPAP cells were maintained in Roswell Park Memorial Institute 1640 (RPMI-1640) medium, and K1 cells were grown in DMEM/F12 (DMEM and Ham’s F-12 nutrient mixture) medium. All cells were grown in media with 10% fetal bovine serum (FBS) and 1% penicillin–streptomycin and incubated at 37°C with 5% CO₂.

For the cell transfection of TPC-1, cell transfection was initiated when the cell density reached approximately 30%. Lipofectamine 2000 (Lip2000) was used as the transfection reagent. First, the optimized minimum essential medium was separately mixed with Lip2000 and siRNA, and each mixture was incubated for 5 min. The two solutions were then combined for an additional 15 minutes. The resulting transfection mixture was added to the complete RPMI-1640 medium without antibiotics. After 6–8 h, the medium was replaced, and SCD knockdown efficiency was evaluated 24 h post-transfection.

According to the manufacturer’s instructions, the TRIzol reagent was employed to acquire total ribonucleic acid). Next, complementary DNA (cDNA) was obtained by performing a reverse transcription assay. Lastly, qPCR was utilized to assess relative messenger ribonucleic acid (mRNA) expression. The 2^−ΔΔCt method was used to calculate relative expression. The primer names and primer sequences are shown in Supplementary Table S1 .

For the CCK8 assay, after successful transfection of TPC-1 cells, 3,000 cells from both the control and knockdown groups were cultured in 96-well plates. Cell viability was checked by incubating the cells for 0, 24, 48, and 72 h, and optical density (OD) at 450 nm was determined at each time point using a microplate reader.

For the scratch assay, after successful transfection of TPC-1 cells, cells were cultured until they reached 100%. A scratch was made using a 200-µL pipette tip, followed by PBS washing to remove the medium. Images were captured immediately after scratching and then again at 24 and 48 h to monitor wound healing.

For the colony formation assay, after successful transfection of TPC-1 cells, 1,000 cells from both the control and knockdown groups were seeded into six-well plates. The medium was replaced every 3–5 days. After 14 days, the plates were collected, washed with PBS, fixed with paraformaldehyde, and stained with crystal violet, and the images were taken after the plates dried using an Olympus microscope.

For the Transwell assay, after successful transfection of TPC-1 cells, 30,000 cells from both the control and knockdown groups were mixed with 260 µL of serum-free medium and transferred to the top chamber of the Transwell plates. The bottom chamber contained 750 µL of complete medium. For the invasion assay, 20% Matrigel was precoated on the bottom of the upper chamber, while for the migration assay, no Matrigel was used. After 24 h, the upper chamber was removed, washed, fixed, stained, and imaged using an Olympus microscope.

RStudio (Version 4.4.1) and GraphPad Prism (Version 10.1.0) were utilized to perform statistical analysis. In this study, a p-value less than 0.05 was considered statistically significant.

FAM-related genes were acquired from the MSigDB database (HALLMARK, KEGG, REACTOME), and 309 FAM-related genes were used for follow-up analyses ( Supplementary Table S2 ). We analyzed DEGs from 309 FAM-related genes based on logFC more than 1 and p-value <0.05, in which 24 DEGs were identified and employed to build the prognostic signature ( Supplementary Table S3 ). We used a heatmap and a volcano plot to display 25 DEGs associated with FAM ( Figures 1B, C ). Furthermore, we selected survival-related genes through univariate Cox regression analysis, and four FAM-related genes were selected ( Figure 1D ). Lastly, multivariate Cox regression analysis was utilized to build a three FAM-related gene prognostic signature ( Figure 1E ).

Through the above analyses, we successfully constructed the prognostic signature based on three FAM-related genes. It was observed through survival analysis that compared with the low-risk group, the high-risk group had a shorter survival probability ( Figure 2A ). The receiver operating characteristic (ROC) curve demonstrates that the prognostic signature exhibits high predictive efficacy, with area under the curve (AUC) values of 0.95 at 1 year, 0.816 at 2 years, and 0.686 at 3 years ( Figure 2B ). Univariate and multivariate regression analyses showed that stage and risk score are independent prognostic factors for PTC patients ( Figures 2C, D ). We depicted the heatmap of the three FAM-related genes in this risk group ( Figure 2E ). In addition, we also displayed the risk score curve ( Figure 2F ) and survival state distribution ( Figure 2G ), which confirmed the reasonable distribution between the two risk groups. Next, we integrated the clinical characteristics of PTC and analyzed their correlation with risk score, presenting the results in a heatmap for clarity ( Figure 2H ). The heatmap revealed that risk score is associated with age, T stage, and N stage. Lastly, a predictive nomogram was generated by integrating multiple clinical factors with risk scores ( Figure 2I ). These scores can be employed to estimate the survival probabilities of PTC patients.

All PTC samples were categorized into high and low groups according to the median expression levels of the hub genes. Next, we performed GSEA to explore the potential functional enrichment of the hub genes. Through GSEA of Gene Ontology (GO), we observed that the high-risk group was related to the AMINO ACID CATABOLIC PROCESS (BP), MONOATOMIC ANION TRANSPORT (BP), FATTY ACID TRANSMEMBRANE TRANSPORT (BP), ABC TYPE TRANSPORTER ACTIVITY (MF), and ATPASE COUPLED TRANSMEMBRANE TRANSPORTER ACTIVITY (MF), and the low-risk group was related to CYTOSOLIC LARGE RIBOSOMAL SUBUNIT (CC), NUCLEAR MEMBRANE REASSEMBLY (BP), BLOC 1 COMPLEX (CC), CYTOSOLIC SMALL RIBOSOMAL SUBUNIT (CC), and CADHERIN BINDING INVOLVED IN CELL-CELL ADHESION (MF) ( Figure 3A ). Through GSEA of KEGG, we observed that the high-risk group was associated with LYSINE DEGRADATION, ABC TRANSPORTERS, FATTY ACID METABOLISM, TYPE II DIABETES MELLITUS, and TRYPTOPHAN METABOLISM, and the low-risk group was associated with RIBOSOME, DNA REPLICATION, PROTEASOME, PATHOGENIC ESCHERICHIA COLI INFECTION, and P53 SIGNALING PATHWAY ( Figure 3B ).

To investigate the relationship between risk score and gene mutations, we analyzed the heterogeneity of gene mutations between the high-risk and low-risk groups. Compared with the high-risk group (64.85%), we uncovered that the low-risk group possessed a higher gene mutation ratio (85.31%) ( Figures 3C, D ). Furthermore, we also observed that BRAF and NRAS gene mutations are the most significant in those genes. The mutation ratio of BRAF is 75% in the low-risk group and 44% in the high-risk group. The mutation ratio of NRAS is 4% in the low-risk group and 12% in the high-risk group.

In order to observe the immune condition of PTC patients and explore the relationship between FAM-related genes and tumor immunity, the CIBERSORT algorithm was employed to visualize the immune landscape of all PTC patients in the TCGA cohorts. The immune cell ratio of PTC patients was displayed in a heatmap ( Figure 4A ). The correlation results provided a comprehensive overview of the associations between risk scores and immune cells and stromal cells. We also employed a heatmap to visualize the distribution of various immune cells across samples with different risk scores ( Figure 4B ). Next, based on seven algorithms, we evaluated the association between immune cells and risk score ( Figure 4C ). We also showed the correlation between immune cells and risk score. The results showed that risk score is closely associated with multiple immune cells ( Figures 4D–O ) based on the xCell algorithm.

Furthermore, we used boxplots to illustrate the changes in immune cell abundance between the two groups, which revealed that 17 immune cell types were present in lower quantities in the high-risk group compared to the low-risk group ( Figure 5A ). Based on the TIMER algorithm, we specifically analyzed the relationship between six immune cells and risk scores, which showed a positive correlation for all six immune cell types ( Figures 5B–E ).

Immune checkpoints and immunogenic cell death (ICD)-related genes showed significant potential in tumor treatment, offering promising avenues to enhance antitumor immunity and improve therapeutic outcomes (27–29). We analyzed the expression differences of ICD-related genes and immune checkpoints between the high-risk and low-risk groups and observed that most immune checkpoints were upregulated in the low-risk group ( Figures 5F, G ). This suggests that, unlike the high-risk group, PTC patients in the low-risk group are more likely to respond favorably to immunotherapy.

As a malignant tumor, DEGs may play an important role in PTC. Firstly, we analyzed the expression of three FAM-related genes across the TCGA, GSE29265, and GSE33630 cohorts. The results showed that, compared to PTC tissues, ACACB and ADH1B were more highly expressed in normal thyroid tissues, whereas SCD exhibited higher expression in PTC tissues ( Figures 6A–C ). Next, we used qPCR experiments to verify the expression of three FAM-related genes in normal thyroid cell lines and PTC cell lines. The results were consistent with the findings from the above bioinformatics analysis ( Figure 6D ). Moreover, using single-cell cohorts of PTC, we validated the relative expression levels of key genes across different cell types ( Figure 7A ). Our analysis revealed that ACACB is highly expressed in fibroblasts and endothelial cells ( Figure 7B ), ADH1B shows elevated expression in fibroblasts ( Figure 7C ), and SCD is predominantly expressed in myeloid cells ( Figure 7D ). Given that SCD is the only highly expressed gene among the three FAM-related genes, we plan to further investigate its role in PTC cells. We used siRNA to downregulate SCD expression in TPC-1 cells. Through the qPCR assay, SCD expression is significantly downregulated after using siRNA ( Figure 8A ). The CCK8 assay revealed that silencing SCD expression decreased the cell proliferation activity in the TPC-1 cell line ( Figure 8B ). The colony formation assay uncovered that colony formation ability was downregulated after silencing SCD expression ( Figure 8C ). Scratch assays observed that silencing SCD expression inhibited TPC-1 cell migration ( Figure 8D ). Transwell assays showed that silencing SCD expression suppressed TPC-1 cell migration and invasion ( Figure 8E ). Therefore, we concluded that SCD is a new PTC biomarker and silencing SCD expression inhibits PTC progression.