ADRB1 was identified as a potential biomarker for breast cancer by the co-analysis of tumor mutational burden and immune infiltration

Breast cancer (BRCA) has traditionally been considered as having poor immunogenicity and is characterized by relatively low tumor mutational burden (TMB). Improving immunogenicity may improve the response to clinical immunotherapy of BRCA. However, the relationship between TMB, immune infiltration, and prognosis in BRCA remains unclear. We aimed to explore their interrelations and potential biomarkers. In this study, based on somatic mutation data of BRCA from The Cancer Genome Atlas (TCGA), patients were categorized into high and low TMB groups utilizing the TMB values. CIBERSOFT algorithm indicated significant infiltration of activated partial immune cells in high TMB group. Besides, ADRB1 had been identified as a prognosis-related immune gene in the mutant genes by the combination of the ImmPort database and the univariate Cox analysis. ADRB1 mutation was associated with lower TMB and manifested a satisfactory clinical prognosis. Various database applications (Gene Set Enrichment Analysis, Tumor IMmune Estimation Resource, Connectivity Map, KnockTF) supported the selection of treatment strategies targeting ADRB1. In conclusion, TMB was not an independent prognostic factor for BRCA and high TMB was more likely to activate a partial immune response. ADRB1 was identified as a potential biomarker and may provide new insights for co-therapy of BRCA.


INTRODUCTION
Programmed death-1 (PD-1) and programmed death ligand-1 (PD-L1) are immune checkpoint inhibitors (ICIs) [1], which is the most studied type of immunotherapy for breast cancer (BRCA) according to relevant statistics [2]. TMB is a novel marker for evaluating the therapeutic effect of PD-1 antibodies, which has been confirmed in the treatment of colorectal cancer with defects in mismatch repair [3,4]. It is worth mentioning that TMB may be a promising tumor biomarker [5], defined as the total number of somatic gene coding errors, base substitutions, gene insertions or deletion errors per megabase (Mb). Higher TMB in tumors was reported to facilitate the formation of more new antigens and enhance tumor immunogenicity, which could improve clinical responses to cancer immunotherapy [6]. For example, patients with high TMB had better responses to ICIs and improved survival rates in melanoma, AGING urothelial carcinoma, non-small-cell carcinoma, and bladder cancer [7][8][9][10].
BRCA has traditionally been considered as having poor immunogenicity and is characterized by relatively low TMB [11]. However, the immune responses vary substantially between BRCA subtypes. Triple negative breast cancer (TNBC) and HER-2 (+) BRCA are generally more immunogenic than hormone-sensitive BRCA, as reflected in a higher proportion of tumor infiltrating lymphocytes [12]. In addition, luminal B subtypes can be more immunogenic than luminal A tumors among hormone-sensitive BRCA [13]. Allison and Vogelstein have reported a large number of new antigens in breast and bowel cancer tissues, and all cancers have the potential to accumulate new antigens that the immune system can recognize during tumorigenesis [14]. These findings suggest that TMB may play a predictive role in BRCA. Studies have demonstrated that the proliferation rate and the intrinsic subtype of BRCA were associated with TMB [15,16], whose role in tumor immunogenicity in BRCA is still unclear.
Meanwhile, certain issues remain with TMB. Previous clinical research found that comparing to patients in the low TMB group, not all patients in the high TMB group benefited from ICIs. Specifically, a subset of patients with mutations in the ERBB family (EGFR /ERBB2) and the deletion of specific 3p segments of the chromosome did not respond to ICIs [17]. Cristescu et al. published a study in Science [18], which has some implications for us: simultaneous detection of T cell activity levels and TMB may be a promising strategy. Indeed, positive correlations between mutations or new antigen loads and immune infiltration have been observed in various cancer types [19,20]. Therefore, we hypothesized that in combination with immune cell groups, TMB as a quantitative indicator of tumor antigenicity may influence the prognosis of BRCA.
In this study, we investigated the association of TMB with gene mutations, immune responses, and prognosis of BRCA in combination with tumor immune infiltration. Using the gene expression profiling data of BRCA from the TCGA database, different gene expressions between high and low TMB groups were compared, and aspects of the clinical characteristics, gene functions and pathways, as well as immune responses were further evaluated. We attempted to elucidate these relationships: different TMB and clinical outcomes, TMB and immune cell populations, immune cells affected by TMB and prognosis. The findings of this study may provide new biomarkers and potential therapy options for BRCA in the future.

