Whole-exome sequencing of 79 xenografts as a potential approach for the identification of genetic variants associated with sensitivity to cytotoxic anticancer drugs

Chemotherapy response remains unpredictable in most patients with cancer. In this study, we performed whole-exome sequencing of 79 cancer xenografts derived from human cancer tissues to identify genetic predictors of chemosensitivity to nine cytotoxic anticancer drugs. Xenografts were harvested from 12 organs with cancer and implanted into nude mice. The mice were exposed to one of nine cytotoxic anticancer drugs (5-fluorouracil, nimustine, adriamycin, cyclophosphamide, cisplatin, mitomycin C, methotrexate, vincristine, and vinblastine) to assess the correlation between chemosensitivity response and variant allele frequency. We found 162 candidate variants that were possibly associated with chemosensitivity to one or more of the nine anticancer drugs (P < 0.01). In a subgroup analysis of breast and gastric cancer xenografts, 78 and 67 variants, respectively, were possibly associated with chemosensitivity. This approach may help to contribute to the development of personalized treatments that may allow for the prescription of optimal chemotherapy regimens among patients with cancer.


Introduction
Cancer is a global health concern, with approximately 18.1 million new cases and 9.6 million deaths in 2018 [1]. Currently, most cancers are treated by surgery, radiation therapy, and/or chemotherapy [2,3]. Chemotherapy remains a gold standard for the treatment of blood cancers, such as leukemia and lymphoma [4][5][6], unresectable or metastatic cancers [7,8], and solid tumors, such as lung, breast and colorectal cancers [9][10][11]. Despite an improved understanding of cancer biology and the development of molecular targeted therapy and immunotherapy for select patients [12,13], chemotherapy still plays a primary role in cancer treatment regimes. Indeed, chemotherapy regimens have improved considerably, now taking into consideration the organ of origin, histological appearance, and stage of progression. Yet, these improvements aside, chemotherapeutic efficacy still varies between individuals [14] and is often complicated by toxic reactions, including nausea, tiredness, diarrhea, and hair loss [14,15], causing physical and mental distress and decreased patient quality of life. As such, it is becoming increasingly important to identify effective treatments with fewer toxic side effects as a first-line therapy for each patient. Several recent studies have sought to establish diagnostic methods for predicting chemosensitivity to cytotoxic anticancer drugs before treatment is undertaken. However, clinically useful genetic markers have yet to be developed [16][17][18][19][20]. Patient-derived xenograft (PDX) models have been established for many types of tumors and have emerged as powerful tools for predicting drug efficacy and for understanding tumor characteristics. With PDX models, fresh human tissue is directly implanted into immunocompromised mice. These models retain the heterogeneity of the original patient tumors and thus allow for tests, predominantly to examine the efficiency of anticancer drugs [21].
Next-generation sequencing technologies have also been developed in recent years, exposing tumor genomic profiles and facilitating the detection of low frequency variants and other genetic mutations that could not otherwise be uncovered by conventional methods [22,23]. Indeed, several studies have reported associations between such genetic mutations in tumors and clinical outcomes [24][25][26]. Guided by these reports, we hypothesized that genetic variants within tumors, including low frequency, rare variants, may underpin patient responses to cytotoxic anticancer drugs, such as chemosensitivity and chemoresistance. To this end, we performed whole-exome sequencing of DNA samples taken from 79 human cancer xenografts prepared from 12 different organs. These xenografts were implanted in mice and treated with one of nine cytotoxic anticancer drugs. We assessed correlations between the chemosensitivities of the xenografts to nine anticancer drugs and the variant allele frequencies (VAFs) as a potential approach to identify variants that may be predictive of drug response.

