Identification of CXCL8 as a potential gene associated with lymph node metastasis in papillary thyroid carcinoma through bioinformatics analysis.

Analysis of pan-cancer data from the UCSC Xena platform revealed elevated CXCL8 expression in tumor samples from multiple cancer types, including Thyroid Cancer (THCA), Bile Duct Cancer (CHOL), Colon Cancer (COAD), Esophageal Cancer (ESCA), Head and Neck Cancer (HNSC), Kidney Papillary Cell Carcinoma (KIRP), Rectal Cance (READ), Stomach Cancer (STAD), and Endometrioid Cancer (UCEC). In contrast, CXCL8 expression was significantly higher in normal samples than in tumor tissues only in Kidney Chromophobe (KICH)(Fig. 1  A). We next focused on PTC by downloading transcriptome data (HTSeq-FPKM) and corresponding clinical information from the TCGA database, comprising 355 PTC and 49 normal thyroid samples. CXCL8 expression was significantly upregulated in PTC compared to normal thyroid tissue ( P  < 0.001, Fig.  1 B). This finding was corroborated by paired differential analysis ( P  < 0.001, Fig. 1  C) and further validated in the independent GSE33630 dataset ( P  < 0.05, Fig.  1 D). To experimentally confirm these observations, we performed RNA sequencing on 7 paired PTC and normal tissues from Gansu Provincial People’s Hospital. Differential expression analysis (threshold: FC > 1.5, |log 2 FC| > 0.585, P  < 0.05) identified 4904 significantly differentially expressed genes (DEGs), including 3035 upregulated and 1869 downregulated genes. Notably, CXCL8 was among the most significantly upregulated genes (log 2 FC = 3.03, P  < 0.001, Fig.  1 E), supporting its potential role as a biomarker in PTC. Unsupervised clustering based on the DEGs clearly separated PTC and normal thyroid samples into two distinct groups, reflecting substantial transcriptome-level heterogeneity (Fig. 1  F). CXCL8 exhibited a consistent high-expression pattern in PTC samples, aligning with the volcano plot results. Minor deviations in gene expression patterns in certain samples (e.g., Cancer-3, Normal-2) suggested underlying molecular heterogeneity within the tumor cohort. Fig. 1 Expression Analysis of CXCL8 in PTC and Various Cancers.( A ) Perform a pan-cancer expression analysis of CXCL8 between normal and tumor samples based on UCSC Xena data.( B ) In the TCGA database, the expression of CXCL8 in thyroid papillary carcinoma tumor samples is higher than in normal samples ( P  < 0.001)( C ) The expression of CXCL8 in PTC samples from the TCGA database is higher than in adjacent normal tissues ( P  < 0.001).( D )In the GSE33630 dataset, the expression of CXCL8 is higher in PTC samples compared to normal samples ( P  < 0.05).( E )The volcano plot shows that CXCL8 is upregulated in tumors based on RNA sequencing data from 7 pairs of PTC and normal thyroid tissues (from Gansu Provincial People’s Hospital)( P  < 0.001).( F ) The heatmap shows some upregulated and downregulated genes, with samples clustered into the PTC group and the normal group (from Gansu Provincial People’s Hospital)( P  < 0.05). Expression Analysis of CXCL8 in PTC and Various Cancers.( A ) Perform a pan-cancer expression analysis of CXCL8 between normal and tumor samples based on UCSC Xena data.( B ) In the TCGA database, the expression of CXCL8 in thyroid papillary carcinoma tumor samples is higher than in normal samples ( P  < 0.001)( C ) The expression of CXCL8 in PTC samples from the TCGA database is higher than in adjacent normal tissues ( P  < 0.001).( D )In the GSE33630 dataset, the expression of CXCL8 is higher in PTC samples compared to normal samples ( P  < 0.05).( E )The volcano plot shows that CXCL8 is upregulated in tumors based on RNA sequencing data from 7 pairs of PTC and normal thyroid tissues (from Gansu Provincial People’s Hospital)( P  < 0.001).( F ) The heatmap shows some upregulated and downregulated genes, with samples clustered into the PTC group and the normal group (from Gansu Provincial People’s Hospital)( P  < 0.05). Through the analysis, it can be seen that CXCL8 and T stage, N stage is statistically significant, and there is no statistical significance between CXCL8 and clinical characteristics including age, gender, clinical stage and M stage(Fig.  2 ). Fig. 2 ( A ) Heatmap illustrating the distribution of key clinicopathological features (including T stage and N stage) between the CXCL8 high- and low-expression groups.( B-G )Boxplots showing the relationship between CXCL8 expression levels and other clinical parameters: age, gender, M stage, and clinical stage. Statistical significance was observed for T stage and N stage ( P  < 0.05), while no significant associations were found for age, gender, M stage, or clinical stage. ( A ) Heatmap illustrating the distribution of key clinicopathological features (including T stage and N stage) between the CXCL8 high- and low-expression groups.( B-G )Boxplots showing the relationship between CXCL8 expression levels and other clinical parameters: age, gender, M stage, and clinical stage. Statistical significance was observed for T stage and N stage ( P  < 0.05), while no significant associations were found for age, gender, M stage, or clinical stage. Based on RNA sequencing data from the TCGA database, a co-expression analysis identified 42 genes that were strongly and positively correlated with CXCL8 expression. Among these, 11 genes with the highest correlation coefficients were selected to construct a gene-gene interaction network (Fig.  3 ). Differential expression analysis between the CXCL8 high-expression and low-expression groups revealed a total of 922 DEGs, including 274 upregulated and 648 downregulated genes (Fig.  3 A-B). GO and KEGG pathway enrichment analyses demonstrated that these DEGs were primarily associated with extracellular matrix organization and cytokine–cytokine receptor interaction pathways (Fig.  3 C–E). Moreover, GSEA further confirmed that the cytokine–cytokine receptor interaction pathway was significantly enriched in patients with high CXCL8 expression (Fig. 3  F). Fig. 3 ( A )Circos plots show the co-expression network of CXCL8 with 11 genes in PTC samples. Each part of the circle represents a gene, and its width indicates the total amount of co-occurrence connecting a gene to another gene. The width of each link represents the total expression time of the linked gene.( B )The heat map shows the top 50 down-regulated genes and top 50 up-regulated genes found in the high expression group.( C-D ) Results from GO analysis.( E )Results from KEGG enrichment analysis.( F )GSEA analysis revealed upregulated and downregulated pathways associated with CXCL8 expression. ( A )Circos plots show the co-expression network of CXCL8 with 11 genes in PTC samples. Each part of the circle represents a gene, and its width indicates the total amount of co-occurrence connecting a gene to another gene. The width of each link represents the total expression time of the linked gene.( B )The heat map shows the top 50 down-regulated genes and top 50 up-regulated genes found in the high expression group.( C-D ) Results from GO analysis.( E )Results from KEGG enrichment analysis.( F )GSEA analysis revealed upregulated and downregulated pathways associated with CXCL8 expression. Based on “CIBERSORT”, the difference of immune cell infiltration levels between the CXCL8 high expression group and the low expression group was analyzed, which confirmed that the expression of CXCL8 was significantly correlated with tumor-infiltrating immune cells(Fig.  4 A-B). We further analyzed the correlation between CXCL8 and immune cells and showed that the infiltration levels of activated dendritic cells, neutrophils, resting dendritic cells(DCs), macrophages M 0 and activated mast cells were positively correlated with CXCL8( Fig.  4 C-G) . However, the infiltration levels of CD8 + T cells, plasma cells, macrophage M 1 , and γδ T cells were negatively correlated with CXCL8 (Fig.  4 H-K).The scores of StromalScore, ImmuneScore, and ESTIMATEScore were significantly different between the CXCL8 high and low expression groups, and the CXCL8 high expression group was significantly higher than the CXCL8 low expression group in all three scores(Fig.  4 L). We then explored the correlation between CXCL8 and immune checkpoints using the expression levels of CXCL8 and 27 immune checkpoint-related genes. CXCL8 was significantly correlated with all 27 immune checkpoint-related genes ( P  < 0.001), and all of them were positively correlated.(Fig.  4 M)(STable 1). We also obtained treatment scores for anti-CTLA4 and anti-PD-1 inhibitors, with higher treatment scores indicating better treatment effect. The analysis of immunotherapy response scores revealed that patients with low CXCL8 expression had significantly higher scores for anti-CTLA4 treatment (median [IQR]: 1.58 [1.21–1.95]) compared to those with high CXCL8 expression (median [IQR]: 0.85 [0.62–1.14]; Mann-Whitney U test, P  = 0.0073).Statistical significance was found only in the CTLA4-positive -PD-1-negative fraction, indicating that anti-CTLA4 therapy is a better choice than anti-PD-1 therapy for patients with low CXCL8 expression(Fig. 4  N). Fig. 4 ( A )Correlation between CXCL8 and immune infiltration.( B )The lollipop plot shows the correlation between immune infiltration and CXCL8 expression.( C-G )Scatter plots show the correlation between CXCL8 expression levels and immune cells: the infiltration levels of activated DCs, neutrophils, resting DCs, macrophages M 0 , and activated mast cells were positively correlated with CXCL8 .( H-K )Scatter plots show the correlation between CXCL8 expression levels and immune cells: infiltration levels of CD8 + T cells, plasma cells, macrophage M 1 , and γδ T cells were negatively correlated with CXCL8 .( L )Immunoscore of CXCL8 in high and low expression groups.( M )Correlation between CXCL8 and immune checkpoint related genes. Twenty-seven immune checkpoints were positively correlated with CXCL8 .( N ) CXCL8 immunotherapy scores of anti-CTLA4 and anti-PD-1 inhibitors in high and low expression groups. ( A )Correlation between CXCL8 and immune infiltration.( B )The lollipop plot shows the correlation between immune infiltration and CXCL8 expression.( C-G )Scatter plots show the correlation between CXCL8 expression levels and immune cells: the infiltration levels of activated DCs, neutrophils, resting DCs, macrophages M 0 , and activated mast cells were positively correlated with CXCL8 .( H-K )Scatter plots show the correlation between CXCL8 expression levels and immune cells: infiltration levels of CD8 + T cells, plasma cells, macrophage M 1 , and γδ T cells were negatively correlated with CXCL8 .( L )Immunoscore of CXCL8 in high and low expression groups.( M )Correlation between CXCL8 and immune checkpoint related genes. Twenty-seven immune checkpoints were positively correlated with CXCL8 .( N ) CXCL8 immunotherapy scores of anti-CTLA4 and anti-PD-1 inhibitors in high and low expression groups. We obtained raw single-cell RNA sequencing data of one PTC sample and one normal thyroid sample from the GEO database GSE241184 . After multiple rounds of cluster analysis and dimensionality reduction, 21 cell clusters were finally identified and labeled as different cell types (Fig.  5 A-C). The data showed that CXCL8 was highly expressed in DCs. We also performed a systematic analysis of CXCL8 expression levels in various cell types (Fig.  5 D-E). Fig. 5 ( A-B )Visualization of reclustered immune cell type annotations in 2D UMAP.( C )Punctate distribution of key genes in the CXCL8 and PI3K-AKT pathways. ( D-E ) shows the differential expression of key genes in the CXCL8 and PI3K-AKT pathways in different cell types. ( A-B )Visualization of reclustered immune cell type annotations in 2D UMAP.( C )Punctate distribution of key genes in the CXCL8 and PI3K-AKT pathways. ( D-E ) shows the differential expression of key genes in the CXCL8 and PI3K-AKT pathways in different cell types. In order to verify the expression of CXCL8 in PTC tissues, we extracted mRNA from the tissues for qRT-PCR analysis.The results showed that compared with the PNTs, the relative expression level of CXCL8 mRNA in PTC tissues was significantly higher than the control tissues, with a significant statistical difference ( P  = 0.0035)(Fig.  6 ). The results are shown in the figure. This result is consistent with the results obtained from the analysis of GEO and TCGA data in bioinformatics, indicating that CXCL8 is abnormally highly expressed in PTC tissues. Fig. 6 Fold change of CXCL8 mRNA expression in PTC and PNT (2^ −ΔΔCt ). Asterisks indicate the level of statistical significance: ** P  < 0.01. Fold change of CXCL8 mRNA expression in PTC and PNT (2^ −ΔΔCt ). Asterisks indicate the level of statistical significance: ** P  < 0.01. To validate CXCL8 expression at the protein level and its clinical significance in PTC, we performed IHC analysis on a tissue microarray constructed from 49 paired PTC and PNT tissues. In PTC tissues (Fig.  7 A-B), HE staining clearly revealed characteristic pathological features of PTC, including papillary architecture, nuclear overlapping, ground-glass nuclei, and nuclear grooves. Corresponding IHC staining (Fig.  7 B) demonstrated strong positive staining for the CXCL8 protein, which was localized in the cytoplasm of PTC cells, appearing intense and widely distributed throughout the tumor regions. In contrast, H&E staining of adjacent normal thyroid tissues (Fig.  7 C-D) showed preserved normal thyroid follicular structure with uniformly sized follicles containing colloid. In these normal tissues, IHC staining for CXCL8 (Fig.  7 D) was negative or showed only very weak positivity, forming a sharp contrast with the PTC tissues. Based on this IHC analysis, we concluded that CXCL8 protein is significantly upregulated in PTC tissues. We subsequently investigated the correlation between CXCL8 protein expression levels and various clinicopathological parameters. Statistical analysis revealed that high CXCL8 expression was significantly associated with lymph node metastasis (χ 2  = 16.034, P  < 0.001). Specifically, all patients with central lymph node metastasis (N 1 a, n  = 15) or lateral (± central) lymph node metastasis (N 1 b, n  = 11) exhibited high CXCL8 expression, whereas only 52.2% (12/23) of patients without lymph node metastasis (N 0 ) showed high CXCL8 expression (Table 2). Furthermore, a significant inverse correlation was observed between CXCL8 expression and the presence of Hashimoto’s thyroiditis (χ 2  = 7.646, P  = 0.006). Among patients with Hashimoto’s thyroiditis ( n  = 10), only 40.0% (4/10) showed high CXCL8 expression, compared to 87.2% (34/39) in patients without Hashimoto’s thyroiditis. No statistically significant correlations were found between CXCL8 expression and other clinicopathological features, including patient gender, age, tumor diameter, tumor location, recurrence risk stratification, capsular invasion, extrathyroidal extension, or multifocality (all P  > 0.05; Table 2). These findings indicate that CXCL8 protein is specifically and highly expressed in PTC tissues, and its elevated expression is closely associated with lymph node metastasis. Furthermore, its expression may be influenced by the presence of Hashimoto’s thyroiditis, providing further evidence for the potential role of CXCL8 in PTC progression and metastasis. Fig. 7 Immunohistochemical analysis reveals significantly elevated CXCL8 protein expression in PTC tissues relative to PNTs.( A ) HE staining of PTC tissues. ( B ) CXCL8 staining was positive in PTC tissues. ( C ) HE staining of PNTs. ( D ) CXCL8 staining was negative in PNTs.All images were acquired at 200× magnification. Scale bars: 100 μm. Immunohistochemical analysis reveals significantly elevated CXCL8 protein expression in PTC tissues relative to PNTs.( A ) HE staining of PTC tissues. ( B ) CXCL8 staining was positive in PTC tissues. ( C ) HE staining of PNTs. ( D ) CXCL8 staining was negative in PNTs.All images were acquired at 200× magnification. Scale bars: 100 μm. Table 2 Relationship between CXCL8 expression and clinicopathological features in PTC. Clinicopathological Feature n CXCL8-High, n (%) CXCL8-Low, n (%) χ² P -value Odds Ratio(95% Confidence Interval) Gender Male 10 8 (80.0) 2 (20.0) – 1.000 1.20 (0.21–6.78) Female 39 30 (76.9) 9 (23.1) 1.00 (reference) Age (years) , Median [IQR] 49 48 [38–57] 52 [45–61] 0.191 - Tumor diameter (cm) 0.416 0.901 ≤ 1 18 14 (77.8) 4 (22.2) - > 1 and ≤ 2 23 17 (73.9) 6 (26.1) - > 2 and ≤ 4 8 7 (87.5) 1 (12.5) - > 4 0 0 (0.0) 0 (0.0) - Tumor location – 0.899 Left lobe 18 13 (72.2) 5 (27.8) 1.00 (reference) Right lobe 17 14 (82.4) 3 (17.6) 1.79 (0.41–7.81) Bilateral 12 9 (75.0) 3 (25.0) 1.15 (0.25–5.26) Isthmus 2 2 (100.0) 0 (0.0) - Recurrence risk stratification 3.523 0.171 Low risk 25 17 (68.0) 8 (32.0) 1.00 (reference) Moderate risk 24 21 (87.5) 3 (12.5) 3.29 (0.83–13.03) High risk 0 0 (0.0) 0 (0.0) - Capsular invasion 2.194 0.138 No 26 18 (69.2) 8 (30.8) 1.00 (reference) Yes 23 20 (87.0) 3 (13.0) 2.96 (0.74–11.82) Extrathyroidal extension 2.667 0.102 No 25 17 (68.0) 8 (32.0) 1.00 (reference) Yes 24 21 (87.5) 3 (12.5) 3.29 (0.83–13.03) Hashimoto’s thyroiditis – 0.006 No 39 34 (87.2) 5 (12.8) 1.00 (reference) Yes 10 4 (40.0) 6 (60.0) 0.10 (0.03–0.36) Multifocality – 0.047 No 47 38 (80.9) 9 (19.1) 1.00 (reference) Yes 2 0 (0.0) 2 (100.0) 0.00 (0.00–0.88) Lymph Node Status N 0 23 12 (52.2) 11 (47.8) – 1.00 (reference) 1.00 (reference) N 1 a (Central only) 15 15 (100.0) 0 (0.0) – < 0.001 – N 1 b (Lateral ± Central) 11 11 (100.0) 0 (0.0) – < 0.001 – Footnotes: (a) Total samples: n  = 49. CXCL8 expression was dichotomized into High and Low groups based on the median IHC score. (b) No missing data were encountered for any variable. (c) For comparisons with any expected cell count < 5, Fisher’s exact test was applied. (d) For 2 × 2 tables with complete separation (e.g., N 1 a, N 1 b: 100% vs. 0%), the two-sided Fisher’s exact test was used. Effect sizes (OR) were estimated using Firth logistic regression (or Haldane–Anscombe 0.5 correction). (e) Extreme P -values are reported as P  < 0.001. (f) The Chi-squared test was applied for other comparisons; Yates’ correction was not used. (g) Age is presented as median [IQR] due to the loss of information from dichotomization. Relationship between CXCL8 expression and clinicopathological features in PTC. Footnotes: (a) Total samples: n  = 49. CXCL8 expression was dichotomized into High and Low groups based on the median IHC score. (b) No missing data were encountered for any variable. (c) For comparisons with any expected cell count < 5, Fisher’s exact test was applied. (d) For 2 × 2 tables with complete separation (e.g., N 1 a, N 1 b: 100% vs. 0%), the two-sided Fisher’s exact test was used. Effect sizes (OR) were estimated using Firth logistic regression (or Haldane–Anscombe 0.5 correction). (e) Extreme P -values are reported as P  < 0.001. (f) The Chi-squared test was applied for other comparisons; Yates’ correction was not used. (g) Age is presented as median [IQR] due to the loss of information from dichotomization.