Somatic mutation landscape in BRCA
Analysis of the 1,044 BRCA mutation samples from TCGA is shown in Figure 1A. Missense mutation was the primary variant classification and all mutations belonged to single nucleotide polymorphisms. C>T was the most common variation in BRCA with the highest number of variations per sample and the median of variation types. In addition, the frequencies of mutations in PIK3CA (29%) and TP53 (27%) were the highest in mutant genes, all of which were missense mutations ( Figure 1B). MUC17, HUWE1, SYNE1, TTN, MUC16, HMCN1 had equally higher co-mutation frequencies, while CDH1 and TP53 showed obvious mutuality of mutual exclusion ( Figure 1C).

Correlation analysis of TMB
The TMB in BRCA ranged from 0.02 to 112.8 per Mb with a median of 0.86 per Mb. With the median TMB value set as the threshold, a total of 986 samples was divided into the high (n=493) and low TMB (n=493) groups. We performed Kaplan-Meier analysis and determined that the 5-year survival rate of was 0.774 for the high TMB group and 0.870 for the low TMB group. Since the high TMB group predicted a better prognosis beyond 10 years, TMB may not be an independent prognostic factor for BRCA ( Figure 2A). In addition, among six clinical characteristics, only age and the N stage were significantly correlated with TMB; specifically, patients over 65 years old or with uninvolved regional lymph nodes had higher TMB ( Figure 2B). The differential expression of 454 mutant genes between groups is shown in Figure 2C.

Relationship between TMB and immune infiltration
The CIBERSOFT algorithm was used to assess the abundance of immune cells in the high and low TMB groups, and to explore the intrinsic relationship between TMB and the survival rate. Compared to those in the low TMB group ( Figure 3A), there were lower levels of B cells and T cells, and higher levels of macrophages in the high TMB group ( Figure 3B). Further comparisons indicated that naive/memory B cells, resting CD4 + memory T cells, follicular helper T cells, gamma delta T cells, resting dendritic cells, and resting mast cells were abundant in the low TMB group ( Figure 3C). For the high TMB group, there were significant infiltration of activated CD4 + memory T cells, M0/M1 macrophages, and activated dendritic cells. Furthermore, there were expressional correlations among the subsets of immune cells in transcriptome, a significant negative correlation between M0 AGING macrophages and resting CD4 + memory T cells, whereas activated CD8 + and CD4 + memory T cells were positively correlated ( Figure 3D). The Venn diagram showed that 44 immune genes in the differentially expressed genes (DEGs) were screened out ( Figure 3E) and ADRB1 was identified as a prognosis-related immune gene by the univariate Cox regression analysis (Table 1).

Functional enrichment analysis
We further examined the functional enrichment of DEGs especially ADRB1. Based on gene ontology categories ( Figure 4A), ADRB1 was significantly enriched in Gprotein coupled receptor binding, neurotransmitter receptor activity, neurotransmitter receptor activity involved in the regulation of postsynaptic membrane potential, and postsynaptic neurotransmitter receptor activity in molecular function, regulation of membrane potential, positive regulation of heart contraction, and heat generation in biological process, synaptic membrane and postsynaptic membrane in cellular component. Gene Set Enrichment Analysis (GSEA) performed with TCGA data indicated that the calcium signaling pathway, dilated cardiomyopathy, endocytosis, and neuroactive ligandreceptor interactions were significantly enriched AGING in samples with ADRB1 ( Figure 4B-4E). The findings also showed that the four pathways ADRB1 located in were all significantly active in the low-TMB group.

CNV of ADRB1, immune cells, and survival in BRCA
Generally, copy number variations (CNVs) refers to the increase or decrease in the copy number of a large segment in the genome whose length exceeds 1 kb. The results were presented in Figure 5A. In B cells and dendritic cells, high amplification of ADRB1 was significantly different compared to other CNVs (p<0.001). In addition, a high level of B cells suggested good prognosis of BRCA, and high expression of ADRB1 may prompt better survival ( Figure 5B).

Various small-molecule drugs of ADRB1 and the transcription factor HIF1A
Among the 43 targeted small-molecule drugs predicted for ADRB1, 17 were identified to be acting on BRCAassociated genes (Table 2), including vascular endothelial growth factor A (VEGFR1), dopamine receptor, prolactin, tumor necrosis factor, and polycyclic aromatic hydrocarbons. In the hypoxia AGING inducible factor A (HIF1A) knocking-down dataset in MCF-7 cells (Table 3), the upregulation of ADRB1 maybe not directly regulated by HIF1A, and it might be combined with the AFF1 factor to cause the knock-on effects, AFF1 was detected to bind to the super enhancer and typical enhancer region of the target gene (ADRB1), the regulatory mechanism within is unclear. However, we accidently found that the PIK3CA gene ( Figure 6) was involved in the VEGFR1-specific signaling pathway that HIF1A participates in, which has the highest proportion of gene mutations in this study. Its role needs to be further studied.

