CD8A as a Prognostic and Immunotherapy Predictive Biomarker Can Be Evaluated by MRI Radiomics Features in Bladder Cancer

Simple Summary It is necessary to explore a reliable predictive biomarker for immunotherapeutic strategies in bladder cancer. As an important member of tumor-infiltrating immune cells, cytotoxic T cells are the ultimate effectors of cancer rejection because of their specificity against cancer cells and cytolytic capabilities. Based on RNA expression data, we found that among T cytotoxic pathway-related genes, CD8A is a novel protective gene in bladder cancer. CD8A expression is highly associated with tumor microenvironment characteristics and tumor mutation burden, which partially explain CD8A as a prognostic and immunotherapy predictive biomarker in bladder cancer. In addition, we performed the radiogenomics in bladder cancer based on the RNA-sequence data and MRI-based radiomics features in our center. In this way, an MRI-based radiomics signature has the potential to preoperatively predict the expression of CD8A, thereby contributing to the preoperative prediction of prognosis and immunotherapeutic susceptibility in bladder cancer. Abstract As an important member of T cytotoxic pathway-related genes, CD8a molecule (CD8A) may be a useful biomarker of immunotherapeutic response and immune cell infiltration. We aimed to investigate the clinical predictive value of CD8A in prognosis and tumor microenvironment (TME) and preoperatively predict the expression of CD8A using radiogenomics in bladder cancer (BCa). Among 12 T cytotoxic pathway-related genes, CD8A was a novel protective gene and had the highest correlations with T cells and Macrophages M1 in BCa. In advanced cancer patients treated with immunotherapy, low CD8A expression was associated with immunotherapeutic failure and poor survival outcomes. CD8A expression was highly related to tumor mutation burden, critical immune checkpoint genes and several types of tumor-infiltrating immune cells, predicting effective response to immunotherapy. The preoperative MRI radiomics features and RNA-sequence data of 111 BCa samples were used to develop a radiomics signature that achieved good performance in the prediction of CD8A expression in both the training (area under curve (AUC): 0.857) and validation sets (AUC: 0.844). CD8A is a novel indicator for predicting the prognosis and immunotherapeutic response in BCa. A radiomics signature has the potential to preoperatively predict the expression of CD8A in BCa patients.


