Genomic analysis of response to bacillus Calmette-Guérin (BCG) treatment in high-grade stage 1 bladder cancer patients.

BACKGROUND: Intravesical bacillus Calmette-Guérin (BCG) therapy is standard treatment for high-risk non-muscle invasive bladder cancer (NMIBC) but overall efficacy is low, and no reliable predictive biomarkers currently exist to refine patient selection. We performed genomic analysis on high-grade (HG) T1 NMIBCs to determine if response to therapy is predicted by certain mutational and/or expressional changes.
METHODS: Patients with HG T1 NMIBC treated with induction BCG were stratified by response into durable and non-durable responders. Baseline tumor samples were subjected to targeted DNA sequencing and whole-exome RNAseq. Genomic variants differing significantly between response groups were analyzed using Ingenuity Pathway Analysis (IPA) software. Variant selection was refined to target potential biomarker candidates for responsiveness to BCG.
RESULTS: Among 42 patients, the median follow-up was 51.7 months and 40.5% (n=17) were durable BCG responders. Deleterious mutations in the RNA sequence of JCHAIN, S100A7, CLEC2B, and ANXA10 were more common in non-durable responders. Mutations in MCL1 and MSH6 detected on targeted sequencing were more commonly found in durable responders. Of all deleterious DNA and RNA mutations identified, only MCL1 was significantly associated with longer recurrence free survival (RFS) (P=0.031).
CONCLUSIONS: Differences in the genomic profiles of HG T1 NMIBC tumors exist between those who show durable response to BCG and those who do not. Using pathway analysis, those differences imply upregulation of several interconnected inflammatory pathways among responders. Specific variants identified here, namely MCL1, are candidates for further study and, if clinically validated, may serve as useful biomarkers in the future. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

The annual incidence of urothelial carcinoma of the bladder (UCB) in the United States is estimated to reach more than 83,000 cases in 2020, making it the 6th most common non-cutaneous malignancy (1). Non-muscle invasive bladder cancer (NMIBC) consists of carcinoma in situ (CIS), Ta, and T1 tumors and accounts for approximately 75% of new diagnoses (2). Overall, NMIBC carries excellent long-term prognosis with a treatment approach focusing on local resection with or without intravesical therapy (3). Tumors with the highest risk for progression (15% at 5 years) to muscle invasive bladder cancer (MIBC) are the subset including high-grade (HG) T1 disease, having demonstrated invasive biology by disrupting the lamina propria. Standard front-line therapy for high-risk NMIBC is visually complete transurethral resection (TUR) plus intravesical bacillus Calmette-Guérin (BCG) administration, but recurrence is common at 40% within 1 year of treatment and there remains a 5-year cancer specific mortality rate of 25% (4,5).

The mechanism of action for BCG therapy is believed to involve immune activation and destruction of local cancer cells but the phenomenon is poorly understood (6). There are no predictors for treatment response outside of classical clinicopathologic features, and the high treatment failure rate highlights the need for additional measures. In the current era of next generation sequencing (NGS) for UCB, much research effort to date has been focused on exploring the prognostic and predictive role for MIBC genomics as it pertains to treatment response, with relatively less attention paid to NMIBC genomics and BCG responsiveness on a post hoc basis (7-10).

There remains ample opportunity and need for NMIBC biomarkers associated with BCG response to improve patient selection for intravesical therapy versus upfront surgery. Furthermore, comparisons between responders and non-durable responders might lead to a better mechanistic understanding of BCG therapy. We used NGS (RNAseq and targeted DNA sequencing) to analyze a cohort of BCG treated NMIBC patients with the highest risk of disease progression (HG T1) stratified by treatment response to compare differences in expressional and mutational profiles in relation to BCG responsiveness. We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/tau-21-158).

Methods

Patient characteristics