DISCUSSION
TMB was calculated based on the BRCA mutation data from TCGA, and the relationship between the survival curve and TMB showed that TMB may not be an independent prognostic factor for BRCA, which is consistent with previous studies on HER2 (-) metastatic BRCA [21]. We speculated that TMB combined with other prognostic factors may have a better predictive effect. To clarify the internal relationship between TMB and immunologic infiltration, we further showed that the low TMB group had abundant levels of B cells, follicular helper T cells, gamma delta T cells, and various resting immune cells. According to a recent study on triple-negative BRCA [22], the research team used a corresponding single anti CD8 + T cells in immune treatment to activate related anti-tumor immune mechanism, while ICIs activated follicular helper T cells that stimulated B cells to produce antibodies. However, the impact on tumor immune responses in inhibiting follicular helper T cells and B cells were more profound than inhibiting CD8 + T cells, which demonstrates that B cells and follicular helper T cells play key roles in tumor immune responses. Moreover, higher levels of gamma delta T cells have been shown to be correlated with better outcomes [23]. On the other hand, tumor immunogenicity was enhanced in the high TMB group, leading to significant infiltration of CD4 + memory T cells, M0/M1 macrophages, and dendritic cells as well as activated immune responses. The relative increase in TMB was also associated with aging and the N stage, consistent with previous literature that mutations of TP53 in lymph node-negative BRCA were higher than those in lymph node-positive BRCA, and mutations in microtubuleassociated proteins may help immune cells recognize tumors and inhibit lymph node metastasis [24].  Furthermore, correlations within the immune cells were determined by analyzing the immune matrix of the entire transcriptome. When investigating the clinical significance of infiltrated immune cells in BRCA, the higher proportion of M0 macrophages indicated a reduced disease-free survival, whereas the increased overall survival was associated with a relatively higher resting CD4+ memory T cells score [25], which corroborates the results of this study. In general, differences in immunogenicity may lead to differences in the activation of immune mechanisms, and the few types of immune cell activation in the high TMB group AGING might indicate that higher TMB suppresses the immune response. It is also worthwhile to note that PIK3CA and TP53 had prominent performance among mutation genes. Previous studies have revealed that mutant allele tumor heterogeneity is positively correlated with TP53 mutation rate, while CDH1 mutation is correlated with a low level of mutant allele tumor heterogeneity [26], confirming the correlation between TP53 and CDH1 observed in this study. The PIK3CA mutation was found in different subtypes such as ER (+), PR (+), HER2 (+), and TNBCs [27,28], but its role in VEGFR1-specific signaling pathway needs to be further  AGING Table 2. Seventeen small-molecule drugs predicted by ADRB1 as well as their effects on BRCA-associated genes.
Name Target  MOA  amiodarone  KCNH2, ADRB1, CACNA1H, CACNA2D2, CHRM3, CYP2C8,  KCNA7, SCN5A   Potassium channel blocker   carvedilol  ADRB1, ADRB2, ADRA1A, ADRA1B, ADRA1D, ADRA2A,  ADRA2B, ADRA2C, ADRB3, CYP2C19, CYP2E1, GJA1, HIF1A  explored. Patients with somatic mutations in TP53 and PIK3CA had reported poor survival [29], and whether the co-mutation of TP53 and PIK3CA can be potential biomarkers for different subtypes of BRCA warrants further investigation. ADRB1 was eventually identified as a prognosis-related immune gene for BRCA, whose functions were further explored. ADRB1, also called β-1 adrenergic receptor (AR), is a member of the G-protein coupled receptor family and an important target in various therapeutic applications. In cardiomyocytes, proteinkinase A activated by ADRB1 phosphorylates troponin I, the Ltype Ca 2+ channel and phospholamban, while increasing cardiac inotropy, chronotropy, and work [30]. In neuroinflammatory diseases, ADRB1 activation may have neuroprotective effects [31]. Furthermore, experiments have revealed that AR signaling can stimulate the transformation of epithelial cells to mesenchymal cells [32], and ADRB1 was observed to be overexpressed in BRCA tissues [33]. High expression levels of ADRB1 can predict better prognosis in this study, possibly because the overexpression of AR enhances the sensitivity of the tumor to β-blockers, although a previous report claimed that there was no correlation [34]. Future research should further clarify this issue. Pharmacoepidemiologic studies have shown that βblockers could reduce disease progression and mortality by inhibiting the metastasizing effect of AR signaling [35], but a retrospective analysis indicated that selective β-blockers alone or in combination were less effective than non-selective β-blockers in reducing cell proliferation in BRCA [33]. Interestingly, long-term deprivation of ovarian sex hormones can induce the upregulation of ADRB1 in the heart of rats [36], which suggested that the expression of ADRB1 is up-regulated when the sex hormone shows negative, and in our study, high expression of ADRB1 predicts better prognosis. Therefore, we speculate that HER2 (+) and triple-negative BRCA may be sensitive to β-blockers comparing to sex hormone types (such as the luminal subtype). That is, these two types of breast cancer may be easier to benefit from co-therapy. Prospective clinical trials of β-blockers on various subtypes of BRCA should be the focus of future research.
We further conducted a series of in-depth analyses on ADRB1. High amplification of ADRB1 in B cells and dendritic cells might indicate that ADRB1 mutation can facilitate two types of antigen presenting cells to efficiently mediate and maintain a normal immune response. The better prognosis in patients with high AGING levels of B cells also supports this observation. In addition to CNV, molecular research has demonstrated that the transcription factor HIF1A drives tumor growth and metastasis, and is associated with poor prognosis in BRCA [37]. The inhibition of HIF1A pathway activation combined with β-blockers may be a promising treatment strategy for BRCA patients. In addition, 17 small-molecule drugs targeting ADRB1 and other cancer-related genes obtained in this study also support this proposed treatment.
In conclusion, our study identified prognosis-related immune genes in BRCA mutations based on a co-analysis of TMB and immune infiltration, and explored the intrinsic correlation between TMB and immune infiltration. ADRB1 was identified as a potential biomarker for BRCA, which may provide new insights for co-therapy.

