Autocrine signaling by receptor tyrosine kinases in urothelial carcinoma of the bladder

Comprehensive characterizations of bladder cancer (BCa) have established molecular phenotype classes with distinct alterations and survival trends. Extending these studies within the tyrosine kinase (TK) family to identify disease drivers could improve our use of TK inhibitors to treat specific patient groups or individuals. We examined the expression distribution of TKs as a class (n = 89) in The Cancer Genome Atlas (TCGA) muscle invasive BCa data set (n >400). Patient profiles of potentially oncogenic alterations (overexpression and/or amplification) clustered TKs into 3 groups; alterations of group 1 and 3 TKs were associated with significantly worse patient survival relative to those without alterations. Many TK pathways induce epithelial-to-mesenchymal transition (EMT), which promotes tumor invasiveness and metastasis. Overexpression and/or amplification among 9 EMT transcriptional activators occurred in 43% of TCGA cases. Co-occurring alterations of TKs and EMT transcriptional activators involved most group 1 TKs; 24% of these events were associated with significantly worse patient survival. Co-occurring alterations of receptor TKs and their cognate ligands occurred in 16% of TCGA cases and several BCa-derived cell lines. Suppression of GAS6, MST1 or CSF1, or their respective receptors (AXL, MST1R and CSF1R), in BCa cell lines was associated with decreased receptor activation, cell migration, cell proliferation and anchorage independent cell growth. These studies reveal the patterns and prevalence of potentially oncogenic TK pathway-related alterations in BCa and identify specific alterations associated with reduced BCa patient survival. Detection of these features in BCa patients could better inform TK inhibitor use and improve clinical outcomes.

oncogenic drivers among cabozantinib's TK targets. Our cited reference 23 also has been updated.
2) Why the authors excluded mutations in their first analysis of the TK genes?
Author Reply: Although we are very interested in the impact of TK mutations in bladder cancer, several prior studies have described the most prevalent mutations and their predicted and/or observed oncogenic impact (our cited references [14][15][16][17][18][19][20][21][22]. A comprehensive assessment of mutations across the entire TK family is also a big task for an individual lab (as opposed to a TCGA or similar large team collaboration), and because a substantial fraction of mutations would have unknown pathogenic impact and thus not bear comment, that portion of the work (often by a trainee) would be unpublishable. On the other hand, TK overexpression as a proxy for oncogenic impact in urothelial carcinoma of the bladder has not yet garnered as much attention from the research community. Because evidence of TK overexpression alone is not as strong as activating genomic alterations in guiding therapeutic strategy, we included additional proxies of pathway activation: coincident ligand overexpression (which might represent autocrine or local paracrine production, in view of how the TCGA samples were processed) and overexpression of EMT transcriptional activators that are known TK pathway effectors.
3) How the TK genes to be analyzed were selected in cell lines?

Author Reply:
We considered determining absolute mRNA copy number for all 52 group 1 and 2 TKs in the 15 BCa cell lines. There were challenges to this task, since we started with reliable primers for much smaller number of TKs. While expanding this primer set, we encountered TKs for which there were no published primers and/or readily available cell lines to serve as positive controls for primer testing. So, we assembled primer sets for group 1 and group 2 TKs in a 3:2 ratio (like the group of 52) to reach a subset that achieved the same correlation with TCGA patient molecular phenotype as the parent 52 TK set (p = 1.00 x10 -4 , q = 8.97 x 10 -5 , F2,405 > 9.42, R 2 > 0.044 for the 3-group comparison F test: LP vs. (LI + L) vs. (BS +N)). This goal was met with 19 group 1 TKs and 12 group 2 TKs, or ~60% of the parent set. We have added text to the Results for clarification (p. 14, lines 369-374). While this strategy limited the cell line data modestly, it was a reasonable and practical way to provide relevant TK coverage and enable TK-based molecular phenotyping of most of the cell lines into the classes described by Robertson et al. (cited reference 22). 4) How their assignment of cell lines to different BCa subtypes and their RNA data are correlated with those previously reported by Earl et al., 2015 (PMID: 25997541), including gene amplifications?

