Low-Coverage Exome Sequencing Screen in Formalin-Fixed Paraffin-Embedded Tumors Reveals Evidence of Exposure to Carcinogenic Aristolochic Acid

The International Agency for Research on Cancer (IARC) classified aristolochic acid (AA) as a group 1 carcinogen (1). Exposure to AA, following intake of Aristolochia herbaceous plants as traditional medicines or due to consumption of bread from flour contaminated by Aristolochia seeds, can lead to AA nephropathy (AAN). AAN is a progressive tubulointerstitial nephropathy with high risk of developing upper tract urothelial carcinoma (UTUC; refs. 2–5). In addition, recent studies proposed AA as a factor contributing to the development of hepatocellular (6–8), renal cell (9, 10) and urinary bladder carcinomas (11), and intrahepatic cholangiocarcinoma (12). Given this growing spectrum of AA-associated tumor types, AA exposure detection methods for screening of disease-risk populations are of key importance.

Following metabolic activation of AA, aristolactam (AL)-DNA adducts accumulate in the proximal tubules of the renal cortex and are used as biomarker of exposure (4, 13, 14). AL-DNA adducts may persist for over 20 years after the exposure had ceased (15, 16) and can be measured by ³²P-postlabeling (14, 17) or by ultra-performance-liquid chromatography-electrospray ionization-multistage scan mass spectrometry (UPLC-ESI-MS/MSⁿ), both applicable to formalin-fixed paraffin-embedded (FFPE) tissues (16, 18, 19). However, the ³²P-postlabeling method lacks specificity, and access to the UPLC-ESI-MS/MSⁿ methodology and its optimization for biomaterial of low quantity are limiting factors.

DNA sequencing established a characteristic AA mutational signature marked by accumulation of A>T transversions within the 5′-Pyr-A-Pur-3′ sequence context (enriched for 5′-CpApG-3′), preferentially located on the nontranscribed strand (8–10, 20, 21). In cancers not associated with AA, such A>T transversions are infrequent (22, 23).

We exploited the unique features of the AA mutational signature to devise a sensitive method for AA exposure detection, based on low-coverage whole-exome sequencing (LC-WES, at approximately 10× in contrast with the conventional 100× coverage), optimized for analysis of tumor-specific DNA of limited quantity and integrity extracted from archived FFPE tissues. The studied urothelial tumor samples originated from a well-characterized population residing in the endemic nephropathy (EN) regions of Croatia and Bosnia and Herzegovina (13), with EN being thus far the only recognized environmental form of AAN (4, 24). For the first time, we report in the urothelial tumors of EN patients the genome-wide signatures of AA, age, and APOBEC cytidine deaminase activity, thereby extending previous mutational analyses of this population based solely on the mutations of the TP53 tumor-suppressor gene (4, 13, 25). In addition, we demonstrate the ability of LC-WES to elucidate the impact of the AA mutation spectra on key homeostatic biologic pathways and to reveal possible mechanisms of tumor dissemination along the urinary tract.

Exposure to AA was investigated in 15 patients with urothelial tumors, diagnosed with EN following established criteria (13, 26). As controls, UTUC samples were obtained from 4 patients from a metropolitan area of the United States, unlikely exposed to AA. All specimens were FFPE-converted in the histopathologic laboratories of the participating centers. The involved anatomical sites were renal pelvis, ureter, and bladder (ICD-10 codes C65, C66, and C67, respectively). Clinicopathologic features and Aristolochia exposure history are listed in Supplementary Table S1. The study protocols included patients' informed consent and were approved by the IARC Ethics Committee and the Institutional Review Boards of the participating institutions.

