Marker combination for thyroid nodule risk assessment and application thereof

Abstract

The invention discloses a marker combination for thyroid nodule risk assessment and application thereof. Based on large-sample multi-omics sequencing data, papillary thyroid carcinoma (PTC) and benign nodules are found to have a common epithelial spectrum origin, epithelial intermediate state cells (EICs) with malignant potential are found in the benign nodules, and typical molecular markers of the EICs are further identified to comprise ZCCHC12, CDLN1 and NPC2. The three genes participate in proliferation and differentiation of cells, the expression in PTC is remarkably up-regulated, and the malignant risk of the benign thyroid nodules can be evaluated by detecting the expression level of the markers in a biological sample. According to the method, a means for recognizing the high-risk benign nodules before malignant transformation occurs is provided for the first time, risk perspectiveness is achieved, and higher specificity and sensitivity are possibly achieved for recognizing the benign nodules with malignant potential.

Description

Technical Field

[0001] This invention relates to the fields of molecular diagnostics and tumor detection, and more specifically to a combination of biomarkers for risk assessment of thyroid nodules and their applications. Background Technology

[0002] Globally, the ultrasound detection rate of thyroid nodules ranges from 19% to 68%, with over 85% being benign. The ATA guidelines recommend annual ultrasound follow-up for thyroid nodules with benign cytological or ultrasound findings to monitor for signs of malignant transformation. This strategy primarily relies on ultrasound imaging features such as the TI-RADS classification, stratifying the risk of nodules based on their morphology, borders, echogenicity, and calcification to guide follow-up frequency and further intervention decisions. However, this method has a degree of subjectivity and struggles to identify high-risk nodules with malignant potential at the molecular level. The large number of thyroid nodules not only creates a heavy follow-up burden but also carries a significant disease burden. The root cause lies in the lack of a clear understanding of whether and how thyroid nodules become cancerous. Therefore, a more precise risk stratification strategy is urgently needed in clinical practice to address the challenges of managing a large number of thyroid nodules.

[0003] The pathogenesis of thyroid cancer remains unclear. Unlike the well-established continuous evolution of colorectal cancer ("normal epithelium-polyp-carcinoma"), and unlike cervical cancer and esophageal squamous cell carcinoma (which originate directly from normal epithelium through dysplasia without undergoing a benign lesion stage), the evolutionary process of thyroid cancer is still not agreed upon. Papillary thyroid carcinoma (PTC), accounting for over 90% of thyroid cancers, has received widespread attention for its origin and evolutionary mechanisms. Previous omics sequencing studies have revealed distinct genomic mutation characteristics between benign nodules and PTC, supporting the independent origin theory of PTC. However, epidemiological studies indicate that some benign thyroid nodules still have the potential to develop into thyroid cancer, although the overall probability is less than 1%. Furthermore, with advancements in molecular diagnostic techniques, researchers have observed epithelial subclonal amplification in some benign nodules, suggesting a potential for further malignant transformation, but direct and compelling evidence to confirm this remains lacking.

[0004] Therefore, the location of benign nodules in the development of thyroid cancer and their relationship with malignant lesions remain not fully understood. Currently, there are no molecular markers to assess the risk of benign thyroid nodules progressing to thyroid cancer.

[0005] It is important to note that PTC is a tumor with a relatively low mutation burden. Apart from identified driver mutations including BRAF V600E (40-60%), RET / PTC gene fusion (10-20%), and RAS mutations (10-20%), it has very few background somatic mutations. This low mutation rate places higher demands on the sensitivity of mutation detection and limits the resolution of traditional DNA sequencing in clonal evolution tracking of PTC.

[0006] In current technologies, imaging methods, including magnetic resonance imaging (MRI), conventional CT, and ultrasound, are commonly used for the detection of thyroid cancer (Gao Jinliang et al., Diagnostic Value of Enhanced CT in Neck Metastasis of Thyroid Cancer, Chinese Journal of CT and MRI). However, imaging detection usually relies on large equipment, resulting in high costs, and various factors can lead to missed detections. Wang Qiqi et al. used serum microRNA-599 and Galectin-3 levels to predict metastasis of papillary thyroid carcinoma (Wang Qiqi et al., Predictive Value of High-Frequency Ultrasound Combined with Serum microRNA-599 and Galectin-3 in Neck Lymph Node Metastasis of Papillary Thyroid Carcinoma, Chinese Journal of Modern Medicine). However, the extraction and separation of miRNA molecules from serum is difficult, usually requiring specific kits for measurement. Ren Yi et al. have provided auxiliary prediction of thyroid cancer metastasis by monitoring the expression of Ki67 and Survivin proteins (Ren Yi et al., Predictive Value of Ki67 and Survivin Protein Expression in Thyroid Cancer Tissue for Postoperative Recurrence and Metastasis, International Journal of Laboratory Medicine). However, commonly used methods for detecting protein levels, such as immunohistochemistry and proteoblotting, are complex to operate. Therefore, there is a need to further explore new risk assessment biomarkers for thyroid cancer to provide a simple and accurate method for detecting thyroid cancer and assessing its prognosis. Summary of the Invention

[0007] The purpose of this invention is to overcome the above-mentioned deficiencies of the prior art and to provide a method for differentiating between benign and malignant thyroid nodules and assessing the risk of benign nodules progressing to malignancy based on a set of novel gene markers (CDLN1, NPC2, and ZCCHC12). This method is achieved by detecting the expression levels of the above markers in biological samples, and is used to assess the nature of thyroid nodules and the risk of malignant transformation of benign thyroid nodules, thereby providing decision support for clinicians.

[0008] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0009] This invention discloses a combination of biomarkers for risk assessment of thyroid nodules, the biomarkers including: CDLN1, NPC2 and ZCCHC12.

[0010] Preferably, the markers consist of the following markers: CDLN1, NPC2, and ZCCHC12.

[0011] This invention discloses a biomarker composition for risk assessment of thyroid nodules, the biomarkers including: CDLN1, NPC2 and ZCCHC12.

[0012] Preferably, the markers consist of the following markers: CDLN1, NPC2, and ZCCHC12.

[0013] This invention discloses the use of CDLN1, NPC2, and ZCCHC12 in the preparation of reagents and / or kits for detecting epithelial intermediate state cells (EIC).

[0014] Preferably, the intermediate state epithelial cells are intermediate state epithelial cells of benign thyroid nodules.

[0015] This invention discloses the use of CDLN1, NPC2, and ZCCHC12 in the preparation of reagents and / or kits for assessing the risk of thyroid nodules.

[0016] Preferably, the risk is assessed by detecting the expression of CDLN1, NPC2, and ZCCHC12.

[0017] This invention discloses the use of CDLN1, NPC2 and ZCCHC12 in the preparation of reagents and / or kits for distinguishing between normal thyroid tissue and thyroid cancer.

[0018] This invention discloses the use of CDLN1, NPC2 and ZCCHC12 in the preparation of reagents and / or kits for predicting the benign and malignant nature of thyroid nodules.

[0019] This invention discloses the use of CDLN1, NPC2, and ZCCHC12 in the preparation of reagents and / or kits for predicting the benign or malignant nature of thyroid tumors.

[0020] Preferably, the assessment of the risk of benign thyroid nodules includes distinguishing between individuals who may progress to thyroid cancer and healthy individuals.

[0021] Preferably, the assessment of thyroid nodule risk includes distinguishing between thyroid cancer and benign thyroid nodules.

[0022] Preferably, the thyroid cancer is papillary thyroid carcinoma.