Author Reply:
This is an interesting question and we thank the reviewer for calling the work of Earl et al. to our attention. The short answer is that there is reasonable agreement between our findings, despite the very different methods that were used.
To recap our method: each cell line TK PCR expression profile (n = 15, 31 TKs) was correlated with each TCGA sample TK RNASeq profile (n = 408, 31 TKs), the correlation coefficients were determined and averaged within each molecular phenotype defined by Robertson et al. (22), and cell lines were assigned to a phenotype by best average coefficient that was distinguished from coefficients for other phenotypes by t-test. (1) "We applied the UBC molecular classifier based on gene expression defined by Sjodahl et al. [12] to identify lines most representative of the taxonomical groups proposed. Figure 6A shows that cell lines could be adscribed to the "Urobasal A", "Urobasal B", and "SCC-like" classifiers (Additional file 1: Table S12). The "Genomically Unstable" group was most commonly represented among the lines." (2) "Rebouissou et al. have recently reported on a 40-gene basal-like signature [25]. We have applied their 40-gene classifier to the cell line dataset and identify 4 major groups: lines with a predominant enrichment in the "Basal-like" signature; lines with a predominant enrichment in the "Non basal-like" signature; lines with enrichment of both signatures; and lines in which none of the signatures is enriched ( Figure 6B and Additional file 1: Table S13)." The most important difference between our respective methods was that we used only TK expression data for developing correlation coefficients, whereas the signature sets used by Earl et al. were developed without restricting genes to any particular functional class.
Eleven of the 15 cell lines that we used were also studied by Earl et al. We adapt below Table S13 from Earl et al. (left side) to compare the classifications of these 11 cell lines (right side): 5) It is not clear whether the differences in clinical outcome between altered and non-altered groups are obtained for each molecular subtype or globally. They should clarify this point, given the intrinsic differences in clinical outcome among groups Author Reply: Thank you for raising this point. We did not limit the Kaplan-Meier analyses to specific molecular phenotypes; in each case a global comparison was made among altered and unaltered gene groups. We have revised the Material and Methods text for clarification (p. 7, lines 157-160).
6) In addition to EMT drivers, they should analyze whether the altered groups also showed enrichment of EMT signatures using GSEA. In the same line, patients harboring overexpression also displayed increased tendency to develop distant metastasis, either at diagnose or in follow up?

Author Reply:
While we agree it would be interesting to determine whether any of the groups with lower survival also display enrichment of EMT signatures by GSEA compared to the unaltered groups, we feel that the amount of work involved is beyond the scope of this manuscript. We hope that an inspired reader of our work might consider the question and use the PLOS ONE online forum to post any findings of interest.
We believe that correlating overexpression with metastasis will be of limited value, since the greatest single factor determining overall survival is metastatic disease. This implies that any cohort that shows reduced overall survival, regardless of the feature(s) used to define that group, is very likely to also show a shorter average time to metastasis. 7) Number of cell lines were they confirm autocrine signaling is too short to generalize autocrine signaling in BC cell lines.

Author Reply:
We did not intend to imply that autocrine signaling is a dominant feature among BCa cell lines in general. We show immunoblot evidence of co-expression of cognate ligand and receptor protein pairs in 9 BCa cell lines (some of which show ligand/receptor pairs for more than one RTK pathway) and functional evidence of autocrine signaling in J82 and FL3. In our summary of these findings at the end of the Introduction (p. 5, lines 107-110) we originally stated: "These findings reveal the prevalence and patterns of autocrine RTK signaling in BCa…" where we used the word "prevalence" to mean "frequency"; to avoid confusion we have substituted the word "frequency" in the revised text. In the Discussion (p. 32, lines 853-854) we originally stated: "Direct evidence of autocrine signaling by AXL, CSF1R and MST1R pathways was found in several BCa-derived cell lines."; we have deleted the word "several" from this sentence in the revised manuscript to avoid any overstatement.
8) The metastatic capacity of J82 cells is increased upon serially passages in mice? Can it be blocked by cabozantinib?

Author Reply:
We did not observe that the metastatic capacity of J82 increased upon serial passages in mice and we do not know whether cabozantinib can block J82 metastasis in mice.

Reviewer #2:
The work is well done but the data is not presented in a manner by which the reader can follow. For example, in Figure 2, authors need to mark on the figure what is the difference between A, B, C, D. Western blots need to be redone -for example, that in Figure 4B upper three panels. In addition, the figure needs a loading control that is the same throughout the blot. Immunoblots in Figure 5D are cut too close to the band, and need to show empty space both before and after. Figure 5 -effect of adding ligands to the receptor or receptor phosphorylation not shown. Effect of genes overexpressing the receptor not shown.