Recognizing the difficulties with clinical heterogeneity in NGS studies, we consciously elected to a priori restrict our analysis to patients with HG T1 UCB treated with complete induction BCG (≥5 of 6 weekly instillations) and compare those who responded to therapy versus those who did not. We retrospectively identified such patients at least 18 years or older at our institution between 2007 and 2019. Patients were excluded from further analysis if they had received any BCG within the previous 1 year. Response to BCG therapy was then assessed for each patient in order to create two discrete categories: BCG non-durable responders and BCG durable responders. Using date of first induction BCG instillation as the reference point, BCG non-durable responders were defined as patients with recurrence of UCB (any stage or grade) within 2 years, and BCG durable responders had no recurrences detected at any point during observable follow-up with a minimum disease-free interval of at least 2 years. BCG non-durable responders included those with a partial response, defined as any recurrence of lower grade or stage than the index lesion. Maintenance BCG was administered in all eligible patients in accordance with accepted guidelines (11). Formalin-fixed paraffin-embedded (FFPE) tissue blocks from the diagnostic TUR immediately preceding induction BCG were available for 45 patients. After microdissection and purification of DNA and RNA, a total of 42 samples were available for analysis.

DNA sequencing and variant calling

DNA samples were sequenced using the Trusight Tumor 170 gene panel (TST170, Illumina Inc., San Diego, CA, USA). Three hundred and fifty paired-end reads in fastq file format were first checked for quality using FastQC, a freely available software. The reads were then trimmed of adapter sequences and low-quality bases using Trimmomatic (12) and checked for quality again following trimming. The reads were aligned to the human genome 38 (hg38) reference genome using the BWA MEM aligner (13). The alignment files were then condensed by lane to a single alignment per sample. These alignments were assessed for quality and depth of coverage using Picard HSMetrics tool (http://broadinstitute.github.io/picard/). The average depth of coverage was 580X, with a median depth of 375X. After being assessed for quality, the alignments were passed to the Genome Analysis Toolkit (GATK) for variant calling (14). A panel of normal references was created using genomes listed in the 1000 Genomes Project. The variants were also quality checked, filtered, and annotated using the GATK package. After variants were called, they were compared between groups (durable responders vs. non-durable responders) using a contingency table and a fisher’s exact test. All variant P values listed are raw P values and any meeting a cutoff P value <0.05 were selected for downstream analysis.

RNA sequencing, gene expression, and variant calling

RNA samples were sequenced using Illumina’s TruSeq RNA whole exome panel. Sequenced reads were quality checked and reads were aligned to the human genome hg38 reference set using STAR aligner (15). The aligned files were used for further analysis of gene expression and variant calling. Counts of mapped reads for the genes were calculated by featureCounts in the R package SubRead (16). These read counts were used as input for DESeq2 to identify significant differentially expressed genes between durable responders and BCG non-durable responder samples. After adjusting for multiple correction using the Benjamini-Hochberg method, those with an adjusted P value of less than 0.05 and log-fold change difference greater than two were considered statistically significant. These significantly differentially expressed genes were inputted into the Ingenuity Pathway Analysis (IPA) software (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) to identify pathways that are potentially affected (17). Clinically annotated DNA sequencing and RNAseq data was uploaded to the National Center for Biotechnology Information Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Statistical analysis

Read counts were normalized and log transformed in R for input into classifyT1BC (Meeks classifier), consensusMIBC, and NMIBC classification (UROMOL classifier) to identify classification subtypes within our dataset (10,18,19). The classifyT1BC and consensusMIBC classifier tools were run using their respective R packages, while the NMIBC classification was run by uploading the log transformed normalized counts at http://nmibc-class.dk. For RNA-seq variant calling, aligned reads were merged per sample and duplicates marked with Picard tools. Sequence recalibration from GATK tools were used as input to Mutect2 for variant calling. Since matched normal samples were not available, a panel of normal references was created with MuTect2 to filter out potential germline variants. Raw variant calls for samples were then used for downstream analyses. After all variants were called, they were passed to bcftools to pull common variants from all samples. Only variants present in at least three samples were included. These common variants were then assessed for significance between groups (durable responders vs. non-durable responders) using fisher’s exact test in R. Variants were filtered to only include those with a raw P value of <0.05 (none had an adjusted P value of <0.05). Further annotation was done on the common variants using the Annotate Variation (Annovar) software to predict amino acid changes (20). The significant common variants were further filtered to only nonsynonymous variants. These filtered variants were further filtered to those that were annotated to genes that were found to be differentially expressed. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by institutional/regional/national ethics/committee/ethics board of the Levine Cancer Institute of Atrium Health (No. LCI-GU-SPEC-PHA-001) and individual consent for this retrospective analysis was waived.

Results

A total of 42 patients were included in the final analysis, of which 17 (40.5%) were BCG durable responders and 25 (59.5%) were BCG non-durable responders. With a median overall follow-up of 51.7 months [interquartile range (IQR), 29.3–69.7 months] baseline characteristics of the cohort () included a median age of 68.9 years old (IQR, 64.7–73.5 years old) comprised of 78.6% men (n=33), 78.6% current or former smokers (n=33), and 12.1% non-white ethnicities (n=5). A markedly higher proportion of females were observed among durable responders (n=7, 41.2%) than non-durable responders (n=2, 8.0%, P=0.019). Additionally, there were no active smokers among the durable responders, as compared to 24.0% (n=6, P=0.057) in the non-durable responder group.

Table 1

Baseline demographics and outcomes in the final cohort of 42 patients with HG T1 UCB

Variables	Total			Durable responders			Non-durable responders		P value
	N	%		N	%		N	%
Age in years, median (IQR)	68.9 (64.7–73.5)			66.1 (64.5–71.1)			69.4 (65.4–74.9)		0.303*
Sex
Male	33	78.6		10	58.8		23	92.0
Female	9	21.4		7	41.2		2	8.0	0.019^
BMI
<25	13	31.0		5	29.4		8	32.0
25–30	17	40.5		8	47.1		9	36.0
>30	12	28.6		4	23.5		8	32.0	0.859†
CCI
0–2	9	21.4		5	29.4		4	16.0
3–9	33	78.6		12	70.6		21	84.0	0.122*
Race/ethnicity
White/not Hispanic or Latino	28	66.7		14	82.4		14	56.0
White/not specified	9	21.4		0	0.0		9	36.0
African American	4	9.5		3	17.6		1	4.0
Asian	1	2.4		0	0.0		1	4.0	0.020‡
Smoking status
Current	6	14.3		0	0.0		6	24.0
Former	27	64.3		14	82.4		13	52.0
Never	9	21.4		3	17.6		6	24.0	0.057‡
Family history
Yes	25	59.5		9	52.9		16	64.0
No	17	40.5		8	47.1		9	36.0	0.534^
BCG responsive
Yes	17	40.5		17	100.0		0	0.0
No	25	59.5		0	0.0		25	100.0	N/A
Survival status
Alive	33	78.6		15	88.2		18	72.0
Dead	9	21.4		2	11.8		7	28.0	0.271^
Follow-up in months, median (IQR)	51.7 (29.3–69.7)			70.4 (56.3–79.7)			32.7 (25.3–51.6)		<0.001*

*, Welch two sample t-test; ^, Fisher’s exact test; †, Wilcoxon rank sum test; ‡, chi-squared test. HG, high-grade; UCB, urothelial carcinoma of the bladder; IQR, interquartile range; BMI, body mass index; CCI, Charlson comorbidity index; BCG, bacillus Calmette-Guérin.

Recurrences in the non-durable responders () occurred at a median of 7.7 months (IQR, 5.4–8.8 months) from initiation of induction BCG and were predominantly HG (n=22, 88.0%). Partial response was observed in 48.0% (n=12). There were 3 patients (12.0%) who progressed to muscle invasion or higher during follow-up, including a single patient who developed distant metastasis in the setting of negative surveillance cystoscopy.

Table 2

Characteristics of recurrences among BCG non-durable responders

Variables	Non-durable responders
	N	%
Grade at recurrence
LG	3	12.0
HG	22	88.0
Stage at recurrence
Ta	11	44.0
CIS	1	4.0
T1	10	40.0
T2 or greater	3	12.0
Partial response*
Yes	12	48.0
No	13	52.0
Time to recurrence months, median (IQR)	7.7 (5.4–8.8)

*, partial response defined as patients with T-stage < T1 at recurrence. BCG, bacillus Calmette-Guérin; LG, low grade; HG, high-grade; CIS, carcinoma in situ; IQR, interquartile range.

Gene expression

In total, 67 genes (Table S1) were found to have statistically significant (P<0.05) differential RNA expression on RNAseq when comparing durable responders to non-durable responders. The top genes found to be mutated across all samples, including both durable responders and non-durable responders, were BRCA2, ATM, and EP300 at 96.73%, 93.55%, and 90.32% respectively. The results of unsupervised hierarchical clustering analysis are shown in along with select clinicopathologic factors. Selected publicly available classification schemes for urothelial carcinoma were applied. To test whether these were involved in mediating inflammatory pathways, the differentially expressed genes were analyzed using the IPA software, using the expression levels for non-durable responders as the reference value. Three of the top pathways with statistically significant differences in gene expression (; ) were directly or indirectly associated with inflammation including macrophage migration inhibitory factor (MIF)-mediated glucocorticoid regulation (Figure S1), MIF regulation of innate immunity (Figure S2), and p38 mitogen activated protein family of kinases (MAPK) signaling (Figure S3). Furthermore, the predicted impact of changes in specific gene expression among durable responders was significant global up-regulation as well as activation of the key inflammatory mediator NF-κB when compared to non-durable responders ().

Figure 1

Gene expression heatmap of 67 statistically significant differentially expressed genes identified. Heatmap shows the count matrix data for 67 genes in all samples. Samples are colored based on their response group, shown in the legend. Additional clinical information is also provided. Publicly available classifiers from UroMol, consensusMIBC, and Meeks have been applied to the RNAseq data. MIBC, muscle invasive bladder cancer; BCG, bacillus Calmette-Guérin.

Table 3

Differentially expressed genes identified within the top 5 pathways as determined by IPA software

Pathway	Gene	Expr log ratio	Expr P value	Expected impact
MIF-mediated glucocorticoid regulation	PLA2G4D	–2.119	0.042	Up
	PLA2G4E	–2.338	0.027	Up
MIF regulation of innate immunity	PLA2G4D	–2.119	0.042	Up
	PLA2G4E	–2.338	0.027	Up
p38 MAPK signaling	PLA2G4D	–2.119	0.042	Up
	PLA2G4E	–2.338	0.027	Up

Expression values included in the table are representative of values among durable responders as compared to non-durable responders as baseline. The expected impact is generated by IPA software to indicate the predicted effect the change in expression in each particular gene would have on the associated pathway. IPA, Ingenuity Pathway Analysis; MIF, macrophage migration inhibitory factor; MAPK, mitogen activated protein family of kinases.

Figure 2

Simplified pathway analysis of genes with statistically significant differential expression between durable responders and non-durable responders. Pathway analysis was performed using IPA software which integrates gene expression results to visualize associated biological pathways. The top affected pathways are summarized here. IPA, Ingenuity Pathway Analysis.

DNA and RNA variant analysis

In addition to gene expression, we analyzed the RNAseq dataset for variants that may be potentially useful biomarker candidates. The pool of RNA variants (Figure S4) was first filtered to variants that appeared in at least three samples (n=126,033 variants). These were further filtered to variants that were annotated to proteins as nonsynonymous mutations (n=44,798 variants), of which 856 variants were significantly enriched (raw P value <0.05). Finally, of these significant nonsynonymous variants, there were six variant genes identified () from among the 228 genes previously identified as being differentially expressed based on BCG responsiveness. Notably, using these criteria, none of the differentially expressed genes from the pathway analysis were also RNA variants considered to be biomarker candidates. We found that JCHAIN variants were observed only in the non-durable responder samples (47.8% vs. 0%, P<0.001), while S100A7 (39.1% vs. 5.9%, P=0.026), CLEC2B (43.5% vs. 11.8%, P=0.041), and S100A7 (52.2% vs. 17.6%, P=0.046) were also significantly enriched in those patients (; Table S2).

Figure 3

Analysis of genetic variation between durable BCG responders and non-durable responders. (A) The heatmap depicts the two statistically significant DNA variants and the post-filtering RNAseq variants. Red indicates the absence of the variant while green indicates the variant is present. Grey squares indicate that either DNA or RNA was unable to be analyzed for that particular sample. The tables numerically outline the presence of each RNA variant (B) and DNA variant (C) as a ratio for the responders and non-durable responders. The associated raw P values for Fisher’s exact test are also listed. BCG, bacillus Calmette-Guérin.

Analysis of DNA variants was conducted with the TST170 gene panel using similar iterative exclusion methodology as above (Figure S4), with 150 DNA variants initially identified in three or more tumor samples, of which 45 were annotated to proteins as nonsynonymous mutations. There was statistically significant (P value <0.05) enrichment of two DNA variants that also passed nonsynonymous filtration (). Each of the identified DNA variants were observed in a higher proportion of durable responders (; Table S3): frameshift deletions in MCL1 (58.8% vs. 15%, P=0.007) and MSH6 (41.2% vs. 10.0%, P=0.034).

All DNA and RNA variants showing a statistically significant association with BCG response were then analyzed for impact in recurrence free survival (RFS) using the Kaplan-Meier survival method. A deleterious mutation in the MCL1 gene () was associated with significantly prolonged RFS as compared to wild type (2-year RFS 77.0% vs. 37.5%; log-rank test P=0.031). Though not reaching predefined statistical significance, mutation of the MSH6 gene () was also associated with longer 2-year RFS (77.8% vs. 42.9%; log-rank test P=0.073).

Figure 4

Kaplan-Meier survival plots. (A) Kaplan-Meier survival plot for RFS based on presence of deleterious mutations in MCL1 (log-rank test P=0.031) on targeted DNA sequencing panel TST170. (B) Kaplan-Meier survival plot for RFS based on presence of deleterious mutations in MSH6 (log-rank test P=0.073) on targeted DNA sequencing panel TST170. RFS, recurrence free survival.

Discussion

Predicting BCG responsiveness in NMIBC has remained an elusive goal, even as NGS techniques have vastly increased availability of genomic data. We sought to identify potential biomarkers for predicting response using a cohort of patients which were a priori selected to have similar baseline pathologic features (HG T1 tumors) who were treated with “adequate” BCG and displayed either durable complete response (median follow-up 70.4 months; IQR, 56.3–79.7) or disease relapse before 2 years of follow-up (median time to recurrence 7.7 months; IQR, 5.4–8.8). In this cohort, deleterious mutations in the MCL1 gene were associated with significantly longer RFS (median 9.3 months vs. not reached; P=0.031). Additional pathway analysis identified notable inflammatory pathways linked to durable responders with predicted increased pathway activation and activation of NF-κB.

MCL1 is an anti-apoptotic protein of the BCL2 family and its down-regulation has been linked to tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) sensitization (6,21). In prior work using immunohistochemical staining for BCL2, expression was shown to be inversely correlated with recurrence of NMIBC following BCG treatment with 47% prevalence at baseline in patients with tumor recurrences following BCG therapy as compared to 79% staining positive among non-recurrent cases (22). TRAIL expression has been shown to be up-regulated in granulocytes following induction of BCG and release of interferon signaling particles by cancer cells, and is therefore believed to be a key death ligand in the anti-neoplastic properties of BCG (23). The location of the frameshift deletion within MCL1 we detected falls in the first exon near the N-terminus within a region believed to be involved in protein regulation (24). We speculate, therefore, that sensitization to TRAIL, induced by a loss of function mutation in MCL1 that is seen more frequently in the durable responders and in association with longer RFS, provides a mechanism by which BCG responsiveness could be enhanced. The combination of a plausible mechanistic explanation for the association between MCL1 and BCG responsiveness makes this a promising target for further biomarker validation.

