Immune Cell Confrontation in the Papillary Thyroid Carcinoma Microenvironment

Background Papillary thyroid cancer has been associated with chronic inflammation. A systematic understanding of immune cell infiltration in PTC is essential for subsequent immune research and new diagnostic and therapeutic strategies. Methods Three different algorithms, single-sample gene set enrichment analysis (ssGSEA), immune cell marker and CIBERSORT, were used to evaluate immune cell infiltration levels (abundance and proportion) in 10 data sets (The Cancer Genome Atlas [TCGA], GSE3467, GSE3678, GSE5364, GSE27155, GSE33630, GSE50901, GSE53157, GSE58545, and GSE60542; a total of 799 PTC and 194 normal thyroid samples). Consensus unsupervised clustering divided PTC patients into low-immunity and high-immunity groups. Weighted gene coexpression network analysis (WGCNA) and gene set enrichment analysis (GSEA) were used to analyze the potential mechanisms causing differences in the immune response. Results Compared with normal tissues, PTC tissues had a higher overall immune level and higher abundance levels and proportions of M2 macrophages, Tregs, monocytes, neutrophils, dendritic cells (DCs), mast cells (MCs), and M0 macrophages. Compared with early PTC, advanced PTC showed higher immune infiltration and higher abundance levels and proportions of M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages. Compared to the low-immunity group, the high-immunity group exhibited more advanced stages, larger tumor sizes, greater lymph node metastases, higher tall-cell PTCs, lower follicular PTC proportions, more BRAF mutations, and fewer RAS mutations. Epstein-Barr virus (EBV) infection was the most significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for key module genes. Conclusions In human PTC, M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages appear to play a tumor-promoting role, while M1 macrophages, CD8+ T cells, B cells, NK cells, and T follicular helper (TFH) cells (including eosinophils, γδ T cells, and Th17 cells with weak supporting evidence) appear to play an antitumor role. During the occurrence and development of PTC, the overall immune level was increased, and the abundance and proportion of tumor-promoting immune cells were significantly increased, indicating that immune escape had been aggravated. Finally, we speculate that EBV may play an important role in changing the immune microenvironment of PTC tumors.


INTRODUCTION
Thyroid cancer (TC) is the fifth most common cancer in women in the USA, with an estimated 52,890 new TC cases nationwide in 2020 (1). The incidence of TC has been increasing in recent decades, and showing a trend of rejuvenation (2)(3)(4)(5)(6)(7). TC also has the highest rate of growth among diagnosed tumors in the USA (8). Papillary thyroid carcinoma (PTC) belongs to differentiated thyroid carcinoma (DTC) and accounts for approximately 80-85% of all TCs (2,9). Although the vast majority of PTCs can be cured by surgery and 131 I treatment, 10% of patients still die from poorly differentiated and advanced tumors (10). A study with a median follow-up period of 27 years also showed that 28% of PTC patients relapsed and 9% died of PTC (11). The reason is partly because the interaction and collective effect of tumors with the surrounding immune microenvironment during the occurrence and development of the tumor promotes immune tolerance and then progresses to immune escape, resulting in the body's inability to completely clear the tumor.
The immune editing hypothesis was proposed by Dunn in 2002 and divided the interaction between the immune system and the tumor into three stages: (1) elimination: tumor cells are cleared by the immune system before a clinical diagnosis is established. This stage is also called immune surveillance; (2) equilibrium: tumor cells with less immunogenicity gradually replace tumor cells with more immunogenicity; and (3) escape: tumor cells eventually escape the body's immune surveillance (12). According to this hypothesis, the immune editing of tumors also sculpts the tumor's immune phenotype while clearing the tumor, enabling the host to generate immune activity against the tumor. When the tumor is clinically diagnosed, it can not only tolerate the immune response but also manipulate the immune system to promote tumor progression. Poschke et al. (13) believe that tumors have two different methods of immune escape: camouflage and sabotage. Camouflage refers to the malformation and loss of the major histocompatibility complex class I (MHC-I) molecules on the surfaces of tumors, allowing tumors to escape the detection of the immune system. Sabotage refers to the ability of some tumors, including PTC, to manipulate part of the immune system to fight against the body's immune response to protect themselves. In this process, tumors attract and even mediate some immune cells, such as tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs), Tregs, and mast cells (MCs). In general, MDSCs, Tregs, and MCs respond to uncontrolled inflammatory reactions in the body, but in tumor progression, they have created an immune microenvironment that allows this process. Existing evidence shows that Tregs, neutrophils, dendritic cells (DCs), MCs, and TAMs play a tumorigenic role in the PTC tumor microenvironment (TME), while CD8+ T cells, B cells, and NK cells play a protective role (14)(15)(16). However, the infiltration and role of some major immune cells in PTC are still controversial. Studies on nonprimary immune cells are scarce, and systematic research on the PTC immune microenvironment remains to be conducted. Therefore, accurate assessment of the specific pattern of immune cells in PTC, including not only its phenotype but also its function, is essential for subsequent immune research and new diagnostic and therapeutic strategies.
The purpose of this study was to systematically understand immune cell infiltration in the PTC microenvironment and to summarize the changes in the immune microenvironment during the occurrence and development of PTC. Our research found that M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages play a tumor-promoting role, while M1 macrophages, CD8+ T cells, B cells, NK cells, and T follicular helper (T FH ) cells play an antitumor role. During the occurrence and development of PTC, the overall immune level increased, and the abundance and proportion of tumor-promoting immune cells increased significantly. Our research will help researchers systematically understand PTC immune cell infiltration and provide a reference for subsequent basic and clinical research related to PTC immunity.