[0023] This invention discloses a kit for detecting epithelial intermediate state cells (EIC), the kit comprising reagents for detecting the expression levels of CDLN1, NPC2 and ZCCHC12.

[0024] This invention discloses a kit for assessing the risk of thyroid nodules, the kit comprising reagents for detecting the expression levels of CDLN1, NPC2 and ZCCHC12.

[0025] Preferably, the reagents include reagents for determining the sequences of CDLN1, NPC2, and ZCCHC12, primers for amplifying CDLN1, NPC2, and ZCCHC12, or reagents for detecting the protein expression levels of CDLN1, NPC2, and ZCCHC12.

[0026] Preferably, the kit includes quality control materials.

[0027] This invention discloses a set of gene markers for assessing the risk of malignant transformation of benign thyroid nodules, detection products containing the markers, and methods for their application.

[0028] This invention discloses the application of gene markers CDLN1, NPC2 and ZCCHC12 in thyroid nodule risk assessment.

[0029] This invention discovers that some benign thyroid nodules can become malignant, and that a special subpopulation of precancerous cells—intermediate epithelial cells (EICs)—exists in their epithelial cells. It provides molecular markers for the malignant transformation of benign thyroid nodules, namely ZCCHC12, CLDN1, and NPC2. The method for quantitatively detecting these markers is applied to products that identify intermediate-state benign thyroid nodules in the early stages of malignant transformation.

[0030] This invention also provides a method for screening the above-mentioned molecular markers for malignant transformation of thyroid nodules, comprising the following steps:

[0031] (1) Single-cell RNA sequencing was performed on samples from patients with papillary thyroid carcinoma, including normal, benign nodules, primary tumors, and lymph node metastases; malignant trajectory construction was performed on epithelial cells to obtain a continuous variation lineage of thyroid cells, and high heterogeneity was found in benign nodule cells;

[0032] (2) The transcriptome similarity between tumor cells and adjacent normal and benign nodule cells was calculated, and it was found that tumor cells have multiple patterns such as benign nodule origin and normal epithelial origin; the somatic mutations identified by deep whole exon sequencing data confirmed that some benign nodules have malignant potential.

[0033] (3) Construct pseudo-temporal trajectories with epithelial single-cell precision. By calculating the transcriptomic distance between benign nodules and tumor cells and normal cells, cells that are closer to tumor cells and farther from adjacent normal cells are screened and defined as EICs.

[0034] (4) By calculating the differential genes between EICs and other benign nodular cells, and further screening, the markers ZCCHC12, CLDN1 and NPC2 were obtained. Their expression levels gradually increased with the malignant progression of the tumor, were negatively correlated with TDS, and were closely related to BRAF phenotype and EMT signal, suggesting their potential role in tumorigenesis.

[0035] (5) Based on candidate biomarkers, a logistic classifier model is constructed to distinguish between papillary thyroid carcinoma and normal tissue.

[0036] This invention, through single-cell transcriptomic analysis of normal tissues, benign nodules, primary tumors, and lymph node metastases from PTC patients, discovered a specific subset of epithelial cells (EICs) within benign nodule epithelial cells. The transcriptomic characteristics of this cell population fall between those of normal epithelial cells and tumor cells, suggesting a potential key role in the transformation from benign to malignant. Further analysis revealed significantly high expression of ZCCHC12, CLDN1, and NPC2 in EICs, with expression levels gradually increasing with tumor malignancy, indicating potential molecular markers for malignant transformation of benign thyroid nodules.

[0037] At the clinical application level, the research findings can be used for:

[0038] (1) By detecting the expression levels of ZCCHC12, CLDN1 and NPC2 genes, the malignant potential of benign thyroid nodules was predicted based on the RiskScore model;

[0039] (2) Guide the development of individualized follow-up strategies, such as shortening the follow-up interval for high-risk patients and avoiding excessive examinations for low-risk patients.

[0040] Based on large-sample multi-omics sequencing data, this invention discovered that PTCs and benign nodules share a common epithelial lineage origin and identified intermediate epithelial cells (EICs) with malignant potential within benign nodules. It systematically depicts the continuous and dynamic malignant evolution of PTCs from "normal tissue (N) - benign nodule (B) - malignant tumor (PT) - metastasis." This provides direct evidence for understanding the origin of PTCs and offers a new explanatory framework for the long-standing theories of "independent origin" and "continuous evolution."

[0041] This invention not only reconstructs the independent evolutionary pathways of benign and malignant thyroid nodules, but also reveals for the first time that approximately one-third of PTC cases have a benign origin, meaning that the tumor and paired benign nodules share somatic mutations and maintain lineage continuity at the transcriptome level. This result overturns the previous view based on bulk sequencing that benign nodules are completely independent of PTC in origin, suggesting that the occurrence of PTC may not be a single pathway, but rather a combination of "independent origin" and "continuous evolution." This dual origin model contrasts with the classic models of "adenoma-carcinoma" in colorectal cancer or "independent origin" in cervical cancer, and provides a reasonable mechanistic basis for explaining the contradictions between previous molecular and epidemiological results.

[0042] This invention employs deeper DNA sequencing than traditional studies (equivalent to three times the depth of traditional studies), thereby capturing rare shared mutation events between benign and malignant nodules. Furthermore, we observed high heterogeneity in benign nodules. Previous studies, due to limited sample size (21 patients with both benign nodules and PTC) and insufficient coverage of the diversity of benign nodules, may have drawn biased conclusions.

[0043] This invention reveals a malignant continuum of thyroid epithelial cells: from morphologically normal tissue, through some benign nodules, a gradual accumulation of transcriptomic alterations, ultimately developing into cancer and lymph node metastasis. In particular, EICs (endothelial cells with thyroid cancer) were found in some benign nodules, exhibiting transitional features between normal and tumor types, including decreased TDS and differentiation, enhanced EMT characteristics, and lipid metabolism reprogramming. These cells stably locate in the middle position in single-cell trajectories, suggesting they may represent a true precancerous state. This result provides a molecular explanation for the epidemiologically suggested low probability of malignant transformation in benign nodules.

[0044] This invention identified typical molecular markers for thyroid nodules (EICs), including ZCCHC12, CDLN1, and NPC2. These three genes are involved in cell proliferation and differentiation and have been reported to be significantly upregulated in various tumor tissues, including PTC. We found that the expression of these genes exhibits high diversity in benign nodules. IHC results showed that their expression levels were significantly increased in tumor-paired benign nodules compared to independent benign nodules, especially in nodules that subsequently became malignant. This result suggests that EIC molecular markers hold promise for risk stratification of benign thyroid nodules, guiding clinical management of nodules with malignant potential.

[0045] This invention has the following significant advantages and positive effects:

[0046] (1) Precision and foresight: Existing technologies (such as ultrasound and cytology) are mainly used to diagnose malignant tumors that have already formed or to perform a one-size-fits-all follow-up on benign nodules. This invention provides for the first time a means to identify high-risk benign nodules before malignant transformation occurs, achieving true prospective risk stratification and changing the clinical practice model.

[0047] (2) High specificity and sensitivity: This biomarker combination is derived from the direct analysis of the malignant evolution trajectory of thyroid epithelial cells and directly targets precancerous lesion cells (EICs). Therefore, compared with detection based on known oncogenes, it may have higher specificity and sensitivity for identifying benign nodules with malignant potential.

[0048] Directly addressing clinical pain points: The products and methods provided by this invention can effectively distinguish between the vast majority of indolent benign nodules and the very few nodules with malignant potential. This is expected to spare most low-risk patients from excessively frequent follow-ups and the resulting psychological burden, while allowing medical resources to be more concentrated on close monitoring and early intervention for high-risk groups, significantly reducing the social medical burden.

