Identification and validation of a combined hypoxia and immune index for triple‐negative breast cancer

The interaction between hypoxia and immune status has been confirmed in various cancer settings, and corresponding treatments have been investigated. However, reliable biomarkers are needed for individual treatment, so we sought to develop a novel scoring system based on hypoxia and immune status. Prognostic hypoxia–immune status‐related signatures of patients with triple‐negative breast cancer (TNBC) were identified in The Cancer Genome Atlas (TCGA) (N = 158), Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) (N = 297), and GSE58812 (N = 107). LASSO Cox regression was used for model construction. Hypoxia and immune status expression profiles were analyzed, and infiltrating immune cells were compared. Quantitative real‐time PCR (qRT‐PCR) was used for validation in the Sun Yat‐sen University Cancer Center (SYSUCC) cohort, and immunofluorescence was applied for the detection of hypoxia and immune markers in cancer tissues. Ten cross‐cohort prognostic hypoxia–immune signatures were included to construct the comprehensive index of hypoxia and immune (CIHI) in the METABRIC cohort. Two subgroups of patients with distinct hypoxia–immune status conditions were identified using CIHI: hypoxiahigh/immunelow and hypoxialow/immunehigh, with a significantly better overall survival (OS) rate in the latter (P < 0.01). The prognostic value of CIHI was further validated in the TCGA, GSE58812, and SYSUCC cohorts (P < 0.01). Hypoxia–immune signatures were significantly differentially expressed between the two groups, and more active immune responses were observed in the hypoxialow/immunehigh group. Cytotoxic lymphocytes were inversely correlated with CIHI in silico. Differentially expressed CA‐IX and stromal PD‐L1 were detected between subgroups of the SYSUCC cohort. A hypoxia–immune‐based cross‐cohort classifier for predicting prognosis was developed and validated, which may guide hypoxia modifier treatment and immunotherapy for TNBC.

Hypoxia-immune signatures were significantly differentially expressed between the two groups, and more active immune responses were observed in the hypoxia low /immune high group. Cytotoxic lymphocytes were inversely correlated with CIHI in silico. Differentially expressed CA-IX and stromal PD-L1 were detected between subgroups of the SYSUCC cohort. A hypoxia-immune-based cross-cohort classifier for predicting prognosis was developed and validated, which may guide hypoxia modifier treatment and immunotherapy for TNBC.

Introduction
Breast cancer is the most common cancer in women worldwide, with an estimated annual death of 41760 cases in women [1,2]. Triple-negative breast cancer (TNBC), a special subtype that represents 10-20% of patients with breast cancer, exhibits the most malignant biological behaviors and worst clinical outcomes [3]. Surgery and chemotherapy are considered the firstline regimens for TNBC. Neither targeted therapy for Her2 (human epidermal growth factor receptor 2) nor endocrine therapy is applied for patients with TNBC in clinical practice [4]. The treatment for TNBC is so limited that it is urgent to develop effective therapies.
In previous studies, TNBC has been clustered into six different molecular subtypes according to genomewide profiling [5]. Lehmann et al. classified patients with TNBC into immunomodulatory (IM), basal-like 1 (BL1), basal-like 2 (BL2), mesenchymal (M), luminal androgen receptor (LAR), and mesenchymal stem-like (MSL) groups. Moreover, Bareche et al. [6] found a higher level of immune features, including checkpoint molecules in the IM subtype. These results indicated that immunotherapies might be effective in particular subtypes of TNBC. Recently, several clinical trials have addressed the effect of immunotherapies on both early-stage and metastatic TNBC, and there have been advances in our understanding of the biological and immunological characteristics of TNBC [7,8]. Furthermore, the addition of immune checkpoint inhibitors could increase the pathologic complete response (pCR) rate of neoadjuvant chemotherapy in early-stage TNBC [9]. Hence, immunotherapies have become novel promising candidates for TNBC treatment, and reliable indicators for treatment effect evaluation are indispensable. Although many different prognostic biomarkers and models have been developed to describe and quantify the immunological characteristics of TNBC, few have considered the effect of the extracellular microenvironment, such as hypoxia and pH, on cancer cells [10,11].
The hypoxia-related mechanism has long been considered one of the hallmarks in the cancer signaling pathway [12][13][14]. Hypoxia is a typical event in solid tumors and is associated with cancer metabolic reprogramming, stem cell signatures, angiogenesis, extracellular matrix (ECM) organization, and cancer cell metastasis [15][16][17]. Recent studies revealed novel roles of hypoxia in tumor progression and metastasis, including the induction of rapid methylation changes to histone and chromatin reprogramming, the development of an reactive oxygen species (ROS)-resistant phenotype, and the induction of epithelial-mesenchymal transition (EMT) through macrophages [18][19][20].
Previous studies have investigated the association between hypoxia and the tumor microenvironment (TME). The hypoxia response in T lymphocytes enhances the expression of CD137 for immunotherapy [21]. The aerobic glycolysis of breast cancer cells could be regulated by lncRNAs from macrophages in the microenvironment [22]. Moreover, under hypoxic conditions, the dampening of NRF1 degradation could impair tumor-associated macrophages (TAM) polarization [23]. However, hypoxia also drives CD8 + Tcell migration and effector function, which indicates that hypoxia plays different roles in immune cells and tumor cells [24]. Although there are explicit links between hypoxia and the immune microenvironment, few studies have focused on the comprehensive status of hypoxia and the immune response in the TME.
Here, we aimed to develop and validate a comprehensive index of hypoxia and immune (CIHI) to reflect and quantify the microenvironment of TNBC based on these results. In this study, hypoxia-related genes (HRGs) and immune-related genes (IRGs) were screened for prognostic relevance and applied to model construction in silico. We confirmed that the CIHI was significantly correlated with hypoxia and immune status using bioinformatics analysis. We further evaluated the CIHI in tissue samples using quantitative reverse transcription polymerase chain reaction (qRT-PCR) and immunofluorescence (IF) to further facilitate clinical application.

