The novel transcriptomic signature of angiogenesis predicts clinical outcome, tumor microenvironment and treatment response for prostate adenocarcinoma

Angiogenesis plays the critical roles in promoting tumor progression, aggressiveness, and metastasis. Although few studies have revealed some angiogenesis-related genes (ARGs) could serve as prognosis-related biomarkers for the prostate cancer (PCa), the integrated role of ARGs has not been systematically studied. The RNA-sequencing data and clinical information of prostate adenocarcinoma (PRAD) were downloaded from The Cancer Genome Atlas (TCGA) as discovery dataset. Twenty-three ARGs in total were identified to be correlated with prognosis of PRAD by the univariate Cox regression analysis, and a 19-ARG signature was further developed with significant correlation with the disease-free survival (DFS) of PRAD by the least absolute shrinkage and selection operator (LASSO) Cox regression with tenfold cross-validation. The signature stratified PRAD patients into high- and low-ARGs signature score groups, and those with high ARGs signature score were associated with significantly poorer outcomes (median DFS: 62.71 months vs unreached, p < 0.0001). The predicting ability of ARGs signature was subsequently validated in two independent cohorts of GSE40272 & PRAD_MSKCC. Notably, the 19-ARG signature outperformed the typical clinical features or each involved ARG in predicting the DFS of PRAD. Furthermore, a prognostic nomogram was constructed with three independent prognostic factors, including the ARGs signature, T stage and Gleason score. The predicted results from the nomogram (C-index = 0.799, 95%CI = 0.744–0.854) matched well with the observed outcomes, which was verified by the calibration curves. The values of area under receiver operating characteristic curve (AUC) for DFS at 1-, 3-, 5-year for the nomogram were 0.82, 0.83, and 0.83, respectively, indicating the performance of nomogram model is of reasonably high accuracy and robustness. Moreover, functional enrichment analysis demonstrated the potential targets of E2F targets, G2M checkpoint pathways, and cell cycle pathways to suppress the PRAD progression. Of note, the high-risk PRAD patients were more sensitive to immune therapies, but Treg might hinder benefits from immunotherapies. Additionally, this established tool also could predict response to neoadjuvant androgen deprivation therapy (ADT) and some chemotherapy drugs, such as cisplatin, paclitaxel, and docetaxel, etc. The novel ARGs signature, with prognostic significance, can further promote the application of targeted therapies in different stratifications of PCa patients.