[0049] (3) Complete chain of evidence: The technical route of this invention includes scRNA-seq, Bulk RNA-seq verification and TCGA data verification. The data foundation is solid and has been verified by immunohistochemistry (IHC) and cell function experiments, proving the existence of EIC cells and the expression pattern and cancer-promoting function of their markers. The theory and technical feasibility are high. Attached Figure Description

[0050] Figure 1a A schematic overview of the single-cell and bulk RNA sequencing cohorts used in this study.

[0051] Figure 1b A UMAP (Uniform Manifold Approximation and Projection) plot containing 547,280 cells from a scRNA-seq dataset. Cells are colored and annotated according to their major cell types, with text labels indicating the corresponding cell subtypes.

[0052] Figure 1c UMAP visualization of epithelial cell subsets, with cells stained according to their assigned clusters. Four panels illustrate the different distributions of epithelial cells in different tissue origins (N, B, T, LNM).

[0053] Figure 1d The left scatter plot depicts the correlation between EMT (epithelial-mesenchymal transition) score and the proportion of C4_SPARC cells, while the right scatter plot illustrates the correlation between IFN-γ response score and the proportion of C5_HLA-DRA cells.

[0054] Figure 1e The violin plot illustrates the distribution of CytoTRACE scores in different epithelial cell clusters, and the Student's test is used to assess statistical significance.

[0055] Figure 1f Malignant trajectories inferred from scRNA-seq data. Each point represents a sample and is stained by the number of differentially expressed genes (DEGs), tissue origin (N, B, T, LNM), or malignancy score, highlighting a continuous spectrum of cellular state transitions.

[0056] Figure 1g Cumulative percentage of malignancy scores for cells of different tissue origins (N, B, T, LNM). Orange shading indicates the difference in malignancy between benign nodules and tumor cells.

[0057] Figure 1h The scatter plot depicts the negative correlation between single-cell malignancy score and TDS score. Each point represents a sample, and the regression line and corresponding R and P values ​​are derived from the Pearson correlation test, indicating the overall trend.

[0058] Figure 1i Malignant trajectories inferred from bulk RNA-seq data from the SYSU cohort. Samples are chromatographic by DEGs, histological type (N, B, NM [non-metastatic], M [metastatic]), and malignancy score. Branch 1 and branch 2 represent different trajectories of tumor progression.

[0059] Figure 1j Cumulative percentage of malignancy scores for samples of different tissue types (N, B, NM, M) along bulk RNA-seq trajectory branches 1 (k) and 2 (l).

[0060] Figure 1k Box plots show the distribution of malignancy scores for IB and MB samples in the single-cell malignancy trajectory. P-values ​​represent statistical significance of differences between groups, assessed using a one-sided Student's t-test.

[0061] Figure 2a UMAP plots from paired scRNA-seq data containing samples of normal tissue (N), benign nodules (B), and primary tumors (T). The top plot shows cells stained by patient ID, and the bottom plot shows cells stained by tissue origin (N, B, T).

[0062] Figure 2b The bar chart shows the proportion of tumor cells in each patient (y-axis) whose transcriptome similarity tends to be normal tissue or benign nodules (x-axis).

[0063] Figure 2c The density UMAP map for each patient illustrates the distribution of cells in different tissue origins (N, B, T). The darker the color, the higher the cell density.

[0064] Figure 2d The scatter plot shows the frequency of somatic mutations in benign nodule and tumor samples from four representative patients (PTC12, PTC32, PTC30, PTC33). Each point represents a mutation, and its coordinates indicate the frequency in both tissue types. Dashed lines represent the frequency thresholds used for selection; red dots indicate common mutations, and gray dots indicate non-common mutations.

[0065] Figure 2e UMAP images of two representative patients (PTC12 and PTC30). Left: Cells stained by tissue origin. Middle: Cells stained by single-cell somatic mutations (scSNVs), distinguishing between shared and non-shared mutations. Right: Cells stained by mitochondrial mutations (mtDNA).

[0066] Figure 3a UMAP plot of all epithelial cells, cells are colored according to pseudo-time inferred from Monocle3, and black arrows indicate the trajectory direction.

[0067] Figure 3b UMAP diagram of all epithelial cells, with cells stained according to tissue origin.

[0068] Figure 3c The ridge map shows the CytoTRACE score distribution of patients grouped by sample type, including benign nodular cells, all normal cells, and all tumor cells.

[0069] Figure 3d UMAP plot of benign nodular cells, with EICs highlighted in orange and all other benign nodular cells shown in gray.

[0070] Figure 3e UMAP plots of all EICs depict their positions along the inferred trajectory of Monocle3.

[0071] Figure 3f --h. Box plots show CytoTRACE scores (f), BRAF scores (g), and TDS scores (h) from normal tissue, other benign nodular cells, EICs, and tumor cells. Data were stratified by patient, and p-values ​​indicate statistical significance between groups.

[0072] Figure 3i The bar chart shows the proportion of EICs in the total epithelial cell population for each patient.

[0073] Figure 3j-- Figure 3l UMAP and pseudotime analysis for three representative patients (PTC12, PTC30, PTC32). For each patient, the first UMAP plot shows cells stained by tissue origin (left), followed by cells stained by pseudotime (middle). The right panel shows the normalized proportion of each cell group (normal, other benign, EICs, tumor) at pseudotime. Highlighted points mark pseudotime intervals where the normalized proportion for a specific cell type exceeds 75%.

[0074] Figure 4a The UMAP plot shows the single-cell expression of three genes (ZCCHC12, CLDN1, and NPC2). Color intensity represents expression level.

[0075] Figure 4b UMAP plots of three representative patients (PTC12, PTC30, PTC32) showing the expression patterns of ZCCHC12, CLDN1, and NPC2.

[0076] Figure 4c Box plots were used to compare the expression of ZCCHC12, CLDN1, and NPC2 from different tissue origins in scRNA-seq data. P-values ​​indicate statistical significance between groups.

[0077] Figure 4d Box plots were used to compare the expression of ZCCHC12, CLDN1, and NPC2 from different tissue origins in SYSU bulk RNA-seq data. P-values ​​indicate statistical significance between groups.

[0078] Figure 4e The scatter plot illustrates the correlation between the average expression of ZCCHC12, CLDN1, and NPC2 and four indicators: malignancy, TDS score, EMT score, and BRAF score. Each point represents a sample, and a Loess curve is fitted to indicate the trend. R and P values ​​represent the correlation strength and statistical significance.

[0079] Figure 4fTop: Growth curves of four thyroid cancer cell lines (Nthy-ori 3-1, KTC-1, IHH-4, TPC-1) overexpressing ZCCHC12, CLDN1, or NPC2. The control group (Vector) was compared with the overexpression group to assess the functional effect on cell proliferation. Error bars represent standard deviation; significance is expressed as *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001. Bottom: Bar graphs show the colony-forming ability of the four thyroid cancer cell lines after overexpression of ZCCHC12, CLDN1, or NPC2. Colony counts quantified the effect on cell proliferation. P-values ​​indicate statistical significance compared to the control group.

[0080] Figure 4g Immunohistochemical (IHC) images show the protein expression of ZCCHC12, CLDN1, and NPC2 in thyroid tissue samples: solitary benign nodule (IB), matched benign nodule (MB), matched tumor (MT), primary benign nodule (OB), and primary tumor (PT). Scale bar: 100 μm (large image), 50 μm (inset).