Animals and tumor xenograft model
Previously, a total of 79 human tumor tissues were obtained aseptically during surgery or at autopsy across 13 hospitals in Japan [27]. The samples included 12 breast cancers, 12 gastric cancers, 10 neuroblastomas, 10 non-small-cell lung cancers, 7 gliomas, 6 pancreatic cancers, 5 colon cancers, 5 choriocarcinomas, 4 small-cell lung cancers, 4 hematopoietic cancers, 3 ovarian cancers, and 1 osteosarcoma. These xenografts were separately transplanted into athymic BALB/c-nu/nu mice (Clea Japan, Inc, Tokyo, Japan) and maintained by serial subcutaneous transplantation of 2×2×2 mm fragments into the flank once a month, as described previously [27]. Microbiological monitoring of the tumor-bearing nude mice was performed for bacteria (e.g., Pasteurella pneumotropica and Mycoplasma pulmonis), viruses (e.g., mouse adenovirus and mouse hepatitis virus), and parasites (e.g., Giardia muris and Spironucleus muris) by culture, serological, or microscopic examinations [27]. Furthermore, histological examination and isozyme testing were carried out to assess for the risk of cross-contamination among the tumor lines, or cross-contamination between human tumor xenografts and a mouse tumor appearing at the inoculation site of the xenograft during passaging [27]. Tumor-bearing mice were euthanized with deep anesthesia followed by cervical dislocation. To minimize discomfort, euthanasia was performed quickly. Tumors were excised from the euthanized mice, and a piece of the tumor tissue was implanted into another mouse using a transplantation needle. All handling of mice was carried out in a gentle to minimize animal suffering and distress. The general conditions of the mice such as appetite and respiratory conditions were monitored patients for the use of their tumor was not obligated at the time. Furthermore, these xenografts are publicly available resources, and chemosensitivity data of them were published in 1996 [27]. Therefore approval from ethics committee is not necessary for this study.

Sample preparation and whole-exome sequencing
Tumor genomic DNA was extracted from 79 xenografts using the QIAmp DNA Mini kit (QIAGEN, Hilden, Germany), according to the manufacturer's protocol and as previously described [28]. Exome enrichment and library preparation were performed using Ion Ampi-Seq Exome RDY Kit PI v3, which targets >97% of consensus coding sequences (CCDS) with 5-bp padding around exons, and Ion Xpress Barcode Adapters (Thermo Fisher Scientific, Inc.). Pooled barcoded libraries were subsequently conjugated with sequencing beads by emulsion PCR and enriched using the Ion PI Hi-Q Chef kit and Ion Chef (Thermo Fisher Scientific, Inc.). Sequencing of templates was performed with 2 samples per Ion PI Chip V3 using the Ion Proton system (Thermo Fisher Scientific, Inc.), according to the manufacturer's protocols.

Variant calling
To avoid false-positive results, we removed reads derived from the mouse genome, as follows. Sequencing reads were aligned to the human genome build 19 (hg19) and two mouse genomes C57BL/6J (mm10, NCBI accession number: GCA_000001635.26) and BALB/c (GCA_001632525.1) using the Torrent Mapping Alignment Program (ver 3.0.1, Thermo Fisher Scientific, Inc.). Reads aligned to the mouse genomes with higher alignment score than to the human genome were considered to be contamination from the host mouse and were removed from subsequent analyses. The Torrent Variant Caller plugin (ver 5.10.1, Thermo Fisher Scientific, Inc.) was used to identify variants. The parameter file, optimized for somatic mutations with low stringency criteria, was obtained from the software vendor. Variants were annotated by ANNOVAR (ver. 2018-04-16) [29] using the following reference databases: RefSeq Gene (refGene); LJB non-synonymous variants annotation (dbnsfp35a); dbSNP version 150 (avsn150); the 1000 Genome Project (1000g2015aug_eas); Clinvar version 20190305 (clinvar_20190305); COSMIC Release v88 (cosmic88); segmental duplication region (geno-micSuperDups); transcription factor binding site (tfbsConsSites); Human Genetic Variation Database version 2.3 [30]; 3.5K Japanese individuals allele frequency panel (3.5KJPNv2) [31]. Variants were filtered and excluded if they: (i) had a quality score < 30; (ii) had segmental duplication or repeat regions identified by Repeat Masker or Tandem Repeats Finder; (iii) were found in homopolymer regions or multi-allelic sites; (iv) were previously detected in 3.5KJPNv2 or Human Genetic Variation Database (HGVD).

Statistical analysis
The correlations between variant allele frequencies (VAFs) and drug sensitivities were assessed using Spearman correlation tests. Significance after Bonferroni correction for multiple testing was P = 1.11 x 10 −6 (P < 0.05; 44,875 variants). Statistical tests were conducted using Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA).