Author Reply:
We thank the reviewer for their positive comments and careful review of our work.
We have revised Figure 2 to identify the groups being compared in each panel's symbol key. The comparisons in each panel are also defined in the figure's legend.
We have replaced the CSF1R immunoblot image for RT4, RT112, J82, UMUC5 and UMUC3 cell lines Figure 4B. We have also added to this panel immunoblot detection results for DDR1, DDR2 and PTK7 proteins across the 9 BCa cell lines, as these are understudied and potentially oncogenic RTKs in urothelial carcinoma of the bladder. The Results text (p. 17, lines 465-473) and Figure 4 legend (p. 18, lines 487-489) have been revised accordingly. The loading control for Figure 4B was omitted in error and this has been corrected in the revised figure (panel B, bottom) and legend.
We wish to emphasize that we did not use immunoblotting here to quantitate protein content, but to simply detect the presence or absence of a protein product under the conditions specified. As shown in Figure 5, we use 2-site electrochemiluminescent immunoassays to quantitatively compare protein-tophosphoprotein ratios under varying conditions. For example, we show that adding Gas6 to serum-deprived J82 cells significantly increased AXL phosphotyrosyl content ( Fig 5A, compare bars for Gas6 "-" and "+" at left). This increase in receptor phosphorylation is reversed by cabozantinib treatment (Fig 5A, 3rd and 4th bars from the left). The remaining panels in Figure 5 show that addition of Gas6, MST1 or CSF1 stimulated significant increases in the phosphorylation of critical downstream effectors Akt and Erk.
Regarding the cropping of bands in Figure 5D, we believe that no information has been obscured by our figure production, and we provide public access to all full-size, uncropped images that are used for figures as required per PLOS ONE policy.

Reviewer #3:
Summary of submission: This work utilizes bladder cancer (BCa) data from The Cancer Genome Atlas (TCGA) database to postulate the role of receptor tyrosine kinases (RTKs) in BCa. Hypotheses are then tested using multiple human BCa cell lines. The results suggest that autocrine signaling by RTKs contributes to the malignant phenotypes of some of the BCa cell lines and that autocrine RTK signaling may be a good target for therapeutic intervention in BCa.
The scope of this work is wide ranging and includes both in silico and wet lab approaches. Overall, the work is well-designed and -performed. Two major and a few minor issues are noted that do not markedly reduce the enthusiasm for this work but that should be addressed prior to publication.

Author Reply:
We thank the reviewer for their positive comments and thorough review of our work.
Major issues 1. The failure to directly measure RTK phosphorylation is noteworthy, particularly in those cases in which autocrine RTK signaling is inferred.

Author Reply:
Please note that direct quantitative measurement of Axl tyrosyl phosphorylation is shown in Fig. 5A and MERTK tyrosyl phosphorylation is shown more qualitatively by immunoblot in Fig. 8C. Any perceived lack of evidence of RTK activation is not due to lack of expertise or diligence (the Bottaro lab has published extensively on TK phosphorylation), but rather by choice. We believe that today, RTK activation through tyrosyl phosphorylation is well-enough established to be assumed from evidence of ligand-stimulated biological effects (by ligand addition or specific pathway inhibition as shown in Figs. 5 -8). Had we failed to find ligand-driven effects in experimental settings where they were expected, we certainly would have measured RTK activation as part of normal troubleshooting. Otherwise, we believe our effort is better used measuring pathway effects that are not yet established.
2. Similarly, the experiments do not appear to explicitly address the possibility that elevated RTK signaling could in some cases arise through RTK overexpression. This is of particular note given the attempt to identify the target(s) for carbozantinib.

Author Reply:
We do appreciate the importance of ligand-independent activation of overexpressed but otherwise unaltered RTKs. We characterize overexpression as a "potentially oncogenic alteration" throughout the manuscript. However, experimental proof of that event is not easy to obtain in cultured cells and much more difficult to prove in tissue samples or animal models. Consistent with our long-term goals, we focused effort on experiments that would convey reasonable evidence of pathway activation and were most likely to enable the development of reliable assays (direct or proxy) for tissue samples.
Minor issues 1. It is unclear how the cell lines, RTKs, and ligands analyzed in Figure 4b were chosen, particularly given the apparent occasional incongruence between transcription data shown in Figure 4a and the protein expression data shown in Figure 4b.