[0081] Figure 4h :based on Figure 4g The resulting image is a bar graph showing the percentage of ZCCHC12, CLDN1, and NPC2 positive cells quantified by IHC analysis.

[0082] Figure 5a This study compares the classification performance of different gene feature combination models across four independent datasets. Receiver operating characteristic (ROC) curves were plotted using single-gene, dual-gene, and triple-gene combined models (ZCCHC12, CLDN1, and NPC2) to distinguish papillary thyroid carcinoma from normal samples, based on the SYSU bulk, TCGA, scBulk, and GSE33630 datasets.

[0083] Figure 5b The bar chart shows the AUC values ​​of single-gene, double-gene, and triple-gene combined models based on ZCCHC12, CLDN1, and NPC2 in distinguishing papillary thyroid carcinoma from normal samples in the SYSU bulk, TCGA, scBulk, and GSE33630 datasets. Detailed Implementation

[0084] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0085] Sample collection

[0086] We collected 96 surgical specimens from the First Affiliated Hospital of Sun Yat-sen University and Guangzhou Women and Children's Medical Center, including 36 patients with percutaneous thyroid carcinoma (PTC) and 5 patients with benign thyroid nodules, for single-cell RNA sequencing. These samples included 39 primary tumors, 29 adjacent normal tissues, 13 lymph nodes, and 15 benign nodules. In-depth WES analysis was also performed on samples from 9 patients with PTC and benign nodules.

[0087] Single-cell RNA sequencing

[0088] Cell preparation: Fresh specimens were processed immediately postoperatively. Tissue was washed with PBS, cut into 1–2 mm pieces, and digested at 37 °C for 15 min using a tumor digestion kit (Miltenyi Biotec, cat# 130-095-929) or trypsin (ThermoFisher, cat# 25300062). The cell suspension was filtered sequentially through 70 μm and 30 μm filters (Miltenyi Biotec, cat# 130-098-462, 130-098-458). Red blood cells were lysed and removed using RBC Lysis Buffer (eBioscience™, cat# 00430054) (centrifuged at 400 × g for 6 min). Cells were washed with PBS, counted using AO / PI, and adjusted to 500–1300 cells / μL for library construction. Cell suspensions were prepared according to the instructions of the Chromium Next GEM Single Cell 5' Reagent Kit (10x Genomics, v2chemistry).

[0089] Library construction and sequencing: GEMs (Gel Bead-In Emulsions) were generated in a droplet system for cell barcoding and reverse transcription, followed by cDNA amplification, quality control, and quantification. 5' gene expression libraries and V(D)J libraries were constructed according to the instructions (10x Genomics, Single Cell 5' Reagent Kits v5.2 User Guide). The libraries were then sequenced on a NovaSeq 6000 platform (Illumina, San Diego, CA).

[0090] bulk RNA sequencing

[0091] This study reanalyzed Bulk RNA-seq data from 216 PTC and 75 BTN patients, data previously generated and published by our team. Only the novel computational analysis methods are described here. Detailed methods for sample processing, library construction, sequencing, and clinical characterization can be found in the original literature.

[0092] Whole exome sequencing (WES)