Deleterious mutations in the DNA mismatch repair protein MSH6 were also closely linked to longer RFS, though this did not meet our prespecified cutoff for statistical significance (P=0.073), a fact that may be explained by the small sample sizes. The MSH6 gene is well known for its role in Lynch syndrome whereby germline mutations lead to accumulation of errors in regions of DNA called microsatellites. The classic genitourinary manifestation of Lynch syndrome is that of upper tract urothelial carcinoma which is related to UCB through a shared histologic origin, but which has significant differences in genomic signatures (25). Beyond the association with Lynch syndrome, MSH6 has been linked specifically to UCB where methylation status of MSH6 has shown an association with progression and recurrence among patients treated with intravesical BCG (26). The deactivating, hypermethylated state of MSH6 corresponded to an increased time to progression and improved survival outcomes following BCG therapy in patients with HG NMIBC, and as such, should be considered for further investigation into its potential for a biomarker of response to therapy.

The pathways identified in our gene expression analysis show that inflammation and immune response are important triggers for BCG response, which we confirm in our formal pathway analysis using IPA software. MIF-mediated glucocorticoid regulation and MIF regulation of innate immunity, two of the top pathways from our analysis, are strongly involved in the immune system response. MIF is released during inflammatory reactions and has been shown to promote the expression of several pro-inflammatory cytokines (27). Upregulation of MIF-related pathways predicted by IPA software (Figures S1,S2) that we observed in our samples would be expected to result in a pro-inflammatory environment conducive to immune-mediated tumor destruction. Interestingly, there is a common linkage between all three of these pathways in the form of decreased expression in the NFKB1 gene leading to upregulation of each associated pathway (). Functionally, the NFKB1 gene is involved in regulating transcription for many pathways associated with inflammation, immunity, and apoptosis (28).

