SMR analysis
We performed SMR analysis between PCD specific cis-eQTL and DTC GWAS from blood. Conservatively statistical significance threshold (P<0.05) along with low heterogeneity test threshold (PHEIDI>0.01) were applied. The results suggested that there were 55 PCD genes in blood which exhibited association with DTC (Supplementary Table 4). To obtain more precise and robust MR analysis results, we further selected thyroid tissue-specific eQTL for SMR analysis, in which we ensured the same set of screening criteria. The results suggested 12 thyroid-tissue-associated genes closely related to outcome (Supplementary Table 5). By combining these results, we gave priority to presumption 6 PCD genes (CTNS, GCC2, NFATC4, RPS3, TM2D1, TPCN2) which could replicate in both blood and thyroid tissues robustly, involving 3 different PCD modes: Lysosome-dependent cell death, Apoptosis and Autophagy (Table 1). Specifically, we observed 3 genes related to Apoptosis had causal effects with DTC: NFATC4 showed risk effect in the development of DTC (βblood= 2.90,Pblood= 2.41E-02;βthyroid= 1.03,Pthyroid= 2.66E-02), while RPS3 (βblood= -3.37,Pblood=4.91E-03;βthyroid= -1.10,Pthyroid= 2.70E-03) and TM2D1(βblood= -0.49,Pblood=3.01E-02; βthyroid=-0.80, Pthyroid=3.69E-02) exhibited protective effects. Additionally, we identified 2 genes associated with the Lysosome-dependent cell death pattern played protective roles, CTNS (βblood= -0.67,Pblood= 2.23E-02;βthyroid= -0.42,Pthyroid= 2.35E-02) and GCC2 (βblood= -0.25,Pblood= 2.52E-02; βthyroid= -0.37,Pthyroid= 2.33E-02). Finally, TPCN2 (βblood= -0.57,Pblood=1.61E-02; βthyroid= -1.22,Pthyroid= 2.55E-02), associated with the Autophagy pattern, also played a protective effect. Locuszoom plots demonstrated consistent genetic effects near these 6 genes in DTC GWAS (Supplementary figure 1-2).
Distribution of DTC associated PCD genes in peripheral blood mononuclear cells (PBMCs)
Based on the results of SMR analysis in blood, we investigated the cell-specific expression in age-dependent single-cell transcriptomics data from peripheral blood. After quality control and batch effect correction, we performed UMAP-based cellular dimensionality reduction and clustering analysis, identifying 12 major cell types: CD4+ T cells, CD8+ T cells, naive T cells, B cells, macrophages, cycling cells, NK cells, monocytes, granulocytes, platelets, plasma cells, and dendritic cells, based on the expression levels of specific classical genes. We visualized the distribution of cell subgroups according to age-dependent accumulation (Figure 2A) and investigated the distribution of identified PCD genes (Figure 2B).
The results suggested that 48 of the 55 PCD genes associated with DTC were enriched in PBMCs. Among these, 38 genes were found to exhibit weak expression (Percent Expressed < 50%) in peripheral blood. Conversely, we observed the uniform distribution of 2 genes involved in cellular energy metabolism (ATP6V1G1, GAPDH) and 1 gene involved in ribosome coding function (RPS3) across the majority of cell subtypes, which is expected given their biological functions. Additionally, it was notable that there was age-specific cell expression of several susceptibility genes.In particular, genes such as FYN, GCC2, HIF1A, MYH9, PLAUR and YWHAQ were decreased in plasma cells among youth group. And TM2D1 expressed higher in dendritic cells of fail group.
Distribution of DTC associated PCD genes in tumor microenvironment
After stringent quality control, a total of 84,822 cells were included. We performed dimensionality reduction clustering and visualization using UMAP and annotated 8 cell subtypes using classical cell type markers, including B cells, T cells, NK cells, endothelial cells, epithelial cells, fibroblasts, myeloid cells, and plasma cells (Figure 3A, B, C). We also explored the expression levels of tissue-specific PCD genes using the “Dot plot” function and visualized their distribution with UMAP (Figure 3D).
The results revealed that all 13 thyroid-specific genes could be identified, with 12 verified in the tumor microenvironment. DDRGK1, GCC2, HEXA, and TM2D1 were mainly expressed in the epithelial cells of DTC. Additionally, we observed notable expression patterns such as SYK in myeloid cells. RPS3 genes were expressed in all cell subgroups and upregulated in B cells, NK cells, and T cells. In contrast, the remaining susceptible genes, such as CTNS, GPSM1, NFATC4, NOC2L, RIPK3, and TPCN2, showed weak expression in the tumor microenvironment.
Clinical correlation between PCD genes and DTC
Considering age played an important role in grading thyroid cancer, we assessed the differential expression of 6 hub genes, which consistently replicated in both blood and thyroid tissues. Among them, the expression of RPS3 in the old (>65 years) was significantly lower (P=0.022) (Figure 4A). Interestingly, the evidence from SMR showed RPS3 acted as protective genes in the pathogenesis of DTC. Therefore, the inactivation mutation of RPS3 should be considered along with carcinogenicity. However, all of these genes did not show significant gender differences (P>0.05) (Figure 4B). Subsequently, based on TNM staging, we found the expression of CTNS, GCC2, TM2D1, and TPCN2 gradually decreased with T stage increasing, which suggested certain inhibitory roles in the development of DTC (Figure 5).
Expression and functional enrichment analysis of PCD genes in DTC
Firstly, based on the TCGA bulk sequencing data, we compared the differential expression in tumor and normal tissues. Remarkably, the expression levels of GCC2, TM2D1 and TPCN2 in tumor tissues were significantly lower, while the expression levels of CTNS and NFATC4 were significantly higher(Figure 6 A).
Subsequently, in order to further explore the biological functions,of these hub genes, we collected 38 kinds of functional pathways related to tumor invasion and metastasis, scored them by ssGSEA algorithm (Figure 6C). We found NFATC4 was significantly correlated with inducing epithelial mesenchymal transformation to promote the metastasis. RPS3 was significantly negatively correlated with metastasis response, VEGFA signaling, HIF1A signaling, EMT to metastasis, cell cycle and blood coagulation, which suggested that RPS3 might limit the potential of tumor growth and spread by inhibiting tumor neovascularization. And the negative correlation between RPS3 and choline metabolism might reflect its potential regulatory effects on tumor energy supplement and cell membrane synthesis. TM2D1 showed significant negative correlation with the regulation of angiogenesis, reflecting the role in inhibiting tumor angiogenesis, epithelial mesenchymal transformation, and tumor cell invasion and metastasis. CTNS, GCC2 and TPCN2 all showed a certain correlation with tumor matrix remodeling and tumor angiogenesis phenotype, suggesting a key role in the progression and metastasis of thyroid cancer.
Immune infltration and immune checkpoints analysis
Subsequently, we evaluated the association of PCD susceptibility genes with nine immune checkpoints and 25 human leukocyte antigen (HLA) molecules to explore their potential role in tumor immunotherapy(Figure 6 B). It is worth noting that TM2D1, GCC2 and TPCN2 are significantly negatively correlated with a variety of HLA molecules, suggesting their potential role in tumor neoantigen recognition and antigen presentation. In terms of immune checkpoint, RPS3, TM2D1, CTNS, GCC2 and TPCN2 genes were significantly positively correlated with OLA1, and NFACT4, RPS3, TM2D1, CTNS, GCC2 were significantly positively correlated with ATIC, suggesting that the targeted screening susceptibility genes may play a potential role in immunocompromised but blocking therapy. We then evaluated the correlation between PCD susceptibility genes and immune infiltration, and the results showed that there was a significant positive correlation between CD56dim natural killer cells except NFATC4. The 6 genes showed significant negative correlation with a variety of immune cells, suggesting that the screened molecules play a potential role in the formation of the "cold tumor" microenvironment. Targeting potential molecules may contribute to a "cold" and "hot" transition in the nature of the tumor microenvironment(Figure 6 D).