Materials
The RNA-seq data (level-three, FPKM) of the thyroid carcinoma (THCA) dataset were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and included 58 normal thyroid samples (N) and 512 PTC samples (T). The clinical data of THCA were downloaded from the University of California at Santa Cruz (UCSC) Xena platform (https://xena.ucsc.edu/).
The other nine PTC microarray datasets were downloaded from the Gene Expression Omnibus (GEO) database (http:// www.ncbi.nlm.nih.gov/geo): GSE3467 (N = 9, T = 9), GSE3678 (N = 7, T = 7), GSE5364 (N = 16, T =  Single-sample gene set enrichment analysis (ssGSEA) The ssGSEA score was used to quantify the enrichment level of 29 immune signatures in each THCA sample (17). To understand the overall situation of immune infiltration and facilitate visualization, we corrected the data as follows. For each ssGSEA score, Xi was replaced with Xi ′ with the equation Xi ′ = (Xi −Xmin)/(Xmax − Xmin), where Xmin and Xmax represent the minimum and maximum values of the ssGSEA score in the THCA sample, respectively (18). The same method was also applied to the other nine GEO datasets.
Additionally, we also measured the infiltration abundance of THCA immune cells through immune cell-specific markers (19). Given that its algorithm and results are highly similar to ssGSEA and are also measurements of immune cell abundance, we present the results in the form of a supplementary graph.

CIBERSORT
CIBERSORT is a deconvolution algorithm based on RNA-seq data to calculate the composition ratio of 22 immune cells in each sample (20). We performed 1,000 permutations and retained samples of p < 0.05 to ensure the reliability of the results, and the sum of various immune cells was 1. Considering the small sample size of each data set except TCGA after deconvolution, the results of nine GEO data sets were integrated into one validation set.

Correlation and Clustering of Immune Cells
The Spearman correlation between each pair of cells was calculated and visualized using the R package "corrplot". Hierarchical agglomerative clustering based on Ward's linkage and Euclidean distance was used to group immune cells in the PTC microenvironment. An unsupervised clustering algorithm (K-means) was applied to identify immune infiltration patterns and classify patients for further analysis. The results were visualized in the form of an immune network (21).

Consensus Clustering of PTC Patients based on Immune Signatures
Based on the ssGSEA scores of 29 immune signatures, the R package "ConsensusClusterPlus" was used to perform unsupervised clustering of the THCA samples (22). The cumulative distribution function (CDF), consensus heat map, and principal component analysis (PCA) (23) were applied to evaluate the optimal K (number of groups). We renamed the two groups to the low-immunity group (L-immunity) and to highimmunity group (H-immunity) based on the differences in the ssGSEA scores between the two groups after clustering.

ESTIMATE
ESTIMATE was used to evaluate the composition of the TME for each THCA sample (24) and included immune score (immune cell infiltration), stromal score (stromal content), ESTIMATE score (stromal-immune comprehensive score) and tumor purity.

WGCNA Reveals Key Modules and Hub
Genes Related to PTC Immunity First, differentially expressed genes (DEGs) between PTC and adjacent normal tissues in TCGA dataset were identified with the R package "limma" (|log 2 FC| > 1 and adjusted p-value < 0.05).
WGCNA aims to identify coexpressed gene modules, explore associations between gene networks and traits of interest, and discover hub genes in networks (25). DEGs were divided into different gene modules by the R package "WGCNA" (minModuleSize = 30).
The module with the highest absolute module significance was defined as the key module. A gene that satisfies the following conditions was defined as a hub gene. Analysis of the main mechanism of PTC immunity

Functional Enrichment Analysis
The genes in the key module were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses with DAVID 6.8 (https://david.ncifcrf.gov/).

Protein-Protein Interaction (PPI) Network Analysis
The hub genes were input into the STRING database (v11.0) (https://string-db.org/) to analyze the interaction network of the hub gene-encoded proteins. The minimum required interaction score was set to 0.4 (medium confidence), and the results were visualized by Cytoscape (v3.7.2).

Statistical Analysis
The t-test and Mann-Whitney U test were used for comparisons between two groups. The chi-square test was used to assess differences in clinical parameters between the L-immunity and H-immunity groups. The Pearson and Spearman method was used for correlation analysis. The log-rank method was used to calculate the significant survival p-values (progression-free survival (PFS) as the ending indicator). R software (version 3.6.1) and SPSS (v25.0) software were used to analyze the data. R, GraphPad Prism (v8.0) software and Excel were used for data visualization. p < 0.05 was considered significant.

RESULTS
The Immune System in PTC Showed Overall Enhancement ssGSEA evaluated 29 immune signatures between PTC and adjacent normal tissues in TCGA ( Figure 1A), and the results showed that Tregs, Th1 cells, T helper cells, plasmacytoid DCs (pDCs),immmature DCs (iDCs), DCs, activated DCs (aDCs), neutrophils, MCs, and macrophages increased in PTC (p < 0.05), while CD8+ T cells decreased in PTC. There were no significant differences in TILs, T FH cells, NK cells, or B cells. Additionally, the type II IFN response, the type I IFN response, parainflammation, MHC class I, HLA, checkpoint, chemokine receptor (CCR) and APC costimulation were elevated in PTC, indicating overall enhancement of the immune system in PTC.

PTPRC (memory T cells), MS4A1 (B cells), IL17A (Th17 cells), CXCR5 (T FH cells), and CD8A (cytotoxic T cells) decreased in PTC (
To summarize the immune cell infiltration results (abundance) of ssGSEA and immune cell markers: 1. Compared with the immune system in normal tissues, the immune system in PTC appears to be enhanced overall. 2. Compared with normal tissues, PTC tissues have an increase in the number of tumor-promoting immune cells, which is particularly significant. 3. Compared with normal tissues, the two algorithms have slightly different evaluations of antitumor immune cells in PTC. The ssGSEA results show that no significant differences in antitumor immune cells, while the immune cell marker results show a downward trend for antitumor immune cell abundance in PTC compared with normal tissues.

Increased Proportion of Tumor-Promoting Immune Cells in PTC
CIBERSORT was used to calculate the proportion of each of the 22 immune cell types based on RNA-seq data. We performed measurements on the TCGA and GEO datasets (Figures 2A-D).
The GEO dataset incorporates valid samples from nine datasets, including GSE3467, GSE3678, GSE5364, GSE27155, GSE33630, GSE50901, GSE53157, GSE58545, and GSE60542. We compared the proportions of immune cells between the PTC and normal tissues in both TCGA and GEO datasets. In TCGA, normal tissues had higher levels of naive B cells, memory B cells, CD8+ T cells, T FH cells, M1 macrophages, activated MCs, and eosinophils than PTC tissues, while PTC tissues had higher levels of M0 macrophages, M2 macrophages, resting DCs, activated DCs, and resting MCs than normal tissues ( Figure 2E). In GEO, normal tissues had higher levels of naive B cells, plasma cells, CD8+ T cells, T FH cells, T cells gamma delta, M1 macrophages, and activated MCs than PTC tissues, while PTC tissues had higher levels of Tregs, activated NK cells, monocytes, M2 macrophages, resting DCs, activated DCs, resting MCs, and neutrophils than normal tissues ( Figure 2F). To assist with the understanding of these results, we listed the results in a balance chart ( Figure 2G). The construction of the balance chart mainly refers to the reviews by Varricchi et al. (14), Galdiero et al. (15) and Ferrari et al. (16) and the literature table that we created ( Table 1).
In summary, compared with normal thyroid tissues, PTC tissues have a higher proportion of tumor-promoting immune cells and a lower proportion of antitumor immune cells.

Immune Cell Infiltration Increased During PTC Progression
ssGSEA was used to analyze the relationship between PTC progression and immune cell infiltration in TCGA. The results showed that patients with high tumor stage, lymph node metastasis, and larger tumor size had an increased number of immune cells, and the increase in tumor-promoting immune cells was particularly significant, including Th1 cells, macrophages, Tregs, monocytes, neutrophils, iDCs and MCs ( Figures 3A-C). Patients with tumor metastases had lower immunity levels than patients without metastases ( Figure 3D).
In addition, the results of immune cell markers also showed that PTC patients with a high tumor stage, lymph node metastasis, and a larger tumor size had an increased number of immune cells, and the increase in tumor-promoting immune cells was particularly significant, including TPSAB1 (MCs), FOXP3 (Tregs), CD1A (iDCs), and neutrophils (CXCR2s) (Figures S1B-D).
In summary, with the progression of PTC, immune cell infiltration increased comprehensively.
The Proportion of Tumor-Promoting Immune Cells Increased During PTC Progression CIBERSORT analyzed the relationship between PTC progression and immune cell infiltration in TCGA. The results show that T FH cells, CD8+ T cells, plasma cells, and M1 macrophages have a higher proportion in stages I and II than in other stages. Resting memory CD4+ T cells, neutrophils, monocytes, resting MCs, M2 macrophages, M0 macrophages, resting DCs, and activated DCs were higher in stages III and IV ( Figure 4A). CD8+ T cells were higher in N0 PTC, and T cells gamma delta, resting memory CD4+ T cells, resting DCs, and activated DCs were higher in N1 PTC ( Figure 4B). T FH cells, CD8+ T cells, plasma cells, M1 macrophages, and memory B cells were higher in T1 + T2 PTC. Resting memory CD4+ T cells, neutrophils, monocytes, M2 macrophages, M0 macrophages, resting DCs, and activated DCs were higher in T3 + T4 PTC ( Figure 4C). These results are further shown in a balance chart ( Figure 4D).
In summary, compared with early PTC, advanced PTC has a higher proportion of tumor-promoting immune cells and a lower proportion of antitumor immune cells.

Correlation and Clustering of Immune Cells
We performed correlation analysis on immune cells based on ssGSEA (abundance) ( Figure 5A) and CIBERSORT (proportion) ( Figure 5B). The ssGSEA results showed a high positive correlation among all immune cells. The CIBERSORT results showed a positive correlation among antitumor immune cells (T FH cells, M1 macrophages, plasma cells, and CD8+ T cells), and these antitumor immune cells were negatively correlated with immune cells with tumor-promoting effects (M0 macrophages, M2 macrophages, Tregs, activated DCs, and resting MCs).
We also clustered immune cells from the perspective of abundance ( Figure 5C) and proportion ( Figure 5D) and examined the effect of each immune cell type on prognosis. CIBERSORT results in immune cells (T FH cells, CD8+ T cells, activated memory CD4+ T cells, plasma cells, M1 macrophages, naive B cells) clustered in cell cluster A. In addition, except for M0 macrophages (proportion) (HR = 11.8, log rank p = 0.032), the effects of other immune cells on PFS did not reach statistical significance.

Consensus Clustering of PTC Patients Based on Immune Signatures
Consensus unsupervised clustering was performed on 502 PTC patients based on 29 immune signatures of ssGSEA. We • Macrophages are present in the immune environment of PTC and are significantly elevated in the BRAF V600E+ group. (27) • The density of macrophages in PTC tissues is significantly higher than that in normal tissues. (28,29), • Compared with thyroid goiter and follicular adenoma, PTC tumors have significantly higher TAM densities.
(30) • The expression of DCs MCs and macrophages in follicular variant of PTC (FVPTC) is higher than that in adenomas and may be involved in tumor development and invasion.
• TAM density is positively correlated with PTC stage, tumor size, lymph node metastasis, and distant metastases.
• Tregs and pDCs together help tumor escape in patients with PTC plus nontoxic goiter (MNG). (38) • The expression of DCs, MCs, and macrophages in FVPTC is higher than that in adenomas and may be involved in tumor development and invasion.
• High CD1a+ DC density is associated with improved disease-free survival (DFS) of PTC. (39) • There is no significant correlation between DC density and DFS.
(36) • DC immunotherapy can be used in patients with papillary or follicular TC without significant side effects.
(40) MCs • MCs are recruited into PTC and promote the proliferation, survival and invasion of cancer cells, thereby promoting the growth and aggressiveness of PTC.
• Increased MC density in TC correlates with increased aggressiveness. (42) Neutrophils • NLR is not significantly increased in patients with PTC.
(44) • TC cells can recruit neutrophils by releasing CXCL8/IL-8. Additionally, in human TC samples, neutrophil density is related to tumor size and has a potential tumor-promoting effect.
• Elevated neutrophils in PTC can predict bilaterality and lymph node metastasis. (46) • Elevated NLR does not appear to be a reliable indicator of DTC progression.
(47) • With the increase in the NLR value of PTC, the histopathological characteristics become worse and the clinical behavior becomes more aggressive. (48-49) • There is no significant association between preoperative NLR and prognostic factors in patients with PTC. (48,50), • Elevated NLR is associated with a high risk of relapse. After treatment, those patients with low stage and good prognosis of DTC were observed to have a significant decrease in NLR. (51-53), • Higher NLR is associated with higher thyroglobulin levels in DTC, which indicates poorer survival. (54) NK cells • Compared with goiter and healthy thyroid, tumor-infiltrating NK cells increase in PTC, while no differences are found in peripheral blood NK cells. (55-57), • In PTC patients, NK cell infiltration is inversely related to the stage of the disease and plays a role of immune surveillance center by killing cancer cells.
• In human PTC, lymphocyte density is associated with improved overall survival and reduced relapse. (65,66), • In terms of CD4+ cell frequency, no difference was found between PTC and MNG patients.
(38) • In DTC, the extent of tumor-infiltrating CD4+ cells does not seem to predict the patient's prognosis.
(64) • Compared with patients with thyroid adenoma, PTC patients have significantly increased CD4+ CD25+ regulatory T cells in the peripheral blood.
• Increased Treg infiltration is positively correlated with advanced PTC. (58-68) • FoxP3 expression in PTC is associated with extrathyroidal invasion and distant metastasis but has nothing to do with overall survival.
• Compared with healthy controls, Th17 levels in the peripheral blood and tissue samples of patients with PTC increased, and Th17 cells in the peripheral blood were negatively correlated with tumor size.
The above table contains the relationship between immune cells and the occurrence, development and prognosis of PTC. Some DTC studies with PTC as the main sample were also included.
considered the relative change in area under the CDF curve at K = 2 to 9 ( Figures 6A, B) and found that the consensus heatmaps showed stable partitions when the samples were clustered into two groups (k = 2) ( Figure 6C). PCA of all genes revealed significant differences between the two groups ( Figure 6D). We renamed the two groups to the H-immunity (n = 293) and L-immunity groups (n = 209) based on the differences in immune signatures ( Figure 6E). The composition of the TME of PTC samples was analyzed by ESTIMATE. The results showed that the H-immunity group had higher immune scores ( Figure 6F), stromal scores ( Figure 6G), and ESTIMATE scores ( Figure 6H) and lower tumor purity ( Figure 6I ) tha n th e Limmunity group.
In summary, we divided PTC patients into L-immunity and H-immunity groups by consensus clustering and verified the excellent discrimination ability of this grouping for all genes, immune signatures, and TME compositions.

Comparison of Clinical Characteristics Between the H-Immunity and L-Immunity Groups
To understand the effects of immune levels on clinical traits in PTC patients, we compared clinical parameters between the H-immunity and L-immunity groups ( Table 2). Compared to patients in the L-immunity group, patients in the H-immunity group had an advanced American Joint Committee on Cancer (AJCC) stage, advanced T classification, lymph node metastasis, higher tall-cell PTC and lower follicular PTC proportions, more BRAF mutations and fewer RAS mutations.

WGCNA Screening Key Modules and Hub Genes Related to PTCImmunity
To construct a weighted coexpression network and screen for modules and genes related to PTC immunity, 6254 DEGs between PTC (n = 512) and normal tissues (n = 58) in TCGA  were selected for subsequent WGCNA (502 patients with complete clinical information were selected) ( Figure 7A). After a series of adjustments to the WGCNA parameters, the DEGs were divided into 19 modules by average linkage hierarchical clustering ( Figures 7B-D). The correlation heatmap showed that the pink module, which contains 276 genes, had the highest correlation with the H-immunity group (Pearson's correlation coefficient = 0.51, p < 0.001) ( Figure 7E). The 34 genes in the pink module with MM > 0.5 and GS > 0.5 were defined as hub genes ( Figure 7F).

Analysis of the Main Mechanism of PTC Immunity
A total of 276 pink module eigengenes were used for functional and pathway enrichment analyses to explore potential mechanisms that cause differences in the immune microenvironment between PTCs. "Response to interferon-gamma", "external side of plasma membrane", and "chemokine receptor binding" were the most common GO terms for biological processes, cellular components, and molecular functions, respectively ( Figure 8A). "Epstein-Barr virus infection" was the most significantly enriched pathway in KEGG analysis (p-adjust < 0.001) ( Figure 8B). Additionally, we constructed a PPI network with the 34 hub genes and found that TYROBP, CTLA4, CSF1R, CTSS, and HLAs play key roles ( Figure 8C). GSEA was also used to identify the mechanism and functional differences between the L-immunity and H-immunity groups. The H-immunity group was enriched in the following terms: cytokinecytokine receptor interaction, autoimmune thyroid disease, cell adhesion molecules (CAMs), natural killer cell-mediated cytotoxicity, chemokine signaling pathway, viral myocarditis, Leishmania infection, intestinal immune network for IgA production, hematopoietic cell lineage, and type I diabetes mellitus ( Figure 8D). Conversely, the L-immunity group was enriched in the following metabolism-related pathways: butanoate metabolism; glycosylphosphatidylinositol (GPI) anchor biosynthesis; glycine, serine, and threonine metabolism; lysine degradation; valine, leucine, and isoleucine degradation;  arginine and proline metabolism; beta alanine metabolism; Huntington's disease; fatty acid metabolism; and glycerolipid metabolism ( Figure 8E).

DISCUSSION
Our findings showed that during the occurrence and development of PTC, the immune system is enhanced; moreover, M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages play a tumor-promoting role, and their abundances and proportions are significantly increased. M1 macrophages, CD8+ T cells, B cells, NK cells, and T FH cells (eosinophils, gd T cells, and Th17 cells with weak supporting evidence) exert an antitumor effect, and the proportions decreased. PTC is subdivided into "inflammatory" tumors (71). Extensive evidence shows that autoimmunity and inflammation are risk factors for TC (72). Increasing evidence also suggests that cancerrelated inflammation may be a useful target for new diagnostic and therapeutic strategies in TCs (73). Therefore, to evaluate the specific pattern of immune cells involved in PTC, including not only its phenotype but also its function, it is essential to understand the immunological characteristics of PTC. Many major advances have been made in the study of immune infiltration in PTC, especially in major immune cells (B cells, TAMs, NK cells, neutrophils, DCs, MCs, CD8+ T cells and Tregs) (14)(15)(16). A considerable part of our research is repeating, verifying, and integrating previously reported results of immune cell phenotypes and functions in PTC. We used ssGSEA, immune cell marker and CIBERSORT to evaluate the infiltration level of immune cells in 10 datasets (a total of 799 PTC and 194 normal samples) from the perspective of abundance and proportion, respectively. Furthermore, it is essential to ensure that the results of the analysis are highly reliable. We also analyzed subtypes of major immune cells and immune cells that lacked attention (such as M1 macrophages, M2 macrophages, Th1, Th2, pDCs, iDCs, eosinophils, monocytes, T FH cells, T cells gamma delta, and Th17 cells) This will help researchers comprehensively understand the PTC immune microenvironment. Our research supports the evidence that M2 macrophages, Tregs, monocytes, neutrophils,  DCs, MCs, and TAMs play a tumor-promoting role, while M1 macrophages, CD8+ T cells, B cells, NK cells, and T FH cells (eosinophils, gd T cells, and Th17 cells with weak supporting evidence) exert an antitumor effect. Among them, related reports of the phenotype and function of monocytes, gd T cells, T FH cells, and eosinophils in PTC patients are still lacking. Inflammation of monocyte recruitment and activation is associated with the development of TC in a mouse model (74). IL-17-producing gdT17 cells play a decisive role in the chemotherapy-induced anticancer immune response (75). T FH cells can regulate antibody production by B cells (76); the frequency of T FH cells increases in patients with autoimmune thyroid disease (AITD), and a significant number of T FH cells are also detected in the thyroid tissues of patients with Hashimoto's thyroiditis (77). The above evidence also supports our conclusions about the function of immune cells.
Our research shows that compared with the immune level in normal tissues, the immune level in PTC tissues is generally higher, and the increase in tumor-promoting immune cells is particularly significant. The proportion of immune cells in the TME is also favoring tumor-promoting immune cells. The elevation of TAMs, DCs, MCs, neutrophils, NK cells, and Tregs in PTC has been reported (28-31, 36, 38, 41, 45, 55, 78). Our study found that these tumor-promoting immune cells were recruited in large numbers in PTC and exhibited a tumorpromoting immune microenvironment phenotype. This indicates that PTCs manipulate partial immune cells such as M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages to fight against the body's immune response to protect themselves, which is also called "sabotage" (13). Eventually, the tumor evaded immune surveillance and escaped the immune system.
In additional research, we found that immune cell infiltration increased during tumor progression, and the increase in tumorpromoting immune cells was particularly significant. The proportion of immune cells in the TME was further tilted towards tumor-promoting immune cells. Previous mainstream studies have shown that TAMs, DCs, MCs, neutrophils, and Tregs are positively correlated with PTC progression (24, 29-34, 38, 41, 42, 46, 48, 49, 51, 58, 61, 68, 69, 79), and B cells, CD8+ T cells, NK cells, and iDCs are negatively correlated with PTC progression (26,37,56,(58)(59)(60)(61)(62). Unlike previous studies, our research found that both the abundance and proportion of tumor-promoting immune cells were positively correlated with PTC progression, while the proportion of antitumor immune cells during the progression of PTC decreased (CIBERSORT), and the abundance increased (ssGSEA) or did not significantly change (immune cell marker). These results show that in escalating immune confrontation, immune escape is aggravated, and the tumor's ability to fight the immune response of the body is further strengthened by recruiting immunosuppressive cells. Eventually, the dynamic balance between antitumor and tumorpromoting immune cells in the original immune system has allowed the irreversible development of tumors. The positive correlations among the H-immunity group, tumor stage, lymph node metastasis, and tumor size also support the above hypothesis.
Our analysis of the correlation of the abundances of immune cells in PTC showed a high positive correlation among all examined cells, indicating that the immune system as a whole has a high degree of consistency. However, antitumor immune cells and tumor-promoting immune cells oppose each other at a proportional level of correlation. Our clustering results on immune cells also show the opposition between immune cells that exert procancer and anticancer effects. We speculate that in the process of tumor progression, the confrontation between tumorkilling immune cells and tumor-controlled immunosuppressive cells in the human immune system continues to escalate, and eventually, the proportion of tumor-promoting immune cells in the TME will have an irreversible advantage. This indicates that the tumor cells escaped the body's immune surveillance by "sabotage" and achieved immune escape. Except for M0 macrophages (proportion) in our prognosis studies, other immune cells did not reach significant levels. Previous studies in PTC have generally shown that TAMs and neutrophils are associated with a worse prognosis (28,(51)(52)(53)(54), while CD8+ T cells and DCs are associated with a better prognosis (26,39,62,63). In fact, the roles of various immune cells in prognosis are still controversial (26,64), and there are many negative results (36,48,50,64,70). Considering the positive correlation between immune cell infiltration and tumor progression, we suspect that the insignificant prognosis of immune cells may be related to the lower number of outcomes, such as death, recurrence and metastasis of PTC.
We divided patients into L-immunity and H-immunity groups by consensus clustering and verified the excellent discrimination ability of this grouping for all genes, immune signatures, and TME composition. Compared to patients in the L-immunity group, patients in the H-immunity group showed a more advanced stage, larger tumor size, more lymph node metastasis, higher tall-cell PTC proportions, lower follicular PTC proportions, more BRAF mutations and fewer RAS mutations. The BRAF V600E mutation is related to the immunosuppressive mechanism of PTC; compared to the ratio in wild-type BRAF PTC, the CD8+/FoxP3+ ratio in BRAF V600E mutant PTC was significantly reduced, and the proportion of M2 type TAMs was increased (27,80,81). We speculate that somatic mutations in genes such as BRAF and RAS in PTC may initiate the changes in the tumor immune microenvironment.
WGCNA and GSEA were used to explore the underlying mechanisms that cause PTC immune differences, and a number   carcinogenesis (86)(87)(88). We speculate that EBV may play an important role in the change in the PTC immune microenvironment. However, this study has some limitations. First, ssGSEA, immune cell markers, and CIBERSORT are algorithms based on RNA-seq data, and the abundant results of ssGSEA and immune cell markers are inconsistent in some antitumor immune cells. We have shown their respective results and common conclusions. Second, the research on the mechanism that causes immune differences was based on the inference of several algorithms and has not been experimentally verified.
In the next study, single-cell analysis of immune cells around and within tumors will be an important direction for PTC  immune research to elucidate the function of immune cells in the PTC TME.
In conclusion, our research shows that M2 macrophages, Tregs, monocytes, neutrophils, DCs, MCs, and M0 macrophages play a tumor-promoting role in PTC, while M1 macrophages, CD8+ T cells, B cells, NK cells, T FH cells (eosinophils, gd T cells, and Th17 cells with weak supporting evidence) play an antitumor role. During the occurrence and development of PTC, the overall immune level increased, and the abundance and proportion of tumor-promoting immune cells significantly increased. We speculate that EBV may play an important role in changing the immune microenvironment of PTC.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
ZX designed the analytical strategies, performed data analyses, and wrote the manuscript. YL, XL, YZH, SW, SYW, JS, and YCH performed data analyses. JZ conceived the research and wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This research was supported by the National Natural Science Foundation of China (81600602).

ACKNOWLEDGMENTS
We would like to thank Mr. Guangqi Li for his support with bioinformatics analysis and pay tribute to the contributions of public databases such as the TCGA and GEO to human medicine. This manuscript has been released as a pre-print at [ResearchSquare] (89).