SMR analysis
We performed SMR analysis between PCD specific cis-eQTL and DTC GWAS from blood. Statistical significance threshold(P < 0.05) along with 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, Autophagy) (Table 1). Specifically, we observed that 3 genes associated with the development of DTC are related to the apoptosis mode. NFATC4(βblood = 2.89, Pblood= 0.024;βthyroid = 1.028, Pthyroid= 0.026)showed a positive correlation with DTC development, while RPS3༈βblood= -3.37253, Pblood= 0.0049༛βthyroid= -1.10, Pthyroid= 0.0026) and TM2D1 (βblood= -0.49, Pblood= 0.030༛βthyroid= -0.80, Pthyroid=0.036༉exhibited a negative correlation. Additionally, we identified 2 genes associated with the lysosome-dependent cell death pattern, CTNS (βblood= -0.67, Pblood= 0.022༛βthyroid= -0.42, Pthyroid= 0.023) and GCC2༈βblood= -0.25, Pblood= 0.025༛βthyroid= -0.37, Pthyroid= 0.023), which played a protective role in DTC. Finally, TPCN2༈βblood= -0.57, Pblood= 0.016༛βthyroid= -1.22, Pthyroid= 0.025), associated with the autophagy pattern, also played a protective role in thyroid carcinogenesis. Locuszoom plots demonstrated consistent genetic effects near these 6 genes in DTC GWAS (Supplementary Fig. 1–2).
Distribution of DTC associated PCD gene in peripheral blood mononuclear cells (PBMCs)
Based on the results of MR analysis above, we have identified DTC associated PCD genes in blood and further investigated its cell specifically expressed distribution in age-dependent accumulation single-cell transcriptomics data among peripheral blood. We obtained 125,035 scRNA-seq data from 17 normal individuals across four age intervals: umbilical cord blood (n = 3) and PBMCs of young (n = 3; age, 30.7 ± 10.0 years), healthy old (n = 6; age, 85.8 ± 11.1 years) and frail (n = 5; age, 88.0 ± 5.8 years) individuals. After quality control and batch effect correction, we performed UMAP-based cellular dimensionality reduction clustering analysis and identified 12 major cell types, including CD4+T cells, CD8+T cells, naive T cells, B cell, macrophages, cycling cells, NK cells, monocytes, granulocytes, platelets, plasma cells, Dendritic cells based on the expression levels of specific classical genes. We visualized the distribution of cell subgroups according to age-dependent accumulation (Fig. 2A). and investigated the distribution of identified PCD genes (Fig. 2B). The results suggested 48 of 55 PCD genes which exhibited association with DTC were enriched in PBMCs. Among them, 38 genes were found to exhibit weak expression (Percent Expressed < 50%) in peripheral blood (AKR1C1, ATP2A3, BBC3, BIRC6, CDK5R1, CHMP7, CLNK, CREB3, CTNS, CREB3, EEF1A2, FZD5, GCLC, ITGAM, KEAP1, LGMN, LIAS, NFATC4, PHLDA3, RHEB, SCARB2, S100A13, SFN, SGMS1, SMPD1, SNAPIN, STAT2, STAT5A, SYNPO2, TF, TFRC, TGFB2, TLR3, TNFRSF25, TNFSF10, TPCN2, TRIM6, USP30). While we observed a uniform distribution of 2 genes involved in cellular energy metabolism (ATP6V1G1, GAPDH) and the ribosome coding function of RPS3 across the majority of cell subtypes, which is expected given their biological functions. In addition, it was worth noting 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 gene in tumor microenvironment
To analyze the distribution of PCD genes susceptible to DTC in thyroid tissues, we analyzed scRNA-seq data from 15 PTC tissues in two datasets. 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 cell, T cell, NK cell, Endothelial, Epithelial, Fibroblast, Myeloid and Plasma cell (Fig. 3A, B, C). Besides, we explored the expression levels of tissue specific PCD genes using the “Dot plot” function. We visualized the distribution of cellular subgroups of PCD genes susceptible to DTC by UMAP (Fig. 3D) The results revealed that 13 PCD related genes were found in thyroid tissue cis-eQTL, 12 of which could be verified in tumor microenvironment. DDRGK1, GCC2, HEXA and TM2D1 were mainly expressed in epithelial cells of PTC. Besides, we still observed remarkably expression patterns such as SYK- myeloid cells. RPS3 genes were still expressed in all cell subgroups and upregulated in B cells, NK cells, and T cells. While the remaining susceptible genes such as CTNS, GPSM1, NFATC4, NOC2L, RIPK3, TPCN2 showed weak expression in tumor microenvironment.
Clinical correlation between PCD genes and PTC
To assess the clinical relevance of thyroid cancer associated PCD genes, we analyzed the clinical information and gene expression difference based on TCGA database. Considering age played an important role in the TNM stage of thyroid cancer, we assessed gender and age differences in the expression of PTC susceptibility PCD genes. Among them, the expression of RRPS3 and TM2D1 genes in elderly (> 65 years) patients was significantly lower than that in younger patients. Interestingly, based on SMR results that RRPS3 and TM2D1 acted as protective causal genes of PTC, therefore their inactivation mutations should be considered along with their carcinogenic susceptibility. However, these PCD associated PTC susceptibility genes did not show significant gender differences(Fig. 4). Subsequently, we further evaluated the susceptibility genes based on TNM (tumor node metastasis) classification(NCCN) staging. The expression of CTNS, GCC2, TM2D1, and TPCN2 genes gradually decreased with increasing T staging, which means certain inhibitory roles in the development of thyroid cancer (Fig. 5).
Expression and functional enrichment analysis of PCD genes in PTC
Firstly, based on the TCGA bulk sequencing data, we compared the expression levels of PTC susceptibility PCD genes with that in corresponding non-tumor normal tissues. it's worth noting that the expression levels of GCC2, TM2D1 and TPCN2 genes in tumor tissues were significantly lower than those in normal tissues, while the expression levels of CTNS and NFATC4 genes were significantly higher than those in normal tissues(Fig. 6A).
Subsequently, in order to further explore the biological functions of PGD-related susceptibility genes, we collected 38 phenotypic gene sets related to tumor invasion and metastasis, scored them by ssGSEA algorithm, and further evaluated their correlation with PGD-related A cancer susceptibility genes(Fig. 6C). Specifically, NFATC4 is significantly correlated with Basal Lamina, Collagen Family, ECM Receptor Interation, EMT to Metastasis and other phenotypes, suggesting that NFATC4 is involved in tumor epithelial mesenchymal transformation. Immune microenvironment remodeling and invasion and metastasis play a key role. RPS3 was significantly negatively correlated with Metastasis Response, VEGFA Signaling, HIF1A Signaling, EMT to Metastasis, Cell Cycle and Blood Coagulation. It is suggested that RPS3 may limit the potential of tumor growth and spread by inhibiting tumor neovascularization. And then inhibit the invasion and metastasis ability of tumor cells. The negative correlation between RPS3 and choline metabolism may reflect its potential regulatory effects on tumor energy supply and cell membrane synthesis. TM2D1 and Regulation of Angiogenesis, MMP Remodeling, Metastasis Response, LOX Remodeling, Hypoxia Response, HIF1A Signaling, Epithelial to Mesenchymal Transition, EMT to Metastasis, ECM Structure, Cellular Differentiation and Blood Coagulation pathway showed significant negative correlation. The association of TM2D1 with these pathways may reflect that it also plays a role in inhibiting tumor angiogenesis, inhibiting epithelial mesenchymal transformation, and limiting 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(Fig. 6B). 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 six 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(Fig. 6D).