Introduction
Bladder cancer (BCa) is one of the most malign neoplasia in the world, with high morbidity and mortality [1]. Approximately 30% of BCa patients have muscle-invasive bladder carcinoma (MIBC) at the time of diagnosis with a high rate of metastasis and mortality [2,3], and the other non-muscle-invasive bladder carcinoma (NMIBC) patients are characterized by a high rate of recurrence and a certain risk of progression to MIBC [4]. Treatment strategies differ greatly between NMIBC and MIBC. NMIBC patients are treated with transurethral resection followed by intravesical therapy with Bacillus Calmette-Guérin (BCG) or intravesical chemotherapy for carcinoma in situ or high-grade T1 [5], while MIBC patients (>T1) are treated with cisplatin-based neoadjuvant chemotherapy and radical cystectomy [6]. Despite the development in adjuvant therapy and surgical techniques, some BCa patients do not benefit from the recommended treatment strategies, with only an 8% 5-year survival rate in BCa patients with metastasis [7]. Therefore, it is urgently needed to investigate new and reliable therapeutic approaches for BCa patients.
Immune checkpoint inhibitor (ICI) is becoming a useful and important immunotherapy in BCa [8]. It is reported that Nivolumab (programmed cell death protein 1 (PD-1) inhibitor) could induce neoplastic cell death in urothelial bladder cancers with metastases [9]. In addition, immunotherapy had been reported to be a useful treatment for both early and late stages of BCa [10,11]. However, the results of clinical trials revealed that only a small fraction of BCa patients benefit from immunotherapy [8,12,13]. The role of inhibitory receptors expression such as programmed death ligand-1 (PD-L1) as a predictor of immunotherapeutic response had not been determined [8]. Thus, it is necessary to explore useful immunotherapy targets to identify subgroups of BCa patients that benefit from immunotherapy.
As members of tumor-infiltrating immune cells (TIICs), cytotoxic T cells (CD8+ T cells) are recruited to the tumor microenvironment (TME) to play an important role in tumor adaptive immunity during tumor immune responses. Cytotoxic T cells are the ultimate effectors of cancer rejection because of their specificity against cancer cells and cytolytic capabilities [14]. The CD8a molecule (CD8A) is a member of T cytotoxic pathway-related genes and encodes the CD8 antigen that is a cell surface glycoprotein found on most cytotoxic T cells. The CD8 antigen acts as a coreceptor with the T-cell receptor on the T cell to recognize antigens displayed by an antigen-presenting cell in the context of class I MHC molecules. CD8A expression may be a useful and measurable predictive marker of immunotherapeutic response and immune cell infiltration [15]. A previous study for pan-cancer has reported that high status of CD8A with high expression of PD-L1 might be a predictive marker of immunotherapeutic response [16]. However, the prognostic value of CD8A and the association between CD8A expression and immunotherapeutic response have not been explored in BCa.
Magnetic resonance imaging (MRI) is emerging as a routine examination for BCa patients in clinical practice. It can intuitively demonstrate the size, number, location and shape of lesions. The Vesical Imaging Reporting and Data System (VI-RADS) is a fivepoint assessment scale based on multiparametric MRI for describing the likelihood of muscle invasion [17]. It has been consistently validated across multiple institutions as a useful tool for local staging in BCa. However, the VI-RADS is generated by the radiologist's subjective evaluation, which limits its application in predicting the biological and molecular characteristics of tumors. In contrast, an emerging technique termed radiogenomics may be a useful tool in clinical oncology. It can extract high-throughput quantitative imaging features not perceptible to the human eye from digital medical images and provide an objective and non-invasive approach in uncovering the TME characteristics in a tumor. A previous study developed a CT-based radiomics signature to predict the CD8 cell's tumor infiltration and immunotherapeutic response in patients with advanced solid tumors [18]. Preoperative evaluation of TME characteristics using radiogenomics may be a useful way to facilitate survival prediction and personalized immunotherapy. At present, there are no studies using radiogenomics to predict the CD8A expression in BCa.
Therefore, the aim of this study was to assess the clinical predictive value of CD8A in prognosis and immunotherapeutic response, investigate the association between CD8A and TME and construct an MRI-based radiomics signature for the prediction of CD8A expression using radiogenomics in BCa.  a useful way to facilitate survival prediction and personalized immunotherapy. At present, there are no studies using radiogenomics to predict the CD8A expression in BCa. Therefore, the aim of this study was to assess the clinical predictive value of CD8A in prognosis and immunotherapeutic response, investigate the association between CD8A and TME and construct an MRI-based radiomics signature for the prediction of CD8A expression using radiogenomics in BCa.   For BCa datasets from the Gene Expression Omnibus (GEO, https://www.ncbi. nlm.nih.gov/, accessed on 6 April 2022) and The Cancer Genome Atlas (TCGA, http: //cancergenome.nih.gov/, accessed on 10 April 2022), patients who met the following inclusion criteria were selected: (1) histologically diagnosed with BCa; (2) available RNA expression data; (3) available prognostic information. In this study, 405 patients from TCGA, 165 patients from GSE13507, 93 patients from GSE31684, 73 patients from GSE48075, 224 patients from GSE32894 and 73 patients from GSE48277 were included. In addition, two datasets (IMvigor210C and GSE93157) with advanced cancer patients treated with immunotherapy were also used to investigate the potential clinical value of CD8A in classifying the immunotherapy susceptibility for patients. IMvigor210 dataset consisted of patients with metastatic urothelial cancer treated with atezolizumab (PD-L1 inhibitor) [19], and the RNA expression data and prognostic information were downloaded from the website http://research-pub.gene.com/IMvigor210CoreBiologies, accessed on 10 April 2022. GSE93157 dataset included 65 patients with melanoma, lung cancer, and head and neck cancer treated anti-PD1 [20]. Supplementary Table S1 summarizes the detailed information on these datasets.

Data Acquisition
In Shanghai Tenth People's Hospital, BCa patients from November 2019 to July 2021 were collected to develop MRI-based radiomics signature for the prediction of CD8A expression. The inclusion criteria in our center included the following: (1) histologically diagnosed with BCa; (2) available RNA-sequence data; (3) available MRI within 20 days before surgery. The exclusion criteria were as followed: (1) poor-quality MRI images; (2) any treatments were performed before MRI examination; (3) tumors for which it was difficult to define the boundaries; (4) lack of baseline clinical factors. The protocols of MRI acquisition, total RNA extraction, paired-end libraries generation and RNA-sequence can be obtained in our previous methods [21,22].

HEmatoxylin and Eosin (H&E) and Immunohistochemistry (IHC) Staining
After fixation, BCa tumor tissues were embedded in paraffin using standardized cassettes, sectioned at 5 µm, and stained with H&E and IHC as previously described in detail [23,24]. Rabbit anti-human monoclonal primary antibodies against CD8A (CST, cat# 85336) were utilized to detect CD8A expression according to the manufacturer's protocol.

Evaluation of TIICs and TME
CIBERSORT analysis was performed to evaluate the relative proportion of 22 TIICs based on the RNA expression data of marker genes in each tumor sample (https://cibersort. stanford.edu, accessed on 10 April 2022). The TCGA dataset that has the maximum number of BCa patients was used to investigate the different TIICs between low and high CD8A expression groups. The relative proportion of 22 TIICs was also evaluated in our center based on the RNA-sequence data. The TME scores (including immune score, stromal score and estimate score) of each BCa patient were calculated using the R package "ESTIMATE".

Association of CD8A with TIICs and Critical Immune Checkpoint Genes
We used the tumor immune estimation resource (TIMER, https://cistrome.shinyapps. io/timer/, accessed on 10 April 2022) to assess the correlation between CD8A and some TIICs including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophil and dendritic cells. In addition, the correlation between CD8A and critical immune checkpoint genes including PDCD1, CTLA4, LAG3, TIGIT, HAVCR2, CD274 and TNFRSF9 was also investigated using TIMER.

Association of CD8A with Tumor Mutation Burden (TMB), Stemness and Microsatellite Instability (MSI)
The UCSCXenaShiny (https://hiplot.com.cn/advance/ucsc-xena-shiny, accessed on 10 April 2022) was used to explore the association of CD8A with TMB, stemness and MSI among 32 cancer types. In addition, the scatter plot showing the correlation between CD8A expression and TMB in TCGA-BCa was also performed.

Association between CD8A and Drug Sensitivity
The drug sensitivity of each BCa sample was predicted based on the Genomics of Drug Sensitivity in Cancer, which is the largest public pharmacogenomics database. Seven drugs for anti-cancer treatment were selected to evaluate the IC50 of each sample using the "pRRophetic" R package.

Survival Analysis and Meta-Analysis
The prognostic value of CD8A expression among 33 types of tumors in TCGA was analyzed using UCSCXenaShiny (https://hiplot.com.cn/advance/ucsc-xena-shiny, accessed on 10 April 2022). The Kaplan-Meier and log-rank tests were used to evaluate the prognostic value of CD8A expression in TCGA, IMvigor210 and six GEO datasets (GSE93157, GSE13507, GSE31684, GSE48075, GSE32894 and GSE48277). Patients were divided into low and high CD8A expression groups based on the optimal cutoff values in each dataset calculated by the "survminer" R package v0.4.8 (Marseille, France). Univariate and multivariate Cox regression analyses were used to obtain the independent prognostic factor in TCGA. In addition, the results of univariate Cox regression analyses (including HR and 95% confidence interval (CI)) in seven BCa datasets (TCGA, IMvigor210, GSE13507, GSE31684, GSE48075, GSE32894 and GSE48277) were used for meta-analysis. I 2 statistic and χ 2 -based Q test were used to evaluated the heterogeneity of meta-analysis. If I 2 value > 25%, the meta-analysis was regarded as high heterogeneity and the random-effect model was used for meta-analysis. The funnel plot, Begg's test and Egger's test were applied to investigate the possible publication bias. A p value > 0.05 was regarded as no publication bias.

Construction of the Nomogram
Factors that were significant in multivariate Cox analysis were used to construct a predictive nomogram. Time-dependent receiver operating characteristic (ROC) curves of the nomogram were used to evaluate its performance in survival prediction. Decision curve analysis (DCA) and calibration plots were performed to investigate the clinical utility and performance of the nomogram, respectively.

Region of Interest (ROI) Segmentation and Radiomics Feature Extraction
A radiologist (F Xu) with over 5 years of experience in bladder MRI reading used the ITK-SNAP software (version 3.6.0; http://itk-snap.org, accessed on 10 April 2022) (Kitware Inc., Clifton Park, NY, USA) to manually outline the ROIs of lesions layer by layer on T2WI images and the delay phase of DCE images. If a BCa patient had multiple lesions, the largest lesion was uniquely outlined for feature extraction. After 30 days, the lesions of 40 randomly selected patients were repeatedly delineated by the same radiologist (F Xu) and another radiologist (T Xu) to evaluate the interobserver reliability.
Radiomics features were extracted from the ROI of each BCa using the Python package Pyradiomics v.  16 GLSZM features and 5 NGTDM features were calculated on the exponential, logarithm, square, square root and gradient images, respectively. In total, 1781 radiomics features were extracted from DCE and T2WI images, respectively. Each radiomics feature was normalized with a Z-score before further analyses.

Feature Selection and Radiomics Signature Construction
BCa patients in Shanghai Tenth People's Hospital were divided into low and high CD8A groups based on the median value of CD8A expression. Then, these patients were randomly allocated into training set and validation set based on a 7:3 ratio.
The intra-and interclass correlation coefficients (ICCs) were calculated between two radiologists to assess the interobserver reliability of radiomics features. The minimum redundancy maximum relevance (mRMR) algorithm was performed to select the relevant and non-redundant features and rank the importance of features [25]. In this way, the top 10 features were obtained to develop radiomics signature for predicting the CD8A expression.
In this study, we used the least absolute shrinkage and selection operator (LASSO) algorithm to construct radiomics signature. The LASSO algorithm can remove redundant features and select features with non-zero coefficients for model construction. The radiomics score was then calculated for each BCa patient by summing the values of radiomics features weighted by their respective coefficients. The radiomics signature was developed via ten cross-validations in the training set to select the optimal values of λ and was assessed in the validation set. The area under the ROC curve (AUC), accuracy, sensitivity, specificity, negative predictive value (NPV) and positive predictive value (PPV) were used to assess the performance of the radiomics signature.

Statistical Analysis
The t-test, Wilcoxon test or one-way ANOVA was properly performed to assess association between different factors. Meta-analysis was performed using the "metafor" R package v2.4 and "meta" R package v4. 16. The "corrplot" R package v0.84 was used to generate correlation matrix. The "glmnet" R package v4.1 was used to develop the LASSO model. The "fmsb" R package v0.7 (Kobe, Japan) was used to plot the comprehensive performance of the LASSO model. SPSS 23.0 (SPSS, Armonk, NY, USA) and R v3.6.1 (https://www.r-project.org/, accessed on 10 April 2022) (R Statistical Foundation, Vienna, Austria) were employed to conduct statistical analyses. A two-sided p value < 0.05 was considered significant. Figure 2A demonstrates the locations of the T cytotoxic pathway-related genes on their respective chromosomes. Figure 2B,C shows the PPI network and correlation network among T cytotoxic pathway-related genes, which suggest that the expressions and functions of T cytotoxic pathway-related genes were highly correlated with each other. Among 12 T cytotoxic pathway-related genes, CD8A had the highest correlations with several TIICs, including T-cell CD8, T-cell CD4 naive, T-cell CD4 memory activated and Macrophages M1 ( Figure 2D,E). Univariate Cox analysis of OS and DFS showed that CD8A was a novel prognostic factor among 12 T cytotoxic pathway-related genes ( Figure 2F,G). In this way, CD8A was selected for further analysis. Patients with low CD8A expression were prone to have shorter overall survival (OS), disease-specific survival (DSS) and progression-free interval (PFI) than patients with high CD8A expression in this cohort of 33 cancer types (Supplementary Figure S1A-C). Then, we divided 8040 cancer patients with different types of cancers in TCGA into low CD8A expression and high CD8A expression groups based on the median value of CD8A gene expression and investigated the prognostic value of CD8A through the Kaplan-Meier method. The result showed that patients with low CD8A expression had significantly shorter OS than patients with high CD8A expression (Supplementary Figure S1D). In TCGA-BCa, patients with low CD8A expression had significantly shorter DFS than patients with high CD8A expression (Supplementary Figure S1E).

The Prognostic Value of CD8A
In three BCa datasets including TCGA, GSE48277 and GSE48075, BCa patients with low CD8A expression still had shorter OS than patients with high CD8A expression (Figure 3A-C). In TCGA, BCa patients in the low CD8A expression group had a significantly higher percentage of lymph node metastasis (p = 0.015, Supplementary Figure S2B) and death (p = 0.048, Supplementary Figure S2C). In two immunotherapy datasets including IMvigor210C and GSE93157, CD8A acted as a protective gene for cancer patients treated with immunotherapy ( Figure 3D,E). In addition, patients in the high CD8A expression Pan-cancer analysis showed that CD8A was highly associated with TMB, stemness and MSI among several tumor types including BCa (Supplementary Figure S1F-I). We further investigated the association between CD8A expression and IC50 values of seven drugs for anti-cancer treatment. The results showed that the IC50 values of Docetaxel, Camptothecin, Paclitaxel, Gemcitabine, FGFR inhibitor, Cisplatin and Tamoxifen were significantly lower in patients with high CD8A expression (Supplementary Figure S2A), which revealed that patients with high CD8A expression were more sensitive to these drugs.

The Prognostic Value of CD8A
In three BCa datasets including TCGA, GSE48277 and GSE48075, BCa patients with low CD8A expression still had shorter OS than patients with high CD8A expression ( Figure 3A-C). In TCGA, BCa patients in the low CD8A expression group had a significantly higher percentage of lymph node metastasis (p = 0.015, Supplementary Figure S2B) and death (p = 0.048, Supplementary Figure S2C). In two immunotherapy datasets including IMvigor210C and GSE93157, CD8A acted as a protective gene for cancer patients treated with immunotherapy ( Figure 3D,E). In addition, patients in the high CD8A expression group had a significantly higher percentage of partial response (PR)/complete response (CR) than those in the low CD8A expression group in IMvigor210 ( Figure 3F, 33.9% and 20.2%, respectively, p = 0.028). In addition, we performed meta-analysis and Cox regression analysis to comprehensively investigate the prognostic value of CD8A. The results of metaanalysis based on seven BCa datasets revealed that CD8A was a protective gene for BCa patients ( Figure 3G, Supplementary Figure S2D) without publication bias (Supplementary Figure S2E, Begg's test: p = 0.453, and Egger's test: p = 0.890). Univariate and multivariate Cox regression analysis showed that CD8A was an independent prognostic factor in BCa ( Figure 3H).
Cancers 2022, 14, 4866 9 of 23 group had a significantly higher percentage of partial response (PR)/complete response (CR) than those in the low CD8A expression group in IMvigor210 ( Figure 3F, 33.9% and 20.2%, respectively, p = 0.028). In addition, we performed meta-analysis and Cox regression analysis to comprehensively investigate the prognostic value of CD8A. The results of meta-analysis based on seven BCa datasets revealed that CD8A was a protective gene for BCa patients ( Figure 3G, Supplementary Figure S2D) without publication bias (Supplementary Figure S2E, Begg's test: p = 0.453, and Egger's test: p = 0.890). Univariate and multivariate Cox regression analysis showed that CD8A was an independent prognostic factor in BCa ( Figure 3H).  A CD8A-based nomogram was developed to predict the 1-, 3-and 5-year OS rates of BCa patients ( Figure 4A). ROC curves revealed that the nomogram had a good performance in survival prediction with the 1-, 3-and 5-year AUC values of 0.679, 0.722 and 0.722, respectively ( Figure 4B). Calibration curves revealed that the prognostic predictions of the CD8A-based nomogram were similar to the actual 1-, 3-and 5-year OS ( Figure 4C-E). DCA demonstrated that the CD8A-based nomogram had the optimal clinical net benefit compared to N stage and T stage ( Figure 4F-H).

Association of CD8A with TME Score, TIICs and Critical Immune Checkpoint Genes
The TME scores of low and high CD8A expression groups were calculated and compared. The results revealed that the high CD8A expression group had higher stromal scores, immune scores and TIDE scores than the low CD8A expression group ( Figure 5C), and CD8A expression was positively related to stromal scores, immune scores and TIDE scores (Supplementary Figure 3A-C). Correlation matrix demonstrated that CD8A was positively related to some TIICs (CD8+ T cell, CD4+ T cell, Macrophages M1) and critical

Association of CD8A with TME Score, TIICs and Critical Immune Checkpoint Genes
The TME scores of low and high CD8A expression groups were calculated and compared. The results revealed that the high CD8A expression group had higher stromal scores, immune scores and TIDE scores than the low CD8A expression group ( Figure 5C), and CD8A expression was positively related to stromal scores, immune scores and TIDE scores (Supplementary Figure S3A-C). Correlation matrix demonstrated that CD8A was positively related to some TIICs (CD8+ T cell, CD4+ T cell, Macrophages M1) and critical immune checkpoint genes (PDCD1, CTLA4, LAG3, TIGIT, HAVCR2, CD274 and TNFRSF9) ( Figure 5D). Scatter plots also showed that CD8A expression was correlated to purity, Figure S3D). In addition, we investigated the association between CD8A and some critical immune checkpoint genes, and the results showed that CD8A expression was positively associated with the expression of critical immune checkpoint genes (Supplementary Figure S3E). Representative immunostainings for high and low CD8A expression in BCa are shown in Figure 5E,F, respectively.

CD8+ T cells, CD4+ T cells, neutrophil and dendritic cells (Supplementary
We then investigated the correlation between CD8A expression and 20 TIICs in a cohort of 30 cancer types. The results revealed that CD8A expression was highly correlated to most of TIICs in pan-cancer (Supplementary Figure S3F). Obviously, CD8A expression was positively correlated to the proportions of CD8+ T cells and Macrophages M1 in pan-cancer (Supplementary Figure S3F).
In BCa, we compared the proportions of 22 TIICs between low and high CD8A expression groups in TCGA ( Figure 6A, Supplementary Figure S3G). There were 13 of 22 TIICs differently distributed between low and high CD8A expression groups. Specifically, the high CD8A expression group had relatively higher percentages of CD8+ T cells and Macrophages M1 than the low CD8A expression group ( Figure 6B, Supplementary Table S2).

Construction and Performance of the Radiomics Signature
In this study, 111 BCa patients with RNA-sequence data and preoperative MRI were collected in our center. There were 77 and 34 BCa patients in the training and validation sets, respectively. Table 1 shows the clinical characteristics of BCa patients. There were no statistical differences between training and validation sets.

Construction and Performance of the Radiomics Signature
In this study, 111 BCa patients with RNA-sequence data and preoperative MRI were collected in our center. There were 77 and 34 BCa patients in the training and validation sets, respectively. Table 1 shows the clinical characteristics of BCa patients. There were no statistical differences between training and validation sets.
The ICCs of radiomics features were between 0.762 and 0.908, suggesting a good inter-and intra-observer reproducibility. We then selected the top ten radiomics features through the mRMR algorithm to construct the radiomics signature. In this way, a LASSO model with nine radiomics features was constructed ( Figure 7A,B). Supplementary Figure  S4A     Patients in the CD8A high expression group had significantly lower radiomics scores than patients in the CD8A low expression group in both training and validation sets (Supplementary Figure S4C,D, both p < 0.001). The waterfall plot revealed that patients with high radiomics scores had a strong tendency for CD8A low expression in the combined training and validation sets, which suggested that the CD8A expression of BCa patients could be correctly predicted based on the cutoff value of the LASSO model ( Figure 7E).

The Potential Biological Mechanisms of the Radiomics-Predicted CD8A Expression
GSEA revealed that several immune-associated pathways, including LYMPHOCYTE_ CHEMOTAXIS, LYMPHOCYTE_MEDIATED_IMMUNITY, NATURAL_KILLER_CELL_AC TIVATION_INVOLVED_IN_IMMUNE_RESPONSE, POSITIVE_REGULATION_OF_CELL_ KILLING, POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY and POSITIVE_REGULATION_OF_T_CELL_PROLIFERATION, were mainly enriched in the radiomics-predicted high CD8A expression group ( Figure 8A). GSVA also showed similar results with the ANTIGEN_PROCESSING_AND_PRESENTATION, NATURAL_KILLER_C ELL_MEDIATED_CYTOTOXICITY, T_CELL_RECEPTOR_SIGNALING_PATHWAY and B_CELL_RECEPTOR_SIGNALING_PATHWAY highly enriched in the radiomics-predicted high CD8A expression group ( Figure 8B). The relative proportions of 22 TIICs based on the RNA-sequence data in our center were evaluated using CIBERSORT analysis (Supplementary Figure S5A). The radiomics score of the LASSO model was negatively related to CD8A expression and degree of infiltration of CD8+ T cells and Macrophages M1 ( Figure 8C-E, Supplementary Figure S5B), and the radiomics-predicted high CD8A expression group had relatively higher percentages of CD8+ T cell and Macrophages M1 than the radiomics-predicted low CD8A expression group ( Figure 8F, Supplementary Table S3). In addition, all of the T cytotoxic pathway-related genes were differentially expressed between radiomics-predicted low and high CD8A expression groups (Supplementary Figure S5C).

Discussion
Previous studies had revealed the important role of the TME in the prognosis and immunotherapeutic response of tumor patients [26,27]. As a crucial component of the TME, TIICs act as a trigger for determining the sensitivity of tumor cells to immune system attack and immunotherapeutic response [28,29]. It has been reported that there were 11 types of TIICs significantly associated with OS in BCa [30]. Among TIICs, cytotoxic T cells (CD8+ T cells) as tumor antigen-specific cytotoxic T cells have a crucial role in the cytolytic killing of tumor cells. The low degree of infiltration of CD8+ T cells at the invasive margin and the tumor center was related to a poor OS in MIBC patients that underwent radical cystectomy [31]. In addition, CD8+ T cells can enhance antitumor activity to inhibit tumor cells through class I MHC [32] and enhance immunotherapy by activating CD4+ T cells through class II MHC [33]. On the cell membrane of CD8+ T cells, there is a cell surface glycoprotein named CD8 antigen that acts as a coreceptor with the T-cell receptor on the T cell to recognize antigens displayed by an antigen-presenting cell in the context of class I MHC molecules. CD8A encodes CD8 antigen, so CD8A expression may be associated with survival outcomes, immune cell infiltration and immunotherapeutic response in cancer patients. Thus, we aimed to explore the clinical predictive value of CD8A in prognosis and immunotherapeutic response and investigate the association between CD8A and TME in BCa.
All the T cytotoxic pathway-related genes were comprehensively analyzed in BCa to select the hub gene. As a result, CD8A was a novel gene that had the highest correlation with several TIICs and had the optimal prognostic value in BCa. Pan-cancer analysis, meta-analysis and multivariate Cox regression analysis were performed to further comprehensively explore the prognostic value of CD8A in cancer patients, especially BCa patients. The results revealed that CD8A was a useful prognostic factor among several types of cancers including BCa, and patients with low CD8A expression were prone to having poor survival outcomes. Previous study also revealed the protective role of CD8A in the prognoses of hepatocellular carcinoma, metastatic melanoma and head and neck squamous cell carcinoma [34][35][36]. However, few studies clarified the association between CD8A and immunotherapeutic response [37], so we further investigated the prognostic value of CD8A in cancer patients treated with immunotherapy in two datasets (IMvigor210 and GSE93157). Our results showed that low CD8A expression was associated with poor survival outcomes among cancer patients treated with immunotherapy, and patients in the high CD8A expression group had a higher percentage of PR/CR than those in the low CD8A expression group. Thus, CD8A can not only be a useful prognostic factor in BCa patients but also a predictive marker of immunotherapeutic response in cancer patients treated with immunotherapy. In addition, we also found that patients with high CD8A expression were more sensitive to several drugs for anti-cancer treatment, which may broaden the application of CD8A in clinical practice.
It is known that immune cell infiltration and activation are associated with the prognosis and immunotherapeutic response in BCa, so we used the bioinformatics analysis to investigate the relationship between CD8A and TIICs. GSEA and GSVA revealed that CD8A played an important role in the positive regulation of immune-associated pathways. In addition, CD8A expression was related to most types of TIICs, especially CD8+ T cells and Macrophages M1. Patients with high CD8A expression are prone to having a higher degree of infiltration of CD8+ T cells and Macrophages M1. As members of TIICs, CD8+ T cells and Macrophages M1 play a crucial role in the enhancement of antitumor activity and immunotherapeutic response, contributing to these two types of TIICs as prognostic and immunotherapy biomarkers in BCa. Jansen et al. reported that BCa patients with a high density of CD8+ T cells had favorite survival outcomes and strong responses to immunotherapy [38]. Zeng et al. reported that Macrophages M1 infiltration is a useful predictive marker for immunotherapeutic response and immunophenotype determination in metastatic urothelial cancer, and patients with high TMB and a high degree of Macrophages M1 infiltration had the best prognoses [39]. Thus, the positive association of CD8A with CD8+ T cells and Macrophages M1 partially explained the CD8A as a prognostic factor and a predictive marker of immunotherapeutic response in BCa.
Previous studies have reported some biomarkers related to cancer immunotherapeutic response such as TMB and the expression of some critical immune checkpoint genes [27,40], so we further explored the association of CD8A with TMB and critical immune checkpoint genes. Pan-cancer analysis revealed that CD8A was related to TMB in several types of tumors. In BCa, patients with higher CD8A expressions were prone to having higher levels of TMB. In addition, CD8A expression was positively related to some critical immune checkpoint genes, including PDCD1, CTLA4, LAG3, TIGIT, HAVCR2, CD274 and TNFRSF9. Therefore, the novel association of CD8A with TIICS, TMB and critical immune checkpoint genes contributes to the crucial role in TME and immunotherapeutic response, which facilitates CD8A as a useful prognostic marker and a predictor of immunotherapeutic response in BCa.
Due to the crucial role of CD8A in predicting the survival outcomes and immunotherapeutic response, we further used the MRI-based radiogenomics to preoperatively predict the expression of CD8A in BCa. The newly developed VI-RADS standardizes the imaging and reporting of BCa on MRI. Many institutions have validated its reliability as an appropriate tool for local staging and determining detrusor muscle invasion in BCa. However, it highly relies on the evaluation of largely morphologic features by radiologists' assessment, such as the size of the lesion, and the enhancement and signal characteristics in digital medical images, which are insufficient and subjective. In addition, it is highly difficult to distinguish the extensive digital characteristics of tumor cells in digital medical images by simple quantitative methods and visual assessment [41]. Radiomics can extract highthroughput quantitative radiomics features that include comprehensive information of tumor biology such as tumor invasive growth patterns, necrosis and neovascularity. This technique is highly sensitive in distinguishing subtle changes in tumor morphology pathophysiology and proteomics, which cannot be deciphered by the human eye. Radiomics had been applied to preoperatively predict the muscle-invasive status, tumor grade, prognosis and treatment response in BCa [42], revealing that the radiomics has potential in predicting the risk groups of BCa patients in training and validation sets, suggesting that some of the MRI-based radiomics features have the potential to reflect the biological behavior on the onset of a tumor. In this study, radiomics features extracted from T2WI and DCE images were used to develop a radiomics signature to preoperatively predict the CD8A expression of BCa patients in our center. A LASSO model based on nine radiomics features achieved good performance in predicting the CD8A expression of BCa patients in both training and validation sets, indicating that the MRI-based radiomics features were associated with the biology of BCa. Interestingly, in this LASSO model, the numbers of radiomics features from T2WI and DCE were four and five, indicating that T2WI-based radiomics features and DCE-based radiomics features have the same importance in the assessment of CD8A expression. In addition, four of nine radiomics features in the LASSO model were wavelet-filtered features, which indicates that the wavelet transform filter has the novel ability to reflect genetic information in BCa. The wavelet transform filter can create eight decompositions per level in each of the three dimensions and generate high-dimensional wavelet-filtered features that cannot be deciphered by visual assessment. Compared with visual assessment by radiologists or low-level radiomics features, wavelet-filtered features include comprehensive information of tumor heterogeneity and tumor biology in several types of tumors, including prostate carcinoma [43], renal cell carcinoma [44], intrahepatic cholangiocarcinoma [45] and BCa [46]. We further investigated the biological behavior of the LASSO model based on the RNA-sequence data in our center. Interestingly, the radiomics score of the LASSO model was negatively related to CD8A expression and the degree of infiltration of CD8+ T cells and Macrophages M1. GSEA and GSVA further showed that several positive regulations of immune-associated pathways were highly enriched in the radiomics-predicted CD8A high expression group. These results revealed distinct TME between radiomics-predicted low and high CD8A expression groups and validated good performance of the LASSO model from an immunological perspective.
This study had limitations. One is that the mechanisms related to the regulation of TME by CD8A in BCa were not investigated, which will require additional mechanistic analyses in the future. Another limitation is that the number of BCa samples in our center is small. A larger number of BCa patients with preoperative MRI and RNA sequence should be collected to further assess the performance of the LASSO model in predicting CD8A expression.

Conclusions
In conclusion, our study illustrated that CD8A can be a useful indicator for predicting the survival outcomes and immunotherapeutic response in BCa. The expression of CD8A was associated with TME characteristics and TMB, which facilitates CD8A as a predictor of immunotherapeutic response. A radiomics signature based on nine MRI-based radiomics features has the potential to preoperatively predict the expression of CD8A in BCa patients, thereby contributing to the prediction of prognosis and immunotherapeutic susceptibility. (B) Percentages of the N stage between high or low CD8A expression groups in TCGA. (C) Percentages of the survival outcomes between high or low CD8A expression groups in TCGA. (D) Sensitivity analysis of the meta-analysis. CI, confidence interval; OS, overall survival. *: p < 0.05., **: p < 0.01, ***: p < 0.001. Figure S3: The association of CD8A with TME characteristics. (A-C) The association of CD8A with ImmuneScore (A), StromalScore (B) and ESTIMATEScore (C). (D) Pan-cancer analysis of the association of CD8A expression with TIICs. (E) The relative infiltration levels of TIICs based on the expression of CD8A. TME, tumor microenvironment; TIICs, tumor-infiltrating immune cells. *: p < 0.05., **: p < 0.01, ***: p < 0.001. Figure S4: The Correlation matrix and coefficients of the radiomics features in LASSO model. (A) The associations between nine radiomics features. (B) The coefficients of nine radiomics features. LASSO, least absolute shrinkage and selection operator. Figure S5: The association of radiomics-predicted CD8A expression with TIICs and T cytotoxic pathway-related genes in our center. (A) The proportions of 22 TIICs in each sample quantified by CIBERSORT algorithm in our center. (B) The relative infiltration levels of TIICs based on the radiomics-predicted CD8A expression. (C) The difference of the expression of T cytotoxic pathwayrelated genes between radiomics-predicted high and low CD8A expression. TIICs, tumor-infiltrating immune cells. *: p < 0.05., **: p < 0.01., ***: p < 0.001. Table S1. Data set information included in this study. Table S2. The relative infiltration levels of TIICs based on the expression of CD8A. Table S3. The relative infiltration levels of TIICs based on the radiomics-predicted CD8A expression status.

Institutional Review Board Statement:
The studies involving human participants were reviewed and approved by the Ethics Committee of Shanghai Tenth People's Hospital (approval number 2021KN108).
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study. Requests to access the datasets should be directed to Shenghua Liu, drfelixliu@163.com.