The RNA sequencing data and related clinical information of patients with PTC were downloaded from the TCGA ( https://cancergenome.nih.gov/ ) database.We downloaded the transcriptome profiling (HTseq-FPKM) gene expression quantitative data of 404 PTC (including 355 tumors and 49 para-cancerous tissues) from TCGA database.The list of immune-related genes was downloaded from the import ( https://www.immport.org/shared/home ) and InnateDB ( https://www.innatedb.com/ ) databases. This study included patients with PTC. The raw data were summarized and processed using Strawberry perl (v5.30.0), and subsequent analyses were conducted using R software (v4.3.3).Our comprehensive from Gene Expression Omnibus (GEO) database ( http://www.ncbi.nlm.nih.gov/geo/ ) for thyroid papillary carcinoma GSE33630 data set, to verify that the related gene expression.The GSE33630 dataset comprises 11 cases of anaplastic thyroid carcinomas (ATC), 49 cases of PTC, and 45 cases of normal thyroid tissues (N). For subsequent analyses in this study, 49 cases of PTC and 45 cases of N were selected.we downloaded pan-cancer data from the UCSC Xena platform ( https://xena.ucsc.edu/ ) to verify the expression levels of CXCL8 between Paired non-cancerous thyroid tissues(PNTs) and tumor samples. TCGA-derived PTC transcriptome and clinical data—including age at diagnosis, sex, survival time, survival status, clinical stage, and pathological TNM staging—were processed with Strawberry Perl and analyzed with R. Associations between CXCL8 expression levels and clinical characteristics (age, clinical stage, T stage, M stage, and N stage) were evaluated using χ² tests or Wilcoxon signed-rank tests. Results were visualized as boxplots and heatmaps. To elucidate the biological role of CXCL8 in PTC, correlation analyses were performed between CXCL8 and other genes, with thresholds set at | Cor | > 0.6 and P   1 and FDR  < 0.05 as cutoffs. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and Gene Set Enrichment Analysis (GSEA) were performed to identify associated molecular pathways and biological processes. R packages used included “ggplot2”, “ggpubr”, “ggExtra”, “circlize”, “corrplot”, and “pheatmap”. To assess the relationship between CXCL8 expression and immune infiltration, CIBERSORT was used to estimate the proportions of 22 immune cell subsets. Infiltration differences between high- and low- CXCL8 groups were displayed as boxplots, while scatter plots and lollipop charts illustrated correlations. To evaluate immunotherapeutic potential, associations between CXCL8 expression and immune checkpoint genes and drug targets were analyzed. Immunotherapy scores for anti-CTLA4 and anti-PD-1 responses were derived from TCIA ( https://tcia.at ), and response differences between CXCL8 expression groups were compared. R packages applied included “reshape2”, “ggpubr”, “vioplot”, and “ggExtra”. Single-cell RNA sequencing data ( GSE241184 ) from lymph node metastases of PTC patients were obtained from the GEO database. Data integration was performed using the “harmony” R package. Dimensionality reduction, clustering, and cell-type annotation were conducted using “Seurat” and “SingleR”. Violin plots were generated to visualize gene expression within annotated cell types. Screening thresholds were set at P   1. To validate bioinformatics findings, we collected clinical tissue samples from 49 patients with PTC who underwent thyroidectomy at Gansu Provincial Hospital (Lanzhou, China), with the following inclusion criteria: (1)pathologically confirmed PTC post-thyroidectomy; (2)no prior anti-tumor therapy (radiotherapy, chemotherapy, or targeted therapy) before surgery; (3)informed consent obtained from patients or their legal guardians. From these 49 patients, two subsets of samples were used for different experiments: 7 paired tissue samples (7 PTC tissues + 7 PNTs): Tissue (0.5–1 cm in diameter) was excised from the center of PTC tissues; PNTs (of similar size) were obtained from areas > 2 cm from the edge of the tumor. Samples were rinsed with normal saline, rapidly cryopreserved in liquid nitrogen, and confirmed by hematoxylin-eosin (HE) staining by a professional pathologist. These samples were used for high-throughput RNA sequencing. 49 paired tissue samples (49 PTC tissues + 49 PNTs): These samples were used for tissue microarray construction (for immunohistochemical analysis) and qRT-PCR validation (a subset of this cohort, as detailed in “RNA extraction and quantitative real-time polymerase chain reaction”). The experimental protocol was approved by the Institutional Review Board of Gansu Provincial People’s Hospital (Approval Number: 2022 − 182), and all procedures complied with the Declaration of Helsinki (1964) and its later amendments. Fourteen tissue samples from 7 patients from Gansu Provincial People’s Hospital were collected in this study. Tissue 0.5 to 1 cm in diameter was removed from the center of the patient’s PTC tissue. PNTs of similar size was obtained from the patient more than 2 cm from the edge of the thyroid papillary carcinoma tissue. Samples were collected, rinsed with normal saline, and rapidly cryopreserved in liquid nitrogen. The tissues were stained with HE and confirmed by a professional pathologist. The patient had not received preoperative radiotherapy or chemotherapy. Total RNA extraction was provided by Cloud Sequencing Biotechnology (Shanghai, China). RNA concentration was measured for each sample using a NanoDrop ND-1000 instrument (Thermo Fisher Scientific), considering an OD 260/280 value of 1.8 to 2.1 to ensure good RNA purity. CloudSeqInc. (Shanghai, China) provided the high-throughput RNA sequencing service. rRNA was removed from total RNA using the GenSeq ® rRNA Removal Kit (GenSeq, Inc.) kit. After RNA removal, RNA sequencing libraries were constructed using the GenSeq ® Low-Input Whole RNA Library Preparation Kit (GenSeq, Inc.). The constructed sequencing libraries were quality-controlled and quantified using a Bioanalyzer 2100 system (Agilent Technologies), and high-throughput sequencing of the libraries was controlled using an Illumina NovaSeq 6000 instrument with a paire-end length of 150 bp. To validate the CXCL8 expression pattern identified by bioinformatics analysis, we selected 10 paired samples (10 PTC tissues + 10 PNTs from the aforementioned cohort of 49 PTC patients (inclusion criteria: same as “Clinical Tissue Sample Collection”; selection method: simple random sampling to avoid selection bias). Total RNA was extracted from these 20 tissues (10 pairs) using the TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions.The concentration and purity of the extracted total RNA were determined using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, USA). The RNA concentrations typically ranged from 178 to 489 ng/µL.Only RNA samples with an A260/A280 ratio between 1.8 and 2.1, indicating high purity without protein or other contamination, were used for subsequent experiments. For reverse transcription, 1 µg of total RNA from each qualified sample was used as input for cDNA synthesis using the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara Bio Co., Ltd., Japan) according to the manufacturer’s instructions. Real-time fluorescence quantitative PCR was conducted using TB Green Premix Ex Taq II (Takara Bio Co., Ltd., Japan) on a Thermo Scientific 7500 Real-Time PCR System (Thermo Fisher Scientific, USA). The reaction conditions were: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 34 s. A melt curve analysis was performed post-amplification to confirm the specificity of the PCR products. The cycle threshold (C t ) values were obtained for both CXCL8 and the reference gene β-actin . The comparative ΔΔC t method was used to analyze the qRT-PCR data. For each sample, the ΔC t = C t(CXCL8) –C t(β−actin) .For statistical comparison between PTC tissues and their PNTs, the non-parametric Wilcoxon signed-rank test was applied directly to the ΔC t values.The relative gene expression fold change for visualization was calculated as 2^ –ΔΔ Ct , where ΔΔC t   = ΔC t(PTC) –ΔC t(PNT) . Each sample was tested in triplicate, and the average C t value was used for calculations. The primers for CXCL8 and β-actin were synthesized by Beijing Qingke Biotechnology Co., Ltd. (Beijing, China), and the primer sequences are shown in Table 1. Table 1 Primer sequences for all genes tested. Gene Sequence(5’−3’) β-actin- F GCACACAAGCTTCTAGGACAAGA β-actin- R GTTCTTTAGCACTCCTTGGCAAA CXCL8 -F AAATCTGGCACCACACCTTCTAC CXCL8 -R CAGCCTGGATAGCAACGTACAT Primer sequences for all genes tested. The primary antibody used was anti - CXCL8 (Clone EPR19824 ;Proteintech Group, Inc., Chicago, IL, USA; Cat#17038-1-AP; RRID: AB_2126217;dilution 1:200).The secondary antibody was a HRP-conjugated goat anti-rabbit IgG antibody(Cell Signaling Technology, Danvers, MA, USA; Cat#7074;RRID: AB_2099233;dilution 1:3000.The chromogenic substrate was 3, 3’ - diaminobenzidine (DAB)(Dako, Agilent Technologies, Cat. #K3468, RRID: AB_2894767)and prepared according to the manufacturer’s instructions.Tissue sections were deparaffinized twice in xylene (20 min each) and rehydrated through a graded ethanol series (100%, 95%, 80%, and 60%; 3 min each). Following three rinses with distilled water (1 min each), antigen retrieval was performed by submerging sections in retrieval solution and heating them in a microwave at medium power for 10 min. The sections were then cooled in buffer for 35 min. Prior to immunostaining, endogenous peroxidase activity was quenched by incubating sections in 3% H₂O₂ (room temperature, 10 min), followed by three washes with 1X TBST (1 min each). Non-specific binding was blocked with blocking buffer (1 h, room temperature), after which sections were incubated with primary antibody (1.5 h, room temperature), washed with TBST, and treated with peroxidase-conjugated secondary antibody (30 min, room temperature). Color development was achieved using DAB (2–5 min), and sections were counterstained with hematoxylin, rinsed with TBST, and immersed in distilled water (5 min). Finally, sections were dehydrated in graded ethanol (60%, 80%, 90%, and 100%; 5 min each), cleared in xylene (twice, 5 min each), mounted with neutral gum, and examined under a microscope.The slides were scanned on the tissue microarray scanner to quantitatively analyze the expression level of CXCL8 .The immunohistochemical staining was evaluated independently by two experienced pathologists who were blinded to the clinical data. The staining intensity was scored as follows: 0 (negative), 1 (weak), 2 (moderate), and 3 (strong). The proportion of positively stained tumor cells was scored as: 0 ( 75%). The final immunoreactivity score (IRS) for each case was calculated by multiplying the intensity score by the proportion score, resulting in a range from 0 to 12. An IRS score ≥ 4 was defined as high expression, and < 4 was defined as low expression for statistical analysis. Statistical analyses were conducted using SPSS version 28.0 (IBM Corp., USA) and GraphPad Prism version 10.0 (GraphPad Software, USA). Data conforming to a normal distribution were expressed as mean ± standard deviation (mean ± SD). For comparisons of CXCL8 mRNA expression levels between paired PTC tissues and PNTs, the non-parametric Wilcoxon signed-rank test was employed due to the sample size and potential deviations from normality assumptions. For the correlation analysis between CXCL8 expression and clinicopathological features, the following specific handling was applied to the age variable: age was treated as a continuous variable and presented as median [interquartile range (IQR)]. This approach was adopted to avoid information loss and potential bias caused by arbitrary dichotomization (e.g., dividing age into ≤ 55 years/>55 years), which could oversimplify the age distribution and obscure its nuanced relationship with CXCL8 expression. For other categorical clinicopathological features (e.g., gender, lymph node status, Hashimoto’s thyroiditis), correlations with CXCL8 expression were analyzed using Pearson’s χ² test; when any expected cell count was < 5, Fisher’s exact test was applied instead. For other categorical clinicopathological features (e.g., gender, lymph node status, Hashimoto’s thyroiditis), correlations with CXCL8 expression were analyzed using Pearson’s χ² test; when any expected cell count was < 5, Fisher’s exact test was applied instead to ensure statistical validity under sparse data conditions. All experiments were independently repeated three times, and only data consistent across two or more repetitions were subjected to statistical evaluation to ensure reproducibility and minimize random error. A P -value < 0.05 was considered statistically significant, and all tests were two-tailed unless otherwise specified.

PTC is the most common pathological type of thyroid malignancy, particularly prevalent in women 13 , 14 . PTC typically grows slowly and has a favorable prognosis, with a key characteristic being its ability to invade adjacent structures, such as lymphatic vessels. Approximately 10% of patients may present with metastatic disease at the time of diagnosis 15 . In recent years, numerous PTC-associated risk genes have been identified, which significantly affect the proliferation, migration, and invasion of PTC cells. Chemokine ligand 8 ( CXCL8 ), also known as interleukin-8 (IL-8), is a member of the CXC family of chemokines 16 , 17 . This widely expressed cytokine plays a critical role in mediating cellular inflammatory responses and regulating autoimmune diseases within the tumor microenvironment. CXCL8 is upregulated in various types of malignant tumors 18 . During tumor progression, CXCL8 promotes the proliferation, migration, invasion, angiogenesis, and epithelial-to-mesenchymal transition (EMT) of tumor cells through autocrine or paracrine signaling mechanisms 7 , 19 . Pan-cancer expression analysis of CXCL8 revealed higher expression levels in tumor samples of several cancers, including CHOL, COAD, ESCA, HNSC, KIRP, READ, STAD, and UCEC. Additionally, data from the TIMER2.0 database showed that CXCL8 expression in thyroid cancer samples was significantly higher than in PNTs. Bioinformatics analysis further confirmed that CXCL8 levels were generally elevated in tumor tissues compared to PNTs. To validate these findings, we performed qRT-PCR experiments, which demonstrated significantly increased CXCL8 mRNA expression in PTC tissues compared to PNTs.To investigate the relationship between CXCL8 expression and clinical-pathological features, we conducted tissue microarray experiments using samples from 49 PTC patients. The results revealed that CXCL8 expression was significantly associated with lymph node metastasis and Hashimoto’s thyroiditis, but not with gender, age, tumor size, tumor location, recurrence risk, capsule invasion, or extrathyroidal invasion. The increased concentration of CXCL8 detected in cancerous tissues contributes to enhanced cancer cell proliferation and migration. A notable association between CXCL8 expression and lymph node metastasis suggests that CXCL8 could serve as a potential biomarker for both tumor progression and patient survival in PTC. Gu et al. showed that CXCL8 expression was significantly correlated with clinical stage, distant metastasis, histological type, and grade in cervical cancer, with high CXCL8 expression serving as an independent poor prognostic marker for cervical cancer patients 20 . Hashimoto’s thyroiditis, a common organ-specific autoimmune disease, is characterized by the presence of various cytokines, including IL-1α, IL-1β, IL-2, IL-3, and IL-8, which amplify the inflammatory response via nitric oxide (NO) and prostaglandins 21 . In our study, the Circos diagram illustrated the co-expression network of CXCL8 and 11 genes in PTC samples. We conducted GO and KEGG pathway analyses of differential proteins, leading us to hypothesize the following mechanisms: CXCL8 exhibits significant co-expression with CXCL1 , CXCL2 , and CXCL3 , binding to CXCR1 or CXCR2 , triggering multiple G-protein-mediated signaling cascades, including the PLC-PKC, PI3K - Akt , Src, and FAK pathways. These pathways promote tumor cell proliferation, survival, angiogenesis, EMT , migration, and invasion 18 , 22 . In our single-cell data, the gene CXCL8 was found to be specifically highly expressed in Dendritic cells. AKT1 in the PI3K-AKT pathway was also highly expressed in Dendritic cells, which provided a basis for the malignant metastasis of CXCL8 in PTC, and our KEGG analysis results also supported CXCL8 as a potential regulator of PI3K-AKT pathway.The PI3K - Akt signaling pathway, crucial for cell membrane receptor signaling, plays an essential role in inflammatory and immune responses 23 . Zhao et al. demonstrated that blocking the PI3K - Akt pathway eliminated IL-8-induced oxaliplatin resistance in gastric cancer derived from cancer-associated fibroblasts (CAFs) 23 . Furthermore, CXCL8 / CXCR1 and/or CXCL8 /CXCR2 signaling induces PI3K expression by activating a cascade that leads to Akt phosphorylation, promoting cell differentiation, proliferation, and survival 24 . Another important pathway is the MAPK cascade, where CXCL8 activates the RAF/ MAPK / ERK pathway 25 . The FAK/IL-8 axis promotes gastric cancer cell proliferation and migration 7 . PELI1, a significant transcriptional target, cooperates with CXCL8 in activating the NF-κB pathway, facilitating tumor metastasis, and driving cell proliferation and migration 26 . High mRNA expression of GPER1 , EGFR , and CXCR1 in PTC correlates significantly with lymph node metastasis (LNM), which may help predict LNM risk 23 . The tumor immune microenvironment becomes active through paracrine, autocrine, and endocrine communication, generating dynamic signals that drive cancer progression 27 , 28 . It is widely accepted that dysfunction in the tumor-immune ecosystem impairs immune surveillance 29 . Using the “CIBERSORT” tool, we performed differential analysis of immune cell infiltration levels between high- and low-expression CXCL8 groups, confirming a significant correlation between CXCL8 and tumor-infiltrating immune cells. The infiltration levels of activated dendritic cells, neutrophils, resting dendritic cells, M 0 macrophages, and activated mast cells positively correlated with CXCL8 , while CD8 + T cells, plasma cells, M1 macrophages, and γδ T cells negatively correlated with CXCL8 . CXCL8 promotes the migration of CD4 + T cells, leading to immune imbalance and mediating pro-inflammatory functions 6 , 30 . In triple-negative breast cancer, CXCL8 inhibits the infiltration of CD4 + and CD8 + T cells and reshapes the tumor immune microenvironment by upregulating CD274 expression 31 . In endometriosis, CXCL8 secretion from endometrial stromal cells impairs the cytotoxic activity of NK cells, contributing to local immune dysfunction 19 . M 0 macrophages promote tumor cell proliferation, migration, and invasion, with CXCL8 being a key regulator of these processes. Zhao et al. found that CXCL8 regulates the proliferation and polarization of M 0 macrophages in cervical cancer progression 32 . DCs initiate adaptive immunity and regulate inflammation by producing inflammatory chemokines 33 , 34 . Infected DCs produce high levels of CXCL8 , further stimulating neutrophil activity, contributing to granuloma formation in cat scratch disease 35 . Mast cells, as multifunctional immune cells, contribute to both innate and acquired immunity. Activated mast cells release various pro-inflammatory mediators, including CXCL8 , which may drive mast cell migration 35 , 36 . CXCL8 and related molecules (CXC family) primarily exert chemotactic effects on neutrophils. By recruiting macrophages, neutrophils, and T cells to the tumor site, CXCL8 plays a pivotal role in tumor invasiveness and can inhibit anti-tumor immune responses 37 , 38 . Based on the immunosuppressive role of CXCL8 in PTC, targeting the CXCL8 pathway 39 , 40 in combination with PD-1 blockade enhances the tumor immune response and inhibits tumor progression 41 . CXCL8 is primarily secreted by macrophages, contributing to the immunosuppressive tumor microenvironment by inducing PD-L1 + macrophages in gastric cancer 42 . Danlos et al. found that mesothelioma patients with primary resistance to PD-1-targeted immunotherapy exhibited elevated CXCL8 levels in both blood and tumors 38 . Blockade of CTLA4 has been shown to reduce immunosuppressive cells such as myeloid-derived suppressor cells (MDSC) and M2 macrophages, enhancing T-cell activation and presenting a new therapeutic target for PTC 43 , 44 .

We detected the expression of CXCL8 in PTC tissues by bioinformatics, qRT-PCR and tissue microarray experiments. The results indicated that CXCL8 plays an important role as a potential biomarker in the biological functions of PTC and can be used as a new target for drug development, providing a possible new direction for the immunotherapy of PTC.

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1