Hematoxylin and eosin preparations of the paraffin block sections were used to identify tumor tissue free of necrotic areas. The tumor cell areas were measured by ImageJ software (27). Ten micromolar sections, cut with the Leica RM 2145 microtome (Leica Microsystems), were used to macrodissect the tumor-enriched areas and isolate genomic DNA yielding 1 to 2 μg (5 to 10 ng/mm²) per sample. Before DNA isolation, slides were deparaffinized for 5 minutes in 100% xylene, kept for 5 minutes in absolute ethanol, 5 minutes in 85% ethanol, 5 minutes in 75% ethanol and stored in milliQ water. DNA isolation was done using the QIAamp DNA FFPE Tissue Kit (Qiagen) following the manufacturer's protocol. DNA yields and concentrations were measured using the Picogreen assay (LifeTechnologies) and Fluoroskan Ascent FL microplate fluorometer (Thermo Fisher Scientific). DNA purity was evaluated by the NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific), and DNA integrity assessed by 0.8% agarose gel electrophoresis.

DNA was isolated from the renal cortex and tumor tissues by standard phenol-chloroform extraction techniques. The level of AL-DNA adducts in the renal cortex DNA (10–20 μg) was determined using ³²P-postlabeling PAGE, as previously described (13). The TP53-specific mutations were identified using the AmpliChip p53 Research Test (Roche Molecular Diagnostics), sensitively detecting all single base-pair substitutions and single-base deletions (13).

Two hundred and fifty (250) ng of genomic DNA was sheared by the adaptive focused acoustics method (Covaris, Inc.) to obtain approximately 300 bp fragments, with water temperature of 4°C, one cycle at 175 Watt peak power, duty factor 10 and 200 cycles per burst. Resulting fragment size was assessed using the 2100 Bioanalyzer and High Sensitivity DNA Kit (Agilent Technologies). The sheared DNA was converted into libraries using the Kapa LTP Library Preparation Kit (Kapa Biosystems). Briefly, the fragmented DNA was subjected to end repair reaction followed by poly-A-tailing and adapter ligation, excess adapters removed by Agencourt AMPure XP beads (Beckman Coulter). Eight cycles of PCR were performed to amplify the libraries with correct adapters on both ends. Four libraries (250 ng each) were pooled per exome capture with the Nimblegen SeqCap EZ Exome reagent. Exome-enriched mixes were PCR-amplified in 10 cycles, post-enrichment libraries pooled in 420 μL of water to a final concentration of 6 pmol/L. This volume was divided and loaded in two lanes of the rapid run mode flow cell for cluster generation and sequencing on the HiSeq2500, in a paired-end 50 bp cycle run. Multiplexing 16 samples per run resulted in the target coverage of approximately 10×.

Four additional EN UTUC samples and two UTUC samples from the metropolitan United States were analyzed in a validation assay using the SOLiD 5500XL sequencer (Life Technologies). See Supplementary Methods for details.

Raw HiSeq2500 sequencing data were deposited to the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) repository (ID SRP042035) to become available from the NCBI's dbGaP database. The annotated list of single-base substitutions (HiSeq2500 data) is provided in Supplementary Table S2.

