- Research
- Open access
- Published:
SLC26A4 correlates with homologous recombination deficiency and patient prognosis in prostate cancer
Journal of Translational Medicine volume 20, Article number: 313 (2022)
Abstract
Background
Homologous recombination deficiency (HRD) is closely associated with patient prognosis and treatment options in prostate cancer (PCa). However, there is a lack of quantitative indicators related to HRD to predict the prognosis of PCa accurately.
Methods
We screened HRD-related genes based on the HRD scores and constructed an HRD cluster system to explore different clinicopathological, genomic, and immunogenomic patterns among the clusters. A risk signature, HRDscore, was established and evaluated by multivariate Cox regression analysis. We noticed that SLC26A4, a model gene, demonstrated unique potential to predict prognosis and HRD in PCa. Multi-omics analysis was conducted to explore its role in PCa, and the results were validated by qRT-PCR and immunohistochemistry.
Results
Three HRD clusters were identified with significant differences in patient prognosis, clinicopathological characteristics, biological pathways, immune infiltration characteristics, and regulation of immunomodulators. Further analyses revealed that the constructed HRDscore system was an independent prognostic factor of PCa patients with good stability. Finally, we identified a single gene, SLC26A4, which significantly correlated with prognosis in three independent cohorts. Importantly, SLC26A4 was confirmed to distinguish PCa (AUC for mRNA 0.845; AUC for immunohistochemistry score 0.769) and HRD (AUC for mRNA 0.911; AUC for immunohistochemistry score 0.689) at both RNA and protein levels in our cohort.
Conclusion
This study introduces HRDscore to quantify the HRD pattern of individual PCa patients. Meanwhile, SLC26A4 is a novel biomarker and can reasonably predict the prognosis and HRD in PCa.
Introduction
Prostate cancer (PCa), a common malignant tumor, is the second leading cause of cancer-related mortality in men worldwide [1]. It can metastasize to bone (80%-100%), lymph nodes, liver, adrenal gland, or lung [2]. Although most early localized prostate cancer can achieve satisfactory results by prostatectomy or radiotherapy with a 5-year survival rate of 98.9%, metastasis is mainly found on initial diagnosis, hampering the avenue to a good prognosis. Therapy for metastatic PCa remains limited, and the current standard therapy is androgen deprivation therapy (ADT) combined with chemotherapy [3, 4]. Although ADT is initially effective, most patients inevitably develop into lethal metastatic castration-resistant prostate cancer (mCRPC) within 2–3 years [5], and the 5-year survival rate of them is only 28.2% [1]. Accordingly, there is significant enthusiasm to improve the stratification of patients with prostate cancer so that high-risk patients can be found earlier and receive active treatment.
The homologous recombination pathway plays a vital role in DNA repair and involves many genes [6], including BRCA (BRCA1/2), ATM, CHEK2, etc. Accumulated evidence has revealed the value of homologous recombination deficiency (HRD) in PCa, representing a high risk of PCa carcinogenesis and aggressiveness. A quarter of patients with recurrent or advanced PCa carry germline or somatic mutations in HRD-related genes [7]. The most commonly altered HRD-related gene in prostate cancer is BRCA2, with a prevalence of 5–6% at the germline level in mCRPC patients [8, 9]. A previous study revealed that BRCA2 mutation carriers have a 5.0 to 8.6-fold increased risk and a 15% absolute risk of developing PCa [10, 11]. Moreover, BRCA2 mutation carriers have higher progression rates from local to systemic disease, higher Gleason scores, shorter metastasis-free survival, and lower overall survival rates when compared to non-carriers [12,13,14]. In general, HRD is closely associated with a worse prognosis in PCa.
By extracting the HRD scores and other information from The Cancer Genome Atlas prostate adenocarcinoma cohort (TCGA-PRAD), we established an HRD signature to distinguish between high-risk and low-risk PCa patients. Through in-depth analysis, we identified and validated the protective effect of Solute Carrier Family 26 Member 4 (SLC26A4) in PCa, which may guide the application of poly(ADP-ribose) polymerase (PARP) inhibitors in PCa complementary to the commonly HRD-related gene mutations.
Material and methods
Prostate cancer datasets and preprocessing
Three open datasets with prostate cancer samples, multi-omics data, and complete clinical information were retrieved from the Cancer Genome Atlas (TCGA), Memorial Sloan Kettering Cancer Center (MSKCC), and Gene-Expression Omnibus (GEO) databases on August 22, 2021, including TCGA-PRAD [15], MSKCC-PRAD [16], and GSE116918 [17] cohorts. Then fragments per kilobase of exon model per million mapped fragments (FPKM) values were transformed into transcripts per kilobase million (TPM) values and log-transformed. The HRD, including loss of heterozygosity (LOH), telomeric allelic imbalance (TAI), and large-scale state transitions (LST), as well as gene-level copy numbers, PARADIGM integrated pathways, immune subtypes, gene-level non-silent mutation, were downloaded from Pan-Cancer (PANCAN) cohort in UCSC Xena (https://xenabrowser.net/) [18, 19]. Patients in the TCGA-PRAD cohorts without specific HRD scores were excluded for further analysis.
Profiling of HRD-related genes
The HRD scores and genome-wide DNA damage footprints were updated on June 13, 2017. Since then, patients in the TCGA-PRAD cohorts without specific HRD scores were excluded for further analysis. We quartered patients in the TCGA-PRAD cohort according to the HRD scores. Quarters 1 and 4 were defined as the bottom HRD group and top HRD group, respectively. Differential analysis was performed based on the transcriptomic data of the two groups using the “limma” R package. Genes with | log2(fold change) |> 0.5 and p value < 0.05 were selected for subsequent univariate Cox analysis, and those significantly correlating with patient progression-free interval (PFI or PFS) were defined as HRD-related genes. Their mutational and expressional profiles were investigated. We also calculated their Spearman’s correlations based on their mRNA expression levels and displayed it as an intra-correlation plot.
Unsupervised clustering for HRD-related genes
Unsupervised clustering analysis was applied to identify distinct HRD patterns based on the expression of the above prognostic HRD-related genes and classify patients for further analysis. The consensus clustering algorithm determined the number of clusters and their stability. We used the ConsensuClusterPlus package to perform the above steps, and 1000 repetitions were conducted to guarantee the stability of classification [20].
The mRNA expression level of each HRD-related gene was depicted among the clusters. Principal Component Analysis (PCA) and Kaplan–Meier survival analysis were performed to assess the power of clustering. The distributions of clinicopathological characteristics, including age at diagnosis, Gleason score, primary outcome, biochemical recurrence (BCR), pathologic T stage, pathologic N stage, original zone of cancer, and immune subtype, were evaluated across the clusters.
Pathway quantification at transcriptomic and proteomic levels
The PARADIGM algorithm integrates pathway, expression, and copy number data to infer activation of pathway features within a superimposed pathway (SuperPathway) network structure. The SuperPathway system comprises 1387 constituent pathways from three pathway databases, NCI-PID, BioCarta, and Reactome (last updated 05/2013), containing 19K pathway features, representing 7369 genes, 9354 complexes, 2092 families, 82 RNAs, 15 miRNAs, and 592 abstract processes. This dataset is ssGSEA scores for 1387 constituent pathways [19].
Reverse-phase protein array (RPPA) data from the PANCAN cohort were used to calculate the pathway activity score of 10 cancer-related pathways. RPPA is a high-throughput antibody-based technique with procedures like Western blots. Proteins are extracted from tumor tissue or cultured cells, denatured by SDS, printed on nitrocellulose-coated slides, followed by an antibody probe. The terms included Apoptosis, Cell Cycle, DNA Damage Response, Epithelial-Mesenchymal Transition (EMT), Hormone a, Hormone b, PI3K/AKT, RTK, and TSC/mTOR pathways. In brief, RBN RPPA data were median-centered and normalized by the standard deviation across all samples for each component to obtain the relative protein level. The pathway activity score is then the sum of the relative protein level of all positive regulatory elements minus that of negative regulatory components in a particular pathway [21].
Estimation of tumor purity and fractions of immune cells
Estimation of stromal and immune components and tumor purity in tumor tissues using expression data was achieved by the “ESTIMATE” R package [22]. Subsequently, the population abundance (fraction) of tissue-infiltrating immune and stromal cell populations was estimated by three well-known algorithms, including MCP counter (10 cell types) [23], ImmuneCellAI (24 cell types) [24], and Cibersort (22 cell types) [25].
Essential molecular characteristics of the tumor
We extracted vital molecular features of malignant tumors from an integrated and in-depth bioinformatics study [26], including proliferation, leukocyte fraction, B cell receptor (BCR) evenness, T cell receptor (TCR) evenness, Th1, Th2, and Th17 cells, aneuploidy score, intratumor heterogeneity (ITH), single nucleotide variant (SNV) neoantigens, insertion-and-deletion (indel) neoantigens, cancer-testis antigen (CTA) score, homologous recombination defects, and fraction of altered genome. The microsatellite instability (MSI) MANTIS score was downloaded from cBioPortol for Cancer Genomics (https://www.cbioportal.org/).
Immunomodulator identification and analysis
A list of 78 immunomodulatory genes was obtained from a previous study that curated them from a literature review performed by immuno-oncology experts within the TCGA immune response working group [26]. Corresponding median mRNA expression levels were used to summarize expression in each cluster. We performed a limma differential analysis across clusters to examine differences in immunomodulatory gene expression and found genes to be significantly differentially expressed. And the immunomodulatory copy number was also outputted from a PANCAN cohort as deep amplifications (2), shallow amplifications (1), non-alterations (0), shallow deletions (− 1), and deep deletions (− 2) of each immunomodulator gene. Proportions of samples with each type of copy number alteration were then compared across HRD clusters.
Profiling of prognostic hub genes and dimensionality reduction
We performed differential expression analysis between pairs in this cohort of HRD clusters and performed Cox survival analysis after taking the intersection of the resulting differentially expressed genes. Those with survival significance were set as prognostic hub genes, whose expression patterns were employed as the basis of subsequent PCA analysis. The risk signature was termed as ‘HRDscore’ and calculated by the following formula:
where “\({expression}_{risk}\)” stood for expression levels of risk genes and “\({expression}_{protective}\)” stood for that of protective genes.
Patients were dichotomized into high HRDscore, and low HRDscore groups based on the best cut-off decided by X-tile software. A Sankey plot was established to investigate the intrinsic relationship among HRD cluster, immune subtype, and HRDscore. Furthermore, we explored the correlations between the HRDscore and clinicopathological features, including survival. For subgroup analysis, TCGA-PRAD patients were divided into different groups based on features as follows: age (≤ 45 years old or > 45 years old) and Gleason score (< 8 or ≥ 8). Finally, multivariate Cox analyses were conducted to test the robustness of the established HRDscore.
Prediction of immunotherapy response and correlation with immune cells
ImmuCellAI was used to predict the response of immune checkpoint blockade (ICB) therapy based on the transcriptomic data [24]. A receiver operating characteristic (ROC) curve was built to illustrate the power of HRDscore in predicting immunotherapy response.
We calculated the correlations between HRDscore and fractions of immune cells and the prognostic value of these cell types. Next, several genes were obtained after the intersection between HRD-related and prognostic hub genes. Their relationship to immune cells was also measured to find critical genes that bridge HRD scores, immune infiltration, and patient prognosis.
Quantitative real-time PCR assay
Quantitative real-time PCR was performed with SYBR Green PCR mixture (Using Roche lightcycler 480 system) according to standard protocols. PCR conditions were: one cycle of 5 min at 95 °C, then 45 cycles of 10 s at 95 °C, 10 s at 60 °C, 10 s at 72 °C. The expression of the SLC26A4 gene was normalized to the expression of the GAPDH gene using the comparative CT method. Primers used were: SLC26A4 (F: 5′-AGGAAATATGCACTGCTCACT- 3′; R: 5′-AGTATTCCCGCAGTTTGCTGA-3′); GAPDH (F: 5′-CAAGGCTGAGAACGGGAAG-3′; R: 5′-TGAAGACGCCAGTGGACTC-3′).
Prostate cancer samples and immunohistochemistry
Prostate cancer samples were acquired from Xiangya Hospital of Central South University. A physician obtained informed consent from the patients. The procedures related to human subjects were approved by the Ethics Committee of Xiangya Hospital, Central South University. Tissues were fixed in 10% buffered Formalin, then transferred to 70% alcohol. These paraffin-embedded tissues were sectioned (4 μm) and stained with antibodies against SLC26A4 (HPA042860, Atlas Antibodies). The following detection and visualization procedures were performed according to the manufacturer's protocol. To quantify the immunohistochemistry (IHC) result of positive staining, five random areas in each tissue sample were microscopically examined and analyzed by an experienced pathologist. The average staining score was calculated by dividing the positive areas by entire regions.
Statistical analyses
The univariate and multivariate Cox analyses were performed to detect the prognostic factors. Kaplan–Meier curves with the log-rank test were used to assess survival differences between groups. Spearman correlation analyses were used to calculate correlations. The cutoff value was determined using the X-tile software (version 3.6.1). All statistical analyses were conducted using R software (version 4.1.2), and most visualization was achieved using the “ggplot2” R package. P < 0.05 was considered statistically significant.
Results
The landscape of genetic variation of HRD-related genes in prostate cancer
The flow chart of our study is summarized in Fig. 1. We first ranked the patients in the TCGA-PRAD cohort (n = 472) in order of HRD scores from high to low, and then subjected their transcriptomic data in the first and fourth quartiles (quartile 4 vs. quartile 1) to differential expression analysis (Additional file 4: Table S1; Fig. 2A, B). Sixty-six differential expressed genes were obtained, and subsequent univariate Cox analysis found 23 of them were prognostic (p < 0.05; Additional file 5: Table S2; Fig. 2C), defined as HRD-related genes.
Of the 23 HRD-related genes, most were protective factors for PFS of patients with prostate cancers, except KRT5 and DDC. Moreover, SLC26A4, KRT16, COL17A1, and AQP9 were genetically unstable. And briefly, the SNV frequencies in MYOT, COL17A1, AQP9, KRT16, and SLC26A4 were high, equal to or exceeding 10% (Fig. 2D). As for the copy number variation (CNV) status, BMPER, SLC26A4, and DDC shared a high amplification frequency, while KRT16, COL17A1, and AQP9 demonstrated a high deletion frequency (Additional file 6: Table S3; Fig. 2E). And most genes were differentially expressed between prostate cancer tissues and adjacent normal controls (Fig. 2F). The intra-correlation among these HRD-related genes was illustrated in Fig. 2G (Additional file 7: Table S4).
Molecular patterns mediated by 23 HRD-related genes
An unsupervised clustering grouped patients from the TCGA-PRAD cohort into three distinct clusters based on the expression patterns of 23 HRD associated genes, including 179, 246, and 47 samples in the HRD clusters 1 to 3 respectively (Fig. 3A, B, E; Additional file 8: Table S5). The expression patterns of HRD-related genes and the distribution of clinicopathological features in the context of HRD clusters are displayed in Fig. 3C. We found that all genes were expressed at the lowest levels in HRD cluster 1 in both the heatmap and boxplot (Fig. 3C, D). Although the HRD cluster failed to distinguish patients’ overall survival (p = 0.53; Fig. 3F), cases in cluster 1 demonstrated the worst progression-free interval (p = 0.0021; Fig. 3G). Moreover, we clustered samples based on the expressions of 23 differentially expressed genes in the MSKCC-PRAD cohort and generated three clusters. Consistently, cluster 1 demonstrated the worst PFI (Additional file 1: Fig. S1).
HRD clusters distinguish clinicopathological characteristics and biological pathways
Chi-square analyses among the three HRD clusters revealed that patients in cluster 1 harbored higher age and higher Gleason scores (Fig. 4A, B). Furthermore, prostate cancer samples in HRD cluster 1 showed more malignant properties as they displayed lower proportions of complete or partial response (CR/PR) but higher ratios of biochemical recurrence (Fig. 4C, D), as well as higher pathologic T and N stages (Fig. 4E, F).
Compared with clusters 2 and 3, HRD cluster 1 showed higher scores of pathways, including Metastasis, cytokine production during the immune response, T cell proliferation during the immune response, but lower scores of PARADIGM pathways such as cell cycle, macrophage activation, DNA damage, and natural killer cell-mediated cytotoxicity (Fig. 4G, H). Importantly, quantification of activated pathways using protein expression data also demonstrated that our HRD cluster could well reflect genetic alterations, as it uniquely distinguished patients regarding the pathway “DNA damage response” (p < 0.001; Fig. 4I).
Immune infiltration characteristics in distinct HRD clusters
HRD cluster 1 distinctly showed the lowest stromal scores, immune scores, and highest tumor purity (Fig. 5A–C; Additional file 9: Table S6). Infiltration abundance estimated by ImmuneCellAI showed a consistent result, as cluster 1 had the lowest infiltration scores (Fig. 5D). Specifically, the fractions of every stromal or immune cell type were lowest in HRD cluster 1 (Fig. 5E–G; Additional file 9: Table S6). But there was a paradoxical situation where we observed that Th1, Th17, central memory cells, and macrophages were most abundant in cluster 1 (based on the ImmuneCellAI algorithm; Fig. 5F). Meanwhile, the fraction of M2 macrophage was highest in cluster 1 (based on the Cibersort algorithm; Fig. 5G).
The HRD clusters showed distinct genetic and immune signatures based on the dominant sample characteristics of their tumor samples (Fig. 5H, I). HRD cluster 1 harbored the highest tumor mutational burden as it had the highest aneuploidy scores and SNV neoantigen counts (Fig. 5H, I). Meanwhile, cluster 1 was genomically unstable since we found that fractions of genomic alteration and MSI Mantis scores were uniquely highest in cluster 1 (Fig. 5H, I). Regarding immune infiltration (signature), HRD cluster 1 showed the highest infiltrating abundance of M2 macrophage but the lowest cytotoxic T cells (Fig. 5H). And it showed the highest Th1/Th2 ratio bias to the adaptive immune infiltrate (Fig. 5H). Additionally, HRD cluster1 had the highest proliferation rate (Fig. 5H, I). The clonal evenness of TCR and BCR, the count of indel neoantigens, and CTA scores were not significantly different among the three HRD clusters (Fig. 5I).
Regulation of immunomodulators
Immunomodulators (IMs) are critical for cancer immunotherapy, with numerous IM agonists and antagonists being evaluated in clinical oncology. To advance this research, understanding their expression and modes of control in different states of the tumor microenvironment (TME) is needed. We examined IM gene expression, CNVs, and SNVs.
Gene expression of IMs (Additional file 10: Table S7, Additional file 2: Fig. S2A) varied across HRD clusters, perhaps indicating their role in shaping the TME. Generally, most genes encoding IMs were at low expression levels in HRD cluster 1. Genes with the most significant differences between clusters (Additional file 2: Fig. S2B) included CX3CL1 (BH-adjusted p < 10–5), most lowly expressed in cluster 1 and TNFSF4 (BH-adjusted p = 0.004), most highly expressed in cluster 1.
Copy-number variations affected multiple IMs and varied across HRD clusters. Cluster 1 showed both frequent amplification and deletion of IM genes, consistent with their greater genomic instability. In particular, BTLA was most frequently amplified in cluster 1, while TIGIT deletion was enriched in cluster 3 (Additional file 2: Fig. S2C). Overall, these marked differences in IM copy number may reflect more direct modulation of the TME by cancer cells. The observed differences in regulation of IMs might have implications for therapeutic development and combination immune therapies, and the multiple mechanisms at play in evoking them further highlight their biological importance.
Dimensionality reduction and construction of the HRD signature
To accurately quantify the prediction power of the HRD clustering, we applied a methodology to establish an HRD signature and calculate HRDscore for all the patients with prostate cancers. We obtained 64 prognostic hub genes when intersecting three group-paired differential expression analyses (DEGs) (Fig. 6A; Additional file 11: Table S8). Using the formula introduced in the “Method” section, we calculate the HRDscore for each patient based on the transcriptomic data. And we found that patients with higher HRDscores (threshold determined by X-tile) had a higher tendency of PFI event (recurrence, metastasis, or death) (55/360 vs. 33/112; p = 0.0013) (Fig. 6B). A Sankey plot was established to understand better the intrinsic relationship between HRD cluster and HRD signature (HRDscore) and other features for individuals (Fig. 6C). It intuitively showed that patients with higher HRDscores mostly came from HRD cluster 1, and they were predominantly classified with immune subtypes C3 (Inflammatory) and C4 (lymphocyte depleted). Subsequent specific analyses found that patients with higher Gleason scores (≥ 8), pathologic T and N stages, and those with PFI events had higher HRDscores (Fig. 6D), suggesting the ability of our HRD signature to represent malignant features.
An important finding of our study was that the established HRDscore could distinguish the survival outcome of patients with prostate cancers. Patients with higher HRDscore had worse prognosis when compared to those with lower scores, in the context of overall survival (OS) (p = 0.049), disease free survival (DSS) (p = 0.025), disease free interval (DFI) (p = 0.0013), and PFI (p < 0.0001) (Fig. 6E). Subgroup analyses stratified by patient’s age and Gleason score confirmed the stability of our HRD signature, as higher HRDscores consistently correlated with unfavorable prognosis (Fig. 6F). Furthermore, multivariate Cox analyses revealed the robustness of our signature since it showed that Gleason score, pathologic T stage, and HRDscore were independent predictors of patient’s DFI and PFI (Fig. 6H, I).
HRDscore correlates with immunotherapy response, genomic instability, and immune infiltration
Patients previously assigned to HRD cluster 1 had the highest level of HRDscore (Fig. 7A). Using the bulk transcriptomic data, we acquired patient responses to ICB by ImmuneCellAI algorithm (Additional file 12: Table S9). High-HRDscore group demonstrated higher response rate to ICB than the low-HRDscore group (9% vs 4%, p = 0.034) (Fig. 7B). Besides, HRDscore outperformed other indicators in ICB response prediction with AUC = 70.82 (Fig. 7C). To better illustrate the characteristics of the HRD signature, we also tested the correlation between the known signatures and the HRDscore (Additional file 13: Table S10). After using the cluster method “ward.D2”, our HRD signature was clustered with recognized signatures related to genomic instability. In brief, it was strongly correlated with Fanconi anemia, cell cycle, DNA damage repair, Nucleotide excision repair, Homologous recombination, DNA replication, and mismatch repair (Fig. 7D).
Next, we explored the relationship between HRD signature and immune infiltration (Fig. 7E, Table). The HRDscore was positively correlated with central memory cells, Th1, Th17, macrophage, and natural killer T-cell (NKT). In contrast, it was negatively associated with B cell, CD8 naïve cell, dendritic cell (DC), and exhausted T cell. Univariate Cox analysis found that CD4 naïve T cell, exhausted T cell, and Th2 were protective cells, whereas Th1 was uniquely related to a poor prognosis in PCa patients (Additional file 14: Table S11). Specifically, patients with higher Th1 abundance had poor survival outcomes (Fig. 7F), and the effects were more significant when combined with the HRDscore group (Fig. 7G).
To better focus on a single gene factor, we intersected HRD-related genes with prognostic hub genes, resulting in four candidate genes: actin alpha 1, skeletal muscle (ACTA1), keratin 6A (KRT6A), orosomucoid 2 (ORM2), and solute carrier family 26 member 4 (SLC26A4). We examined their links with 22 immune cells (Fig. 7H) and found that only SLC26A4 was significantly but not strongly correlated with the four prognostic immune cell types (CD4 naïve T cell, exhausted T cell, Th1, and Th2). The scatter plots illustrating correlations are displayed in Fig. 8A.
SLC26A4 serves as a critical gene and correlates with immune infiltration and clinical prognosis
SLC26A4 was positively correlated with protective immune cells (CD4 naïve T cell, exhausted T cell, and Th2), while was negatively correlated with risk cell type (Th1) (Fig. 8A). Importantly, SLC26A4 consistently acted as a protective factor for PFI in TCGA-PRAD (p = 0.0003), MSKCC (p < 0.0001), and GSE116918 (p = 0.0014) cohorts (Fig. 8B). A meta-analysis revealed that, compared to low-SLC26A4 group, high-SLC26A4 group had a higher risk for PFI events (risk ratio = 0.39, 95% CI 0.29–0.53). And the pooled estimate showed no heterogeneity with I2 = 0% and tau-square = 0 (Fig. 8C).
Pan-cancer analysis focusing on SLC26A4 found that it was differentially expressed in most cancer types (Fig. 8D). Elevated in cholangiocarcinoma (CHOL) and liver hepatocellular carcinoma (LIHC) compared to corresponding normal tissues. Whereas SLC26A4 expression level was decreased in various tumor tissues in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), rectum adenocarcinoma (READ), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC). Moreover, SLC26A4 was a risk for uveal melanoma (UVM) and LGG. It was protective for UCEC, THCA, SKCM, PRAD, and KIRC (Fig. 8E).
To better understand the function of this gene in cancer, we next obtained expression data and functional state scores at the single-cell level in the CancerSEA database (Additional file 15: Table S12). In all available cohorts, SLC26A4 has negatively correlated “DNA damage” and “invasion” functions (Fig. 8F). Unfortunately, the prostate cancer single-cell cohort there didn’t provide the expression data for SLC26A4, hampering our understanding of the correlations in prostate cancer. Finally, DNA methylation of SLC26A4 (by Illumina HumanMethylation450 BeadChip) was increased in prostate cancer tissues than the normal controls (Additional file 3: Fig. S3). The methylation was inversely correlated with its mRNA expression in TCGA-PRAD cohort (R = -0.60, p = 2.85e−50; Fig. 8G), suggesting epigenetic silencing. The methylation probe ID was cg15320854, and the methylation site was cpg 107660494.
SLC26A4 was down-regulated in prostate cancer samples with HRD in independent external validation
Among the 62 included patients enrolled in Xiangya Hospital, nine (14.5%) were diagnosed with benign prostatic hyperplasia (BPH). Of the remaining 53 patients with prostate cancer, 33 carry HRD mutation (Additional file 16: Table S13).
Compared to BPH tissues, SLC26A4 mRNA expression levels and IHC scores were decreased in prostate cancer samples (Fig. 9A). SLC26A4 mRNA surpassed baseline total prostate specific antigen (PSA) value in predicting prostate cancer (AUC = 0.845; Fig. 9B). Furthermore, SLC26A4 was significantly down-regulated in prostate cancer tissues with HRD than those without HRD at mRNA and protein levels (Fig. 9C). And SLC26A4 demonstrated excellent performance in predicting HRD in the context of prostate cancer (AUC = 0.911; Fig. 9D). Representative IHC results in three tissue types are shown in Fig. 9E–G.
Discussion
For many years, people have been exploring the initiation, development, and treatment of PCa. Gleason score and serum PSA level are still the most important prognostic factors of PCa. Recently, increasing evidence has suggested that HRD plays a key role in the biological process and therapeutic response in various tumors, one of the most influential factors for the prognosis of tumors [27]. In addition, many studies have found that mCRPC patients with HRD-related gene mutations show impressive responses to PARP inhibitors, even in very advanced disease settings [28,29,30,31,32]. The common HRD-related gene mutations in PCa are BRCA2, ATM, and CHEK2, all of which are included in the molecular eligibility criteria of virtually all PARP inhibitor trials involving mCRPC patients [33]. However, their germline mutations were found in 5–6%, 1–2% and 1–2% of mCRPC patients, respectively [8, 9].Therefore, new biomarkers need to be developed for molecular typing of PCa patients. In this study, we deeply analyzed the molecular characteristics of PCa patients with different HRD scores and identified a biomarker that could be complementary to the HRD scores.
HRD score integrates three indicators focusing on DNA-based genomic instability, which has been less explored in prostate cancer. The previous study has found that patients with primary prostate cancer have lower HRD scores, while patients with germline BRCA2 mutations have higher HRD scores [34]. Since BRCA2 mutation is the indication of two PARP inhibitors recently approved by the Food and Drug Administration (FDA) for the treatment of PCa, HRD score analysis may help improve treatment options. In this study, by analyzing the HRD scores of the TCGA-PRAD cohort, we obtained 23 genes associated with HRD scores, defined as HRD-related genes. We identified three molecular patterns with distinct clinicopathological characteristics based on these genes, and HRD cluster 1 was particularly correlated with worse clinicopathological types and poor prognosis.
The HRD clusters demonstrated distinct immune landscapes. In general, T cells, CD8 T cells, cytotoxic lymphocytes, and natural killer (NK) cells were less enriched in HRD cluster 1 indicating that the inhibited immune response may explain the poor outcome of patients in cluster 1. The HRD clusters also showed different immunogenomic characteristics. Specifically, HRD cluster 1 harbored the highest mutational burden, highest proliferation potentials, and lowest genomic stability, indicating an absolute potential to derive mutation and subsequent carcinogenesis. Besides, HRD cluster 1 demonstrated the lowest leukocyte abundance but the highest M2 macrophage infiltration. Macrophage infiltration in solid tumors is associated with poor prognosis [35]. The previous study has found that macrophages infiltrating PCa were mainly M2 type and associated with invasiveness and unfavorable outcome. We also noticed that Th1 to Th2 ratio was lowest in cluster 1. Cellular immunity mediated by Th1 mainly plays an anti-tumor role. Once it shifts from Th1 to Th2, resulting in immunosuppression. Thus, the anti-tumor immunity of the body will be seriously disturbed. Yamamura et al. and Kharkevitch et al. first found that Th2 cells were dominant in tumor patients [36, 37], and then found that Th2 shift occurred in many types of tumors such as non-small cell lung cancer, choriocarcinoma, glioma, gastric cancer, ovarian cancer, melanoma, colorectal cancer, and lymphoma. The above results define an immunosuppressive microenvironment phenotype of prostate cancer and an unstable genomic condition in HRD cluster 1.
Given the importance of IMs in cancer immunotherapy, we compared the differences in IM gene expression between these three clusters. The genes with the most obvious difference among clusters. In general, most of the IMs were in a relatively low expression level in HRD cluster 1 than those in clusters 2 and 3, suggesting that immune responses regulated by membrane checkpoints were less common there. Consistently, copy number variations of IMs were more frequent in HRD cluster 1 in amplification and deletion, confirming the unstable genomic phenotype. Although such a trend was not evident in SNV, several untypical checkpoints still had higher variation frequencies like GZMA, PRF1, ENTPD1, and ARG1. Paradoxically, TNFSF4 was significantly up-regulated in HRD cluster 1. Recent studies have shown that stimulation of OX40, the ligand of TNFSF4, is helpful for therapeutic immunization strategies for cancer [38]. It has been found that TNFSF4 is enriched in bone metastatic PCa [39]. Combined with our results, it may serve as a new therapeutic target in PCa, especially for those patients with high expression.
Furthermore, we established a signature (termed HRDscore) with excellent power to predict prognosis with stability. Based on proteomic data, the HRDscore was tightly correlated with existing signatures related to genomic instability, including homologous recombination, DNA damage repair, and Fanconi anemia (correlation coefficient ≥ 0.5, p < 0.001). This result suggested that the HRD-derived risk system could represent the signature of genomic defects. Besides, the HRDscore was positively related to macrophages, the unfavorable cell type, which was consistent with the above suppose. A recent article has explored HRD scores in PCa, which focused on the correlation between HRD scores and mutations of BRCA2 and ATM [34]. However, these mutations are not common in PCa, especially in non-mCRPC, so the HRD score is of little value to numerous PCa patients without these mutations. In comparison, our HRDscore has excellent value for predicting the prognosis and even guiding treatment in PCa.
To further explore valuable biomarkers, we finally focused on a single gene, SLC26A4, correlated with immune infiltration and clinical diagnosis. It showed protective effects in several independent PRAD cohorts (RR 0.39, 95% CI 0.29–0.93, I2 = 0). Functional single-cell analysis suggested that SLC26A4 was negatively correlated with "DNA damage" and "invasion" functions. Nevertheless, the lack of prostate cancer single-cell cohort with SLC26A4 expression data hampered our understanding of its functions in PCa. Previous studies have mostly believed that SLC26A4 plays a vital role in maintaining normal hearing and never explored its significance in malignancies [40]. Our study revealed its potential value in tumorigenesis and development for the first time, which is worthy of in-depth exploration in future research.
SLC26A4 encodes a membrane protein called pendrin that permits the anion exchange between the cytosol and extracellular space, maintaining the proper function of auditory sensory cells. It is mainly expressed in the inner ear and thyroid gland, and its mutation is related to dyshormonogenic goiter and Pendred syndrome [41, 42]. Hypermethylation of SLC26A4 often occurs in cancers such as thyroid cancer and acute myoid leukemia [43, 44], consistent with our results. All the above findings indicated that the epigenetic changes of SLC26A4 may be involved in tumorigenesis. Our study uniquely found that SLC26A4 was highly associated with HRD in prostate cancer.
In our own Xiangya cohort, the SLC26A4 expression in PCa samples was lower than that in benign prostatic hyperplasia tissues at both mRNA and protein levels, which was inconsistent with the results of the TCGA-PRAD cohort. This may be due to the insufficient sample size of our cohort. Therefore, it needs to be further confirmed. Importantly, we found that SLC26A4 performed well in predicting HRD in patients with PCa. Patients with HRD-related gene mutations are often sensitive to PARP inhibitors, so we proposed that SLC26A4 may be a novel biomarker to screen patients sensitive to PARP inhibitors.
Consequently, we herein provided a potential biomarker for the treatment of PCa with PARP inhibitors. However, several limitations should be addressed in our study. First, there is a lack of SLC26A4 expression data in the prostate cancer single-cell cohort, which has been mentioned above. Secondly, our analyses were also limited by the relatively small sample size. Finally, due to the lack of prognostic and treatment information in the cohort, we failed to thoroughly verify the value of SLC26A4 in suggesting prognosis and guiding treatment. Therefore, further validation based on a large cohort is warranted.
Conclusion
We introduced HRDscore to quantify the HRD pattern of individual PCa patients, which can predict the prognosis of PCa with stability and universality and has a specific value in guiding treatment. A new biomarker, SLC26A4, plays a protective role in PCa and can screen patients suitable for PARP inhibitor treatment.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
References
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68(1):7–30.
Bubendorf L, Schopfer A, Wagner U, Sauter G, Moch H, Willi N, et al. Metastatic patterns of prostate cancer: an autopsy study of 1,589 patients. Hum Pathol. 2000;31(5):578–83.
Mohler JL, Antonarakis ES, Armstrong AJ, D’Amico AV, Davis BJ, Dorff T, et al. Prostate cancer, Version 2.2019, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw. 2019;17(5):479–505.
Cornford P, van den Bergh RCN, Briers E, Van den Broeck T, Cumberbatch MG, De Santis M, et al. EAU-EANM-ESTRO-ESUR-SIOG Guidelines on Prostate Cancer. Part II-2020 update: treatment of relapsing and metastatic prostate cancer. Eur Urol. 2021;79(2):263–82.
Zong Y, Goldstein AS. Adaptation or selection–mechanisms of castration-resistant prostate cancer. Nat Rev Urol. 2013;10(2):90–8.
Hoppe MM, Sundar R, Tan DSP, Jeyasekharan AD. Biomarkers for homologous recombination deficiency in cancer. J Natl Cancer Inst. 2018;110(7):704–13.
Chung JH, Dewal N, Sokol E, Mathew P, Whitehead R, Millis SZ, et al. Prospective comprehensive genomic profiling of primary and metastatic prostate tumors. JCO Precis Oncol. 2019. https://doi.org/10.1200/PO.18.00283.
Pritchard CC, Mateo J, Walsh MF, De Sarkar N, Abida W, Beltran H, et al. Inherited DNA-repair gene mutations in men with metastatic prostate cancer. N Engl J Med. 2016;375(5):443–53.
Nicolosi P, Ledet E, Yang S, Michalski S, Freschi B, O’Leary E, et al. Prevalence of germline variants in prostate cancer and implications for current genetic testing guidelines. JAMA Oncol. 2019;5(4):523–8.
Mersch J, Jackson MA, Park M, Nebgen D, Peterson SK, Singletary C, et al. Cancers associated with BRCA1 and BRCA2 mutations other than breast and ovarian. Cancer. 2015;121(2):269–75.
De Felice F, Tombolini V, Marampon F, Musella A, Marchetti C. Defective DNA repair mechanisms in prostate cancer: impact of olaparib. Drug Des Devel Ther. 2017;11:547–52.
Castro E, Goh C, Leongamornlert D, Saunders E, Tymrakiewicz M, Dadaev T, et al. Effect of BRCA mutations on metastatic relapse and cause-specific survival after radical treatment for localised prostate cancer. Eur Urol. 2015;68(2):186–93.
Castro E, Goh C, Olmos D, Saunders E, Leongamornlert D, Tymrakiewicz M, et al. Germline BRCA mutations are associated with higher risk of nodal involvement, distant metastasis, and poor survival outcomes in prostate cancer. J Clin Oncol. 2013;31(14):1748–57.
Gallagher DJ, Gaudet MM, Pal P, Kirchhoff T, Balistreri L, Vora K, et al. Germline BRCA mutations denote a clinicopathologic subset of prostate cancer. Clin Cancer Res. 2010;16(7):2115–21.
Cancer Genome Atlas Research N. The molecular taxonomy of primary prostate cancer. Cell. 2015;163(4):1011–25.
Taylor BS, Schultz N, Hieronymus H, Gopalan A, Xiao Y, Carver BS, et al. Integrative genomic profiling of human prostate cancer. Cancer Cell. 2010;18(1):11–22.
Jain S, Lyons CA, Walker SM, McQuaid S, Hynes SO, Mitchell DM, et al. Validation of a metastatic assay using biopsies to improve risk stratification in patients with prostate cancer treated with radical radiation therapy. Ann Oncol. 2018;29(1):215–22.
Goldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38(6):675–8.
Hoadley KA, Yau C, Hinoue T, Wolf DM, Lazar AJ, Drill E, et al. Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer. Cell. 2018;173(2):291–304.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
Akbani R, Ng PK, Werner HM, Shahmoradgoli M, Zhang F, Ju Z, et al. A pan-cancer proteomic perspective on The Cancer Genome Atlas. Nat Commun. 2014;5:3887.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17(1):218.
Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, et al. ImmuCellAI: a unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy. Adv Sci (Weinh). 2020;7(7):1902880.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.
Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity. 2018;48(4):812-30.e14.
Nguyen L, Martens JWM, Van Hoeck A, Cuppen E. Pan-cancer landscape of homologous recombination deficiency. Nat Commun. 2020;11(1):5584.
Mutetwa T, Foulkes WD, Polak P. Olaparib for metastatic castration-resistant prostate cancer. N Engl J Med. 2020;383(9):890.
Abida W, Campbell D, Patnaik A, Shapiro JD, Sautois B, Vogelzang NJ, et al. Non-BRCA DNA damage repair gene alterations and response to the parp inhibitor rucaparib in metastatic castration-resistant prostate cancer: analysis from the Phase II TRITON2 Study. Clin Cancer Res. 2020;26(11):2487–96.
Smith MR, Sandhu SK, Kelly WK, Scher HI, Efstathiou E, Lara P, et al. Phase II study of niraparib in patients with metastatic castration-resistant prostate cancer (mCRPC) and biallelic DNA-repair gene defects (DRD): preliminary results of GALAHAD. J Clin Oncol. 2019;37(7_suppl):202.
Bono JSD, Mehra N, Higano CS, Saad F, Buttigliero C, Mata M, et al. TALAPRO-1: A phase II study of talazoparib (TALA) in men with DNA damage repair mutations (DDRmut) and metastatic castration-resistant prostate cancer (mCRPC)—First interim analysis (IA). J Clin Oncol. 2020;38(6_suppl):119.
Abida W, Patnaik A, Campbell D, Shapiro J, Bryce AH, McDermott R, et al. Rucaparib in men with metastatic castration-resistant prostate cancer harboring a BRCA1 or BRCA2 gene alteration. J Clin Oncol. 2020;38(32):3763–72.
Marshall CH, Antonarakis ES. Therapeutic targeting of the DNA damage response in prostate cancer. Curr Opin Oncol. 2020;32(3):216–22.
Lotan TL, Kaur HB, Salles DC, Murali S, Schaeffer EM, Lanchbury JS, et al. Homologous recombination deficiency (HRD) score in germline BRCA2- versus ATM-altered prostate cancer. Mod Pathol. 2021;34(6):1185–93.
Cassetta L, Pollard JW. Targeting macrophages: therapeutic approaches in cancer. Nat Rev Drug Discov. 2018;17(12):887–904.
Yamamura M, Modlin RL, Ohmen JD, Moy RL. Local expression of antiinflammatory cytokines in cancer. J Clin Invest. 1993;91(3):1005–10.
Kharkevitch DD, Seito D, Balch GC, Maeda T, Balch CM, Itoh K. Characterization of autologous tumor-specific T-helper 2 cells in tumor-infiltrating lymphocytes from a patient with metastatic melanoma. Int J Cancer. 1994;58(3):317–23.
Lu X. OX40 and OX40L interaction in cancer. Curr Med Chem. 2021;28(28):5659–73.
Ihle CL, Provera MD, Straign DM, Smith EE, Edgerton SM, Van Bokhoven A, et al. Distinct tumor microenvironments of lytic and blastic bone metastases in prostate cancer patients. J Immunother Cancer. 2019;7(1):293.
Honda K, Griffith AJ. Genetic architecture and phenotypic landscape of SLC26A4-related hearing loss. Hum Genet. 2021;141(3–4):455–64.
Dossena S, Rodighiero S, Vezzoli V, Nofziger C, Salvioni E, Boccazzi M, et al. Functional characterization of wild-type and mutated pendrin (SLC26A4), the anion transporter involved in Pendred syndrome. J Mol Endocrinol. 2009;43(3):93–103.
Cangul H, Liao XH, Schoenmakers E, Kero J, Barone S, Srichomkwun P, et al. Homozygous loss-of-function mutations in SLC26A7 cause goitrous congenital hypothyroidism. JCI Insight. 2018;3(20):e99631.
Xing M, Tokumaru Y, Wu G, Westra WB, Ladenson PW, Sidransky D. Hypermethylation of the Pendred syndrome gene SLC26A4 is an early event in thyroid tumorigenesis. Cancer Res. 2003;63(9):2312–5.
Kroeger H, Jelinek J, Estécio MR, He R, Kondo K, Chung W, et al. Aberrant CpG island methylation in acute myeloid leukemia is accentuated at relapse. Blood. 2008;112(4):1366–73.
Acknowledgements
Not applicable.
Funding
This research was supported by the Key Research and Development program of Hunan Province (2021SK2014), the National Natural Science Foundations of China (81902606), and National Natural Science Foundation of Hunan Province (2020JJ5882, 2020JJ5891).
Author information
Authors and Affiliations
Contributions
Acquisition of data: CL, YZ, and WY. Analysis and interpretation of data: CL, ZL, and YG, Conception and design: CL, YZ, WY, and YC. Data curation: XG and XZ. Development of methodology: WY and YC. Writing the manuscript: CL and WY. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
The procedures related to human subjects have been performed in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Xiangya Hospital, Central South University. We have obtained informed consent to participate in the study from all participants.
Consent for publication
Not applicable.
Competing interests
No potential conflicts of interest were disclosed.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1: Figure S1.
Consensus clustering of HRD-related genes in the MSKCC-PRAD cohort. (Left) The consensus matrix when k = 3. (Right) The survival differences in three clusters regarding PFI.
Additional file 2: Figure S2.
The immune checkpoint landscape in distinct HRD clusters. (A) Regulation of immunomodulators in distinct HRD clusters. From left to right: mRNA expression (median normalized expression levels); amplification frequency (the difference between the fraction of samples in which an immunomodulator is amplified in a particular subtype and the amplification fraction in all samples); and the deletion frequency (as amplifications) for 75 immunomodulatory genes by HRD clusters. (B) Representative mRNA expression levels (CX3CL1 and TNFSF4) in distinct HRD clusters. (C) Representative CNV frequency (of BTLA and TIGIT) in distinct HRD clusters. The statistical significance is indicated as asterisks (*), *p < 0.05, ** p < 0.01, ***p < 0.001, ns, not significant.
Additional file 3: Figure S3.
The DNA methylation level of SLC26A4 in prostate cancer tissues and paired normal tissues in the TCGA-PRAD cohort.
Additional file 4: Table S1.1.
HRD scores of PCa patients from the PANCAN cohort. Table S1.2. Differentially expressed genes.
Additional file 5: Table S2.
Univariate Cox analysis results of differentially expressed genes.
Additional file 6: Table S3.
The detailed CNV percentages of HRD-related genes in PCa.
Additional file 7: Table S4.
Correlations among HRD-related genes.
Additional file 8: Table S5.
Information of HRD clusters and patient survival.
Additional file 9: Table S6.1.
ESTIMATE score. Table S6.2 Immune cell fractions based on ImmuCellAI algorithm. Table S6.3 Immune cell fractions based on ImmuCellAI algorithm. Table S6.4 Immune cell fractions based on CIBERSORT algorithm.
Additional file 10: Table S7.
The mRNA expression of IMs in PCa.
Additional file 11: Table S8.1.
Differentially expressed genes between HRD cluster 1 and cluster 2. Table S8.2. Differentially expressed genes between HRD cluster 1 and cluster 3. Table S8.3 Differentially expressed genes between HRD cluster 2 and cluster 3.
Additional file 12: Table S9.
Estimated response to immune checkpoint blockades.
Additional file 13: Table S10.
HRDscore and existing famous biological pathways based on proteomic data.
Additional file 14: Table S11.
The correlation between immune cell type and patient survival based on Cox analyses.
Additional file 15: Table S12.
Correlation between SLC26A4 expression and cancer functional states at single-cell level in multiple cancer types.
Additional file 16: Table S13.
The clinical characteristics, pathological diagnosis, HRD gene mutation, and SLC26A4 expressions of the samples in Xiangya Cohort.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Luo, C., Liu, Z., Gan, Y. et al. SLC26A4 correlates with homologous recombination deficiency and patient prognosis in prostate cancer. J Transl Med 20, 313 (2022). https://doi.org/10.1186/s12967-022-03513-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12967-022-03513-5