Study population and data acquisition
The workflow is displayed in Fig. 1 and Fig. S1. For TNBC datasets from the Molecular Taxonomy of Breast Cancer International Consortium (METAB-RIC), The Cancer Genome Atlas (TCGA), and the Gene Expression Omnibus (GEO), patients who met the following selection criteria were included: (a) histologically diagnosed with malignant breast cancer; (b) available RNA expression data; and (c) available OS (overall survival) data. Patients without active followup were excluded. In this study, 297 patients from METABRIC, 158 patients from TCGA, and 107 patients from GSE58812 were included. For CIBER-SORT analysis, 293 patients from METABRIC, 137 Fig. 1. The work flow of this study. A panel of integrated prognostically relevant hypoxia-related and immune-related genes was identified in the METABRIC, TCGA, and GSE58812 datasets. A novel quantified index, the CIHI, was estimated using the LASSO Cox regression model, and its crucial roles in hypoxia and immune status were further validated in multiple cohorts. patients from TCGA, and 103 patients from GSE58812 were included after screening. For TNBC RNA sequencing data from FUSCCTNBC cohort, 360 patients were included without overall survival information.
The METABRIC data were obtained from cBioportal (http://www.cbioportal.org/) [25]. The FUSCCTNBC data were downloaded from the Sequence Read Archive (SRA) (SRP157974, https:// www.ncbi.nlm.nih.gov/sra) [26]. In addition, the expression profiles of the TCGA cohort were downloaded from the TCGA data portal (https://portal. gdc.cancer.gov/repository). The ensemble IDs were mapped to gene symbols according to the annotation of Homo_sapiens.GRCh38.91.chr.gtf from the ENSEMBLE website. The 'limma' package in R was used for gene expression normalization using the scale method [27]. The average RNA expression was calculated for duplicates, and genes with low abundance were discarded.
The GSE58812 expression profile was based on the GPL570 platform, and the expression matrix was obtained from the GEO database (https://www.ncbi. nlm.nih.gov/geo/) [28]. The probes were mapped according to the GPL570 annotation file (https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570), and the average RNA expression was calculated for duplicates.
We complied with the access policies of the TCGA, METABRIC, SRA, and GEO databases in this study.
For the Sun Yat-sen University Cancer Center (SYSUCC) cohort, the inclusion criteria were as follows: (a) histologically diagnosed as malignant breast cancer; (b) samples were available after surgery; (c) the molecular subtypes were confirmed by immunohistochemistry (IHC), and Her2 status was further validated using fluorescence in situ hybridization (FISH) if indecipherable in IHC; (d) AJCC 7 th TNM stage; (e) free from any other malignant tumors; and (f) no immunotherapies administered. Patients without active follow-up records were excluded. Thirty patients with triple-negative breast cancer between 2008 and 2011 were included in the SYSUCC study cohort. Overall survival was calculated from the date of diagnosis to the date of the last follow-up or death. Disease-free survival (DFS) was defined as the time from the date of diagnosis to the date of the first recurrence, including local recurrence and distant metastasis. Written informed consent was required for all patients, and this study was approved by the institutional review board of Sun Yat-sen University Cancer Center and followed the guidance of the Declaration of Helsinki.

Generation of hypoxia and immune profiling
The HRGs were extracted from the hallmark gene sets in the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/), and IRGs were obtained from ImmPort (https://www.immport. org). The available HRGs and IRGs in METABRIC, TCGA, and GSE58812 were included in this study. The functional protein-protein interaction network analysis was conducted by STRING (v11.0) (https:// string-db.org). The position in chromosomes, expression, and interactions of these genes were exhibited in Circos plots by the 'circlize' package in R [29]. In total, 200 HRGs and 1811 IRGs were identified in all these cohorts and analyzed by univariate Cox regression for prognostic relevance. For analysis of previously reported hypoxia scores, the scores of the TCGA cohort were obtained from previous studies [30][31][32][33][34].

Development of the CIHI
The prognostically relevant HRGs and IRGs were identified using univariate Cox regression. The least absolute shrinkage and selection operator (LASSO) Cox regression model was further applied to determine the crucial signatures and the corresponding coefficients for model construction [35]. LASSO Cox regression was conducted using the 'glmnet' package in R software, and the ideal coefficients were estimated according to the partial likelihood deviance with tenfold cross validation [36]. The optimal log λ was determined at −3.37. To quantify the comprehensive effects of hypoxia and immune status, a novel score was calculated from the signatures selected by the LASSO model. The following formula was used: For clinical utility, the comparative value (2 −Ct ) was calculated from qRT-PCR results and used for score calculation in the SYSUCC cohort. The score was further standardized and simplified to generate a comprehensive index of hypoxia and immune. The score was subsequently mapped by subtracting the minimum and dividing by the maximum. Mapping was conducted to facilitate the interpretation of results from different platforms. The CIHI was calculated as follows:

Tumor microenvironment analysis
For Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis, the expression matrices were uploaded to the online analytical platform (https://cibersort.stanford.edu), and the proportions of infiltrating immune cells were estimated according to LM22 signatures with 1000 permutations [37]. The qualified samples were then selected using a criterion of P < 0.05. xCell analysis was performed according to the guidance of the website https://xcell.ucsf.edu [38]. Immune and stromal scores were further estimated to quantify the immune and stromal components by the ESTIMATE algorithm using the 'estimate' package in R [39]. MCP-counter scores regarding immune-related activity and fibroblasts were evaluated using the 'MCPcounter' package in R [40]. ssGSEA was performed to calculate the enrichment score of specific immune signatures in samples using the 'GSVA' package in R according to a previous study [41]. For comparison of the hypoxia and immune status between the low-risk and high-risk groups, the key expression profiles were compared.

Functional study and heterogeneity
The GSEA was conducted to explore the pathway enrichment between low-risk and high-risk groups using GSEA 3.0. The annotated gene set was downloaded for reference (http://software.broadinstitue.org). The tumor mutation burden (TMB) and the mutant-allele tumor heterogeneity score (MATH) for the TCGA cohort were estimated by the 'maftools' package in R software [42]. PCA was used to examine the clustering efficacy of the selected signatures in the LASSO model.

RNA isolation and quantitative real-time PCR
The total RNA of breast cancer tissues was extracted using TRIzol reagent (Invitrogen, Carlsbad, California, USA). Reverse transcription was performed according to the manufacturer's instructions (Takara, Kusatsu, Japan). The expression levels of the target genes were further examined in triplicate using the SYBR Green method (Takara). The primers involved in this study are provided in Table S1. The expression levels were normalized to that of β-actin with the comparative Ct method.

Immunofluorescence
Immunofluorescence (IF) was performed as previously described [43]. The formalin-fixed paraffin-embedded (FFPE) sections were deparaffinized, and we performed the antigen-retrieval procedure at 98°C in citrate buffer (pH 6.0) for 10 min. The endogenous peroxidase blocking procedure was then carried out using 3% hydrogen peroxide for 10 min at room temperature. The sections were further incubated with anti-carbonic anhydrase 9 (CA9) (Proteintech, Chicago, Illinois, USA) or programmed death-ligand 1 (PD-L1) (Proteintech) at 4°C overnight. IgG (CST, Danvers, Massachusetts, USA) was used as a negative control. Then, the samples were washed with PBS three times and incubated with anti-mouse or anti-rabbit secondary antibodies, namely, Alexa Fluor-594 or Alexa Fluor-488 (1 : 1000, Invitrogen). A Zeiss LSM 880 confocal microscope was used to observe the results, and images were acquired. The reagents are shown in Table S2.

Statistical analysis
Univariate Cox regression was used to identify the prognostically relevant HRGs and IRGs in the METABRIC, TCGA, and GSE58812 cohorts with a cutoff value of P < 0.1. Crucial signatures involved in hypoxia and immune status were identified by the LASSO Cox regression model. The multivariate Cox regression model was constructed using the 'survival' package to include the CIHI and clinical predictors. The optimal cutoff value of survival analysis was determined using the 'survminer' package in R, and the OS and DFS of different subgroups were compared using the Kaplan-Meier method with the logrank test. Time-dependent receiver operator characteristic (ROC) analyses were performed using the 'time-ROC' package in R [44]. Decision curve analysis (DCA) was performed in R software using a previously reported method [45]. The significance of the difference between immune cell fractions was assessed by the Wilcoxon test. Spearman's correlation test was used for CIHI-related analysis. Statistical analyses were performed using R software (Version 3.6.0). A P value of < 0.05 was considered statistically significant, and all P values were two-tailed.

Overview of hypoxia and immune signatures
A population of 297 patients from METABRIC, 158 patients from TCGA, and 107 patients from GSE58812 were identified and included in this study. The transcriptome data were used to construct a comprehensive indicator from hypoxia and immune profiling ( Fig. 2A,B). Two hundred hypoxia-related genes and 1811 immune-related genes were used for model construction. A total of 715 genes were not identified in all datasets and were excluded from the survival analysis (Fig. 2C). A univariate Cox regression model was applied to determine the prognostic relevance of these signatures, and 32 prognostically relevant genes in the METABRIC, TCGA, and GSE58812 cohorts were identified (Fig. 2D,E).

Construction of a comprehensive index of hypoxia and immune status in TNBC
Given that a hypoxic microenvironment might affect the activation status of infiltrating immune cells and the immune response of tumor cells, an integrated analysis of both hypoxia and immune response might have potential prognostic value and quantify the TME. Hence, 32 prognostically relevant signature genes were applied to the LASSO Cox regression model to construct a predictive model for the overall survival of patients with TNBC in the METABRIC dataset (N = 297). Ten genes were selected according to the partial likelihood deviance method, and the corresponding coefficients were generated at the optimal log λ of −3.37. The results are shown in Fig. 3A,B. Among the 10 signature genes, PFKL, SLC25A1, and SERPINE1 were hypoxia-related genes, while the others were immune-related genes (TAPBPL, CXCL11, HLA-A, TCF7L2, TANK, IL12B, and IL18RAP). Kaplan-Meier analysis further confirmed the prognostic value of the individual genes (Fig. S2). The selected signature genes were applied to the formula above, and the CIHI was calculated. Moreover, the CIHI was also estimated in the TCGA and GSE58812 cohorts using this method.
To further facilitate the application of the CIHI in TNBC, the patients were divided into low-risk and high-risk groups according to the median value of the CIHI. PCA showed that patients in different groups could be distinctively clustered according to the selected signatures in all datasets (Fig. 3C-E). In addition, Spearman's correlation test indicated that the standardized CIHI was significantly correlated with the selected genes. The correlations among these signatures are shown in Fig. 3F.
Previous studies suggested that a high level of hypoxia could lead to immune suppression in the TME [46,47]. These results demonstrated that hypoxia was associated with the immune response in the microenvironment. Thus, this novel scoring method, the CIHI, might be indicative of the hypoxia-immune status of patients and has importance in clinical guidance.

Validation of hypoxia profiling in CIHI
The analysis above was used to construct a comprehensive indicator and identify two distinctive subgroups of patients with TNBC. Correlation analysis in the METABRIC cohort showed a correlation between the CIHI and hypoxia-related genes, indicating that the CIHI might reflect hypoxia in the TME (Fig. 3G). Hence, we further sought to validate the correlation of the CIHI and hypoxia. A previous study addressed the key hypoxia-related expression profiles in cancer [48]. We first compared the key hypoxia-related signatures in the high-risk and low-risk groups. Subsequent analyses mapped the expression level of genes in the two phenotypes of the METABRIC cohort to reflect the distinct hypoxia status. Sixteen genes showed a statistically significant difference between the subgroups (low-risk group < high-risk group). The expression of ALDOA, ANGPTL4, CA9, DCBLD1, ENO1, FOSL1, HK1, KCTD11, LDHA, P4HA1, PDK1, PFKL, PGAM1, SDC1, and VEGFA was significantly higher in the high-risk group (Fig. 4  A). The results showed that hypoxia-induced metabolic rearrangement and angiogenesis were more common in the high-risk group.
To further validate the results, qRT-PCR was performed in the SYSUCC cohort (N = 30) to facilitate the clinical application of the CIHI. Based on the relative Ct values of the ten selected signature genes described above, the CIHI was estimated using our previous formula. Previous studies addressed the indicative function of Carbonic anhydrase IX (CA-IX) in breast cancer and made it a reliable reflection of hypoxia [49][50][51]. Therefore, the protein level of CA-IX was examined in these patients using IF (Fig. 4B,C). Patients with high CIHI values were prone to exhibit a higher level of CA-IX protein in tumor cells. We also detected the PD-L1 protein because PD-L1 is a crucial target for immunotherapy and reflects a potential immune response in tumors. Interestingly, patients with low CIHI values exhibited a lower level of hypoxia but a higher level of PD-L1, which was consistent with the following results regarding immune status. Moreover, a stromal PD-L1 expression was more prevalent. Hence, the low-risk and high-risk groups were more likely to exhibit hypoxia low /immune high and hypoxia high /immune low phenotypes, respectively.

Immune profiling
Further results showed a continuity feature of the CIHI and its correlation with the TME and clinical prognosis  . (A, B) The LASSO Cox regression model was constructed from 32 signature genes, and the tuning parameter (λ) was calculated based on the partial likelihood deviance with ten-fold cross validation. An optimal log λ value was defined at −3.37, as shown by the vertical black lines in the plots. The ten signature genes were identified according to the best fit profile. (C-E) PCA based on the expression profile of the 10 selected signature genes according to different risk groups. (F) Correlation between the CIHI and the selected signature genes in the METABRIC cohort. (G) Correlation network between the CIHI and its correlated hypoxia-related genes in the METABRIC cohort. Data were analyzed using Spearman's rank correlation analysis.
in the METABRIC, TCGA, and GSE58812 cohorts (Fig. 5A,D and G). The value of the CIHI was mapped into a continuous scale, which indicated potential clinical relevance in predicting the prognosis of patients. For overall survival, patients with higher CIHI values seemed to have a higher rate of death. Moreover, the ESTIMATE method was employed to explore the overall TME status. A low immune score and stromal score were more commonly observed in patients with a higher CIHI. Further analysis indicated a significant inverse correlation between the CIHI and immune score in the METABRIC (r = −0.69, P < 0.01), TCGA (r = −0.68, P < 0.01) and GSE58812 (r = −0.80, P < 0.01) cohorts (Fig. 5B,E and H). Interestingly, the CIHI was also negatively correlated with the stromal score, indicating that the CIHI might also be associated with stromal nonimmune components, such as fibroblasts (Fig. 5C,F and I).
To further facilitate the indicative function of the CIHI in the immune response, a panel of immune-response signatures was also mapped into the hypoxia high /immune low and hypoxia low /immune high phenotypes (Fig. 6A). Immune response-related signatures were mostly differentially expressed: hypoxia high / immune low < hypoxia low /immune high .
Our results demonstrated a different pattern of immune response between the phenotypes. Overall, patients with the hypoxia high /immune low phenotype exhibited a suppressive immune microenvironment compared with patients with the hypoxia low /immune high phenotype. This finding was consistent with the fact that these signatures were responsible for mediating pro-or antitumoral activity. The phenotypic and functional markers of T cells, namely, CD3E, CD4, CD8B, GZMB, PRF1 and TBX21, were expressed at higher levels in the hypoxia low /immune high phenotype, although FOXP3 did not show a significant difference. For the myeloid lineage phenotypic and functional markers, CD14 and CD163 were highly expressed in the hypoxia low /immune high group, indicating a higher percentage of monocytes and M2 macrophages. However, CD33, a previously reported marker for M-MDSCs, was highly expressed in the hypoxia high /immune low group. Other markers, including IFNγ signatures (CXCL10, CXCL9, IDO1, IFNG, and STAT1) and immune modulators (ENTPD1), also exhibited a higher expression level in the low-risk group (hypoxia low /immune high ). Interestingly, a higher level of inhibitory immune receptors or ligands (PD-L1, CTLA4, LAG3, and PD-1) and activating immune receptors (CD27, CD40, CD80, ICOS, TNFRSF4, and TNFRSF9) were both observed in the hypoxia low /immune high phenotype, indicating a complex immune response in the low-risk group. Then, we analyzed the tumor mutation burden (TMB) and tumor heterogeneity based on genomewide methods. Although a distinct immune response was observed in different phenotypes, the TMB did not show a consistent result (Fig. 6B). However, the MATH score, a quantified scale for genome heterogeneity, demonstrated a higher level of tumor heterogeneity in the high-risk group (P = 0.091) (Fig. 6C). This result was consistent with that of a previous study showing that the hypoxic microenvironment might exert an effect on the genome.

The relevance of the CIHI and hypoxia-immune markers
As described above, qRT-PCR of the ten-gene signature was performed, and the CIHI was calculated for patients in the SYSUCC cohort. The hypoxia and immune markers were then examined in tumor samples by IF. CA-IX was previously reported to be a reliable marker for hypoxia in breast cancer and was mainly expressed in cancer cells instead of stromal cells. A  high level of CA-IX often indicates a hypoxic microenvironment in tissues. PD-L1 is widely recognized as a potential immune marker and therapeutic target for cancers. Patients with a higher CIHI were more likely to exhibit positive CA-IX and negative stromal PD-L1 expression (Fig. 4B,C). However, due to the limited number of patients, we did not find any convincing evidence for the correlation between CA-IX and tumoral PD-L1 expression.
To further facilitate our hypothesis in larger cohorts, we first estimated the immune cell fractions in the METABRIC cohort using the CIBERSORT algorithm (Fig. 6D). Patients in the low-risk group exhibited a higher percentage of antitumoral immune cells, including naive B cells (P < 0.01), CD8 + T cells (P < 0.01), activated memory CD4 + T cells (P < 0.01), activated NK cells (P < 0.01), and activated dendritic cells (P < 0.01), while patients in the high-risk group showed a higher level of regulatory T cells (P = 0.08) and M2 macrophages (P < 0.01). Moreover, the CIBERSORT results for the TCGA and GSE58812 cohorts are shown in Fig. S3. The high-risk group presented with a higher percentage of M1 macrophages and a lower percentage of M2 macrophages in both cohorts. Infiltrating cell scores in TME were also estimated by xCell, and similar results were shown in Table S3.
We then applied GSEA to examine the relevant signaling pathways involved in patients with a low CIHI. Our results demonstrated that antigen processing and presentation, B-cell activation, and T-cell activation were significantly enriched in the hypoxia low /immune high group (Fig. 7A-C). Furthermore, more algorithms were employed to validate the association between the CIHI and infiltrating immune cells. Using a previously reported ssGSEA method, the crucial signals regarding immune cells were estimated and scored in the METABRIC, TCGA, and GSE58812 cohorts. Significant reverse correlations between the CIHI and activated B cells (r = −0.70, P < 0.01), activated CD4 + T cells (r = −0.75, P < 0.01), and activated CD8 + T cells (r = −0.74, P < 0.01) were observed in the METABRIC cohort (Fig. 7D). The results were further validated in the TCGA and GSE58812 cohorts (Fig. 7E,F). The MCP-counter score demonstrated that the CIHI was inversely correlated with cytotoxic lymphocytes and NK cells in these cohorts (Fig. 7G--I).

Prognostic value of the CIHI
The development of convenient tools for the early diagnosis of diseases and treatment guidance remains a crucial clinical issue. Previous studies have demonstrated that hypoxia and immune status are indicators of malignant tumors. In our study, as shown above, the CIHI exhibited a continuity feature, and remarkably heterogeneous characteristics of hypoxia and immune status were observed between the hypoxia high / immune low and hypoxia low /immune high phenotypes. More death events were observed in patients with higher CIHI values. To further elucidate the predictive and prognostic value of the CIHI in TNBC, ROC analysis and the Kaplan-Meier method were used to examine and validate prognosis in the cohorts.
Time-dependent ROC analysis was performed, and the area under the curve (AUC) was calculated at different time points according to data availability ( Fig. 8A-C). Generally, the ROC analysis of the 5year follow-up in the METABRIC (AUC = 0.690), TCGA (AUC = 0.757), and GSE58812 (AUC = 0.762) cohorts indicated a favorable predictive value of the CIHI in long-term follow-up. Moreover, the CIHI seemed to present with an evidently better predictive ability for OS in the TCGA cohort during short-term follow-up (AUC = 0.887). These results indicated that the CIHI could serve as a clinical biomarker. The decision curve analysis also confirmed the predictive value of the CIHI in cohorts (Fig. S4).
Overall survival was compared in patients with different CIHI values. The patients in the hypoxia low /immune high group exhibited a better OS than those in the high-risk group in the METABRIC cohort, the TCGA cohort, and the GSE58812 cohort (Fig. 8D-F). Moreover, patients with higher CIHI values also showed a poorer outcome for DFS in the TCGA cohort (Fig. 8  G). The online database KM plotter facilitated the construction of predictive models and evaluation of prognosis in breast cancer patients. The risk score was calculated with the KM plotter, and progression-free survival (PFS) and distant metastasis-free survival (DMFS) were calculated as described above (Fig. 8H, I). Patients with a low risk score had a better PFS [hazard ratio (HR), 2.68; 95% confidence interval (95% CI), 1.7-4.23; P < 0.01] and DMFS (HR, 3.91; 95% CI, 0.79-19.43; P = 0.072). The CIHI was calculated for patients in the SYSUCC cohort, and the hypoxia low /immune high group showed significantly better OS and DFS compared with the hypoxia high /immune low group (Fig. S5).
Furthermore, we sought to validate the indicative role of the CIHI for hypoxia and immune characteristics in a larger cohort. The data from FUSCCTNBC cohort were analyzed, and patients were divided into groups according to the median value of CIHI. The hypoxia markers were highly expressed in the high-CIHI group (Fig. S8a). The ESTIMATE immune and stromal scores were also inversely correlated with the CIHI (Fig. S8b,c). MCP-counter and ssGSEA scores indicated that patients with lower CIHI value were characteristics of more tumor-infiltrating cytotoxic immune cells (Fig. S8d,e). The CIHI value was compared among the previously reported subtypes, and patients with the IM (immunomodulatory) subtype showed a significantly lower CIHI value, which further supported our findings (Fig. S8f). Previously reported hypoxia scores were significantly higher in CIHI-high group (Fig. S9). Therefore, a lower CIHI might predict relatively active immune response, and a relatively higher CIHI was more likely to predict a hypoxia high / immune low status (Fig. 9).

Discussion
Previous studies have addressed the important role of hypoxia and immune status features in breast cancer [52,53]. Recent studies have advanced our understanding of the efficacy of immunotherapies in TNBC, and promising results have emerged in several clinical trials [7][8][9]. Researchers have focused on the effect of targeted therapies on immune checkpoints in triple-negative breast cancer, and the effect of the combination of PD-L1/PD-1 blockade with routine therapies has been investigated.
As shown in previous studies, the hypoxic microenvironment drives the suppression of immune status in breast cancer. The tumor microenvironment is reprogrammed under the pressure of hypoxia, resulting from the modulation of CD8 + cytotoxic T cells and macrophages [21,23,24,46,54]. Moreover, the metabolic reprogramming of tumor cells provides an adaptation under hypoxia by increasing lactate production and glucose uptake, contributing to the development of an immunosuppressive microenvironment [55,56]. Interestingly, crosstalk from TAMs to breast cancer cells by exosomes regulates HIF-1α-mediated aerobic glycolysis metabolism and contributes to tumor progression in breast cancer [22].
Although clinical indicators regarding hypoxia or immune status have been developed, few studies have focused on their comprehensive effects and their potential roles in clinically relevant classification and therapy selection. Generally, patients with similar clinical characteristics still present with great heterogeneity in clinical outcomes, and integrated predictors from single biomarkers could effectively improve the prognostic value.
Given that hypoxia modifiers and immune checkpoint inhibitors have been shown to exhibit a potential effect on breast cancer, we explored the potential role of a combined hypoxia and immune status classifier for TNBC in this study [57,58]. The use of a combined immune-hypoxia signature in a cross-cohort manner helped to develop a continuous index for comprehensive TME assessment. Subgroup classification divided the population into the hypoxia high /immune low and hypoxia low /immune high groups, correlating with distinct clinical prognosis, transcriptional hypoxia-immune patterns, and activated pathways that could serve as therapeutic targets. Therefore, we addressed the importance of separating them.
The key markers of hypoxia (CA9, ALDOA, ENO1, LDHA, etc.) are highly expressed in the hypoxia high /immune low group compared with the hypoxia low /immune high group, indicating heterogeneous hypoxia statuses within the population. In addition, immune cell infiltration was also investigated to further validate our model. Infiltrating immune cells, including cytotoxic T cells, macrophages and NK cells, are widely recognized as crucial immune indicators in Fig. 9. Proposed relatively indicative roles of the CIHI in the relative hypoxia and immune status of the tumor microenvironment. A lower CIHI predicted relatively active immune response, and a relatively higher CIHI was more likely to predict a hypoxia high /immune low status.
cancers [59][60][61]. In our study, a higher percentage of cytotoxic T cells, M2 macrophages, and NK cells were observed in the hypoxia low /immune high group. PD-L1 expression, especially within stromal cells in the TME, is correlated with clinical prognosis in breast cancer and may reflect an adaptive immune response. Previous studies have revealed the distinct role of PD-L1 expression in stromal immune cells and tumor cells [62,63]. The PD-L1 protein is always positive in tumor cells when that in immune cells is positive [64]. Moreover, patients with positive PD-L1 expression are more likely to exhibit better prognosis. Hence, it is more reasonable and convenient to use the stromal PD-L1 as an indicator for detection [8]. In our study, immunofluorescence showed a higher level of CA-IX, a marker for hypoxia, and correspondingly a lower level of PD-L1 in the hypoxia high /immune low group.
Interestingly, both the activating and inhibitory immune markers were highly expressed in the hypoxia low /immune high group, which seems contradictory in common sense. This result may indicate a consuming anticancer immune response in the hypoxia low /immune high group with an increased anti-immune reaction in tumor cells. The active anticancer immune response may be restored through immune checkpoint blockade. As described above, myeloid-derived suppressor cells (MDSCs) and TAMs may be regulated by the hypoxia status and contribute to immunosuppression.
Overall, hypoxia exerts its function on tumor cells and the TME, tipping the balance of the immune response of chemokines, immune effector production, and immune cell infiltration. For patients with a relatively hypoxia-high/immune-low feature, hypoxic modifier combined with subsequent immunotherapies is potentially applicable.
There are certain limitations to our study. First, this is a cross-cohort and retrospective study, and further prospective analysis in multicenter cohorts is necessary due to the observed heterogeneity between different populations. Second, although a prognostic continuous index has been developed, further inclusion of clinical factors (TNM stage, histological grade, age at diagnosis, and anatomic site) may minimize the potential of bias. Third, this is a retrospective study, and the size of the SYSUCC cohort was relatively small; only stromal PD-L1 was analyzed. Finally, qRT-PCR analysis of all selected signatures might be costly in clinical practice, and more convenient methods could be further applied, such as IHC using CA-IX and relevant immune markers.
The development of a novel comprehensive hypoxia-immune status classifier in our study implies that individual treatment should be applied in different subgroups. For example, immunotherapies might be more favorable in patients with hypoxia low /immune high features. A higher CIHI was indicative of an active immune response and potential benefit from immune checkpoint blockade, such as anti-PD1/PD-L1 treatment. For patients in the hypoxia high /immune low group with a poorer prognosis, a lower CIHI might imply a hypoxic microenvironment and lower efficacy for immunotherapy. However, hypoxia modifiers might help to sensitize patients to further radiotherapy and chemotherapy after surgery. The results demonstrated that subsequent immunotherapy might work after hypoxia modification. Further subsequent immunotherapy might prevent the immunosuppressive response after treatment. For those with an intermediate CIHI exhibiting a mixed status of hypoxia and immune response, both anti-hypoxia and immunomodulatory drugs may harbor certain effects on patients.

Conclusions
In conclusion, a novel prognostic index and classifier based on hypoxia and immune expression profiles were developed and cross-cohort validated. This classifier could be used for prognostic prediction and selecting patients for hypoxia modification and immunotherapies, such as PD-L1/PD-1 and CTLA-4 blockade. We recommend further validating the efficacy of the CIHI in prospective studies and developing a simplified version using more convenient methods, such as IHC.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Fig S1. Detailed flowchart of this study. Fig S2. Kaplan-Meier analysis of the genes included for CIHI construction (TAPBPL, CXCL11, HLA-A, TCF7L2, TANK, IL12B, IL18RAP, PFKL, SLC25A1, and SERPINE1). Fig S3. (a, b) Comparison of 22 infiltrating immune cells in different subgroups of the TCGA and GSE58812 cohorts. Fig S4. Decision curve analysis of the CIHI in the METABRIC, TCGA, and GSE58812 cohorts. Fig S5. (a, b) Validation of the prognostic value of the CIHI for OS and DFS in the SYSUCC cohort. Fig S6. (a, b) Univariate and multivariate Cox regression analyses of the CIHI and clinical indicators. Fig S7. (a, b) Correlation between the CIHI and tumor size in the METABRIC cohort or AJCC T stage in the TCGA cohort. Fig S8. (a) Box and whisker plots showing the expression of the selected hypoxia-related genes in the FUSCCTNBC cohort. Fig S9. Correlation between the CIHI and previously reported hypoxia scores in the TCGA-TNBC cohort. Table S1. Primers used in this study. Table S2. Reagents used in this study. Table S3. Correlation between CIHI and xCell results. Appendix S1. Cell cycle correlation. Appendix S2. Cell proliferation correlation.