[0093] WES was performed on primary tumors, benign nodules, and adjacent normal tissues from 9 patients with PTC. Genomic DNA was extracted from single-cell suspensions or frozen tissues stored at –80°C using the DNeasy Blood & Tissue Kit (Qiagen, cat# 69506). DNA concentration was determined using the Qubit dsDNA BR Assay Kit (Thermo Fisher, cat# Q32853), and integrity was verified by 1% agarose gel electrophoresis. 1–3 μg of DNA from each sample was fragmented to approximately 250 bp using a Biouptor Pico (Diagenode) sonication. Fragment size was determined using an Agilent 2100 Bioanalyzer. Library preparation was performed using the SureSelectXT Target Enrichment System (Agilent, cat# G9661B), and exon capture was performed using the SureSelectXT Human All Exon V6 kit (Agilent, cat# 5190-8863). After library construction, quantification was performed using Qubit, Agilent 2100, and qPCR (KAPA Library Quantification Kit, Roche, cat# KK4824). Finally, the library was sequenced at 150 bp paired ends on the NovaSeq 6000 (Illumina, San Diego, USA) platform (performed by Hangzhou Lianchuan Company).

[0094] Sequencing data analysis

[0095] Data preprocessing and quality control: Raw scRNA-seq sequencing data were processed using Cell Ranger software (v5.0.0) provided by 10X Genomics and aligned to the human reference genome GRCh38 (hg38). Detailed quality control metrics were then generated and evaluated, and high-quality cells were obtained through rigorous filtering for downstream analysis. The R package Seurat (v4.4.0) was primarily used for quality control and subsequent bioinformatics analysis. The filtering criteria were: removing cells with <200 detected genes and cells with a mitochondrial gene ratio >25%; cells with <3 detected genes were also removed.

[0096] Two-cell detection and batch effect correction: To ensure high-quality single-cell expression profiles, a multi-step strategy was employed to remove potential two-cell or multi-cell mixtures. First, cells with extremely low or high library complexity were removed: low-quality / dead cells with <200 genes or potential multi-cell mixtures with >6000 genes were detected. Second, based on cluster distribution and classic marker gene expression patterns, cells with abnormal co-expression (e.g., the epithelial marker gene EPCAM in T cell clusters or the myeloid marker gene CD68 in B cell clusters) were eliminated. These cells typically represent two-cell mixtures with mixed expression profiles and abnormally high gene counts. Further, the DoubletFinder algorithm was used to systematically identify and remove two-cell / multi-cell candidates. After two-cell removal, the filtered gene-cell matrix was processed using Seurat. UMI counts were normalized using NormalizeData (LogNormalize method, scale.factor=10,000) and transformed using log1p. To correct for technical and batch-to-batch differences, batch-to-batch correction was performed using Harmony in the PCA embedding space, and was performed independently in the major immune cell populations (B cells, myeloid cells, T cells, and endothelial cells) to ensure integration consistency between donors and experimental conditions.

[0097] Unsupervised clustering and dimensionality reduction analysis: After quality control and double cell removal, approximately 50,000 single cells were retained for downstream analysis. Gene expression values ​​were normalized using NormalizeData (LogNormalize method, scale.factor=10,000), and then high-variable genes (HVGs) were screened using FindVariableFeatures (selection.method="vst", top 2000). After centering and normalizing HVGs using ScaleData, principal component analysis (PCA) was performed using RunPCA. The number of principal components was determined using elbow plot and JackStraw analysis. Based on the selected PCs, a shared nearest neighbor (SNN) graph was constructed using FindNeighbors, and cell clustering was performed using FindClusters. The resolution parameter was optimized based on cluster stability and biological interpretability. Dimensionality reduction visualization was performed using UMAP (RunUMAP), combined with t-SNE for comparing cluster separation. To further characterize cellular heterogeneity, this workflow employs multiple iterative clustering methods to identify and annotate refined cell subpopulations. For the analysis of adjacent tumor epithelial cells, the Seurat data integration framework is used to reduce inter-patient heterogeneity. Specifically, SplitObject is first used to split and normalize single-cell data from different patients (NormalizeData), and then HVGs (FindVariableFeatures, top 2000) are identified separately. FindIntegrationAnchors (top 30 PCs) are used to find integration anchors, and IntegrateData is used to generate batch-corrected expression matrices (integrated assay). Subsequently, the data is normalized using ScaleData, dimensionality reduced using RunPCA (30 PCs), and clustered using UMAP (dims=1:30) and FindNeighbors / FindClusters (resolution=0.4) in the PCA space. This workflow ensures integration across patients while preserving biological differences. For epithelial cells from other tissues, Seurat's projection-based mapping method is used for annotation transfer. Epithelial cells from adjacent cancerous tissues were used as a reference dataset. The reference and query data were normalized and HVGs (top 2000) were identified.We used FindTransferAnchors (PCA dimensionality reduction of the reference data, using the first 30 dimensions, combined with k.filter to improve matching) to find transfer anchors. Then, we used TransferData (k.weight=50) to transfer the reference cell type labels to the query dataset and added them to the metadata. Finally, we used MapQuery to project the query cells into the reference UMAP space to ensure annotation consistency. The reference and query data were merged into a unified UMAP expression space, and cell identities were assigned according to the transferred labels, thus achieving unified annotation and comparison of epithelial cells from different sources.

[0098] Cell type annotation: Marker genes for each cluster were identified using Seurat's FindAllMarkers function (Wilcoxon rank-sum test, minimum percentage = 0.25, logfc threshold = 0.25). Cluster-specific markers were obtained using FindMarkers. The putative cell identity was determined by cross-referencing significant marker genes with established cell type markers in select databases (CellMarker, PanglaoDB) and relevant literature (especially pan-cancer scRNA sequencing atlases). Annotations were refined based on the expression of classic lineage-specific marker genes (e.g., EPCAM for epithelial cells, PTPRC (CD45) for immune cells, PECAM1 (CD31) for endothelial cells, ACTA2 (αSMA) for fibroblasts, and typical immune cell markers). To comprehensively define the functional characteristics of the epithelial cell procedure, three complementary approaches were employed:

[0099] (1) Metabolic activity inference: We applied the METAFlux framework to assess cellular metabolic status. First, response scores were calculated for each procedure, and then metabolic flux was calculated using the compute_sc_flux function under human blood culture conditions. This analysis was able to interpret the metabolic functions associated with each epithelial cell procedure.

[0100] (2) Functional enrichment analysis: Differentially expressed genes (DEGs) for each procedure were identified using Seurat's marker detection workflow. Subsequently, gene ontology (GO) enrichment analysis was performed on these procedure-specific gene sets using the clusterProfiler software package, focusing on multiple GO categories (biological processes, cellular components, and molecular functions). The enriched terms provided functional annotations for individual epithelial cell procedures.

[0101] (3) Knowledge-based reference dataset verification: To verify the program identity based on known epithelial cell states, we compiled a list of published reference genes annotating the function of epithelial cell programs. The differentially expressed genes (DEGs) identified in step (2) were compared with this set of reference genes using a hypergeometric test. The overlap between program-specific differentially expressed genes (DEGs) and the reference gene set allowed us to establish a correspondence between the identified programs and previously reported epithelial functional states.

[0102] Cross-validation for annotation consistency: Integrating the above analyses provides a robust and multi-faceted definition of epithelial procedures. Metabolic activity, enriched pathways, and reference-based matching are comprehensively considered to ensure a reliable assignment of biological function to each epithelial procedure.

[0103] Malignant lineage inference:

[0104] (1) scRNA-seq analysis: To investigate malignant progression within epithelial cells, we compared the gene expression profiles of malignant epithelial cells (derived from primary tumors, benign lesions, and lymph node metastases) with their matched normal epithelial cells. Differential expression analysis was performed on each sample using the Wilcoxon rank-sum test. Genes were retained if there was a significant difference between at least two independent samples (FDR ≤ 0.05 and |log2FC| ≥ 0.5). The union of these significant genes was then used to construct the log2FC matrix for all samples. The pseudo-temporal trajectory was reconstructed using the DDRTree algorithm implemented in Monocle3 (v1.0). This method projects the epithelial procedure onto a low-dimensional manifold, thereby capturing the underlying structure of malignant progression. For each sample, the nearest point on the trajectory was determined by minimizing the Euclidean distance between the sample and the fitted trajectory curve. The samples were then sorted along the inferred pseudo-temporal axis according to their trajectory positions, providing a relative measure of progression. By embedding principal components and smoothed trajectory curves into the DDRTree for visualization, it is possible not only to depict the overall progression trend, but also to locate the distribution of individual samples in the malignant lineage.

[0105] (2) Batch RNA Sequencing Analysis: For batch tissue samples, salient features identified at the single-cell level were used as input. Log2FC values ​​were calculated by comparing benign nodule or tumor samples with neighboring non-malignant tissues. Log2FC matrices were constructed for all batch samples, and the same DDRTree-based pseudo-time inference method as used for scRNA sequencing was applied to these batch lineages.

[0106] (3) Definition of trajectory root: To determine the starting point of the trajectory, we quantified the number of differentially expressed genes (DEGs) on the pseudo-time path. Samples with a higher number of differentially expressed genes (DEGs) were interpreted as having a greater difference from the normal epithelial state. Therefore, the sample with the lowest number of differentially expressed genes (DEGs) relative to normal tissue was selected as the trajectory root, representing the most "normal" state.

[0107] WES data analysis: Matched neighboring non-malignant tissues were used as references. Raw whole-exome sequencing (WES) reads (FASTQ) were aligned to the human reference genome GRCh38 using BWA-MEM (v0.7.17). Somatic single nucleotide variant (SNV) detection employed a consensus approach integrating results from VarScan2 (v2.3.9) and Strelka2 (v2.9.10) to improve the reliability of variant detection. Variants were screened to retain high-confidence somatic detection results. Specifically, for Strelka2-derived variants, only those marked "PASS" in the FILTER field were retained; and for VarScan2-derived variants, only those annotated "Somatic" were retained. At this stage, variants were not limited to exon regions, therefore both coding and non-coding SNVs were included in downstream analysis. Furthermore, only variants meeting the minimum sequencing depth (DP ≥ 10) and quality threshold (QUAL ≥ 20) were considered, and a variant required support from at least two independent detectors to be preserved. High-confidence SNVs were then functionally annotated using ANNOVAR (v105) to obtain gene-level annotations and predicted functional consequences, and variant recall was further visually validated using IntegrativeGenomics Viewer (IGV, v2.16.2).

[0108] Single-cell data single nucleotide variant detection: SComatic was used to identify somatic mutations in tumor and benign nodule cells at the single-cell level, with a minimum base quality threshold (min_bq) set to 30. To improve confidence, SNVs were filtered, retaining only high-confidence somatic variants, requiring somatic state = "somatic", p-value < 0.01, tumor variant allele frequency > 20%, and normal variant allele frequency ≤ 5%. Since sequencing employed a 5′-UTR targeting approach, downstream analysis considered only SNVs located within the 5′ untranslated region. This filtering strategy improved the reliability of mutation detection and reduced false positives, thus ensuring reliable identification of cell type-specific variants.

[0109] Mitochondrial DNA variant detection in single-cell data: Mitochondrial DNA (mtDNA) variants were identified from scRNA-seq data using the 10x Genomics mode of the mgatk software. To ensure high-confidence results, variants were further screened based on the following criteria: variants detected in at least two cells (n_cells_conf_detected ≥ 2), strand correlation ≥ 0.25, and log10 variance-to-mean ratio (vmr) > -2. This screening strategy robustly identifies mtDNA mutations while minimizing false positives due to technical noise.

[0110] Tumor epithelial cell origin inference: Patient-specific analyses were performed on datasets containing paired samples (normal adjacent tissue, benign nodules, and tumors). For each patient, cells from normal adjacent tissue and benign nodules were first merged and subjected to dimensionality reduction and clustering (see Section 1.3) to define a reference cell state. Tumor cells were then mapped onto a reference atlas constructed from paired normal and benign cells. Anchor points between tumor cells and reference cells were identified to guide label transfer, and tumor cells were projected into a reference low-dimensional space to assess their similarity to the reference population. Finally, tumor cells were assigned to categories corresponding to normal adjacent tissue or benign nodules based on the transferred labels. The proportion of tumor cells mapped to each reference category was quantified to measure relative similarity and infer the potential origin of tumor epithelial cells.

[0111] Identification of EICs: To explain the intrinsic heterogeneity of epithelial cells in benign nodules, we first used the SC3 clustering algorithm to subdivide benign nodules into subclones. For each patient, the number of SC3 clusters was set between 2 and 6, and the optimal number of clusters was selected based on the maximum profile index to ensure robust subclonal classification. Next, using the previously defined tumor mapping pattern, we assessed which benign nodule subclones exhibited potential for malignant transformation. Two complementary indices were calculated for each subclone: ​​(1) Similarity to tumor cells: For each benign subclone, the proportion of cells mapped to the tumor cell class was calculated. A higher proportion indicates a higher transcriptional similarity to tumor cells; (2) Distance to normal tissue: The Euclidean distance between the UMAP center of each benign subclone and the center of cells in the normal adjacent tissue was calculated. A larger distance indicates a greater deviation from normal epithelial features. Finally, subclones with high similarity to tumor cells and low similarity to normal tissue were classified as epithelial intermediate state cells (EICs), representing a transitional state that may be in the process of malignant transformation.

[0112] Identification of upregulated marker genes in EIC: We used a three-step strategy to identify significantly upregulated marker genes in epithelial intermediate state cells (EIC): (1) Differential expression in EIC and other benign nodule cells: Using the FindMarkers function, EIC was compared with non-EIC benign nodule cells to identify differentially expressed genes (DEGs). Genes with adjusted p < 0.001 and log2 fold change > 0.5 were defined as significantly upregulated, while genes with log2 fold change < -0.5 were considered downregulated. (2) Screening for epithelial-specific expression: To ensure epithelial specificity, we performed differential expression analysis on epithelial cells and other cell types. Genes were sorted according to fold change and adjusted p value, and genes with log2 fold change > 0.5 and adjusted p < 0.001 were retained as epithelial-specific genes. (3) Screening for tumor-enriched expression: Finally, we compared tumor epithelial cells (“T”) with normal adjacent epithelial cells (“P”). Genes with a log2 fold change > 0.5 and adjusted p < 0.001 were selected as tumor-enriched genes. Genes meeting all three criteria (EIC upregulation, epithelial specificity, and tumor enrichment) were defined as robust EIC marker genes.

[0113] Monocle3: To investigate potential lineage relationships among epithelial cells, Monocle3 (v1.4.26) was used to reconstruct pseudo-temporal trajectories for each patient, including paired samples of normal adjacent normal tissue, benign nodules, and tumors. For each patient, a Cell Dataset (CDS) object was constructed and preprocessed using 40 principal components. The `align_cds` function (k = 10) was used to mitigate batch effects between different tissue types. Dimensionality reduction and unsupervised clustering were then performed, and finally, the master graph was constructed using the `learn_graph` function (with default parameters). To define the trajectory root, a node in the normal cell population was interactively selected as the starting point. The `order_cells` function was then used to sort the cells along the trajectory. Pseudo-temporal values ​​were extracted and normalized to the [0,1] range for subsequent comparative analysis, quantitatively representing the inferred differentiation continuum among benign, normal neighboring cells, and tumor epithelial cells.

[0114] Statistical analysis

[0115] Unless otherwise specified, all statistical analyses were performed in the R statistical environment (v4.2.1). Between-group comparisons were performed using nonparametric or parametric methods, depending on the data distribution. Specifically, for single-cell data, the Wilcoxon rank-sum test was used to compare continuous variables (e.g., gene expression levels, label scores, cell proportions) between groups, while Student's t-test was used for normally distributed clump or pseudo-clump data. For comparisons of three or more groups, the Kruskal-Wallis H test was used. For paired analyses of matched patient samples, the Wilcoxon signed-rank test was used.

[0116] Correlation analysis: For nonparametric data, Spearman's rank correlation coefficient was used; for normally distributed variables, Pearson's correlation coefficient was used. Fisher's exact probability method was used to assess differences in classification frequencies between groups. Ordinal logistic regression analysis was performed using the `polr` function in the R package `MASS`.

[0117] Survival analyses were performed to assess the prognostic significance of molecular subtypes and gene expression characteristics in the TCGA-THCA sample. Clinical and survival data from 505 patients were downloaded from the TCGA Pan-Cancer Atlas via cBioPortal (https: / / www.cbioportal.org, Study ID: thca_tcga_pan_can_atlas_2018). Overall survival (OS) and disease-free survival (DFS) were assessed. Kaplan-Meier curves were generated to visualize survival probabilities, and the log-rank test was used to determine statistical differences between groups. All reported P-values ​​are two-sided. To correct for multiple hypothesis testing, the Benjamini-Hochberg method was used to control for false discovery rate (FDR), and significance was defined as P < 0.05 or FDR < 0.05.

[0118] For trend visualization and interaction analysis, local weighted scatter plot smoothing regression is used to fit smooth curves to the data.

[0119] For all high-dimensional analyses, such as differential gene expression, the Benjamini-Hochberg procedure was used to adjust p-values ​​to control for the false discovery rate (FDR). Unless otherwise specified, statistical significance was defined as p < 0.05 or FDR < 0.05. For extremely small p-values, results were reported as p < 2.2 × 10⁻¹. 6 All tests were two-tailed tests.

[0120] Immunohistochemical staining (IHC)

[0121] Formalin-fixed, paraffin-embedded (FFPE) tissue sections were dewaxed, hydrated, and subjected to antigen retrieval. Sections were blocked sequentially with 3% hydrogen peroxide (20 min, 20–25°C) and 5% BSA (30 min, 37°C). Primary antibodies were then added: anti-Claudin1 antibody (Cell Signaling Technology, 1:200, 13255S), anti-NPC2 antibody (Proteintech, 1:600, 19888-1-AP), and anti-ZCCHC12 antibody (Atlas, 1:60, ATL-HPA034940), and incubated overnight at 4°C. After washing, horseradish peroxidase (HRP)-conjugated secondary antibody (Long Island Antibody, 2406271) was added, and incubation was performed at room temperature for 30 min. Signal was developed using DAB and counterstained with hematoxylin, followed by dehydration, clearing, and mounting for imaging.

[0122] Cell experiments

[0123] Cell lines and cell culture: Thyroid cancer cell lines Nthy-ori 3-1, KTC-1, IHH-4, and TPC-1 were purchased from Procell™ (Wuhan, China). Nthy-ori 3-1, KTC-1, and IHH-4 cells were cultured in RPMI 1640 medium (Pricella, PM150110) with 10% fetal bovine serum (FBS, Pricella, 164210) and 1% penicillin-streptomycin (Gibco, 15140122); TPC-1 cells were cultured in DMEM medium (Pricella, PM150210) with the same supplements. All cell lines were identified by STR and cultured in a 37°C, 5% CO2 incubator. All cells were negative for mycoplasma according to routine tests.

[0124] Lentiviral transfection and cell line construction: Nthy-ori 3-1, KTC-1, IHH-4, and TPC-1 cells were infected with lentiviruses overexpressing Claudin1, NPC2, and ZCCHC12 (LV-Claudin1, LV-NPC2, LV-ZCCHC12) and an empty vector control (LV-control), respectively. Viruses were purchased from GeneChem (Shanghai, China). After selection with puromycin, the expression of the target genes was verified using qRT-PCR and Western blot.

[0125] RNA extraction, RT-PCR, and real-time quantitative qRT-PCR: Total RNA was extracted using TRIzol reagent (Invitrogen) and reverse transcribed into cDNA using PrimeScript™ RT Master Mix (Takara, RR036A). Amplification and detection of the target gene were performed using real-time quantitative PCR (qRT-PCR) with Taq Pro UniversalSYBR qPCR Master Mix (Vazyme, Q712-02), using β-actin as an internal control. Relative expression levels were calculated using the 2^(-ΔΔCt) method and normalized to ACTB.

[0126] Protein extraction and Western blot: Cells were lysed in RIPA lysis buffer containing 1% protease and phosphatase inhibitors on ice for 15 min, then centrifuged at 13,000 rpm, 4°C for 15 min to collect the supernatant. Protein concentration was determined using the BCA protein quantification kit (EpiZyme). Proteins were separated by SDS-PAGE and transferred to PVDF membranes, which were then blocked with 5% BSA. The membranes were incubated overnight at 4°C with primary antibodies: anti-Claudin1 (1:1000, CST, 13255S), anti-NPC2 (1:1000, Proteintech, 19888-1-AP), anti-ZCCHC12 (1:1000, Invitrogen, PA5-56615), and anti-β-actin (1:5000, CST, 4967S). HRP-labeled secondary antibody (CST, 7074S) was then added, and the cells were incubated at room temperature for 1 h. Chemiluminescence signals were detected using the Omni-ECL Femto Light Chemiluminescence Kit (EpiZyme, SQ201), with β-actin as an internal reference.

[0127] Cell proliferation assay: Cell proliferation was detected using CCK-8 reagent (MCE, HY-K0301). Thyroid cancer cells were seeded at 1×10³ cells / well in 96-well plates, and CCK-8 solution (1:10) was added every 24 hours for 2 hours of incubation. The absorbance at 450 nm was measured on days 1, 2, 3, 4, and 5 post-transfection.

[0128] Colony formation assay: 1 × 10³ cells were seeded into 6-well plates and cultured for 2 weeks. After colony formation, the cells were washed with PBS, fixed with 4% paraformaldehyde for 15 minutes, stained with 0.1% violet, rinsed with running water, and then observed.

[0129] Example 1: The malignant continuum of thyroid epithelial cells reveals the precancerous evolution trajectory in heterogeneous benign nodules.

[0130] To understand the potential for transformation from benign to malignant thyroid nodules, we performed single-cell RNA sequencing (scRNA-seq) on 96 samples from 41 patients, including normal tissue (N), benign nodules (B), primary tumors (PT) of papillary thyroid carcinoma (PTC), and lymph node metastases (LNM). Some samples also underwent deep whole-exome sequencing (WES). Furthermore, validation was performed using a separate bulk RNA-seq dataset containing 577 samples from 291 patients in another cohort. Figure 1a ).