Author Reply:
The nine BCa cell lines analyzed in Figure 4B are a subset of the 15 cell lines analyzed in Figure 4A; the number was reduced from 15 to 9 for practicality, but the subset includes each of the molecular phenotype classifications shown in Table 1. The subset is weighted toward the more aggressive phenotypes (basal squamous and neuronal) out of an interest in the potential coincidence of multiple ligand/RTK pair co-occurrences that were suggested by our analysis of the TCGA dataset [22] ( Table 2A and Table S2A). Table S3A shows significant co-occurrence of group 1 TKs (i.e. TKs whose overexpression was associated with lower survival among the TCGA patients [22]), so the potential for multiple group 1 ligand/RTK pair co-occurrences in BCa cell lines was plausible.
Regarding the pathways chosen, as shown in Table 2A, the Gas6/AXL, CSF1/CSF1R and MST1/MST1R pathways have the highest Odds Ratio (Log2) values for cooccurrence, so these were prioritized in the immunoblot analysis in Figure 4B.
We regret that journal word limits that we encounter frequently (though not for PLOS ONE) have had a lasting negative impact on our tendency to provide detail. We have revised the Results text (p. 17, lines 465-469) to include our rationale for these choices.
2. Is there strong congruence between the cell line gene expression data reported in this work (Fig 3A and Fig 4A) and data available from the Broad CCLE? Similarly, if other biomarkers (from ref 21 cited in this work) characteristic of the N, BS, L, LI, and LP subtypes are applied to Broad CCLE data for the 15 BCa cell lines, are the subtype assignments the same as those reported in Table 1?

Author Reply:
We thank the reviewer for raising this interesting question. To address it, we downloaded the file "CCLE_expression.csv" from the DepMap portal (Source: Broad Institute; https://depmap.org/portal/download/). Per the DepMap web page, the file contains "RNAseq TPM gene expression data for just protein coding genes using RSEM. Data is Log2 transformed, using a pseudo-count of 1. Rows correspond to cell lines (Broad IDs) and columns correspond to genes (HGNC symbol and Entrez ID)". Of the 15 BCa-derived cell lines used in our study, SLT3, FL3, T24T and T24M2 were absent from the CCLE dataset. Of the TKs analyzed in Fig. 3A, the CCLE dataset omits DDR1 and LYN. Using all other TKs in Fig. 3A, we composed the scatter plot shown below, with mRNA expression as determined by RT-PCR (Y.H. Lee et al., x-axis) vs. mRNA expression as determined by RNAseq (CCLE, y-axis). Displaying both axes using log2 scaling allows resolution of all data points. The data are reasonably concordant, considering the differences in method, date of analysis and established variability among the same cell lines cultured in different laboratories at different times: Regarding whether the biomarkers used by Robertson et al. (our cited reference 22) characteristic of the N, BS, L, LI, and LP subtypes, if applied to Broad CCLE data for the 15 BCa cell lines, would result in the same subtype assignments as those reported in Table 1: while this is also an interesting question, it would take considerable time to address. Robertson used >3000 genes to generate their phenotypic groupings. Considering the work involved, we would be compelled to include all BCa cell lines for which data was available, an effort that we feel is beyond the scope of the current manuscript. Even if we were to consider it for an addendum, we would caution readers that real biological differences between cultured cell lines and tumor tissue samples affect the reliability and potential utility of this type of characterization, regardless of how definitive the approach might be.
We asked whether a limited set of differentially expressed TKs could be used to classify BCa cell lines into established molecular phenotypes and the results show that this can be done. The utility of this classification may be limited to TK-related studies, but in view of the importance of TKs as oncology targets and TK inhibitors as a drug class, it is still appreciable.
We appreciate the possible confusion here and we apologize for being overly brief. We have added text to the manuscript to provide greater detail regarding our goals and rationale and/or better clarify our findings (Introduction,Materials and Methods,p. 6,Results,p. 14,. Regarding the number of cabozantinib targets: 16 targets were known prior to our work and we identified DDR1, DDR2 and NTRK3 that had not been reported previously, bringing the total to 19. We have revised the Results text to clarify this point (pp. 15-16, lines 422-429).
The 14 RTK cabozantinib targets (ctRTKs) associated with significantly worse OS and DFS are group 1 and group 3 TKs only. Group 2 TKs are associated with improved survival relative to the other groups ( Fig. 2B and Table S2C) and three group 2 TKs are also cabozantinib targets (DDR1, MERTK and MST1R), so the combined groups 1-3 ctRTK set is not associated with significantly worse outcome. We have revised the Results text to clarify this point (p. 16, lines 444-448 ).
4. Finally, the manuscript appears to occasionally use TK instead of RTK.
Author Reply: This is intentional: RTKs are a subset of TKs, so only in TK sets that are exclusively RTK do we use that abbreviation. Since abbreviations used in the paper are listed separately in a footnote, we did not elaborate further on this in the text.