We first downloaded the RNA-seq expression data and the survival data of TCGA-HNSC cohort from UCSC Xena website (https://xenabrowser.net/datapages/). FPKM values of expression matrix were then transformed into TPM values, and Ensembl IDs were annotated by the gencode annotation file (version 22) which was also downloaded in UCSC Xena. The pan-cancer expression data of CD73 protein was collected and analyzed in the Human Protein Atlas (HPA) database (9, 10) by using the “HPAanalyze” (11) R package. Somatic mutation data processed with Mutect software were downloaded in the TCGA GDC website (https://portal.gdc.cancer.gov/). Single-cell RNA-seq data were obtained and visualized by TISCH (http://tisch.comp-genomics.org/) (12) and CancerSEA (http://biocc.hrbmu.edu.cn/CancerSEA/) (13) database.

Formalin-fixed, paraffin-embedded tumor tissues and their paired peritumor or normal tissues were collected from the Tumor Sample Repository of Sun Yat-Sen University Cancer Center (SYSUCC), with a total of 15 HNSCC tissues and 14 peritumor or normal tissues. These tissues were cut to a thickness of 4 μm sequential sections and sliced for immunohistochemical staining and analysis. This study complies with the Declaration of Helsinki and was approved by the ethics committee of Sun Yat-Sen University Cancer Center (Approval Number: B2022-507-01). All patients obtained informed consent.

The slides were sequentially baked at 60°C for 3 hours, deparaffinized with xylene, and rehydrated with graded ethanol. Antigen retrieval was performed after blocking endogenous peroxidase with 3% hydrogen peroxide. Next, slides were blocked with 10% fetal bovine serum at room temperature for 30 minutes, and then incubated with primary antibody including CD73 Rabbit mAb (Cell Signaling Technology, #13160, 1:400) or CD8α Rabbit mAb (Cell Signaling Technology, #85336, 1:400) overnight at 4°C. The next day, slides were incubated with anti-mouse/rabbit HRP-labeled secondary antibody at 37°C for 30 minutes in the dark environment. Subsequently, diaminobenzidine and hematoxylin were used for staining and counterstaining, respectively, and then the slides were quickly immersed in 1% hydrochloric acid alcohol for 1 second. Finally, the slides were dehydrated in reverse order through graded alcohol and xylene, and sealed with resin.

Interpretation of immunohistochemical results was conducted by Dr. Qinian Wu, Department of Pathology. Briefly, immunohistochemical scores of CD73 were analyzed by semi-quantitative method. That is: the staining intensity is 0 (negative), 1 (weakly positive), 2 (moderately positive) and 3 (strongly positive), separately; and the staining area is 1 (<25%), 2 (26%-50%), 3 (51%-75%), 4 (>75%), separately. Multiplying the staining intensity score and the staining area score to obtain the final score for each slide. The score of CD8 is to evaluate the proportion of CD8-positive cells to lymphocytes in the tumor stroma.

First, we divided the samples into CD73-high group and CD73-low group according to the median expression value of CD73, and then we used the “limma” (14) R package for differential expression gene (DEG) analysis. All genes were ranked according to the logFoldChange (logFC) values ​​obtained from the differential expression analysis, and GSEA analysis was performed using the GSEA function in the “clusterProfiler” (15, 16) R package and the HALLMARK gene set in MSigDB (17) database. Gene set permutations were performed 1000 times. GSEA results were visualized using the ridgeplot function in the “clusterProfiler” (15, 16) R package and the gseaplot2 function in the “enrichplot” R package.

We assessed immune infiltration in HNSCC samples using three different methods. First, we employed the “Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data” (ESTIMATE) algorithm (18), which scores each tumor sample for immune, stromal and tumor purity. The R package “estimate” was used to perform this analysis. Next, in order to further evaluate the specific immune cell infiltration in each tumor sample, we further adopted the CIBERSORT algorithm (19) to quantitatively evaluate the infiltration proportions of 22 types of immune cells in HNSCC tumor samples. Furthermore, Microenvironment Cell Populations (MCP)-counter algorithm (20) was also performed by “MCP-counter” R package.

As described previously, the somatic mutation data of HNSCC patients called by Mutect software were downloaded from TCGA GDC website. According to the median expression of CD73, we divided the group into high expression group and low expression group. We used the subsetMaf function in the “maftools” (21) R package to extract the somatic mutation subsets of the two groups, and used the tmb function to calculate the tumor mutation burden (TMB) of the two groups of patients for comparison. In addition, we also used the mafCompare function to compare the mutant genes of the two groups, and the parameter minMut (which means considering only genes with minimum this number of samples mutated in at least one of the cohorts for analysis) was set as 20.

We performed mutational signature analysis through the pipeline of the “maftools” (21) R package. First, we used trinucleotideMatrix function to extract single 5’ and 3’ bases flanking the mutated site of each subset maf files for de-nove signature analysis, and the reference genome was set as “BSgenome.Hsapiens.UCSC.hg38”. Next, we used default parameter of estimateSignatures function and NMF algorithm (22) to estimate the number of signatures, and used compareSignatures function to compare these signatures with known updated/refined 65 COSMIC signatures (23).

First, we used the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (24) to predict the immune efficacy response of HNSCC samples. According to the instructions on the TIDE website, we first performed the median centering of the expression matrix in two directions, and submitted the matrix to the TIDE website for calculation. We then extracted the TIDE score, MSI score, T cell exclusion score and T cell dysfunction score from the output results to compare the CD73 high expression group and the CD73 low expression group. Higher TIDE score means the increased potential to develop immune escape. Next, to further confirm the predicted results of TIDE, we performed validation analysis using the IMvigor210 cohort, a locally advanced or metastatic urothelial carcinoma (mUC) cohort receiving anti-PD-L1 therapy (25). The expression data and clinical information of this cohort are integrated in the “IMvigor210CoreBiology” R package, from which we extracted the expression of CD73 and the efficacy evaluation information of immunotherapy for analysis and verification. Besides, we calculated the correlation between the expression of CD73 and the expression of immune checkpoints (PDCD1, CD274, PDCD1LG2, HAVCR2, TIGIT, LAG3 and CTLA4) or CD8⁺ T cell marker (CD8A) in TCGA-HNSC cohort, using the “ggstatsplot” R package for statistical analysis and visualization.

All statistical analyses were conducted by the R software V4.1.2 (The R Project for Statistical Computing, https://www.r-project.org/). The comparisons between two groups were performed by using the Wilcoxon test for statistical analysis. The correlations of gene expression were analyzed using the Pearson correlation coefficient. Kaplan-Meier survival analysis used the “survival” and “survminer” R packages for statistics and visualization, and the log-rank test was used to compare survival differences between the two groups. The receiver operating characteristic (ROC) curve was used to analyze the diagnostic value of CD73 in populations, and this was calculated and plotted using the “pROC” (26) R package. P < 0.05 was considered as significant.

We first analyzed the mRNA expression difference of CD73 between the normal samples (n = 44) and tumor samples (n = 502) in TCGA-HNSC cohort. Compared to normal samples, the mRNA expression of NT5E (encoding CD73 protein) is significantly elevated in tumor samples (P = 9.7×10^-10) ( Figure 1A ). To further validate the expression of CD73 at protein level, we performed the immunohistochemistry assay in SYSUCC cohorts and collected the IHC data in Human Protein Atlas (HPA) database, respectively. In consistent with the mRNA expression result, we also found that CD73 protein is highly expressed in HNSCC samples (P = 7.6×10^-5), especially with high/medium expression ( Figures 1B, C and Supplementary Figure 1 ). Next, the prognostic and diagnostic value of CD73 in HNSCC were also determined. Kaplan-Meier analysis revealed that patients with higher CD73 expression had worse overall survival (log-rank P = 0.0094) ( Figure 1D ). Besides, we evaluated the diagnostic role of CD73 in HNSCC through ROC curve, and the area under the curve was 0.778 ( Figure 1E ). In summary, all these results suggest that CD73 is elevated in tumors, and its overexpression indicates a dismal survival for HNSCC patients.

To elucidate the concrete role and mechanism how CD73 exerts in HNSCC tumors, we employed the GSEA analysis to calculate the enrichment score of hallmarks of cancer pathways ( Figure 2A ). We observed many oncogenic pathways, such as KRAS signaling pathway, angiogenesis pathway, TGF-β pathway and hypoxia pathway; and immune-related pathways, such as IFN-α pathway, inflammatory response and complement pathway, were more enriched in CD73-high group ( Figure 2A ). Meanwhile, the top 2 pathways are epithelial-mesenchymal transition (EMT) pathway ( Figure 2B ) and interferon-gamma (IFN-γ) response ( Figure 2C ), which further illustrate the oncogenic and immune-regulatory role of CD73. Therefore, these results suggest that CD73 is an immune-related oncogene in HNSCC.

As mentioned above, we speculate that CD73 could regulate tumor immune microenvironment. To confirm this speculation, we applied multiple methods to evaluate the immune infiltration status of tumor samples. As depicted in Figures 2D-G , we found that the stromal scores ( Figure 2D ) and the ESTIMATE scores ( Figure 2E ) are both higher in CD73-high group than in CD73-low group (P = 1.6×10^-7, P = 0.00046, respectively). Accordingly, the tumor purity is lower in CD73-high group when compared to CD73-low group (P = 0.00046) ( Figure 2F ). No significant difference of immune score was found between two groups (P = 0.31) ( Figure 2G ). In addition, we also assessed the cell infiltration status in the tumor microenvironment by two different algorithms, CIBERSORT and MCP-counter ( Figures 2H, I ). The CIBERSORT algorithm analysis found that compared with the low CD73 expression group, both the B cell lineages and T cell lineages had less infiltration in the samples with high CD73 expression (Fighter 2H). Further analysis with the MCP-counter algorithm also confirmed this phenomenon (Fighre 2I). Interestingly, we found that in both algorithms, CD8⁺ T cells decreased in the CD73-high group ( Figures 2H, I ). CD8⁺ T cells are the cells that mainly exert anti-tumor function in the tumor microenvironment. When the infiltration of CD8+ T cells is reduced, the anti-tumor immune response will be correspondingly weakened. Therefore, it may partially explain the inhibitory mechanisms of CD73 on tumor progress. Furthermore, we also observed a higher proportion of fibroblasts and endothelial cells in CD73-high group ( Figure 2I ). Collectively, all these data suggest the immunosuppressive role of CD73 in regulating the tumor microenvironment.

To exactly specify which cell types express CD73 in HNSCC tumor microenvironment, we performed a single-cell RNA-seq analysis by using GSE103322 data. After dimension reduction and cell annotation, the cells can be divided into 11 clusters ( Figure 2J ). NT5E, encoding the CD73 protein, is mainly expressed in the fibroblasts, endothelial and malignant cells ( Figure 2K ), which was consistent with the result of higher stromal scores and higher proportions of these cells in CD73-high group ( Figures 2D, I ). Furthermore, we also evaluate the correlation between CD73 expression and tumor signatures. The result showed that CD73 expression is considerably positive correlate to EMT, metastasis and invasion signatures ( Figure 2L ), which was same as the results from bulk RNA-seq ( Figure 2B ). Therefore, scRNA-seq data further verified the expression and oncogenic role of CD73.

To further explore the possible mechanism of CD73, we divided the somatic mutation data of TCGA-HNSC into two groups according to the median expression of CD73. Figure 3A depicted the landscape of somatic mutations in two groups. We further compared differences in somatic mutations between the two groups ( Figure 3B ). TP53, HRAS and CDKN2A had more mutations in CD73-high group, and genes like FBXW7 and MUC16 had more mutations in CD73-low group ( Figure 3B ). In addition, we found that the presence of somatic mutations co-occurring is considerably lower in CD73-high group compared to the CD73-low group ( Figures 3C, D ). TMB was also lower in CD73-high group (P = 0.0055) ( Figure 3E ), with a median of 1.73/MB ( Supplementary Figure 2A ), while the TMB in the CD73-low group was a median of 2/MB ( Supplementary Figure 2B ). Furthermore, we also deconvolute the mutational signature to investigate the potential risk factors corresponding to the two groups ( Figures 3F–I ). In the CD73-high expression group, we extracted a total of 5 mutational signatures ( Supplementary Figure 3A ). According to the cosine similarity, these 5 mutational signatures are corresponding to SBS13 (APOBEC cytidine deaminase), SBS 7b (UV exposure), SBS16 (Unknown), SBS4 (exposure to tobacco/smoking mutagens) and SBS1 (spontaneous or enzymatic deamination of 5-methylcytosine) of COSMIC database, respectively ( Figures 3F, G ). Meanwhile, we also extracted 2 signatures in CD73-low expression group ( Supplementary Figure 3B ), and these two signatures were identified as SBS7a (UV exposure) and SBS3 (Defects in DNA-DSB repair by HR), respectively ( Figures 3H, I ).

As mentioned above, we found that CD73-high group has a considerably lower TMB than CD73-low group ( Figure 3E ). Given that high TMB is more effective for immunotherapy, we therefore speculate that patients with high CD73 expression are more likely to be unresponsive to immunotherapy. To validate this hypothesis, we calculated several metrics related to the prediction of immunotherapy efficacy to reflect the potential clinical response of CD73. According to the prediction scores of TIDE algorithm, CD73-high group had a higher TIDE score (P = 2.6×10^-11) ( Figure 4A ), lower MSI score (P = 0.00034) ( Figure 4B ) and high T cell exclusion score (P = 2.3×10^-11) ( Figure 4C ), indicating that higher CD73 expression was more likely to generate immune escape and thus be resistant to immunotherapy. However, we did not observe any correlations between CD73 and T cell dysfunction score, which suggests that CD73 may not influence the T cell exhaustion (P = 0.48) ( Figure 4D ). We performed immunohistochemistry assay to further validate the relationship between the expression of CD73 and the infiltration of CD8⁺ T cells. The results showed that the infiltration of CD8⁺ T cells was relatively less in tumor tissues with high CD73 expression, while the infiltration of CD8⁺ T cells was more in tumor tissues with low CD73 expression ( Figures 4E, F ). Although the relationship between CD73 and CD8⁺ T cells did not reach statistical significance (P = 0.07), this may be due to the small sample size ( Figure 4F ). In addition, we used a cohort of metastatic urothelial carcinomas (IMvigor210 cohort) that had actually received anti-PD-L1 immunotherapy to validate this result, and CD73 expression was indeed lower in the group that responded to immunotherapy (P = 0.0025) ( Figure 4G ). We also performed a correlation analysis between the expression of immune checkpoints (PDCD1, CD274, PDCD1LG2, HAVCR2, TIGIT, LAG3 and CTLA4) or CD8⁺ T cell marker (CD8A) and NT5E, and the results showed that these checkpoints were not strongly correlated with NT5E expression levels ( Supplementary Figures 4A–H ). To sum up, as illustrated in Figure 5 , CD73 could be act as a predictor for immunotherapy, and increased CD73 expression indicates less clinical benefit for immunotherapy.