[0131] After quality control, a total of 547,280 high-quality cells were retained, among which epithelial cells showed significant intertumor heterogeneity. After batch effect correction, we identified six epithelial subtypes ( Figure 1c Subtypes C4 and C5 were significantly enriched in PT and LNM samples. Figure 1c ), and were characterized by epithelial-mesenchymal transition (EMT) and interferon response features, respectively. Figure 1d ), exhibiting a low level of differentiation ( Figure 1e ).

[0132] Next, we compared the transcriptomic profiles of epithelial cells from B, PT, and LNM samples with those of normal epithelial cells to identify differentially expressed genes (DEGs). We aggregated epithelial cells from the same sample and performed principal component analysis (PCA) based on DEGs, ranking the samples according to their position relative to the spline curve fitted in the PCA space. This continuous trajectory reflects the cumulative transcriptomic alterations and the transition from normal to benign, to malignant, and then to metastatic states. Figure 1f And g). This malignant trajectory is associated with thyroid differentiation score (TDS), BRFA phenotypic characteristics, transcription factor activity, and metabolic pathways related to lipid metabolism reprogramming ( Figure 1h ).

[0133] Inference of the thyroid malignancy continuum based on bulk RNA-seq revealed two branches ( Figure 1i The findings from branch 1 were highly consistent with those from the scRNA-seq cohort, and correlated with the transition from normal to benign to cancerous states. Figure 1j However, branch 2 is enriched in normal and benign samples. Figure 1j This highlights the heterogeneity within benign nodules. The different locations of independent benign nodules (IB) and tumor-matched benign nodules (MB) along the malignant trajectory in the scRNA-seq cohort further support this heterogeneity. Figure 1k In summary, these results highlight the transcriptional diversity and malignant potential of heterogeneous benign nodules.

[0134] Example 2: PTC has a benign nodule origin.

[0135] To further investigate the cellular origin of PTC, we used paired normal and benign samples as references to perform etiological analysis on individual tumor epithelial cells (Figure 2a). The results revealed significant inter-patient heterogeneity (Figure 2b): one-third (3 / 9) of patients showed a benign-dominant origin, with some tumor cells traceable to benign precursors; another third showed a normal-dominant origin, with all tumor cells originating from normal epithelium; the remainder showed a mixed origin. UMAP visualization further validated this pattern: in benign-dominant cases, tumor cells and benign nodules highly overlapped, while in normal-dominant cases they were relatively independent (Figure 2c). Somatic mutation detection based on deep WES showed shared mutations between benign and tumor tissues in benign-dominant patients (Figure 2d). In two patients, these shared mutations were also confirmed at the single-cell level (Figure 2e), supporting the model that some PTCs can gradually evolve from benign nodules.

[0136] Example 3: Identification of precancerous epithelial intermediate-state cells in benign nodules

[0137] Based on scRNA-seq data, we reconstructed the state transition trajectory of epithelial cells (Fig. 3a, b). Epithelial cells of benign nodules were distributed in different branches, showing a high degree of heterogeneity in differentiation state (Fig. 3c).

[0138] In tumor-matched benign nodules (MBs), a special subset—interepithelial intermediate state cells (EICs)—was identified, whose transcriptomic characteristics were more closely similar to those of the matched primary tumor than to normal tissue (Fig. 3d). EICs stably oscillated between normal and tumor cells on their trajectory (Fig. 3e), exhibiting poor differentiation, higher EMT scores, and lower TDS scores (Fig. 3f–h).

[0139] EICs were identified in three patients with benign dominant or mixed origins (Fig. 3i). Personalized trajectory analysis further confirmed that EICs were stably distributed between normal and tumor cells, supporting their heterozygous phenotype and transformation potential (Fig. 3j–l).

[0140] Example 4: Pro-tumor gene characteristics of EIC

[0141] ZCCHC12, CLDN1, and NPC2 were significantly overexpressed in EICs (Fig. 4a), and their expression gradually increased from non-malignant cells to EICs and then to malignant cells in patients (Fig. 4b). They were also significantly upregulated in tumor tissues (Fig. 4c–d).

[0142] The expression levels of the three genes were positively correlated with malignant progression, negatively correlated with TDS, and closely related to BRAF phenotype and EMT signaling (Figure 4e), suggesting their potential role in tumorigenesis.

[0143] In in vitro experiments, CLDN1, NPC2, and ZCCHC12 were overexpressed in Nthy-ori 3-1, TPC1, KTC1, and IHH4 cells via lentiviral transfection, and their increased expression at both mRNA and protein levels was verified. The results showed that overexpression of these genes significantly enhanced cell proliferation and colony formation (Fig. 4f). Immunohistochemical analysis of FFPE samples from 100 PTC patients further confirmed that these three genes were significantly upregulated in PT and MB nodules compared to IB nodules (Fig. 4g, Fig. 4h). Furthermore, in three patients who initially underwent surgical resection of benign nodules, subsequently developed PTC, and underwent secondary surgery, paired benign and malignant tissue analysis showed that the expression of these three genes in their benign nodules was also higher than in IB nodules (Fig. 4g, Fig. 4h). These results suggest that the expression levels of ZCCHC12, CLDN1, and NPC2 are highly correlated with the malignant progression of thyroid epithelial cell lines.

[0144] Example 5: Biomarker-based classification model for papillary thyroid carcinoma

[0145] To evaluate whether the expression patterns of the three genes CLDN1, ZCCHC12, and NPC2 can be used to distinguish between positive and negative samples, we constructed a supervised classification model based on existing expression data and ground truth labels. Specifically, we used the Random Forest algorithm, with the expression levels of the three genes as input features and the PTC sample (positive) / normal thyroid sample (negative) status of the sample as the output label, thereby establishing a mapping relationship between gene expression levels and positive diagnosis.

[0146] The basic form of model learning can be represented as:

[0147]

[0148] in This represents the predicted probability that a sample belongs to the positive category. The model captures the nonlinear association between trigene expression and disease state through ensemble voting of a large number of decision trees.

[0149] We employed 10-fold cross-validation, repeated three times, to reduce bias from random sample partitioning. Specifically, the data was randomly divided into 10 subsets, with one subset used as the validation set and the rest as the training set each time, and this process was repeated three times. Within this validation framework, the random forest model learned the expression patterns of three genes based on the labeled raw data (including independent single-cell pseudobulk samples and 75% of the PTC bulk RNA-seq samples (SYSU bulk)), and formed stable and reliable classification rules under the constraint of multiple cross-validations.

[0150] The constructed model can not only establish the association between expression and diagnosis based on training set samples, but can also be directly applied to new patient samples. For any newly measured sample of three gene expression levels:

[0151]

[0152] When the predicted probability exceeds the preset threshold ( When the expression levels of the three genes are reached, the sample can be classified as highly probable positive. Therefore, this model constitutes a quantifiable and reproducible association framework, making positive predictions based on the expression levels of the three genes possible.

[0153] To compare with the three-gene joint model, we further constructed classification models using single genes and arbitrary two-gene combinations as features. Specifically, for each feature combination (single-gene, two-gene, and three-gene), we adopted the same modeling process: using labeled samples as training data, setting 10-fold cross-validation and repeating it three times to ensure the stability and comparability of the model evaluation. Within this unified framework, random forest models with different gene combinations all selected the optimal model parameters through the same hyperparameter search strategy and evaluated their classification performance under the same data partitioning and validation conditions. This strategy allows us to systematically evaluate the relative contributions of different gene combinations in positive discrimination and determine whether the three-gene joint model is superior to the single-gene or two-gene model, thereby clarifying the gain effect of the joint biomarker.

[0154] We systematically compared the classification performance of single-gene, dual-gene, and tri-gene combination models on different training sets and external validation sets. The training sets used were: (1) the SYSU bulk dataset mentioned above, which contains 137 normal thyroid samples and 137 PTC samples, of which 103 normal thyroid samples and 103 PTC samples were randomly obtained for training; (2) the pseudobulk dataset mentioned above, which contains 6 normal samples and 8 PTC samples. The validation sets included: (1) the TCGA thyroid cancer cohort, which contains 67 normal samples and 505 PTC samples; (2) another independent single-cell dataset, scBulk, which contains 23 normal samples and 31 PTC samples; (3) the GSE33630 dataset, which contains gene chip data of 49 PTC samples and 45 normal samples. Overall, the tri-gene combination model (NPC2 + ZCCHC12 + CLDN1) performed best in all datasets, with the highest AUC values ​​among all combinations (Figures 5a-b).

[0155] In the SYSU bulk dataset, the three-gene model achieved an AUC of 0.9906, significantly outperforming arbitrary two-gene combinations (AUC 0.9826–0.9901) and single-gene models (AUC 0.9639–0.9854). Figure 5a -b).

[0156] On the independent scBulk dataset, a consistent trend was observed: the AUC for the three-gene model was 0.9947, while that for the two-gene combination ranged from 0.9721 to 0.9912, and that for the single-gene model it was slightly lower (AUC 0.9478–0.9819). Figure 5a -b).