Identification of variants associated with chemosensitivity
Whole-exome sequencing was used to identify genetic variants associated with chemosensitivity to one or more of nine cytotoxic anticancer drugs (MMC, CPM, ACNU, DDP, ADR, VER, VLB, 5FU, and MTX). Drugs were administered to 79 PDX models of cancers prepared from 12 different human tissues. Between 30,012 and 45,638 variants were detected for each xenograft, with a coverage of 62 to 249 (mean; 37,921 variants, with depth of 138×). Variants were filtered using an in-house program (see Materials and methods), leaving a total of 44,875 variants for correlation analysis. Chemosensitivity was calculated as T/C and the variants whose allele frequency was higher in xenografts with lower T/C as were defined as 'chemosensitive variants' and variants whose allele frequency were higher in xenografts with higher T/C as 'chemoresistant variants'.

Subgroup analysis
We further performed a subgroup analysis based on cancer type to identify tissue-specific chemosensitivity-related variants. Subgroups of breast and gastric cancers were analyzed because more than 10 xenografts of these cancer types were available. In breast and gastric cancer xenografts, 78 and 67 variants, respectively, were possibly associated with chemosensitivity to one or more drugs, with P values < 0.01 (Tables 11 and 12). rs386792906 (chr16:g.81253642 AG>TC) in polycystin 1 like 2 (PKD1L2), which was associated with resistance to MTX, showed the strongest association of the nine tested anti-cancer drugs among the breast cancer   Table 11). rs73302038 (chr17:g.21215682 G>A) in mitogen-activated protein kinase kinase 3 (MAP2K3), which was associated with resistance to VCR, showed the strongest association for the gastric cancer subgroup (P = 8.32 × 10 −6 , r s = 0.935, Table 12). However, of the variants with P < 0.01 in the subgroup analyses, only three (chr15:g. 22960698 insG in MTX, rs59277111 in VCR, and rs3217238 in DDP) were significant at P < 0.01 in the whole-group analysis (Tables 11 and 12).