Other variants identified in our cohort, CLEC2B and S100A7, can also be linked to NF-κB activity but the correlation is less direct. CLEC2B is a member of the C-type lectin superfamily of proteins, carbohydrate binding proteins that can also self-bind and modulate many physiological activities (29). Increased TRAIL levels have been linked to C-type lectin-like molecule binding, but no specific work has shown direct correlation with expression of the CLEC2B gene. S100A7 is a less well described pro-survival protein present in other cancer models (30). It has also been linked to NF-κB signaling through suppression of the Yes-associated protein (YAP) oncogene cascade (31) and can be evaluated with immunohistochemical staining, though a previous study failed to find statistically significant association with BCG response (32). Though speculative until confirmed by further study, these observations lend additional support to support the hypothesis that alterations in NF-κB and its associated inflammatory regulatory networks may impact the response to BCG in patients with HG T1 UCB. ANXA10 is from the annexin family of proteins that are calcium dependent membrane binding proteins (33). While ANXA10’s specific function is still unknown, high expression levels have previously been linked to lower rates of NMIBC progression and improved survival outcomes (33,34).

The present study is limited by the small sample size available for genomic analysis, though we believe a strength of our data lies in the homogeneity of primary tumor samples. By narrowing our criteria to only HG T1 tumors that had clearly either responded to BCG, with a median disease-free follow-up of nearly 6 years, versus patients with relapse at a median of less than 8 months, we have intentionally maximized the difference in clinical outcomes. The use of a targeted 170 gene panel (Illumina’s TST170) is another potential limitation. While there may be additional interesting variants to identify through whole-exome sequencing, the targeted panel increases confidence as it provides higher sequencing coverage depth of the selected genes. Due to the retrospective nature of our dataset, we were not able to collect normal samples for matched analysis. Matching tumor/normal sample pairs are important in filtering out potential germline variants from the true somatic variants and represents an important future step in this work.

Conclusions

Intravesical BCG remains the preferred frontline treatment for most patients with HG T1 bladder cancer who have the highest risk of progression to muscle invasive disease. In this vulnerable patient population there is an ongoing need to identify reliable biomarkers that can predict patient responsiveness to BCG immunotherapy. We have identified several differentially mutated and expressed genes that vary based on responsiveness to therapy that can serve as candidate biomarkers to be further interrogated with validation studies, specifically MCL1 and MSH6. Furthermore, based on our work we hypothesize that NF-κB signaling may be a common thread linking these variants and pathways we analyzed. This serves as the foundation for further investigation into these alterations as possible rational biomarkers for high-risk UCB treatment.

The article’s supplementary files as