[0157] In the TCGA dataset, although the overall AUC is slightly lower than the aforementioned datasets, the three-gene model is still the best (AUC 0.8674), higher than the two-gene combination (AUC 0.8515–0.8638) and the single-gene model (AUC 0.8067–0.8419). Figure 5a -b).

[0158] In the external cohort GSE33630, the three-gene model also achieved the highest AUC (0.9739), further validating its robustness on different platforms. Figure 5a -b).

[0159] Overall, the three-gene combined model outperformed the single-gene and two-gene models in both the training set and multiple independent validation sets, indicating that the combined expression of the three genes can more comprehensively reflect the transcriptional differences between tumor tissues and normal tissues, thus providing a more stable and accurate classification capability. It can serve as a potential molecular tool for clinical thyroid nodule risk assessment and auxiliary diagnosis.

[0160] The present invention has been illustrated through the above embodiments, but the present invention is not limited to the above process steps, that is, it does not mean that the present invention must rely on the above process steps to be implemented. Those skilled in the art should understand that any improvements to the present invention, equivalent substitutions of the raw materials used in the present invention, additions of auxiliary components, and selection of specific methods, etc., all fall within the protection scope and disclosure scope of the present invention.

Claims

1. A combination of biomarkers for risk assessment of benign thyroid nodules, characterized in that, The markers include: CDLN1, NPC2 and ZCCHC12.