FastQ reads were aligned to the human genome (hg19) using Burrows–Wheeler Aligner. Realignment and base quality score recalibration was done by the Genome Analysis Toolkit (GATK) and the duplicate-read removal by Picard. GATK HaplotypeCaller was used to call variants subsequently annotated on the RefSeq Gene transcript contents by ANNOVAR (28). Polymorphisms present in normal population and removed from our data originated from these collections: 1,000 genomes (1,000 g, http://www.1000genomes.org/), Exome Sequencing Project (ESP, http://exome.gs.washington.edu/) and the SNP database build 137 (dbSNP, http://www.ncbi.nlm.nih.gov/SNP/). We removed variants with frequency above 0.1% in either the 1,000 g or ESP databases, or annotated in the dbSNP database, or present in a custom germline variant catalog built from 560 cases from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/). Variants mapping to repetitive sequences contained in the genomic segmental duplication database (29) alongside variants with ≥90% homology with multiple regions were excluded. R functions were developed to compute the mutation type distributions and strand bias. The strand bias significance was determined by the Pearson χ² test. These tertiary analysis parameters were computed in two separate coverage ranges, ≥3× with no defined maximum, and between 3 and 9× to emulate ultra-low coverage.

Nonnegative matrix factorization (NMF) decomposes mutational patterns based on factorization of one matrix (n × m) in two matrices W (n × r) and H (r × m) with the constraint that all three matrices must be composed of nonnegative elements (30). The r is the rank of factors to be extracted from the input matrix, corresponding to the number of signatures. The input matrix contained one column per patient (only HiSeq2500 data considered) and in rows the frequency of mutations types in 96 possible two-base sequence contexts. The R package NMF (31) was used to extract mutational signatures. The correlation between the extracted signatures and previously published ones (7, 21, 22, 32) and/or available in COSMIC (23) was computed as the inner product of the two signatures (vectors) divided by the product of their norms.

To examine the biologic impact of the gene mutants in the AA signature-positive samples, analysis was performed using the DAVID tool (33), with two input gene lists: (i) genes harboring nonsynonymous (missense, stop-gain, or stop-loss) SBS and (ii) genes nonsynonymously mutated in the EN dataset and in AA signature-positive samples from at least one of the two published datasets on UTUC in Taiwanese patients (8, 21). The list was further narrowed by classifying the mutated genes as established oncogenes or tumor suppressors listed by the Gene Set Enrichment Analysis (GSEA) database (34) and/or a cancer driver genes defined by recent seminal studies (35–39).

We applied HiSeq2500 LC-WES to genomic DNAs isolated from FFPE urothelial tumors from 11 EN cases, two of whom had concurrent UTUC and bladder carcinoma, and from two U.S. patients providing non-EN control samples (13 patients and 15 tumors in total, see Supplementary Table S1). Features of AA signature had been described earlier, as follows: mutational load of ≥40 SBS or ≥10 A>T in exonic positions, high proportion of A>T (>35% of all SBS types or as the predominant type) with a strand bias of ≥1.25, and ≥33% enrichment of A>T in the 5′-(C/T)pApG-3′ sequence context (8, 21). We used analogous criteria for the AA signature (≥50 SBS/sample of which ≥15% are A>T SBS, of which ≥20% are in the 5′-CpApG-3′ context), applying more stringent statistical analysis of the strand bias ratio combined with a cutoff of ≥1.5 (9, 11). Under these criteria the AA signature was readily observed in 10 of the 13 analyzed EN tumor samples, with 33% to 77% of A>T transversions per sample (Fig. 1), a nontranscribed strand bias of 2.0 to 3.3 and the 5′-C_G-3′ context enrichment above 19% (mean 24.6%, SD = 4.9, range of 19.1%–27.4%; Table 1). In contrast, A>T mutations and their enrichment in the 5′-CpApG-3′context are generally low in cancers of non-AA etiology, based on our analysis of 7,160 tumors of 52 cancer types in the COSMIC database (average 5.8% A>T, range 0%–12.1%, of which 10% are in the 5′-CpApG-3′ context). Similarly, the average percentage of A>T in the 5′-CpApG-3′ context in TCGA urothelial carcinoma data (only bladder data available) is 10.8% (0%–50%), whereas the mean percentage of all A>T mutations is low (average of 3.9%, range 0.8%–8.3%).

A weaker signature marked by 18.7% A>T, strand bias of 2.1 and the 5′-CpApG-3′ context proportion of 12.5% was observed in the bladder tumor sample (EN-01-B) of a patient with a concurrent AA signature-positive UTUC (EN-01-RP, see Table 1 and Supplementary Table S1). Two EN samples (EN-06 and EN-07) and the two non-EN controls were found negative for the AA signature, with A>T transversions present at 4% to 8%. In the case of EN-07 (bladder carcinoma with no history of UTUC), despite the presence of AL-DNA adducts in the patient's renal cortex, the mutation profile (Fig. 1) suggested AA-unrelated etiology. Among the AA signature-positive samples, we detected an average of 1,142 (range, 349–2,707) mutations per tumor [∼18 SBS/sample per exome megabase (Mb)] whereas the mutation rate in the control and negative samples (including the weaker AA-signature bladder cancer) was on average 357 (range, 258–440) mutations per sample (∼6 per sample/exome Mb). As shown in Table 1, the predominant A>T transversions substantially contributed to the high SBS counts.

Thus, LC-WES analysis of the EN UTUC generates results consistent with previous reports on the highly mutagenic potential of AA (8, 21, 40), and our results justify the use of exome sequencing for reliable detection of exposure to AA in archived FFPE material.

We next investigated whether the AA signature can be identified at ultra-low coverage. Upon considering 3 to 9 non-duplicate per-base reads, mutation counts in the AA-associated samples decreased to 233 per tumor on average (∼4 per sample/exome Mb) and to average 67 per tumor (1/sample/exome Mb) in the negative samples and the weakly positive bladder tumor (EN-01-B). The 10 tumors shown in Fig. 1 (two top rows) still exhibited the AA signature at ultra-low coverage (Supplementary Fig. S1), with the strand bias ratios between 1.7 and 4.7, and retained prominent enrichment of the 5′-CpApG-3′context (>25%). Thus, the specific and unique features of the AA signature can be reliably detected in FFPE tumor samples by superficial coverage sequencing.

These results open an attractive opportunity for retrospective analyses of archived pathologic specimens from the regions of AA exposure risk. In comparison with the ³²P-postlabeling and mass spectrometry adduct detection techniques, the LC-WES approach is based on a commodity technology that generates genome-wide information. LC-WES is also very sensitive, using low input DNA amounts (250 ng compared with 5–10 μg required for adduct analysis). Finally, it can indicate exposure to AA when neither AL-DNA adducts nor mutations in TP53 are detected, as we demonstrate for the AA signature-positive cases EN-01, EN-03, EN-04, and EN-11 (Fig. 1; Supplementary Table S1).

NMF extracts individual mutational signatures from complex alteration patterns observed in primary tumors, reflecting thus the specific effects of etiologic factors (7, 9, 22, 41). NMF was used to describe the AA signature in human UTUC, bladder, liver, and renal carcinomas (7, 9, 11) and in experimental in vitro system designed to model mutational signatures of carcinogens (32). Here, in the EN urothelial tumors, the NMF approach identified three distinct signatures, the AA-specific signature (Signature 22; ref. 42), the signature related to age (C:G>T:A in the 5′-XpCpG-3′ context, Signature 1; ref. 22), and the Signatures 2 and/or 13 associated with the cytidine deaminase activity of the APOBEC enzymes (Fig. 2; ref. 22). All three signatures are currently listed in the COSMIC database (23). Furthermore, NMF aided in classifying the EN-01-B bladder tumor as positive due to nonnegligible sample contribution to the AA signature (18%), in contrast with the negative samples (EN-06-RP and EN-07-B, with contributions of 0% for both) and non-EN controls (each 0% contribution, see Table 1 and Fig. 2B). The identified EN UTUC AA signature correlated highly (>90%) with the COSMIC Signature 22 (23), derived from AA-associated primary UTUC tumors from Taiwanese patients (8, 21), and with the AA signature modeled in vitro (Fig. 2C; ref. 32). The other EN tumor signatures matched their COSMIC counterparts with 72% similarity (age) and 70% and 64% similarity (APOBEC, Signature 2 and 13, respectively; Fig. 2C).

To validate the LC-WES performance using another sequencing chemistry and platform, we analyzed four additional EN UTUCs positive for AL-DNA adducts and p53 A>T mutations (EN-12, EN-13, EN-14, and EN-15), and two control UTUCs from U.S. patients (Non–EN-03 and Non–EN-04), on the SOLiD 5500xl sequencer. At the average 14.5× coverage, we observed the AA signature in all EN samples, although in samples EN-13-RP and EN-14-RP, the A>T transversion was the second most abundant mutation type following C>T (see Supplementary Fig. S2A). The signature remained detectable at ultra-low coverage (∼4.6×), when considering only the 3 to 9 read interval (Supplementary Fig. S2B).

In the A>T enriched samples, A>T transversions were randomly distributed along the sequenced regions, with linear correlation between A>T SBS counts and chromosome size (R² = 0.9; Fig. 3A). Similar correlation was maintained in the minimum coverage interval of 3 to 9× (data not shown). This result was confirmed by the analysis of the Taiwan UTUCs (8, 21) in which a similar, although less linear trend was observed (R² = 0.61–0.64). These findings suggest a stochastic A>T mutation distribution within the gene/transcription units represented by the exome.

Despite this apparently random pattern, we identified 83 cancer driver genes carrying protein sequence-altering A>T SBS, that were recurrently mutated across the three datasets of the AA signature-positive tumor samples (this study, n = 10, and the two previously reported Taiwanese sets of n = 18 (21) and n = 9 (8). These findings are summarized in Fig. 3B and in Supplementary Table S4. The recurrently mutated genes included numerous known drivers and chromatin-associated factors such as TP53, ARID1B, ATRX, CREBBP, CHD2, CHD5, CHD8, FAT1, KDM6A, MLL2 (KMT2D), SETBP1, TRRAP. TP53 was the most frequently mutated gene [17/37 (46%) samples] with all its mutations being A>T transversions. Fifteen samples exhibited mutations in the histone methyl-transferase KMT2D (MLL2), with varying SBS types, suggesting that secondary mutation processes possibly linked due to high mutational loads and increased genomic instability. Further systematic investigations should be undertaken to establish possible recurrent alterations in particular genes and pathways in UTUC across studies of different populations/geographical areas. For instance, data in Fig. 3 and in Supplementary Table S3 indicate that TP53, CREBBP, and LRRK2 are mutated mostly in the Taiwanese samples whereas mutations in the AHNAK, ATRX, SMCHD1, and XIRP2 genes are enriched in the EN UTUC samples. Other factors contributing to these differences merit further investigations, including varying modes of AA exposure (low-dose chronic intake in the EN regions as compared with higher-dose, (sub)acute exposures resulting from the use of traditional herbal medicines in Asia) and disease susceptibility due to the patients' genetic background.

Using NIH DAVID, we performed Gene Ontology (GO) and KEGG pathway analyses of the genes harboring nonsynonymous A>T mutations in the AA-signature–positive samples analyzed by HiSeq2500 (n = 10). We identified gene targets from the functional classes of cell adhesion, cell-matrix contact, cell migration, cell-cycle, cell signaling (MAPKKK/RAS and PI3K cascades, mTOR pathway), pathways of WNT, insulin and ERBB signaling, nucleotide excision repair, and the DNA-dependent ATPase and helicase activity, chromatin modification and histone binding related to gene-expression regulation, with dozens to hundreds of mutated genes per category (Supplementary Table S3). This observation suggests massive deregulation and/or destabilization of key homeostatic pathways by the high A>T mutation loads.

Next, for the 83 recurrently mutated cancer genes (Fig. 3B; Supplementary Table S4), we observed enrichment of GO and KEGG categories related to regulation of transcription, chromatin/histone modification, and categories of DNA damage response and DNA repair (Supplementary Table S5). These included numerous previously established cancer driver genes (TP53, AHNAK, ARID1B, ATRX, BLM, CHD2, CHD5, CHD8, CHD9, CHEK2, CLTC, ERBB4, FN1, HUWE1 IARS2, KALRN, LRRK2, MLL2, NEB, RXRA, SMCHD1, SPEG, STAG2, SYNE1, TRIO; refs. 35–39). Thus, the LC-WES analysis of AA-exposed urothelial tumors and associated data mining can reveal biologic information contents, particularly upon meta-analysis with data from different populations characterized by identical etiology and tumor types.

Two EN patients had synchronous urothelial tumors in distinct anatomical sites (renal pelvis and bladder, samples EN-01-RP and EN-01-B; and renal pelvis and ureter, samples EN-02-RP and EN-02-U). By using LC-WES, we investigated the common genetic origins of these synchronous tumor pairs. In patient EN-01, the overlapping SBS were enriched for C>T mutations (42%) followed by A>G (20%), and only 7.7% of the overlapping SBS were A>T transversions affecting the coding sequence of mere 3 non-cancer genes (VWA3B, KDM3B, and ACIN1). However, the A>T SBS were enriched among the mutations unique to the renal pelvis and to the bladder tumor (77% and 28%, respectively, Supplementary Fig. S3A), suggestive of a common precursor carrying mainly non–A>T driver mutations, giving rise to two tumor progenies subsequently accumulating distinct patterns of A>T alterations in either anatomical site. The distinct AA signature in the bladder tumor is in keeping with a recent study of Asian bladder cancer patients in whose tumors the AA signature manifested without the involvement of upper tract or a history of renal disease (11). In contrast, the tumors in the renal pelvis and ureter of patient EN-02 shared the majority of mutations contributing to a prominent AA signature, suggesting a common precursor carrying mostly A>T alterations (Supplementary Fig. S3B). This genetic relationship between same-patient tumors suggests cell seeding along the tract as the basis for tumor dissemination. However, further investigations of a larger multiple-tumor case series and with the use of deep sequencing is needed to further elucidate the exact mechanisms of multifocal and recurrent tumorigenesis in the urinary tract of AA-exposed patients.

In summary, we report successful detection of the genome-wide AA signature in urothelial tumors of EN patients, using archived FFPE specimens and a customized low-coverage exome sequencing. The described technique is a cost-effective screening tool potentially applicable to molecular epidemiology studies aiming at identifying cancers associated with AA exposure. This ability of the LC-WES and its applicability to archived biomaterial may be exploited in future systematic studies on AAN and associated cancers, in support of established or future disease prevention programs.

S.F. Shariat has received honoraria from speakers bureau from Astellas, Takeda, Ipsen, Janssen, Sanofi, Olympus, Wolff, Pierre Fabre, and Sanochenia, has ownership interest (including patents) in prostate and bladder cancer biomarkers, and is a consultant/advisory board member for Astellas, Ipsen, Sanochemia, Olympus, Wolff, and Janssen. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.P. Grollman, B. Jelaković, J. Zavadil

Development of methodology: X. Castells, M. Ardin, G. Durand, N. Forey, F. Le Calvez-Kelm, I. Dolgalev, J. McKay, A. Heguy, J. Zavadil

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Castells, S. Karanović, K. Tomić, E. Xylinas, G. Durand, N. Forey, F. Le Calvez-Kelm, D. Dittrich, S.F. Shariat, V.S. Sidorenko, A. Heguy, K.G. Dickman, A.P. Grollman, B. Jelaković, M. Mišić

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): X. Castells, M. Ardin, K. Tomić, C. Voegele, I. Dolgalev, J. McKay, S.F. Shariat, M. Olivier, J. Zavadil

Writing, review, and/or revision of the manuscript: X. Castells, S. Karanović, M. Ardin, K. Tomić, E. Xylinas, F. Le Calvez-Kelm, C. Voegele, S.F. Shariat, K.G. Dickman, M. Olivier, A.P. Grollman, B. Jelaković, J. Zavadil

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): X. Castells, S. Karanović, K. Tomić, S. Villar, K. Karlović, D. Dittrich, I. Dolgalev, A. Fernandes, K.G. Dickman, M. Ardin, M. Mišić

Study supervision: B. Jelaković, J. Zavadil

The authors thank the Genetic Platform, IARC, and the Genome Technology Center, NYU Langone Medical Center for expert assistance with exome sequencing; the authors thank Christine Carreira, Dr. Behnoush Abedi-Ardekani, and Dr. Elisabetta Kuhn for assistance with the biospecimen evaluation and processing. The authors thank the staff of the General Hospital Slavonski Brod, the University Hospital Centre Zagreb and of the Weill Medical College of Cornell University, New York, for expert assistance and support. The authors are grateful to Drs. Monica Hollstein and Michael Korenjak for critical comments on the article.

IARC Regular Budget; grant P01 ES04068 (to A.P. Grollman) from the NIH/NIEHS; R03-TW007042 grant (to A.P. Grollman) from the NIH Fogarty International Center; grant (to A.P. Grollman) from Henry and Marsha Laufer; grant no. 108-0000000-0329 (to B. Jelaković) from the Croatian Ministry of Science, and grant 04/38 (to B. Jelaković) from the Croatian Foundation for Science; and grant NIH/NCI P30 CA016087-33, which partially supports the NYU Genome Technology Center.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.