Introduction
Prostate adenocarcinoma (PRAD), accounting for 95% of all prostate cancers (PCa), is a complicated disease threatening the health of men worldwide (Rebello et al. 2021). As reported, PCa has become the second most common male malignancy, with a incidence of 29.3/10 5 and morality of 7.6/10 5 worldwide (Feng et al. 2019). Similarly, the incidence of PCa in China has been rising to over 10/10 5 and now ranks the seventh most common male cancer and tenth leading cause of cancer death (Feng et al. 2019;Zhu and Ye 2020). With the development of screening techniques and modern medicine, prostate cancer has a favorable 5-year survival rate among the common types of cancers. However, patients with PCa in China are often diagnosed at locally advanced stage or metastasis stage, resulting in a high mortality-to-incidence ratio, nearly 50% (Zhu and Ye 2020). Though the 5 year-survival of patients with localized PCa is over 95%, about half patients who received radical prostatectomy will experience biochemical recurrences (Dell'Oglio et al. 1990). Thus, reliable prognostic biomarkers should be proposed to improve the clinical management of PCa.
Some risk stratification tools based on gene expression have been developed, such as Decipher (Spratt et al. 2017), Oncotype Dx Genomic Prostate Score (Eure et al. 2017) and Prolaris (Shore et al. 2016). These tools have been proved to be associated with disease recurrence after surgery and/or adjuvant therapy. However, some of these models did not make contributions to clinical practices (Cucchiara et al. 2018). To our knowledge, a number of molecular biomarkers have been discovered as well, which are capable of identifying aggressive subtypes of indolent PCa with higher specificity (Bertoli et al. 2016;Lam et al. 2020;Faisal and Lotan 2020;Zhuo et al. 2018). For instance, miRNA prognostic biomarkers, such as let-7a, miR-141, miR-375 and miR-182, are correlated with lymph node metastasis and clinical staging in PCa (Bertoli et al. 2016). Moreover, epigenetic biomarkers, such as GSTP1 and APC, were found to be correlated with the relapse of PCa (Lam et al. 2020). Next generation sequencing (NGS) technologies enabled us to have a deep understanding of the genomic information of PCa patients, and some genomic alteration signatures have been applied in diagnosis, prognosis, and clinical therapeutics (Faisal and Lotan 2020). Furthermore, increasing transcriptomic data of PCa patients provided deep insights into molecular mechanisms of tumorigenesis and progression, and also promoted the development of molecular biomarkers. For example, the upregulation of SRPK2 was shown to affect cancer cell proliferation, migration, cell cycle, and was found to be correlated with a higher Gleason score, advanced clinical stage, tumor metastasis and poor prognosis (Zhuo et al. 2018). However, the prognosis and management of PCa are still complicated due to the molecular, cellular, and clinical heterogeneity in PCa . Therefore, more rigorous and reliable prognosis prediction models are urgently needed to stratify the PCa patients and improve the clinical outcomes.
Angiogenesis is well-known to play critical roles in the formation of nascent blood vessels, which are mainly involved in the vascular development and repair of damaged blood vessels (Rajabi and Mousa 2017). Whereas, pathological blood vessels, unlike normal blood vessels, are immature and can influence the tumor microenvironment (TME) and promote proliferation and migration of cancer cells and support tumor progression, invasion and distant metastasis (Katayama et al. 2019;Carmeliet and Jain 2000). Thus, pathological angiogenesis is considered as one of the tumor hallmarks, correlated with poor prognosis in many cancer types (Ramjiawan et al. 2017). Nowadays, plenty of ARGs had been identified to be involved in the development of tumors (Pavlakovic et al. 2001). The overexpression of some certain ARGs, such as vascular endothelial growth factor A (VEGFA), neuropilin-1 (NRP1) and PlexinA2 (PLXNA2), were found in PCa and their overexpression correlated with tumor metastasis and poorer prognosis (Melegh and Oltean 2019;Yin et al. 2021). In addition, some studies revealed that PCa progression and metastasis were promoted by other angiogenesis-enhancing factors (Zhao et al. 2021;Strohmeyer et al. 2000). These studies focused on the prognostic value of a single ARG or angiogenesis-related processes, but the roles of ARGs related to the whole angiogenesis activity in predicting the prognosis of PCa are not clearly clarified.
In our study, we investigated the RNA-seq and clinical features of PRAD samples downloaded from the TCGA, and 23 angiogenesis-related genes correlated with worse DFS were identified from the Hallmark Angiogenesis gene set (36 ARGs in total). LASSO Cox regression with tenfold cross-validation further confirmed a novel 19-ARG signature, which was a predominant prognostic factor for the DFS of PRAD. Furthermore, a nomogram based on the ARGs signature was generated to predict the prognosis of PRAD patients. Finally, the role of the signature in predicting treatment response was investigated to explore the potential targeted therapies for PRAD patients.

Data resources and selection of angiogenesis-related genes
The original genomic alternations, mRNA expression and clinical features of 501 PRAD samples (4/501 PRAD samples lacking one of these information were excluded, only 497/501 PRAD samples were included in the following analyses) in PRAD cohort (TCGA, Firehose Legacy) were downloaded via the cbioportal (https:// www. cbiop ortal. org/ study/ summa ry? id= prad_ tcga) as the discovery dataset (Table 1). GSE40272 dataset, which contains 89 PRAD samples, was downloaded from the GEO (https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE40 272) as the first validation dataset (Additional file 7: Table S1), meanwhile, the PRAD_MSKCC cohort, including a total of 156 PRAD samples with mRNA expression data and there were 140/156 PRAD samples with DFS information), retrieved via the cbioportal was regarded as another validation dataset (Additional file 8: Table S2).

Differential expression analysis between PRAD patients and normal people, and expression and mutation profiles of ARGs in PRAD
The mRNA expression data in the cohort (containing 500 tumor samples & 51 normal samples), namely GDC TCGA Prostate Cancer (PRAD) that was collected in the website of UCSC xena (https:// xenab rowser. net/ datap ages/) and could be directly downloaded from the following: https:// gdc-hub. s3. us-east-1. amazo naws. com/ downl oad/ TCGA-PRAD. htseq_ counts. tsv. gz, were retrieved to explore the differentially expressed ARGs using the "DeSeq2" package, while |log 2 (Fold Change) |> 1 and p < 0.05 were considered significantly differentially expressed.
Angiogenesis-related genes data were extracted using cBioPortal to show the relationship between genetic alterations (missense, splice, truncating, amplification, and deep deletion alterations, etc.) and patients' survival or between expression alterations and patients' survival among these patients from the TCGA-PRAD cohort. The relatively high or low expression was defined by the custom parameter (z-score) of cBioPortal, while "z-score > 2" represented the relatively high mRNA expression and "z-score < − 2" represented the relatively low mRNA expression.

Construction of the angiogenesis-related prognostic signature and nomogram
The ARGs related to DFS (P < 0.05) were identified via a univariate Cox regression analysis using optimized algorithm. Following least absolute shrinkage and selection operator (LASSO) Cox regression with tenfold crossvalidation was performed using the "glmnet" package in R studio (Sauerbrei et al. 2007;Wang et al. 2019). Multivariate Cox regression analysis was then conducted for the construction of this ARG signature, and the ARGs signature score formula was established. ARGs signature score = Expression level of gene 1 × ɑ 1 + Expression level of gene 2 × ɑ 2 + Expression level of gene 3 × ɑ 3 + … + Expression level of gene n × ɑ n , where ɑ n denotes the coefficient for the corresponding gene in this model. The cutoff value was drawn from median ARGs signature score to sort the patient samples into high-and low-ARGs signature score groups. The survival curves were drawn using the "survival" package in R studio. The analysis of receiver operating characteristic curve (ROC) was demonstrated by the "timeROC" package in R studio, and the values of area under ROC curve (AUC) for DFS were used to evaluate the ability of ARGs signature or clinical feature in prognosis prediction. Underlying ARGs signature, clinical association analysis was subsequently conducted. Moreover, ARGs signature and clinical features, such as diagnosis age, PSA (prostate specific antigen) level, Gleason score and T, N, M stage, were further enrolled in the multivariate Cox regression analysis to determine the independent prognostic factors. Based on the independent prognostic factors, the nomogram was constructed via the Cox regression model by using the "rms" package in R studio (Chun et al. 2006). Also, the calibration plots were drawn by the "rms" package in R studio to estimate the accuracy of predicted DFS by the nomogram.

Genomic alteration and functional enrichment analysis
Then, the genomic alteration profiling between high-and low-ARGs signature score groups was compared. Initially, the comparison analysis of mutation count between these two groups was conducted by using the "violinplot" package in R studio. Subsequently, the oncoplots were drawn by using the "maftools" package in R studio to exhibit genomic alteration landscapes of high-and low-ARGs signature score groups. The oncogenic genes involved in the DDR, PI3K, and Wnt signaling pathways were highly correlated with the development and progression of PCa (Armenia et al. 2018), thus, genomic alterations in these three pathways were further investigated in the present study.
To explore the biological processes involved in the high-and low-ARGs signature score groups, we used the "clusterProfiler" package (Yu et al. 2012) in R studio to perform Hallmark gene set enrichment analysis (http:// www. gsea-msigdb. org/ gsea/ msigdb/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https:// www. kegg. jp/) pathway enrichment analysis, and all of these analyses were visualized by "ggplot2" and "enrichplot" in R studio. The threshold was set at an adjusted p-value < 0.05.

Evaluation of TME between high-and low-ARGs signature score groups
Transcriptomic data from the TCGA-PRAD cohort were used to compare the proportions of infiltration among 22 types of tumor-infiltrating immune cell by a deconvolution algorithm (CIBERSORT: Cell type Identification By Estimating Relative Subsets Of known RNA Transcripts) between the high-and low-ARGs signature score groups. Furthermore, immune and stromal scores illustrating the infiltration status of TME cells were calculated using ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumours using Expression data) (Yoshihara et al. 2013a). In addition, the expression levels of immune checkpoint genes, including PD-L1, PD-1, CTLA4, and TIM3, were compared between the high-and low-ARGs signature score groups. Moreover, the relationship between the expression levels of immune checkpoint genes and ARGs signature score (or signature-related ARGs) was further investigated by the correlation analysis.

The analyses of endocrine therapy and chemotherapy responses
The GSE102124 dataset containing 19 primary prostate cancer samples with the treatment of neoadjuvant ADT and 3 primary prostate cancer control samples (Sowalsky et al. 2018) and the GSE74685 dataset containing 117 castration-resistant prostate cancer (CRPC) samples with the treatment of luteinizing hormone-releasing hormone (LHRH) agonist, 19 CRPC samples receiving the orchiectomy treatment, and 13 CRPC control samples (Kumar et al. 2016) were retrieved from the GEO for the analysis of the ADT treatment responses.
Meanwhile, data were retrieved from the Genomics of Drug Sensitivity in Cancer (GDSC, https:// www. cance rrxge ne. org/) to determine the difference in the response to chemotherapy between high or low-ARGs signature score groups (Yang et al. 2013a;Iorio et al. 2016;Garnett et al. 2012). In the present study, the half-maximal inhibitory concentration (IC50) was utilized to evaluate the drug response. Meanwhile, the IC50 of each sample was estimated between the high-and low-ARGs signature score groups. Furthermore, mRNA expression data of docetaxel sensitive and resistant variants of PC3 PCa cell line from GSE140440 (Schnepp et al. 2020) was further employed to testify the predictive role of ARGs in chemotherapy response.

Immunohistochemistry (IHC) analysis
The protein expression of 19 ARGs involved in the signature, among TCGA-PRAD samples, was analyzed by IHC using the available scanned tumor staining from Human Protein Atlas (HPA) database (https:// www. prote inatl as. org/). The information of IHC staining was determined and manually adjusted by experts from the HPA database, and ARGs IHC staining was defined and exhibited as high, medium, low staining or not detected. A total of 9-12 PRAD samples were enrolled in the analysis of ARGs IHC staining. In addition, there existed the IHC staining of 16 ARGs protein expression, except CXCL6, OLR1, and LPL, retrieved from the HPA database.

Statistical analyses
Chi-square, Fisher test, and Wilcoxon rank test statistical analyses were performed using R software (v. 3.4.3, https:// rstud io. com/). A (adjust) p-value < 0.05 was considered statistically significant. A log-rank test was used to estimate the survival curves between the high-and low-ARGs signature score groups. A schematic representation for the design of the present study was shown in Fig. 1.

Genomic and expression alterations of ARGs in PRAD patients
Initially, we analyzed the genomic alteration profiling of ARGs gene set in PRAD patient samples from the TCGA database, and it was showed that 43% (214/497) of PRAD patient samples had at least one genomic alteration, including missense, splice, truncating, amplification, and/or deep deletion alterations in ARGs ( Fig. 2A). Among those altered genes, STC1 was the most prevalent (17%), while others were altered in 0.4%-16% of patient samples. However, there was no significant difference (p = 0.791) in the DFS between groups with or without genetic mutations in ARGs (Fig. 2B).
Then, we explored the expression data of the ARGs in PRAD patient samples from TCGA. Almost 68% (336/497) of PRAD patient samples had altered expression of ARGs (Fig. 2C). Notably, patient samples with the altered expression of ARGs had significantly poorer DFS than that of patient samples without changes (median DFS: 77.33 months vs unreached, p < 0.0001, Fig. 2D). While the comparison of expression data between tumor and normal samples showed that 8 out of 36 ARGs were significantly differentially expressed, including CCND2, CXCL6, KCNJ8, LPL, SERPINA5, SLCO2A1, VTN, and APOH (Additional file 9: Table S3). Interestingly, only APOH was downregulated but other 7 ARGs were upregulated in the tumor samples, compared to the normal samples.

Construction and evaluation of a novel prognostic signature for PRAD patients based on the ARGs expression profiling
The univariate Cox regression analysis was used to identify the specific association between each individual ARG gene from the ARGs gene set and the DFS of PRAD patient samples from the TCGA cohort, and we found a total of 23 ARGs were significantly correlated with the DFS of PRAD patients (p < 0.05, Fig. 3A). Among these identified ARGs, the overexpression of 17 ARGs were associated with poorer DFS of PRAD patients; on the By the LASSO Cox regression analysis, 19 ARGs were further determined to construct a novel prognostic signature for PRAD patients (Fig. 3B). The coefficients of these 19 ARGs were calculated via the multivariate Cox regression analysis (Table 2), and a median ARGs signature score was defined as a cutoff value which was used to divide PRAD patient samples into the high-and low-ARGs signature score groups. Notably, patients in the high ARGs signature score group had the significantly worse prognosis (median DFS: 62.71 months vs unreached, p < 0.0001, Fig. 3C), compared to those in the low ARGs signature score group. Moreover, the AUC values for DFS at 1, 3, and 5 years were 0.77, 0.75 and 0.78, respectively (Fig. 3D).
In addition, the constructed signature was further verified in two independent validation cohorts: GSE40272 dataset & PRAD_MSKCC cohort. The Kaplan-Meier curve for PRAD patients in GSE40272 showed that patients with high ARGs signature score had the worse survival (median DFS: 60.30 months vs unreached, p < 0.05, Fig. 3E), in contrast to those with low ARGs signature score. Besides, the AUC values for DFS at 1, 2, and 3 years were 0.51, 0.61 and 0.65, respectively (Fig. 3F). In addition, the survival analysis of PRAD patients from PRAD_MSKCC cohort also demonstrated that patients in high ARGs signature score group had the poorer outcomes (median DFS: 55.39 months vs unreached, p < 0.0001, Fig. 3G), and the AUC values for DFS at 1, 3, and 5 years were 0.76, 0.68 and 0.63, respectively (Fig. 3H).

Association analysis between the ARGs signature and clinical features
Next, we investigated the clinical features related to ARGs signature in the TCGA-PRAD cohort. There was no significant difference in the age at diagnosis between two groups (median age at diagnosis: 61.02 vs 60.75 years old, p > 0.05, Fig. 4A Gleason score is an extensively used marker for histological grading in PCa, thus, we investigated the difference in the distribution of Gleason score between patients with different ARGs signature score. It was found that more PRAD patients in the high ARGs signature score group had a high Gleason score (16.94% vs 8.87% in Gleason 8, 42.34% vs 12.91% in Gleason 9, 1.61% vs 0% in Gleason 10), compared to those in low ARGs signature score group (p < 0.01, Fig. 4C). Further comprehensive analyses of clinical stages revealed that patients with high ARGs signature score were significantly associated with advanced tumor stages and lymph node invasion. Among these patients with high ARGs signature score, 14.92% and 25.41% of them had disease of T3 and N1, respectively (p < 0.01, Fig. 4D, E). In contrast, only 6.4% of patients with low ARGs signature score were at T3 stage and/or N1 stage. No significant difference in distant metastasis between two groups was observed due to the limited patients in M1 stage involved in the TCGA dataset (Fig. 4F). Of note, it was observed that the expression of common tumor angiogenesis markers of CD34 & Fig. 3 Construction and evaluation of a prognostic 19-ARG signature. A Hazard ratio in high and low expressions of ARGs in PRAD cohort, and corresponding 95% confidence intervals were calculated by the univariate regression analysis. B A 19-ARG signature was confirmed by the LASSO regression analysis with tenfold cross-validation. C Kaplan-Meier DFS curves for patients from the TCGA cohort in the high-or low-ARGs signature score group. D The ROC curves for DFS at 1, 3, and 5 years were analyzed in the TCGA cohort. E Kaplan-Meier DFS curves for patient samples assigned to high-and low-ARGs signature score in the validation group, GSE40272 from the GEO database. F The ROC curves for DFS at 1, 2, and 3 years were analyzed in GSE40272. G Kaplan-Meier DFS curves for patient samples assigned to high-and low-ARGs signature score in another validation group of the PRAD_MSKCC cohort. H The ROC curves for DFS at 1, 3, and 5 years were analyzed in the PRAD_MSKCC cohort CD105 was higher in high ARGs signature score group, moreover, ARGs signature score was positively with the expression of CD34 & CD105 (p < 0.001, Additional file 1: Fig. S1).

Construction of a nomogram
Compared to the representative clinical features, the ARGs signature score, Gleason score, and T stage were shown to be the predominant prognostic factors for PRAD patients by the multivariate Cox regression analysis and the ROC curve analysis (Table 3). The AUC values for DFS at 1, 3, and 5 years showed that the specificity and sensitivity of the 19-ARG signature were higher than those of clinical features ( Fig. 5A-C). Besides, the ARGs signature outperformed a single ARG involved in this signature to predict the 1, 3, and 5 years DFS of PRAD patients (Additional file 2: Fig. S2). Subsequently, a nomogram based on the ARGs signature score was established, together with the clinical parameters: T stage and "Gleason primary score + secondary score" (Fig. 5D). While the concordance index (C-index) of the nomogram was 0.799 (95%CI = 0.744-0.854). As indicated in the nomogram, every PRAD patient would have a total score and the patients with higher total score could have worse prognosis. In addition, the calibration plots revealed that the predicted outcomes were basically in agreement with the observed DFS in the TCGA-PRAD cohort (Fig. 5E). Meanwhile, the AUC values of 1-, 3-, 5-year DFS for the nomogram were 0.82, 0.83, and 0.83, respectively, showing that the nomogram had better and more stable predicting ability (Fig. 5F). Of note, the AUC values of 1-, 2-, 3-year DFS for the nomogram in GSE40272 were 0.59, 0.72, and 0.76, respectively, furthermore, the AUC values of 1-, 3-, 5-year DFS for the nomogram in the PRAD_ MSKCC cohort also increased to 0.79, 0.74, and 0.68, by comparison with the performance of ARGs signature in prognosis prediction (Additional file 3: Fig. S3).

Fig. 4
The analyses of clinical features between high-and low-ARGs signature score groups in PRAD cohort. A Boxplots representing age of diagnosis for high-and low-ARGs signature score groups (p = 0.092). B Boxplots of the PSA level between high-and low ARGs signature score groups (p = 0.005). C Percentage-staked bar plot representing patient distribution with Gleason primary + secondary score between high-and low-ARGs signature score groups (p < 0.001). D Percentage-staked bar plot representing tumor stage (T1/T2/T3/T4/NA) between high-and low-ARGs signature score groups (p < 0.001). E Percentage-staked bar plot representing lymph node status (N0/N1//NA) (p < 0.001) between high-and low-ARGs signature score groups (p < 0.001). F Percentage-staked bar plot representing tumor metastasis (M0/M1//NA) between high-and low-ARGs signature score groups (p = 0.12). (NA not applicable)

Enrichment of genomic alterations
As is shown in Fig. 6A, more genomic mutations (mutation count: 35 [4, 6524] vs 34 [2, 104], p < 0.01) were presented in patients from high ARGs signature score group. Then we deeply investigated the mutation profiles of the top ten most frequently altered genes in high-and low-ARGs signature score groups. Obviously, the most prevalently altered genes in high ARGs signature score group were significantly different from those in the low ARGs signature score group (p < 0.05, Fig. 6B, C, Additional file 10: Table S4). Consequently, the systematic investigation was further conducted to explore which events were correlated with the high-risk PRAD patients. Genomic alterations happened more frequently in the high ARGs signature score group, and there was a statistically significant difference of genomic alterations in a total of 47 altered genes between the high-and low-ARGs signature score groups (Additional file 11: Table S5). By statistical analysis, 44 of 47 altered genes, including TP53, PIK3CA and ALMS1, etc., were significantly enriched in the high ARGs signature score group (p < 0.05, Fig. 6D, Additional file 11: Table S5), reversely, the low ARGs signature score group had a significant enrichment of only 3 altered genes, including FLG2, KHDC4 (KIAA0907), and KRTAP4-6 (p < 0.05, Fig. 6D, Additional file 11: Table S5). For the level of signaling pathway, we found that more frequently altered pathways of DNA Damage Response (DDR), PI3K and Wnt were enriched in the high ARGs signature score group, meanwhile, genes PIK3CA in PI3K signaling pathway and POLE in DDR signaling pathway were more frequently altered in the high ARGs signature score group (p < 0.05, Fig. 6E), and the altered genes were excluded as of the frequency < 1%. In addition, the DDR, PI3K and Wnt pathways were altered in nearly 10.62%, 6.01%, and 6.01% of the PRAD patient samples from the TCGA-PRAD cohort (Fig. 6F-H).

Functional enrichment analysis
Initially, Hallmark gene set (Additional file 12: Table S6) and KEGG (Additional file 13: Table S7) enrichment analysis underlying the expression of signature-related ARGs revealed that the aberrant regulation of ARGs expression in PRAD samples was highly correlated with biological activities of E2F targets, G2M checkpoint, Myc targets V1, oxidative phosphorylation, cell cycle, angiogenesis, Androgen or estrogene response, extracellular matrix, and immune-related signaling pathways, etc. Based on the integrative ARGs signature, the Hallmark gene set enrichment analysis showed that E2F targets, G2M checkpoint, Myc targets V1, mitotic spindle, allograft rejection, DNA repair, Myc targets V2, interferon γ response, angiogenesis, epithelial-mesenchymal transition, inflammatory response, and interferon α response were highly enriched in the high ARGs signature score group (adjusted p-value < 0.01, Fig. 7A; Additional file 14: Table S8). And KEGG pathway enrichment analysis indicated that pathways related to cell cycle, DNA replication, and spliceosome etc., were abundantly enriched in the high ARGs signature score group (adjusted p-value < 0.01, Fig. 7B; Additional file 15: Table S9).

The estimate of tumor microenvironment and ARGs signature-related immune checkpoint gene expression analysis
The stromal score, immune score and ESTIMATE score were significantly higher in high ARGs signature score group (p < 0.05, Fig. 8A). Meanwhile, infiltrating Fig. 6 Genomic alteration enrichment analysis. A Vioplots of gene mutation counts in high-and low-ARGs signature score groups. B The alteration profile of the top 10 most frequently altered genes in high-ARGs signature score group. C The alteration profile of the top 10 most frequently altered genes in low-ARGs signature score group. D The distribution of the altered genes between the high-and low-ARGs signature score groups. E The distribution of the altered signaling pathways between the high-and low-ARGs signature score groups. F The alteration profile of the frequently altered genes of DDR, PI3K, and Wnt pathways between the high-and low-ARGs signature score groups immune cells were analyzed to describe the immune profiles in the high-and low-ARGs signature score groups (Fig. 8B, C). Naive B cells, plasma B cells and CD4 + memory resting T cells are significantly enriched in the low ARGs signature score group. On the contrary, the high ARGs signature score group had a significantly higher abundance of memory B cells, Treg cells, activated NK cells, and Macrophage M2 (p < 0.05, Fig. 8C). Additionally, the expression of immune checkpoint genes was investigated, and which demonstrated that the expression of PD-1, CTLA4, and TIM3 was significantly higher in the high ARGs signature score group (p < 0.05, Fig. 8D), however, the expression of PD-L1 was nearly equivalent between two groups (p = 0.73, Fig. 8D) and was not associated with ARGs signature score (p = 0.96, Fig. 8E). Of note, it was found that the ARGs signature score had the extremely positive correlation with the expression of PD-1 (p < 0.01), CTLA4 (p < 0.0001), and TIM3 (p < 0.0001) (Fig. 8E). Correlation analysis between each signature-related ARG and the selected immune checkpoint genes further demonstrated that OLR1, COL3A1, S100A4, POSTN, LPL, CXCL6, FGFR1, JAG2, and SERPINA5 were all positively correlated with PD-1, CTLA4, and TIM3 (p < 0.05, Fig. 8F), while APP was negatively correlated with CTLA4 and TIM3 but positively correlated with PD-1 (p < 0.05, Fig. 8F).

The evaluation of endocrine therapy and chemotherapy responses
The endocrine therapy has been one of the standard treatment strategies to improve the clinical outcomes of PCa patients. Remarkably, it was found that the ARGs signature score of the primary prostate cancer patients from the GSE102124 dataset was significantly downregulated (p = 0.0013, Fig. 9) after receiving the conjunction treatment of leuprolide, abiraterone acetate, and prednisone, which was a neoadjuvant ADT. However, it was observed that LHRH agonist and orchiectomy treatments had no influence on the regulation of the ARGs signature score of CRPC patients from the GSE74685 dataset (p > 0.05, Additional file 4: Fig. S4). The data from the GDSC database were utilized to predict the possible responses to some traditional chemotherapy drugs in the cancer. It could be observed that there was a statistically significant difference of IC50 among six chemotherapy drugs between the high-and low-ARGs signature score groups. Moreover, PRAD patients with the high ARGs signature score seemed to be more sensitive to five identified chemotherapy drugs (cisplatin, paclitaxel, docetaxel, vinblastine and doxorubicin, p < 0.05, Fig. 10). Nevertheless, PRAD patients with low ARGs signature score seemed to be more sensitive to vinorelbine (Fig. 10). Of note, it was further found in PC3 PCa cell line from the GSE140440 dataset that high ARGs signature score group was significantly correlated with the sensitivity of docetaxel (p < 0.0001, Additional file 5: Fig. S5). In addition, correlation analysis revealed that in the high ARGs signature score group not all signature-related ARGs were associated with sensitivity to the five identified chemotherapy drugs, while only LRPAP1 was significantly correlated with sensitivity to all these five chemotherapy drugs and positively correlated with cisplatin and doxorubicin but negatively correlated with Fig. 8 Tumor microenvironment (TME) in high-and low-ARGs signature score groups in PRAD cohort. A Comparison of the immune score, stromal score, and ESTIMATE score between the high-and low-ARGs signature score groups. B The immune profiles of high-and low-ARGs signature score groups. C Comparison of infiltrating immune cells between high-and low-ARGs signature score groups. D Differential expression of immune checkpoint genes (PD-L1, PD-1, CTLA4, TIM3) between high-and low-ARGs signature score groups. E Correlation analysis between the expression level of immune checkpoint genes (PD-L1, PD-1, CTLA4, TIM3) and the ARGs signature score. F Correlation analysis between the signature-related ARG and immune checkpoint genes. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, cor.sp represented "Spearman Correlation Coefficient") paclitaxel, docetaxel, and vinblastine (p < 0.05, Additional file 6: Fig. S6). Whereas, it could be obviously observed that most ARGs in signature were significantly correlated with sensitivity to the six identified chemotherapy drugs (p < 0.05, Additional file 6: Fig. S6).

Discussions
Although emerging strategy of anti-angiogenesis has been validated as an effective treatment for multiple solid tumors including renal carcinoma, lung cancer and stomach cancer (Schmidt and Carmeliet 2011;Sarkar et al. 2015), unfortunately, it did not exhibit promising clinical efficacy as expected in the PCa. On the other side, the integrated role of ARGs in predicting the prognosis of PCa has been neglected and scarcely studied. The present study offered a deep insight of the integrated function of angiogenesis in the prognosis of PCa. The angiogenesis-related gene set, namely the Hall-mark_Angiogenesis gene set from Molecular Signatures Database, was selected to establish the signature. The hallmark gene set has reduced noise and redundancy, as well as controlled false discoveries by identifying gene set overlaps and analyzing coordinate expression (Liberzon et al. 2015). We investigated the specific association between these ARGs and DFS of PRAD and found that patients with altered expression of ARG(s) were significantly correlated with poorer DFS, compared to those without any alterations. However, no significant difference was observed in the DFS of groups with or without ARG genetic mutations. Then, we identified a total of 23 ARGs that significantly correlated with PCa patients' survival. In concordance with previous studies which investigate the ARG solely in prostate cancer, genes including NRP1 (Tse et al. 2017), JAG2 (Kwon et al. 2016), COL5A2 (Ren et al. 2021), POSTN (Cattrini et al. 2018), FSTL1 (Zhao et al. 2019), PF4 (Baselga et al. 2008), JAG1 (Terada et al. 2014), COL3A1 (Angel et al. 2020), VCAN (Asano et al. 2017), VEGFA (Zhan et al. 2013), SPP1 (Pang et al. 2019), VTN (Niu et al. 2016), and S100A4 (Ganaie et al. 2020) were identified as unfavorable prognosis-related biomarkers in prostate cancer. Meanwhile, the prognosis roles of APP (Takayama et al. 2009) and FGFR1 (Yang et al. 2013b) were controversial with previous studies, which needed further validation and study. Notably, genes including PRG2, VAV2, LPL, OLR1, SERPINA5, CXCL6, LRPA1, and KCNJ8 have not been identified to be correlated with the prognosis of PCa before, meanwhile, it was also found that CXCL6, KCNJ8, LPL, and SERPINA5 were upregulated in the tumor samples by comparison with normal samples. Thus, it is worthy to investigate the specific biologial function of these genes in the carcinogenesis and development of PCa in the further study.
Based on the expression profiling of the ARGs gene set, the signature (a 19-ARG signature) was developed, which had the ability to predict the prognosis of PRAD patients. Meanwhile, the result of clinical association analysis demonstrated that the higher ARGs signature score indicated the advanced clinical status as well as the higher expression of clinically used tumor angiogenesis markers (CD34 & CD105), which could further help clinicians promote the PRAD management. In addition, a prognostic nomogram based on the combination of ARGs signature score and two clinical features (Gleason score & T stage) was established with the high AUC values of 0.82, 0.83, and 0.83 for the DFS at 1, 3, and 5 years, respectively, demonstrating more reliable and stable capacity of nomogram to predict prognosis of PRAD patients and which was also verified in other two independent PRAD patient cohorts, by comparison with the performance of the ARGs signature solely in prognosis prediction. Thus, the ARGs signature-based prognostic nomogram was highly recommended for future clinical applications. Previous studies had demonstrated the influence of angiogenesis on tumor growth and progression (Katayama et al. 2019;Carmeliet and Jain 2000;Ramjiawan et al. 2017;Pavlakovic et al. 2001;Schmidt and Carmeliet 2011;Baselga et al. 2008;Goos et al. 2016;Sbiera et al. 2017;Wei et al. 2020;Dunne et al. 2006;Townsend et al. 2019;Javaherian and Lee 2011;Lin et al. 2003;Maurer et al. 2006;Yang et al. 2020). Moreover, inhibiting angiogenesis has emerged as an effective therapeutic strategy for most solid tumors. But past researches of ARGs or angiogenesis-related pathways in PCa are limited, most studies are intrigued by the VEGF-related pathways and develop new medicines for the VEGF-related targets in PCa (Eisermann and Fraizer 2017; Sarkar et al. 2020). However, very few novel therapies showed significant improvement on the survival of PCa patients (McKay et al. 2016;Michaelson et al. 2014;Horti et al. 2009;Tannock et al. 2013). Endocrine therapy is still the common treatment for advanced PCa patients by far, including ADT and novel hormone treatments, etc. (Zhu and Ye 2020; Tietz and Dehm 2020; Mollica et al. 2021). Nevertheless, endocrine therapies give rise to CRPC (Eisermann and Fraizer 2017). Whereas, the integrated angiogenesis activity has been found to be correlated with PRAD progression, thus the identified 19 ARGs instead of one single ARG could become the molecularly therapeutic targets for PRAD.
Hallmark gene set analysis and KEGG pathway enrichment analysis showed significant differences in the biological processes of patients between the high and low ARGs signature score group. Interestingly, the E2F targets, G2M checkpoint pathways, and cell cycle pathways positively correlated with patients in high ARGs signature score group, and these signaling pathway were all vital to the tumor progression and metastasis (Löbrich and Jeggo 2007). For the targeted therapy, cell cycle pathway inhibitors could be used to improve the prognosis of PRAD. Besides, the functional analyses further showed the significant differences of cell metabolism, including fatty acid metabolism, bile acid metabolism and peroxide metabolism/peroxisome, between the high and low ARGs signature score group. As was reported previously, the aberrant regulation of fatty acid metabolism, bile acid metabolism and peroxide metabolism would cause the accumulation of reactive oxygen species (ROS) and/or peroxide, leading to the oxidative stress (Sies and Jones 2020). Generally, oxidative eustress promoted the activities of angiogenesis and cell proliferation, differentiation, and migration, but oxidative distress contributed the inflammation, fibrogenesis, and tumor metastasis (Sies and Jones 2020;Hayes et al. 2020). Therefore, specific targeted drugs in the regulation of redox signaling pathways could be suggested for the treatment of PRAD, while which needs to be further explored. We also investigated the infiltration of immune cell, and found that patients in the high ARGs signature score group had significantly higher ESTIMATE score correlated with the lower tumor purity (Yoshihara et al. 2013b), but who had significantly higher immune score as well as higher expression of immune checkpoint genes, suggesting that the immune therapies could be more effective for the high-risk PRAD patients, however, the higher presence of Tregs in patients with high ARGs signature may hinder them from benefits from immune checkpoint inhibitors.
As demonstrated, neoadjuvant ADT combining leuprolide, abiraterone acetate, and prednisone had the capacity of diminishing the ARGs signature score, which predicted the relationship between ARGs and PCa, suggesting this signature had the potential to serve as a response-related biomarker (Sowalsky et al. 2018). However, it was also found that LHRH agonists and orchiectomy exerted no significant difference on the ARGs signature score in CRPC patients. These differences may attribute from the involved patients were encountered with different status of disease and different treatments. On the other hand, the ARGs signature score could also predict the efficacy of treatments, which would be of great merit to the management of PCa. In addition, the chemotherapy was often used in treating advanced PCa as well. Findings in the evaluation of chemotherapy responses demonstrated that relatively advanced patients who had the high ARGs signature score were more suitable to take the treatment with cisplatin, paclitaxel, docetaxel, vinblastine and doxorubicin, all of which had the ability to inhibit angiogenesis in prostate tumors (Shankar et al. 2005;Quesada et al. 2005;Karashima et al. 2002;Kumar et al. 2005;Reiner 2006). Among them, docetaxel is the standard care for treating the mCRPC patients, and as suggested in previous study, it was also recommended to be early used for metastatic hormone-sensitive PCa (mHSPC) (Markowski and Carducci 2017). However, there were few experimental studies demonstrating whether these 19 ARGs were collectively associated with sensitivity to the five identified chemo drugs except studies that docetaxel could markedly suppress the expression of VEGFA (Zhu et al. 2014), and the VEGFA-VEGFR2 signaling pathway could be regulated when receiving the chemo drugs of cisplatin or paclitaxel in prostate cancer (Yun et al. 2021). Therefore, in the future work it is necessary to figure out if these 19 ARGs were collectively associated with the sensitivity to the identified chemotherapy drugs in clinical trials. The established ARGs signature may have the potential to distinct the mHSPC patients who may have more clinical benefit to ADT with docetaxel than ADT alone or combination with novel hormone treatments, which merits further validation. Overall, the high ARGs signature score can indicate the relatively advanced PRAD, which is of great value for clinicians to make better decisions for the management of PRAD.
Though we have established a robust and meaningful tool in PCa based on the ARGs, there were some limitations. Firstly, our study was based on the public database, thus it would be necessary to further verify these findings in prospective cohort from our hospital. In addition, the roles of angiogenesis involved in the progression of PCa are also needed to be explored in the mouse model and/ or the cell lines.

Conclusions
This is the first reported angiogenesis-related multigene signature in PCa. The present study determined a 19-ARG signature scoring model and a prognostic nomogram, which were very robust, stable and reliable for predicting the clinic outcomes of PRAD patients. In addition, this defined 19 ARGs expression could also be promoted as a decision-making tool for selecting treatment strategy to improve the outcomes of PRAD patients.


s2fieldsofstudy: Biology Medicine

text:
A prognosis marker MUC1 correlates with metabolism and drug resistance in bladder cancer: a bioinformatics research

Background MUC1 is a type I transmembrane protein that plays an important role in tumor cell signal transduction. Although current studies have shown that MUC1 is upregulated in bladder cancer (BC), the specific mechanism is still unclear. Methods We performed expression analysis, gene set enrichment analysis, survival analysis, immune infiltration analysis, drug sensitivity analysis, and metabolism-related gene expression analysis on TCGA-BLCA, GES31684 and GSE13507. Results The expression of MUC1 in the tumor and lymphatic metastasis positive samples was significantly increased. Genes related to MUC1 expression were significantly enriched in immune response, ribosomes, exosomes, and energy metabolism. The results of the immune infiltration analysis showed that M1 macrophages in BC with high MUC1 expression were significantly decreased. Expression of MUC1 increases drug resistance in BC patients. In addition, MUC1 increases glycolysis, glucose uptake, and lactate production by inducing metabolic reprogramming. Conclusion MUC1 has a significant effect on the metabolism and immune cell infiltration of BC, which may be the cause of increased drug resistance, and can be used as a molecular target for the diagnosis and treatment of BC.

Introduction
Bladder cancer (BC) is a malignant tumor that seriously endangers human health. There are approximately 573,000 new cases of BC each year worldwide, and the number of deaths caused by BC is approximately 213,000 ranking tenth among all tumors [1]. BC is divided into muscle invasive bladder cancer (MIBC) and non-muscle invasive bladder cancer (NMIBC) according to the depth of invasion [2]. BC mainly includes transitional cell carcinoma, squamous cell carcinoma, and adenocarcinoma on the basis of the pathological type, among which transitional cell carcinoma accounts for more than 90% [3]. There are plenty of risk factors in the development of BC, including gender, smoking prevalence, occupational carcinogen exposure, Schistosoma haematobium infection, etc. [4]. Recently, researchers have considered the interaction of hypoxia and tumor, immune cell infiltration and immune microenvironment as the key step in the occurrence and development of BC [5,6]. Currently, early diagnosis of BC can rely on serum markers, urine markers, and cystoscopy. However, the lack of clear markers and poor patients' compliance make the early diagnosis of BC very difficult [7]. Despite the 5-year survival rate for BC is as high as 77%, the recurrence rate of BC is still high, and the 5-year survival rate for metastatic BC is less than 5% [8]. Therefore, it is particularly urgent to find highly sensitive biomarkers for BC.
MUC1 encodes a type I transmembrane protein with a polar distribution, called Mucin 1. N-terminal alpha subunit of Mucin1 functions in cell adhesion and the C-terminal beta subunit is involved in cell signaling [9]. Recent studies have found that MUC1 plays an important role in tumor hypoxia microenvironment and tumor metabolism [10]. MUC1 can not only mediate hypoxia-driven angiogenesis by regulating a variety of angiogenic factors; it can also stabilize and activate hypoxia inducible factor-1α (HIF-1α) to promote the metabolic reprogramming of cancer cells [11,12]. Recently, studies have developed antibody-drug conjugated (ADC) based on MUC1 employed in the treatment trastuzumab-resistant breast cancer patients [13]. Besides, a multitargeted recombinant Ad5 PSA/ MUC-1/brachyury-based immunotherapy vaccine developed by Bilusic et al. was recently introduced into phase I trials [14]. The types and abundance of immune cell infiltration in the tumor microenvironment have proved to be very different from normal tissues [14]. Studies have found that the expression of MUC1 induces immune suppression in colon cancer, and this suppression can be reversed by blocking the PD1/PDL1 pathway [15]. MUC1 has shown a strong potential in the diagnosis and treatment of tumors.
Although some studies claim that MUC1 has significance in BC, the specific mechanism remains unclear. In this study, we used gene expression data obtained from the the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database to explore the expression of MUC1 in human BC samples. We used R (version 4.0.3) to analyze the correlation between MUC1 expression and clinical characteristics. To better understand the mechanism of MUC1 in the occurrence and development of BC, we conducted Gene Set Enrichment Analysis (GSEA) and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis [16]. Next, we used Tumor Immune Estimation Resource (TIMER) and CIBERSORT to analyze the relationship between MUC1 expression and tumor immune cell infiltration. Cancer Cell Line Encyclopedia (CCLE) and Genomics of Drug Sensitivity in Cancer (GDSC) was used to analyze the relationship between MUC1 expression and drug sensitivity. Finally, we use Gene Expression Profile Interaction Analysis (GEPIA), Kaplan-Meier survival analysis (KM), and Human Protein Atlas (HPA) to analyze the relationship between MUC1 and patient prognosis.

Data retrieval and download
To explore the expression of MUC1 in BC and its relationship with clinical characteristics, gene expression data and clinical data of BC in TCGA (TCGA-BLCA cohort) were downloaded from the Genomic Data Commons Data Portal (https:// portal. gdc. cancer. gov). Besides, we downloaded the gene expression data and survival data of GSE31684 and GSE13507 as a validation dataset from the GEO (https:// www. ncbi. nlm. nih. gov/ geo) [17,18].

Gene expression analysis
We compared the expression of MUC1 in cancer and adjacent normal tissues, as well as in positive and negative lymphatic metastasis tissues. With the GSE13507 dataset, we compared MUC1 expression in normal, paraneoplastic, primary and recurrent tumour tissues. In addition, we used GEPIA to show the expression of MUC1 in bladder tumor tissues and adjacent tissues (http:// gepia. cancer-pku. cn/) [19]. Next, we removed samples with incomplete clinical data (age, gender, American Joint Committee on Cancer (AJCC) pathologic stage, AJCC pathologic N, AJCC pathologic M, AJCC pathologic T). The protein expression of MUC1 between normal and cancer tissues was retrieved and compared from the HPA database (www. prote inatl as. org).

Prognostic analysis
Univariate and multivariate analyses were conducted using the Cox Proportional Hazards Regression model on the remaining samples to obtain risk scores. To further understand the relationship between MUC1 expression and clinical characteristics, we conducted receiver operating characteristic (ROC) analysis and logistic regression analysis. After dividing the samples into high and low expression group according to the median expression of MUC1, KM survival analysis was conducted to analyze the overall survival (OS) rate of BC patients related to MUC1 expression. Furthermore, we conducted a similar analysis in GSE31684 and GES13507 to verify the results obtained previously.

Gene set enrichment analysis
To explore the mechanism of MUC1 in BC, we divided the samples into MUC1-H and MUC1-L groups according to the level of MUC1 expression. Next, the gene differential expression analysis was performed on two sets of samples, and the genes that were differentially expressed were obtained (logFC > 2, p value < 0.05). Incorporate differentially expressed genes in GSEA analysis. The number of permutations was set to 1000. To investigate the possible biological functions of MUC1, we used GSEA to analyze the GO pathway and the KEGG pathway. The absolute value of normalized enrichment score (|NES|) is greater than 2, nominal p value (NOM p value) is less than 0.05, and false discovery rate (FDR) is less than 0.05 of all enrichment pathways.

Immune infiltrate analysis
Li et al. developed a deconvolution algorithm (TIMER) to calculate the type and abundance of immune cell infiltration based on the gene expression profiles (http:// timer. cistr ome. org/) [20]. We evaluated MUC1 expression in BC and its correlation with the type and abundance of immune cell infiltration, including B cells, CD4 + T cells, CD8 + T cells, Dendritic cells, Neutrophils and Macrophages. Another deconvolution algorithm based on gene expression called CIBERSORT (http:// ciber sort. stanf ord. edu/) was used to further assess the relationship between MUC1 and immune cell [21,22]. And we performed a nonparametric Wilcox test to assess whether the infiltration of immune cells differed between MUC1-H and MUC-L groups. In addition, we also assessed the correlation between different immune cells.

Metabolism analysis
The samples in the TCGA, BLCA, and GSE13507 datasets were divided into two groups according to the expression of MUC1, and the expression levels of genes related to glucose metabolism (HK2, GLUT1, PKM, RWDD3, SLC16A3, SLC5A12, ENO1) in different groups were compared. Data were analyzed using t-test and p value < 0.05 was considered statistically significant.

Expression of MUC1 at the mRNA and protein levels and its relationship to clinical characteristics
We compared MUC1 mRNA expression between BC and normal tissues in TCGA database, and also between lymphphatic metastasis negative and positive samples. A total of 411 tumor tissues and 19 normal tissues were enrolled in the expression analysis, including 233 lymphatic metastasis negative and 81 lymphatic metastasis positive samples. The expression of MUC1 in different groups was represented by Scatterplot (Fig. 1A). The results showed that the expression of MUC1 in cancer tissues was significantly higher than that in normal tissues (p value = 2.2e−2), and it was significantly higher in lymph node metastasis positive samples than in lymph node metastasis negative (p value = 4.2e−5). The GSE13507 dataset contains normal, paraneoplastic, primary and recurrent tumour samples. We analysed the expression of MUC1 in different subgroups of GSE13507. The results showed significant differences in MUC1 between primary tumours and paraneoplastic tissues (p value < 0.002). Similarly, MUC1 was significantly different between primary tumours and normal tissues (p value = 0.0298) (Fig. 1B). Furthermore, we used GEPIA to compare the expression of MUC1 in normal and tumor tissues and found that the mRNA expression of MUC1 significantly increased in tumor tissues (p value < 0.01) (Fig. 1C).
We evaluated the connection between clinical characteristics and MUC1 expression levels of BC patients. There are 382, 382, 192, 350 samples with complete data of AJCC pathologic Stage, AJCC pathologic T stage, AJCC pathologic M stage, and AJCC pathologic N stage, respectively. The results showed that MUC1 has little correlation with AJCC pathological T stage (p value > 0.05, Fig. 1D) and AJCC pathological Stage (p value > 0.05, Fig. 1F), but is significantly correlated with AJCC pathological N stage (p value < 0.05, Fig. 1E) and AJCC pathological M stage (p value < 0.05, Fig. 1G). In addition, we conducted logistic regression to analyze the correlation between MUC1 expression and clinical characteristics ( Table 1). The increase of MUC1 expression level in BC was significantly correlated with lymphatic metastasis (N1 vs. N0, p value = 0.007). These indicate that compared with BC patients with low levels of MUC1 expression, BC patients with high levels of MUC1 expression are more likely to develop lymphatic metastasis. Next, to understand the difference in protein levels of MUC1, we searched HPA for the immunohistochemical data of MUC1 in BC. The results showed that MUC1 has higher levels of expression in tumor compared to normal tissues (Fig. 1H).

The role of MUC1 in the prognosis of bladder cancer
We used univariate Cox proportional hazards regression analysis to explore the relationship between MUC1 expression, clinical characteristics, and OS. The results showed that Age (HR  (Table 2). Next, multivariate Cox proportional hazards regression analysis showed that MUC1 expression is an independent prognostic factor ( Fig. 2A, Table 3).
The distribution of MUC1 expression, survival status of patients with BC, and expression profiles of MUC1 are shown in Fig. 2B. ROC curve analysis showed that the AUC values of 1, 3, and 5 years were 0.6, 0.634, and 0.634, respectively, indicating that the expression of MUC1 has a potential prognostic ability (Fig. 2C).
Next, we performed a survival analysis in TCGA based on the expression of MUC1. The results of TCGA show that MUC1 expression is significantly related to OS (p value = 1.47e−03) (Fig. 2D). To further verify the influence of MUC1 on the OS of BC patients, we further conducted similar analysis in GSE13507 and GSE31684. As expected, similar results appeared in GSE31684 (p value = 1.21e−03) and GSE13507 (p value = 2.66e−02) (Fig. 2E, F).

MUC1 related enrichment pathway
To explore the potential biological functions of MUC1, we performed GO term and KEGG pathway enrichment analyses of differentially expressed genes. Signal pathways that meet the following conditions (|NES| > 2, NOM p value < 0.05, FDR < 0.05) are considered enriched. As shown in Table 4, there are five positive GO Biological Process terms related to MUC1 (nuclear-transcribed mRNA catabolic process, nonsense-mediated decay; cotranslational protein targeting to the membrane; SRP-dependent cotranslational protein targeting to the membrane; protein targeting to ER; establishment of protein localization to endoplasmic reticulum) and five negative GO Biological Process terms related to MUC1 (cornification; keratinization; immune responseregulating signaling pathway; immune response-regulating cell surface receptor signaling pathway; defense response to other organism) (Fig. 3A). There are five positive GO Cellular Component terms related to MUC1 (cytosolic large ribosomal subunit; cytosolic ribosome;   (Fig. 3C). As shown in Table 5, there are two positive KEGG pathways related to MUC1 (Coronavirus disease, COVID-19 and Ribosome) and two negative KEGG pathways related to MUC1 (Staphylococcus aureus infection and Intestinal immune network for IgA production) (Fig. 3D). These results indicate that the pathways that regulate ribosomal and lipid metabolism are critical to BC patients and are closely related to the expression of MUC1.

Relationship between MUC1 expression and immune cell infiltration
We used TIMER to clarify the correlation between the expression of MUC1 and the level of immune infiltration of BC. The results have shown that MUC1 expression is significantly correlated with B cells (p value = 1.50e−03), macrophages (p value = 3.593e−02) and Neutrophils (p value = 3.94e−04), indicating that MUC1 plays a key role in the immune infiltration of BC (Fig. 4A). To explore the impact of MUC1 expression on the immune microenvironment of BC patients, we divided the samples into MUC1-H and MUC1-L according to MUC1 expression. CIBERSORT was used to calculate the infiltration abundance of 22 immune cells in each sample to assess the difference between different groups. The results indicate that M1 macrophages are affected by the expression of MUC1 (p value = 0.033) (Fig. 4B). We also assessed the potential correlation between immune cell infiltration and clinical features. The results showed a strong correlation between immune cell infiltration and clinical features (Fig. 4C).

Discussion
Recently, the role of MUC1 in tumors has gradually been emphasized. In pancreatic cancer, MUC1 reduces radiation-induced pancreatic cancer cell toxicity and DNA damage by enhancing glycolysis, the pentose phosphate pathway, and nucleotide biosynthesis [23]. In triple-negative breast cancer, MUC1-C increases PD-L1 transcription and promotes tumor immune escape by recruiting MYC and NF-κB p65 to the PD-L1 promoter region [24]. In prostate cancer, MUC1-C inhibits androgen receptor (AR) axis signaling and induces the expression of the neural BRN2 transcription factor, increasing the self-renewal capacity and tumorigenicity of prostate cancer cells [25].
In addition, MUC1 also plays an important role in ovarian cancer [26], non-small cell lung cancer [27] and pancreatic cancer [28]. One study found that the positive rates of MUC1 and MUC2 in urothelial carcinoma were 89.7% and 44.3%, respectively, and MUC1 expression was significantly correlated with tumor grade, while MUC2 was the opposite. Different staining patterns of MUC1 and MUC2 have important prognostic implications in patients with difficult-to-recognize low-grade papillary urothelial lesions and precancerous lesions [29]. Urothelial carcinoma patients expressing high levels of MUC1-C have poor survival and are strongly associated with cisplatin resistance [30]. However, the specific mechanism of MUC1 in BC remains unclear. We analyzed the potential mechanism of MUC1 in BC (Fig. 5D). In our current study, to explore the possibility of MUC1 as a prognostic biomarker for BC, we used TCGA BLCA cohort to evaluate the prognostic value of MUC1. First, we analyzed the expression of MUC1 (protein level and mRNA level), and the results showed that the expression of MUC1 in tumor tissues was significantly higher than that in adjacent normal tissues, and the expression of MUC1 in lymphatic metastasis positive samples was significantly higher than that in lymphatic metastasis negative samples. BC patients with high expression of MUC1 showed advanced N stage, M stage, and tumor status. We analyzed the role of MUC1 and clinical features in prognosis, and the results showed that high MUC1 expression was strongly associated with worse prognosis. These results suggest that MUC1 may serve as an independent prognostic factor for overall survival in BC.
MUC1 is closely related to metabolism and resulting drug resistance. Recent studies have found that there is a crosstalk between the MUC1 and HIF-1 signal pathways. The HIF-1 signaling pathway has been shown to be abnormally activated in numerous cancers, especially solid tumors [31]. HIF-1 consists of two subunits, HIF-1α and HIF-1β, of which the α subunit is the active unit. In oxygenated cells, HIF-1α will be hydroxylated by proline hydroxylase (PHD1), which will then be recognized and degraded by the ubiquitin-protease complex; while in hypoxic cells, the activity of PHD1 decreases and the lack of a substrate for the hydroxylation reaction leads to the accumulation of HIF-1α [32]. The transcriptional targets of HIF-1α include genes that play key roles in angiogenesis, erythropoiesis, metabolic reprogramming, vasomotor function, and apoptosis proliferation response [33]. In the hypoxic condition, MUC1 promotes the recruitment of HIF-1α and p300 to the promoter region of glycolysis genes, upregulates the expression of glycolysis-related genes, increases glucose uptake and lactate production, and assists tumor cells to survive in hypoxic environments [12]. In addition, MUC1 can increase the carbon flux of the non-oxidative pentose phosphate pathway (PPP) and pyrimidine synthesis pathway of pancreatic cancer cells by promoting the expression of HIF-1α, leading to gemcitabine resistance [34].
In this study, GLUT1, PKM, and SLC16A3 were significantly upregulated in the MUC1-H group, while RWDD3 was significantly downregulated. GLUT1 encodes the major glucose transporter in cell [35]. Pyruvate kinase is encoded by the PKM and is the main ratelimiting enzyme in the glycolysis process [36]. SLC16A3 removes lactate produced by glycolysis from the plasma membrane and plays a key role in the growth and proliferation of hypoxic cancer cells [37]. RWDD3 negatively regulates the HIF-1α signaling pathway by increasing the sumoylation of HIF-1α, promoting its stabilization, transcriptional activity, and the expression of its target gene VEGFA during hypoxia [38]. Therefore, we can draw the hypothesis: MUC1 increases the uptake of glucose and the production of lactate through the HIF-1α signaling pathway, so that tumor cells can survive under hypoxic conditions. In addition, GO and KEGG pathway enrichment analysis also showed that up-regulated MUC1 was mainly closely related to metabolism. Through the CCLE and GDSC databases, we found that the expression of MUC1 in BC was associated with drug resistance, which may be caused by increased pyrimidine metabolism, and the specific mechanism remains to be further studied. We used the TIMER database to explore the relationship between the expression of MUC1 and the level   [40]. Almost simultaneously, Maureen et al. successfully constructed a BCG-human interleukin-2 (hIL2)-MUC1 breast cancer vaccine, which significantly reduced the incidence and rate of tumorigenesis in mice, suggesting that co-expression of MUC1 and IL-2 plays a key role in the induction of tumor-specific immune responses [41]. Yuan et al. constructed a BCG-MUC1-based breast cancer vaccine (rBCG-MVNTR4-CD80), which significantly suppressed tumors, induced interferon production, and stimulated both CD4+ and CD8+ positive lymphocyte production [42]. In addition, the authors constructed two other rBCG-MUC1 vaccines (rBCG-MVNTR4-CSF and rBCG-MVNTR8-CSF, with similar effects [43]. In 2016, an anti-tumor vaccine based on MUC1 (MUC1maltose-binding protein (MBP)/BCG) entered preclinical trials and showed no significant toxicity in mice, rats or crab-eating monkeys.MUC1-MBP/BCG acts as an antitumor agent by inducing Th1 cell activation and MUC1specific IgG antibody secretion [44]. At present, most of the research on MUC1 and BCG is targeting MUC1 and building tumor vaccines based on BCG to better play the role of BCG. Therefore, high MUC1 expression can theoretically improve the responsiveness of tumors to such rBCG drugs, but there is no research on the relationship between traditional BCG and MUC1 expression. This needs further study.
MUC1 has shown great potential as a widely expressed membrane protein in the diagnosis and treatment of tumors. However, the mechanism of action of MUC1 in BC has not been reported in the literature. The role of MUC1 in BC was investigated by expression analysis, prognostic analysis, immunological analysis, metabolic analysis and drug sensitivity analysis. The potential mechanisms of MUC1-mediated drug resistance in bladder cancer were uncovered. However, there are some shortcomings in this study: no experimental study was conducted to verify the function of MUC1, and the role of MUC1 needs to be verified by cellular and animal experiments subsequently.

Conclusion
In short, this is the first report to explore the possibility of MUC1 as a marker of BC at the mechanistic level. This work helps to elucidate the role of cellular metabolism, immune cell infiltration, and genes regulating these processes in BC development and drug resistance. In addition, MUC1 can be used as a target for monoclonal antibodies, antibody-drug conjugate, cancer vaccine, Aptamer, and radiotherapy. Understanding the mechanism of MUC1 will help to develop corresponding treatment and diagnosis methods.


extfieldsofstudy: Medicine Biology

text:
Genetic background and mistranslation frequency determine the impact of mistranslating tRNASerUGG

Abstract Transfer RNA variants increase the frequency of mistranslation, the misincorporation of an amino acid not specified by the “standard” genetic code, to frequencies approaching 10% in yeast and bacteria. Cells cope with these variants by having multiple copies of each tRNA isodecoder and through pathways that deal with proteotoxic stress. In this study, we define the genetic interactions of the gene encoding tRNASerUGG,G26A, which mistranslates serine at proline codons. Using a collection of yeast temperature-sensitive alleles, we identify negative synthetic genetic interactions between the mistranslating tRNA and 109 alleles representing 91 genes, with nearly half of the genes having roles in RNA processing or protein folding and turnover. By regulating tRNA expression, we then compare the strength of the negative genetic interaction for a subset of identified alleles under differing amounts of mistranslation. The frequency of mistranslation correlated with the impact on cell growth for all strains analyzed; however, there were notable differences in the extent of the synthetic interaction at different frequencies of mistranslation depending on the genetic background. For many of the strains, the extent of the negative interaction with tRNASerUGG,G26A was proportional to the frequency of mistranslation or only observed at intermediate or high frequencies. For others, the synthetic interaction was approximately equivalent at all frequencies of mistranslation. As humans contain similar mistranslating tRNAs, these results are important when analyzing the impact of tRNA variants on disease, where both the individual’s genetic background and the expression of the mistranslating tRNA variant need to be considered.

Introduction
Mistranslation occurs when an amino acid that differs from what is specified by the standard genetic code is incorporated into a growing polypeptide chain during translation. Mistranslation occurs in all cells but can be enhanced by environmental conditions or mutations in the translational machinery (Lee et al. 2006;Netzer et al. 2009;Ling and Sö ll 2010;Jones et al. 2011;Wiltrout et al. 2012;Reverendo et al. 2014;Schwartz and Pan 2017;Lant et al. 2018). Mutations in tRNAs that cause mistranslation were initially identified as intergenic suppressors that change the meaning of the genetic code Stadler and Yanofsky 1959;Yanofsky and Crawford 1959;Benzer and Champe 1962;Gorini and Beckwith 1966). tRNA Ser variants are particularly prone to mistranslate because the anticodon is not a major identity element for aminoacylation by the cognate serine aminoacyl-tRNA synthetase (Gieg e et al. 1998). Rather, specificity for aminoacylation comes from the long variable arm positioned 3 0 of the anticodon stem (Asahara et al. 1994;Biou et al. 1994;Himeno et al. 1997).
The toxic effects of mistranslating tRNAs are buffered through multiple copies of each tRNA isodecoder (e.g. there are 275 tRNAencoding genes in Saccharomyces cerevisiae; Chan and Lowe 2016) and through protein quality control mechanisms that deal with misfolded protein and protein aggregates (reviewed in . When mistranslation reaches a threshold, protein quality control mechanisms no longer protect the cell and growth is impaired (Berg, Zhu, et al. 2019). The extent of growth impairment is inversely related to the frequency of mistranslation in a linear fashion with yeast growth being arrested when mistranslation approaches approximately 12% .
We previously demonstrated that the negative genetic interactions with mistranslating tRNAs depend on the amino acid substitution (Berg, Zhu, et al. 2021). At similar frequencies of mistranslation, a tRNA variant substituting serine at arginine codons has more genetic interactions than one substituting alanine at proline codons. In this report, we identify the genetic interactions of a mistranslating serine tRNA variant that incorporates serine at proline codons. Using a regulated tRNA expression system, we show that although there is a general correlation between the frequency of mistranslation and impact on growth, the impact of different mistranslation frequencies depends on a strain's specific genetic background. As similar mistranslating tRNAs are found in the human population, these results suggest that genetic background contributes to the impact of tRNA variants on health and disease.

Materials and methods
Yeast strains and growth BY4742 (MATa his3D0 leu2D0 lys2D0 ura3D0; Brachmann et al. 1998) and Y7092 (MATa can1D::STE2pr-SpHIS5 lyp1D his3D1 leu2D0 ura3D0 met15D0) strains are derivatives of S288c. Y7092 was a kind gift from Dr. Brenda Andrews (University of Toronto). Strains from the temperature-sensitive collection are derived from the wild-type MATa haploid yeast strain BY4741 and described in Costanzo et al. (2016). The strains containing the gene expressing tRNA Ser UGG,G26A (CY8613) were made by integrating modified SUP17 and flanking sequence into Y7092 at the HO locus and selecting for the natMX marker, as previously described in Zhu et al. (2020), using the construct described below. The control strain (CY8611) was made by integrating only the natNT2 marker at the HO locus. Transformants were selected on 100 mg/mL nourseothricin-dihydrogen sulfate (clonNAT) and integration was verified by PCR.
Yeast strains were grown at 30 C in yeast peptone media containing 2% glucose (YPD) or synthetic media supplemented with nitrogen base and amino acids, unless otherwise indicated. Growth curves were generated by diluting saturated cultures to OD 600 ¼ 0.1 in synthetic complete media and incubating at 30 C. OD 600 was measured every 15 min for 24 h using a BioTek Epoch 2 microplate spectrophotometer. Doubling time was calculated using the R package "growthcurver" (Sprouffske and Wagner 2016).

DNA constructs
The construct to integrate the gene encoding tRNA Ser UGG,G26A at the HO locus was created using a synthetic DNA containing 200 bp up and downstream of the HO translational start, previously described in Zhu et al. (2020). The construct was cloned into pGEM-T Easy (Promega Corp.) as a NotI fragment to create pCB4386. The natNT2 marker from pFA6-natNT2 was PCR amplified using primers UK9789/UK9790 (Supplementary Table 1) and cloned into pCB4386 as an EcoRI fragment to generate the control SGA integrating vector (pCB4394). The gene encoding tRNA Ser UGG,G26A was PCR amplified from pCB4023 (Berg et al. 2017) using primers UG5953/VB2609 and cloned as HindIII fragments into pCB4394 to create pCB4397.

Synthetic genetic array analysis and validation
The SGA assay was performed as described by Tong et al. (2001) with minor modifications. Strains CY8611 (HO::natMX) and CY8613 (HO::tRNA Ser UGG,G26A -natMX) were crossed to a yeast temperature-sensitive collection (Ben-Aroya et al. 2008;Li et al. 2011;Kofoed et al. 2015;Costanzo et al. 2016) in quadruplicate 1,536 colony array format using a BioMatrix (S&P Robotics Inc.) automated pinning system. In this format, each allele of the temperaturesensitive collection is present in technical quadruplicate on the plate. Double mutants were selected on YPD plates containing 200 mg/L G418 and 100 mg/L clonNAT. Diploids were sporulated on enriched sporulation media and MATa haploid double mutants selected using standard SGA media. The entire SGA procedure was carried out at room temperature, except for the colony scoring in order to minimize growth defect of the temperature-sensitive strains. To identify genetic interactions, double mutants were pinned onto double mutant selection SGA medium and grown at 30 C for 5 days. Images were taken every 24 h to determine colony size computationally. SGATools (Wagih et al. 2013) was used to determine genetic interaction scores using a multiplicative model (e ¼ W AB ÀW A * W B ; where e is the interaction score, W AB is the fitness of the double mutant and W A and W B are the fitness values of each single mutant). Double mutant strains with an average interaction score less than À0.2 and Benjamini-Hochberg corrected P-value less than 0.05 were considered as potential negative genetic interactions.
Double mutants that were identified as negative genetic interactions from the screen were validated by recreating the double mutant strain, starting from the single mutant haploid strains, using the SGA approach. Double mutant strains were grown in liquid media to saturation, cell densities were normalized, and cultures were spotted on SGA media. The temperature-sensitive mutant crossed with the control strain CY8611 and the mistranslating tRNA strain crossed with a control his3D strain were also spotted to determine fitness of the single mutants. Intensity of each spot was measured with ImageJ (Schneider et al. 2012). Expected double mutant growth was calculated based on the growth of the single mutants and compared to the experimental growth of the double mutant. Double mutants that grew more slowly than expected were considered validated negative genetic interactions. Raw and validated data can be found in Supplementary File 2.
Synthetic interactions with various frequencies of mistranslation were assessed by transforming the relevant temperature-sensitive strains with URA3-containing centromeric plasmids expressing tRNA Ser (pCB3076), tRNA Ser UGG,G26A (pCB4023), tRNA Ser UGG,G26A with 5 0 GAL1pr (pCB4568), or tRNA Ser UGG,G26A with 3 0 GAL1pr (pCB4566). At least 3 individual transformants for each plasmid and strain were grown in synthetic complete medium lacking uracil and containing 2% galactose as the carbon source. Cells were grown to confluency, diluted 1:32 in 1Â yeast nitrogen base and 5 lL was spotted onto solid media lacking uracil and containing 2% galactose. Cells were grown at 30 C for 32-56 h, depending on the strain, to achieve a level of growth (for the strain without mistranslation) consistent with the wild type. Mean density of growth was determined with ImageJ (Schneider et al. 2012) and normalized growth for each mistranslating tRNA was calculated as a percent of the wild-type tRNA Ser containing strain.

Heat shock assay
Yeast strains containing the HSE-GFP reporter and a mistranslating tRNA variant were grown to stationary phase in medium lacking uracil and containing 0.6% casamino acids, diluted 1:100 in the same medium and grown for 18 h at 30 C. Cell densities were normalized to OD 600 before measuring fluorescence with a BioTek Synergy H1 microplate reader at an excitation wavelength of 488 nm and emission wavelength of 528 nm. The mean relative fluorescence units were calculated from 3 technical replicates for each biological replicate.

Mass spectrometry
Liquid chromatography-tandem mass spectrometry to identify mistranslation was performed on 5 biological replicates of each strain. Starter cultures of each strain were grown overnight in YPD before being diluted to an OD 600 of 0.1 in the same media and grown to an OD 600 of approximately 1.0. Cells were lysed in a urea lysis buffer (8 M Urea, 50 mM Tris pH 8.2, 75 mM NaCl) and proteins were reduced with dithiothreitol and alkylated with iodoacetamide. Robotic purification and digestion of proteins into peptides were performed on the KingFisher Flex using LysC and the R2-P1 method described in Leutert et al. (2019).
MS/MS spectra were searched against the S. cerevisiae protein sequence database (downloaded from the Saccharomyces Genome Database resource in 2014) using Comet (release 2015.01; Eng et al. 2013). The precursor mass tolerance was set to 50 ppm. Constant modification of cysteine carbamidomethylation (57.0215 Da) and variable modification of methionine oxidation (15.9949 Da) and proline to serine substitution (À10.0207 Da) were used for all searches. A maximum of 2 of each variable modification was allowed per peptide. Search results were filtered to a 1% false discovery rate at the peptide spectrum match level using Percolator (Kä ll et al. 2007). The mistranslation frequency was calculated using the unique mistranslated peptides for which the non-mistranslated sibling peptide was also observed. The frequency is defined as the counts of mistranslated peptides, where serine was inserted for proline, divided by the counts of all peptides containing proline, and expressed as a percentage.

Results and discussion
Synthetic genetic interactions with tRNA Ser

UGG,G26A
Yeast cells expressing tRNA Ser UGG,G26A , which contains a proline anticodon, mistranslate serine at proline codons (Berg, Zhu, et al. 2019;. The G26A mutation is required in combination with the anticodon change to dampen otherwise lethal levels of mistranslation (Berg et al. 2017). To perform the SGA analysis, the gene encoding tRNA Ser UGG,G26A (Fig. 1a), including approximately 300 base pairs of 5 0 and 3 0 flanking sequence and a clonNAT resistance marker, was integrated at the HO locus. A control strain was created with only the clonNAT resistance marker integrated at the HO locus. Mass spectrometrybased analysis of the cellular proteome identified 4.9% proline to serine substitution in the strain expressing tRNA Ser UGG,G26A compared to only 0.6% substitution in the control strain (Fig. 1b). As shown in Fig. 1, c and d, tRNA Ser UGG,G26A reduces cell growth (doubling time of 98 min vs 84 for the control strain) and results in a heat shock response (6.4-times greater than the control strain). We then performed an SGA analysis to identify genetic interactions with tRNA Ser UGG,G26A using a collection of 1,016 temperature-sensitive alleles. The robotic screen identified 125 alleles with negative genetic interactions with tRNA Ser UGG,G26A . Genetic interactions were validated by remaking the double mutant strains, spotting normalized densities of the double mutants and their control strain on selective plates, and measuring growth after 2 days (Supplementary File 2). After validation, 109 alleles representing 91 genes were classified as having a negative genetic interaction with tRNA Ser UGG,G26A (Fig. 2a). To further analyze the network of genes associated with the mistranslating tRNAs, we identified areas of the yeast genetic interaction map (Costanzo et al. 2016) that were enriched for is statistically different compared to the control strain (Welch's t-test; Bonferroni corrected P-value < 0.05). c) Doubling times for the strains described in B were determined from growth curves of the strains diluted to an OD 600 of approximately 0.1 in synthetic complete media containing clonNAT and grown for 24 h. Doubling time was calculated with the R package "growthcurver" (Sprouffske and Wagner 2016). Each point represents 1 biological replicate (n ¼ 4). Doubling time is statistically different between the strain expressing tRNA Ser UGG,G26A and the control strain (Welch's t-test; Bonferroni corrected P-value < 0.05). d) Strains described in B were transformed with a GFP reporter transcribed from a promoter containing heat shock response elements, grown to saturation in media lacking uracil, diluted 1:100 in the same media and grown for 18 h at 30 C. Cell densities were normalized and fluorescence measured. Each point represents 1 biological replicate (n ¼ 5). Relative heat shock induction is statistically different in the strain expressing tRNA Ser UGG,G26A compared to the control strain (Welch's t-test; Bonferroni corrected P-value < 0.05).
negative genetic interactions with the mistranslating tRNA Ser UGG,G26A using SAFE (Baryshnikova 2016;Fig. 2b). Areas of the yeast genetic interaction network annotated with roles in protein turnover, cell polarity, mitosis, DNA replication and repair, transcription, and mRNA processing were significantly enriched.
In a previous screen looking at negative synthetic genetic interactions with tRNA variants that mistranslate alanine at proline codons (tRNA Pro G3:U70 ; Hoffman, Berg, et al. 2017) and serine at arginine codons (tRNA Ser UGG,G26A ; Berg, Zhu, et al. 2019), we identified 10 and 47 negative genetic interactions, respectively (Berg, Zhu, et al. 2021). While these tRNA variants mistranslated at lower frequency ($3%) than tRNA Ser UGG,G26A making specific comparison of genetic interactions difficult, it is possible to compare enriched pathways with a SAFE analysis (Fig. 2c) (Pro!Ala) and tRNA Ser UCU,G26A (Arg!Ser), respectively. The unique genetic interactions identified for Pro!Ala and Pro!Ser mistranslating tRNAs suggest that it is not simply the proline content of the encoded protein that predicts genetic interactions with these tRNA variants. This agrees with our previous result that genes with negative genetic interactions are not enriched for the amino acid being mistranslated (Berg, Zhu, et al. 2021).
Only one strain had a positive genetic interaction with tRNA Ser UGG,G26A . The strain contains a temperature-sensitive allele of eco1, an acetyltransferase required in sister chromatid cohesion. As we demonstrated previously (Zhu et al. 2020), the positive interaction results from tRNA Ser UGG,G26A restoring serine at the S213P mutation of eco1-1.  Fig. 2. Negative genetic interaction network of the mistranslating tRNA Ser UGG,G26A . a) Genes validated as having a negative genetic interactions with tRNA Ser UGG,G26A are arranged according to their predicted function based on gene descriptions in the yeast genome database (www.yeastgenome.org; last accessed February 2022). b) SAFE analysis of genes that have a negative genetic interaction with tRNA Ser UGG,G26A were mapped onto the yeast genetic interaction profile map (Costanzo et al. 2016) using TheCellMap (Usaj et al. 2017). Blue dots represent genes within the local neighborhood of genes validated to have negative genetic interactions with tRNA Ser UGG,G26A . Terms in blue boxes are network regions that are significantly enriched (Bonferroni corrected P-value < 0.05). c) SAFE analysis as performed in B with genetic interactions for tRNA Pro G3:U70 (yellow) and tRNA Ser UGG,G26A (red) from Berg, Zhu, et al. (2021). Terms in boxes represent network regions that are significantly enriched for the, respectively, mistranslating tRNA (Bonferroni corrected P-value < 0.05). et al. 2019; Lant et al. 2019). Due to their ability to generate proteotoxic stress, we and others have suggested that these variants may be genetic modifiers of disease (Reverendo et al. 2014;Berg et al. 2017). In a previous analysis, we demonstrated that, when comparing the same amino acid substitution, there is a near-linear negative correlation between mistranslation frequency and cell growth in a wild-type S. cerevisiae background (Berg, Zhu, et al. 2019). Our goal was to determine if the genetic background of strains having a synthetic interaction with tRNA Ser UGG,G26A changes the impact of different mistranslation frequencies.
To regulate the frequency of mistranslation, we took advantage of our finding that placing a GAL1 promoter (GAL1pr) sequence up or downstream of a tRNA gene represses tRNA expression when cells are grown in galactose . Previously, using mass spectrometry, we determined the frequency of proline to serine mistranslation in strains containing centromeric plasmids expressing wild-type tRNA Ser , tRNA Ser UGG,G26A , tRNA Ser UGG,G26A with 5 0 GAL1pr, and tRNA Ser UGG,G26A with 3 0 GAL1pr in galactose containing medium to be 0.3%, 5.6%, 3.5%, and 0.9%, respectively ( Fig. 3a; . Three frequencies of mistranslation above background are thus achieved: the greatest mistranslation with no flanking GAL1pr, intermediate mistranslation with the GAL1pr upstream of the tRNA, and the least mistranslation with the GAL1pr downstream of the tRNA. To determine the impact of mistranslation frequency, we selected 16 of the strains identified in the SGA analysis, representing different pathways, that had a robust negative genetic interaction but not so severe that there was no dynamic range to investigate the effects of different levels of mistranslation. These and the control strain were transformed with centromeric plasmids containing wild-type tRNA Ser , tRNA Ser UGG,G26A 5 0 GAL1pr-tRNA Ser UGG,G26A , and 3 0 tRNA Ser UGG,G26A -GAL1pr. Triplicate cultures of independent transformants for each strain were grown to confluency, diluted 1:32, spotted onto minimal medium with galactose as the carbon source, and grown at 30 C. Spot plates are shown in Supplementary Fig. 1. The density of the spotted cultures was measured and then expressed as a percentage of the spot density for the same strain background containing wildtype tRNA Ser . For example, tRNA Ser UGG,G26A reduces the growth of the wild-type BY4742 strain to 37% 6 5% of growth seen in BY4742 containing wild-type tRNA Ser . The 5 0 GAL1pr- , and tRNA Ser UGG,G26A [tS(UGG)] used to regulate proline to serine mistranslation frequency. Mistranslation frequencies were measured by mass spectrometry in . b) Wild-type BY4742 or the indicated strains from the temperature-sensitive collection were transformed with the constructs described in A. Strains were grown to confluency in media lacking uracil and diluted 1:32 and spotted on media lacking uracil with galactose as the carbon source. The spot intensity of the strain containing the mistranslating tRNA was divided by the intensity of the strain containing the wild-type tRNA Ser to determine normalized growth. Each point represents 1 biological replicate. tRNA Ser UGG,G26A and 3 0 tRNA Ser UGG,G26A -GAL1pr reduce growth to 55% 6 9% and 71% 6 7%, respectively.
As shown in Fig. 3b and Supplementary Fig. 2, the wild-type strain displays a near-linear decrease with increasing mistranslation, consistent with our previous observations that increasing mistranslation frequency is negatively correlated with effects on growth for the same amino acid substitution (Berg, Zhu, et al. 2019). All the synthetic strains showed a graded response where increased mistranslation results in more severe loss of growth, but interestingly, the pattern of decreased growth in response to changing the frequency of mistranslation differed amongst the strains. This difference suggests that the genetic background influences the impact of mistranslation.
To look at the patterns in more detail, we plotted the normalized growth of each temperature-sensitive strain as a percentage of the growth of the wild-type strain (BY4742) for each of the three tRNA constructs that result in low, medium, and high mistranslation frequency (Supplementary Fig. 3). We note that in these plots, 100% indicates a lack of a negative synthetic effect, not the absence of an impact of mistranslation. Although the patterns appear to represent a continuum, we will focus the analysis on representative examples of four categories (Fig. 4). In the first category are strains where the synthetic interaction with tRNA Ser UGG,G26A increases proportionately with the frequency of mistranslation. In the representative examples ctf8-62 and cdc32-e72g (Fig. 4a), little synthetic interaction is seen at the lowest frequency of mistranslation. Many of the strains, including arc35-6 and las17-14 (Fig. 4b), are in the second category. These show a modest synthetic effect at a low mistranslation frequency, but have a strong negative synthetic interaction at both moderate and high frequencies of mistranslation. The third category, represented best by act1-4, has a nearly equivalent synthetic effect at all three levels of mistranslation (Fig. 4c). The last group includes cdc20-1 and is related to group 3 but appears relatively more impacted by the lowest frequency of mistranslation (Fig. 4d).
Many factors determine the impact of a mistranslating tRNA. Factors intrinsic to the tRNA include the anticodon sequence and its resulting amino acid substitution (Berg, Zhu, et al. 2021), the level of tRNA expression , and the presence of secondary mutations that alter the stability of the tRNA (Berg et al. 2017;Berg, Zhu, et al. 2019). Other factors are extrinsic to the tRNA. These include the number of competing tRNAs that buffer the mistranslating tRNA (Zimmerman et al. 2018), the environment in which cells expressing the mistranslating tRNA is found (Berg, Zhu, et al. 2021), and the genetic background of the organism.
Genetic background contributes to the impact of a mistranslating tRNA in numerous ways. In yeast (this work and Hoffman, Berg, et al. 2017) and Escherichia coli (Ruan et al. 2008), loss of genes that regulate proteotoxic stress increase the severity of mistranslating tRNA variants, most likely because the tRNA variants increase the load of mismade proteins. The impact of an extragenic mutation will depend on the extent to which it disrupts protein quality control or otherwise contributes to proteotoxic stress (Redler et al. 2016). Complexities arise since there are multiple quality control pathways that act independently but ultimately overlap to regulate proteostasis (reviewed in Chen et al. 2011). As such, impairing each pathway has the potential to show a different response to both changing levels of mistranslation and type of amino acid substitution. The nonlinear nature of the response likely arises because growth effects are not observed until a threshold of proteotoxicity is reached; different mutations will approach or exceed this threshold to different extents.
Other  and High: tRNA Ser UGG,G26A [tS(UGG)]g for temperature-sensitive strains expressing ctf8-162 or cdc33-e72g (a), arc35-6 or las17-14 (b), act1-4 (c), and cdc20-1 (d). The growth of the wild-type strain, 100%, is plotted as gray dots. Each point is the average of at least 3 biological replicates as in Fig. 3. function caused by mutation will be compounded by the reduced level of functional protein caused by mistranslation. The closer a genetic mutation brings the protein to the critical level of expression, the more impact the mistranslating tRNA will show when combined with that genetic mutant. Similar to the example above, until a critical threshold is exceeded, mistranslation may have little effect. At a more global level, genes that impact translation or mRNA processing can further limit expression of proteins already reduced by mistranslation.
Factors that alter the gene expression profile will also influence the impact of a mistranslating tRNA. Synonymous codon usage varies across the open reading frames in a genome (reviewed in Liu et al. 2021). Specific tRNA variants will mistranslate at a subset of these synonymous codons as determined by wobble rules and base modifications. tRNA variants only affect genes that are translated, and those that are more highly expressed will lead to greater proteotoxic stress when mistranslated. For yeast, we have shown that one such factor is the environment in which the cells are grown (Berg, Zhu, et al. 2021). As cell type determines gene expression profile in multicellular eukaryotes, different cell types are expected to be impacted differently by a tRNA variant. This argument is relevant to genetic background because mutations alter the internal and potentially external environment of the cell and often invoke a transcriptomic response (Hughes et al. 2000). The altered gene expression will in part determine the impact of a mistranslating tRNA. In addition, the impact of different anticodons decoding the same amino acid should affect the distribution of genetic interactions because anticodon sequence determines the extent and placement of mistranslated residues.
The genetic background could also directly or indirectly alter the expression of the mistranslating tRNA or that of competing endogenous tRNAs, altering the frequency of mistranslation. Genetic mutations that alter tRNA expression could occur in genes involved in regulating tRNA transcription (e.g. by altering chromatin structure) or in genes with roles in tRNA processing, modification, nuclear import, or stability. In our experiments, changes in expression of the mistranslating tRNA could occur if a secondary mutation occurs in a gene regulating transcription from the GAL promoter, though we note that none of the 16 temperature-sensitive alleles we analyzed were in genes with annotated GAL regulatory roles.

Conclusion
The impact of a mistranslating tRNA depends on a cell's genetic background and the expression of the variant as it relates to the frequency of mistranslation. Approximately 20% of humans contain a tRNA variant that potentially mistranslates (Berg, Giguere, et al. 2019). The variants are found within different tRNA isodecoders at different loci, which is highly relevant since tRNA genes in the human genome are expressed at different levels (Torres et al. 2019). Our findings demonstrate that genetic background must be taken into account when determining the contribution of tRNA variants to disease. Indeed, due to genetic or epigenetic differences, some individuals may be particularly sensitive to even low levels of mistranslation.

Data availability
Strains and plasmids are available upon request. The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and supplemental material. Supplementary File 1 contains all supplemental figures.
Supplementary File 2 contains raw and validated SGA data. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (Perez-Riverol et al. 2019) partner repository with the dataset identifiers PDX025934 and PXD032063.
Supplemental material is available at G3 online.