2. The marker combination according to claim 1, characterized in that, The markers consist of the following markers: CDLN1, NPC2, and ZCCHC12.

3. Use of CDLN1, NPC2 and ZCCHC12 in the preparation of reagents and / or kits for detecting epithelial intermediate state cells (EIC).

4. Use of CDLN1, NPC2 and ZCCHC12 in the preparation of reagents and / or kits for assessing the risk of thyroid nodules.

5. The use according to claim 4, characterized in that, The assessment of the risk of benign thyroid nodules includes distinguishing between thyroid cancer and healthy individuals.

6. The use according to claim 4, characterized in that, The assessment of thyroid nodule risk includes distinguishing between thyroid cancer and benign thyroid nodules.

7. The use according to any one of claims 5-6, characterized in that, The thyroid cancer mentioned is papillary thyroid carcinoma.

8. A kit for detecting epithelial intermediate state cells (EIC), characterized in that, The kit includes reagents for detecting the expression levels of CDLN1, NPC2, and ZCCHC12.

9. A kit for assessing the risk of benign thyroid nodules, characterized in that, The kit includes reagents for detecting the expression levels of CDLN1, NPC2, and ZCCHC12.

10. The kit according to any one of claims 8-9, characterized in that, The reagents include reagents for determining the sequences of CDLN1, NPC2, and ZCCHC12, primers for amplifying CDLN1, NPC2, and ZCCHC12, or reagents for detecting the protein expression levels of CDLN1, NPC2, and ZCCHC12.

Images