Data collection
Gene expression profiling for BRCA tissue samples (n=109, t=1109) and patients' clinical data (n=1097) were downloaded from the TCGA portal (https://portal.gdc.cancer.gov/) (Data Release 24.0 -May 07, 2020). In addition, tumor mutation data (n=1044) of BRCA including the names of the mutation genes, the mutation types, and the mutation locations were obtained from the "SomaticSniper variant aggregation and masking" platform.

Analysis of mutation genes
BRCA samples of TCGA were assessed using the R package "BiocManager" MAF files containing somatic variants and visualized with the maftools package. TMB was obtained by calculating the number of tumor mutations per Mb in each sample. The survival curve was plotted to present the survival rate in relation to TMB. A p-value < 0.05 was considered significant. The limma package was performed to assess the relationship between TMB and clinical characteristics including age, sex, and the stages T, N, and M (p<0.05).

TMB grouping and differential expression analysis
Normal samples in the tumor mutation data were deleted and the remaining tumor samples were cross-analyzed with the transcriptome samples. The median value of TMB was used as the threshold to divide samples into high and low TMB groups. The DEGs between the two groups were identified using the Wilcoxon rank test. The p value was adjusted by the false discovery rate (FDR) to improve the accuracy of the results, and the thresholds were set as FDR < 0.05 and logFC (fold change) > 1.0.

Co-analyses of TMB and immune infiltration
The deconvolution algorithm CIBERSORT [38] was used to evaluate the relative abundance of immune cells and the gene expression of tissue samples utilizing the gene expression characterization system of 22 different tumor-infiltrating lymphocyte subsets. The number of permutations was set to 1000, and a p-value <0.05 was regarded as successful. The immune cell matrix was obtained for each sample in the transcriptome data by the CIBERSORT R script v1.03. Similarly, the intrinsic differences in the abundance of immune cells between the high and low TMB groups were further explored and visualized by bar plots. The differences of immune cell infiltration between the high and low TMB groups were visualized by violin plots. Additionally, the list of immunologically relevant genes was downloaded from the ImmPort database (https://www.immport.org/) (Data Release 34, April 2020). Immune genes in DEGs were screened out using the Venn diagram. The univariate cox regression analysis was performed to identify prognosis-related immune genes (p<0.05).

Functional enrichment analysis
The gene ontology categories including biological processes, molecular functions, and cellular components were assessed for DEGs. Moreover, to determine whether ADRB1-related pathways were statistically and consistently different between the high and low TMB groups, we performed pathway enrichment analysis using the GSEA software (version 4.0.3) with FDR<0.05 considered statistically significant.

CNV and immune cells
Tumor IMmune Estimation Resource (TIMER v2.0, https://cistrome.shinyapps.io/timer/), a web server for comprehensive analysis of tumor-infiltrating immune cells, was used to estimate the abundances of six immune infiltrates (B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells) [39]. Changes in CNV were observed in prognosis-related immune genes, and the correlations between CNV and immune cell abundance, and between immune cells and survival were further assessed.