Discussion
Precision medicine demands the development of biomarkers to detect patient chemosensitivity to anti-cancer drugs. Here, we sought to identify clinically useful genetic markers for chemosensitivity to one or more of nine cytotoxic anticancer drugs by whole-exome sequencing for 79 xenografts. Although none of the genetic variants achieved a significance level after Bonferroni correction for multiple testing (P = 1.11 × 10 −6 ), numerous variants showed possible associations with chemosensitivity to each of the nine tested drugs. Moreover, the subgroup analysis indicated chemosensitivity markers specific for breast and gastric cancers. We propose that our method could contribute to the development and optimization of personalized chemotherapy regimens among patients with cancer.
In the whole-exome sequencing analysis of 79 xenografts, we found that, the variant chr8: g.22960701insC, located in TNFRSF10C and LOC254896, had the most significant (i.e., lowest) P value for its associated chemosensitivity to ADR (P = 7.15 × 10 −5 , r s = 0.437, Table 3, Fig 1). Although the function of LOC254896 remains to be clarified, the down-regulated expression  [32,33] and hypermethylation [34,35] of TNFRSF10C in colorectal, prostate, and breast cancers has been reported previously. Additionally, in vitro experiments have suggested that an upregulation in TNFRSF10C in response to ADR treatment may induce resistance to ADR [36]. TNFRSF10C is reported to protect cells from TRAIL-induced apoptosis [37], and thus may be associated with resistance to ADR through these pathways. A variant (chr15:g.63673951 insG) located in the 5'UTR of carbonic anhydrase 12 (CA12) was also associated with resistance to ADR (Table 3). CA12 is a membrane carbonic anhydrase and plays important roles in several physiological functions, such as acid-base balance and calcification [38]. A recent in vitro study showed CA12 overexpression in chemoresistant colon cancer cells expressing the drug efflux transporter P-glycoprotein (Pgp). Moreover, ADR chemosensitivity in tumors overexpressing both CA12 and Pgp can be increased using CA12 inhibitors [39]. Therefore, the chr15:g.63673951 insG variant may increase resistance to ADR by altering CA12 expression; further functional analyses would be required to verify this hypothesis.
Moreover, a variant (chr2:g.189916175 T>G) associated with resistance to MTX was located in exon 42 of collagen type V alpha 2 chain (COL5A2) ( Table 7). COL5A2 is upregulated in colorectal and breast cancers [40,41] and is associated with poor clinical outcome and poor survival rates in bladder cancer [42]. Studies have suggested that collagen expression increases tumor drug resistance by inhibiting drug penetration into the cancer tissue and increasing cellular resistance to apoptosis [43]. As shown in Table 10, we identified genetic variants that could be associated with multi-drug resistance or sensitivity. Some of these genes may be involved in the proliferation and invasion of tumor cells; for example, NYAP2 is reported to activate PI3K, Akt and Rac1, and mediates remodeling of the actin cytoskeleton   [44], whereas ZNF277 regulates cell migration and invasion through phosphatase and tensin homolog (PTEN) [45].
In the subgroup analysis using breast and gastric cancer xenografts, we identified possible tissue-specific biomarkers in the response to anticancer drugs; however, most of these variants showed weak or no association in the whole-group analysis. These results suggest a degree of tissue specificity in sensitivity to cytotoxic anticancer drugs. rs386792906 (chr16:g.81253642 AG>TC), which showed the strongest association with MTX chemosensitivity in breast cancer xenografts, was located in intron 1 of PKD1L2. PKD1L2 is a member of the polycystin protein family, and may function as a component of cationic channel pores [46]. According to a previous study using The Cancer Genome Atlas (TCGA) dataset, overexpression of PKD1L2 mRNA is associated with improved prognosis in patients with breast cancer [47]. Although the functional association between PKD1L2 and MTX is unknown, this variant may be a useful marker for predicting sensitivity to MTX, and may act as an indicator of prognosis for breast cancer in the clinical setting. We investigated the functional consequences of the associations between the top variants and the response to chemotherapy by interrogating the expression quantitative trait loci (eQTL) information in the Genotype-Tissue Expression (GTEx) database [48]. rs3830265, which showed the strongest association with sensitivity to VCR, was associated with the expression of PRKCD in the skin (P = 7.3 × 10 −7 ) and esophagus (P = 1.7 × 10 −5 ). Moreover, rs3842515, which showed the strongest association with sensitivity to ACNU in breast cancer xenograft, displayed a cis-regulatory effect on CCDC82 expression in several tissues, including esophagus, thyroid, skin, and nerve (P min = 2.3 × 10 −13 ). However, the functional associations between these genes (PRKCD and CCDC82) and sensitivities to the aforementioned drugs or mechanisms of drug metabolism remain unknown and require further investigation.
There were several strengths and limitations in our study. The main strength of our study is that we sought to identify tissue-agnostic predictive markers for chemosensitivity to nine cytotoxic anticancer drugs. As we have entered a new era of precision medicine, tissue-agnostic cancer therapy will continue to grow and expand treatment options for patients with cancer [49]. In addition to our tissue-agnostic approach, we also performed subgroup analyses of breast and gastric cancers as a deeper understanding of the genomic profiles of specific tumor types is also important. There were several limitations in our study. First, the total number of xenografts and the total number of each tumor type are small, and there were differences in the numbers of tumor types. Therefore, our study is likely to be underpowered to detect

Fig 2. Correlation between variant rs1292467126 and chemosensitivity to CPM (A), ADR (B), VCR (C), and ACNU (D).
Chemosensitivity to each drug is represented by relative tumor volume of treated mice (T) with respect to that of the control mice (C). rs1292467126 was commonly associated with increased sensitivity to four of the nine tested anticancer drugs.
https://doi.org/10.1371/journal.pone.0239614.g002 Table 11. Variants associated with chemosensitivity to each drug (P < 0.01) in breast cancer xenografts. statistically significant variants or perform a subgroup analysis for all tumor types. Second, the results need to be confirmed using a larger number of samples, along with a functional analysis of the identified genes.

Drug Chr Position
In conclusion, using whole-exome sequencing and a PDX model, we identified 162 genetic variants as possible susceptibility factors for sensitivity to one or more of the nine tested cytotoxic anticancer drugs. This method and the results presented herein may contribute to the development of personalized treatments for the prescription of optimal chemotherapy regimens. Although the underlying mechanisms should be further investigated using a larger number of clinical samples and molecular analysis, we propose that our findings may help to contribute to understanding the mechanisms of chemoresistance and chemosensitivity, and aid in the improved prognosis and quality of life for patients with cancer.