Germline variant burden in multidrug resistance transporters is a therapy‐specific predictor of survival in breast cancer patients

Multidrug resistance due to facilitated drug efflux mediated by ATP‐binding cassette (ABC) transporters is a main cause for failure of cancer therapy. Genetic polymorphisms in ABC genes affect the disposition of chemotherapeutics and constitute important biomarkers for therapeutic response and toxicity. Here we correlated germline variability in ABC transporters with disease‐specific survival (DSS) in 960 breast cancer (BRCA), 314 clear cell renal cell carcinoma and 325 hepatocellular carcinoma patients. We find that variant burden in ABCC1 is a strong predictor of DSS in BRCA patients, whereas candidate polymorphisms are not associated with DSS. This association is highly drug‐specific for subgroups treated with the MRP1 substrates cyclophosphamide (log‐rank p = 0.0011) and doxorubicin (log‐rank p = 0.0088) independent of age and tumor stage, whereas no association was found in individuals treated with tamoxifen (log‐rank p = 0.13). Structural mapping of significant variants revealed multiple variants at residues involved in protein stability, cofactor stabilization or substrate binding. Our results demonstrate that BRCA patients with high variant burden in ABCC1 are less prone to respond appropriately to pharmacological therapy with MRP1 substrates, thus incentivizing the consideration of genomic germline data for precision cancer medicine.


Introduction
Chemotherapy constitutes an essential cornerstone of oncological therapy, especially for metastatic cancers. While the number of newly approved molecules for the treatment of cancer has seen a drastic increase in recent years, the longterm survival for many cancers has remained low. 1 Importantly, common estimates are that 90% of therapeutic failures are attributed to chemotherapy resistance. 2 While a variety of mechanisms are implicated in drug resistance, facilitated drug efflux constitutes the most prevalent escape pathway for a variety of the clinically most important chemotherapeutics, such as taxanes, topoisomerase inhibitors, vinca alkaloids and antimetabolites. [3][4][5] Drug efflux and chemotherapy resistance is primarily attributed to transporters of the ATP-binding cassette (ABC) transporter superfamily, whose members mediate drug export from the cytosol into the extracellular space or designated vesicles against their electrochemical gradient. In total, the human genome encodes 48 functional ABC transporters, of which ABCB1 (encoding MDR1 or P-gp), ABCC1 (MRP1) and ABCG2 (BCRP) have been most extensively characterized. Besides their role in drug efflux, ABC transporters mediate the transport of a variety of endogenous substrates, such as amino acids, sugars, nucleosides, vitamins, lipids and bile acids, as well as a multitude of natural toxins. 6 Genetic germline variations in various ABC transporters can impact transport activity with important implications for drug toxicity and response. For instance, polymorphisms in ABCB1 have been associated with pharmacokinetics, response and toxicity of imatinib, 7,8 clearance of methotrexate, 9 as well as outcomes of doxorubicin, vincristine and prednisolone induction therapy. 10 Similarly, variants in ABCG2 were reproducibly associated with significant gefitinib toxicity 11 and outcomes to antracycline 12 and tyrosine kinase inhibitor (TKI) therapy [13][14][15] and ABCC1 polymorphisms predicted hematological toxicity in patients treated with fluorouracil, epirubicin and cyclophosphamide chemotherapy. 16 While the architecture of somatic alterations has been extensively investigated in a multitude of cancers, studies that investigate the effects of genetic factors in the germline genome on drug resistance have mostly interrogated only few selected candidate variants in relatively small, heterogeneous cohorts. Importantly, the technological progress in Next Generation Sequencing (NGS) methods and the rapid increase in available sequencing data allows for the first time to systematically profile the genetic landscape in multidrug resistance transporters and to relate this genetic complexity to clinical outcomes in patients undergoing chemotherapy. Specifically, we analyzed whole exome and whole genome sequencing (WES and WGS, respectively) data from 138,632 individuals and provide a consolidated overview of the genetic variability in key drug transporters on an unprecedented scale. Importantly, 98.4% of the 3,476 variants identified in ABCB1, ABCC1 and ABCG2 were rare with minor allele frequencies (MAF) <1%. Strikingly, analyses of survival data from 960 breast cancer (BRCA), 314 clear cell renal cell carcinoma (ccRCC) and 325 hepatocellular carcinoma (HCC) patients revealed that variant burden in ABCC1 and ABCG2 were cancer and treatment-specific biomarkers of disease-specific survival (DSS), whereas previously reported candidate variants did not predict DSS, thus indicating that sequencing-based profiling can add valuable information for the prediction of clinical outcomes. The presented study constitutes the most comprehensive analysis of genetic variation in the multidrug resistance transporters published to date and demonstrates the effects of this genetic complexity on drug resistance in cancer therapy.

Population-scale sequencing data
Single nucleotide variant (SNV) and indel data of 48 human ABC transporters were assembled based on WES and WGS data from 138,632 individuals acquired from the Genome Aggregation Database. 17 Information about copy number variations (CNVs) from 59,451 individuals were obtained from the Exome Aggregation Consortium and analyzed as previously described. 18 Linkage disequilibria were computed by leveraging linkage from the 1000 Genomes Project using LDLink. 19

Cancer cohorts
Analysis-ready bam files storing aligned (genome build hg38) whole-exome sequencing data from blood samples as well as tumor-adjacent normal tissue from the three investigated cohorts (BRCA: breast invasive carcinoma; KIRC: kidney renal clear cell carcinoma; LIHC: liver hepatocellular carcinoma) from The Cancer Genome Atlas (TCGA) were downloaded from the Genomic Data Commons Data Portal. Germline variants were called using the GATK software (version 4.1). Variant calling was confined to exome regions as defined in the Agilent SureSelect Version 7 interval list. HaplotypeCaller, variant quality score recalibration and genotype refinement were applied as suggested in GATK Best Practices. Clinical and survival data of TCGA samples were retrieved from Liu et al. 20 Demographic and clinical data of the three analyzed TCGA cohorts are summarized in Table 1.
In all subsequent analysis, only blood samples were considered for the BRCA and LIHC cohort. For the KIRC cohort, we only considered samples classified as ccRCC. 21 As less blood samples were available for this tumor type, we considered both blood samples and tumor-adjacent normal tissue for genomic analyses. In case that genotype data were available for more than one sample of a patient, blood samples were selected. If for a patient several blood samples were genotyped or no blood sample but several tissue samples were available, differing genotypes were set to missing. In addition, What's new? Chemotherapy resistance is due in large part to ATP-binding cassette (ABC) transporters, which actively remove drugs from cells. Multiple genetic variants, however, can alter ABC transporter activity, often modifying therapeutic responses. In this study, germline mutational burden in the multidrug resistance transporter ABCC1 (MRP1) strongly predicted survival in breast cancer patients treated with the MRP1 substrates cyclophosphamide and doxorubicin. Structural analyses of mutated transporters revealed changes in protein stability and substrate binding. The findings delineate a role for high ABCC1 variant burden in driving chemoresistance in breast cancer patients and highlight the value of germline sequencing for predicting therapeutic response.
the genotypes of variants with genotype quality <20 were set to missing. For each patient and each gene, we computed the number of variants carrying the respective alternative allele.

Computational functionality predictions
The functional consequences of missense variants were analyzed using computational algorithms that assess sequence information as well as consequences of genetic variants on the respective protein structure. Specifically, we selected SIFT, Poly-phen2, MutationAssessor, VEST3 and Eigen, as these tools were among the best performing models in four independent benchmarking data sets. 22,23 Variants were categorized as deleterious when the majority of algorithms predicted functional consequences. In addition, we classified all variants as deleterious that resulted in frameshifts, inframe deletions or insertions, start-lost, stop-gained or that affected canonical splice sites.

Structural analyses
The crystal structure of the human BCRP transporter (6ETI) was obtained from RCSB PDB. As no crystal structure for human MRP1 was available, we modeled its structure using Phyre 2 onto the resolved bovine MRP1 structure (6BHU) with 100% prediction confidence. 24 PyMOL (version 2.1.1) was used to map the genetic variability onto transporter structures.

Statistical analyses
Optimal variant burden cut-offs for cohort stratification with regard to DSS were determined using conditional inference trees implemented in R using the library partykit_1.2-2 [25][26][27] with the minimal required group size set to 20 and the p-value criterion relaxed until a split was identified. Only the first binary split was used in subsequent analyses and corresponding permutation p values, properly accounting for the plurality of investigated splits are stated. To avoid small group sizes, cut-offs of the respective overall cohort were transferred for drug subgroup analyses. Since permutation p values are not applicable for these investigations, log-rank test p values were used instead. To facilitate comparisons, we also provide log-rank p values for overall cohort analyses. These analyses were complemented by linear tail-restricted cubic spline modeling (with five knots) using the R package rms (version 5.1-2) in order to detect potential effects of dichotomization, such as the masking of nonmonotone trends. 28 Log-rank tests were also applied to compare DSS between carriers and noncarriers of the single variants.  (8), pazopanib (7), sorafenib (7), gemcitabine (6) Sorafenib (26) Associations between single variants and variant burden were investigated by asymptotic linear-by-linear association tests implemented in the R-package coin 1.3-0. 29 Here, alternative alleles were merged if different variants were present for a given genomic coordinate. For each variant, the linear-bylinear association test investigates the independence of variant burden and genotype grouping, coded as ordinal variables.
The underlying model assumes that deviation from independence is linear in variant burden for a fixed genotype group and vice versa. All tests were two-sided and statistical significance was defined as p < 0.05. Where indicated, Bonferroni correction was used to correct for multiple testing.

Data availability
TCGA bam files used in the analyses were downloaded from the Genomic Data Commons Data Portal (https://portal.gdc. cancer.gov). Clinical and survival data of TCGA samples were retrieved from Reference 20. Genetic variability data of ABC transporters in the general population were retrieved from the Genome Aggregation Database (http://gnomad.broadinstitute. org/). Data retrieved from the Genome Aggregation Database is available under the Fort Lauderdale agreement, supporting free and unrestricted use of the retrieved information. Germline data from the TCGA are managed by the National Cancer Institute (NCI) and the National Human Genome The vast majority of these variants were rare with 98.4% occurring in <1% of alleles worldwide. In addition, 52% of all variants were only found in a single individual. (c) Besides single nucleotide variations and indels, the analyzed ABC transporters harbor rare copy number variations (CNVs) in which at least one exon is deleted or duplicated. (d) Computational predictions using five partly orthogonal algorithms suggest that ABCB1, ABCC1 and ABCG2 harbor hundreds of variants with putative functional effects. Error bars indicate standard error of the mean (SEM) across five computational algorithms (see methods for details). (e) On average each individual (diploidy considered) was predicted to have between 0.15 variants in ABCC1 and 1.2 variants in ABCB1 that modulate the functionality of the encoded transporter. Rare variants were predicted to account for <5% of the genetically encoded functional differences in ABCB1 and ABCG2, whereas they were estimated to contribute up to 40% for functional ABCC1 variability.
Research Institute (NHGRI) and access authorization was granted by the National Institute of Health (NIH) Data Access Committee. All data was retrieved and analyzed in compliance with NIH genomic data sharing policy.
To estimate the functional importance of the observed ABC variability, we used five partly orthogonal computational algorithms that predict the functional consequences of genetic variations based on a variety of diverse features, including sequence information, evolutionary conservation, structural features as well as functional genomic data, indicating that of the 1908 missense variants, 920 (48.3%) were predicted to Most variants with putative functional consequences were found in ABCC1 (n = 453), followed by ABCB1 (n = 344) and ABCG2 (n = 315; Fig. 1d). Notably, only 12.5% (1 of 8) of common missense variants with MAF >1% were putatively deleterious, compared to 48.5% (921 of 1900) for rare variations.
A diploid human genome was found to harbor on average 1.2 variants with putative functional effects in ABCB1, whereas ABCG2 and ABCC1 contained only 0.48 and 0.15 variants, respectively (Fig. 1e). Interestingly, in ABCB1 and ABCG2, rare variations are estimated to account for 2.2 and 3.2% of the genetically encoded functional variability in these genes, respectively, whereas rare variants explained 38.4% of the functional variability in ABCC1 (Fig. 1e). Combined, each individual was found to harbor 1.8 variants that entail functional alterations in major multidrug resistance transporter. These data demonstrate that the genetic variability in these transporter genes is complex and might contribute to interindividual variability in the response to chemotherapeutics.

The burden of ABC variability predicts clinical outcomes
To analyze whether ABCB1, ABCC1 and ABCG2 variability could predict chemotherapeutic outcomes, we analyzed disease-specific survival (DSS) of 960 BRCA, 314 ccRCC and 325 HCC patients. Strikingly, we found that the overall germline variant burden constituted a drug-and cancer typespecific predictor for DSS ( Fig. 2 and Table 2), whereas common ABC polymorphisms were not associated with altered sensitivity to chemotherapy (Supporting Information Table S1). Specifically, variant load in ABCC1 was significantly associated with DSS in BRCA patients (HR of high vs. low variant load = 3.22; 95% CI = [1.62-6.4]; log-rank p = 0.00043; permutation p = 0.031) and this association remained significant after adjusting for patient age and tumor stage (Wald test p = 0.02). Similar associations were not detected for ABCB1 (log-rank p = 0.08; permutation p = 0.076) and ABCG2 (logrank p = 0.061; permutation p = 0.81).
To explore the mechanisms underlying this observation, we performed subgroup analyses in which we stratified the cohort by prescribed medicines. Doxorubicin (n = 329), cyclophosphamide (n = 238) and tamoxifen (n = 238) were most widely prescribed. Interestingly, when applying the conditional inference tree cut-offs from overall BRCA analysis, associations were confirmed for ABCC1 with doxorubicin (HR = 4.57; 95% CI = [1.31-15.93]; log-rank p = 0.0088; Wald test p = 0.015) and cyclophosphamide (HR = 9.22; 95% CI = [1.83-46.36]; log-rank p = 0.0011; Wald test p = 0.0085). In contrast, ABCC1 variability was not found to correlate with outcomes in patients treated with tamoxifen (log-rank p = 0.13; Wald test p = 0.081). In addition, we identified a drug-specific association between variant load in ABCG2 and DSS in patients treated with doxorubicin (log-rank p = 0.018). Variant burden did not correlate with DSS in ccRCC and  1 Hazard ratio and log-rank test comparing patients with variant burden above versus below the cutoff determined by conditional inference trees. 2 Permutation p-value from conditional inference tree analysis. 3 For drugs subgroup analyses, conditional inference trees cutoffs for the respective overall cohort were applied.
HCC patients of which only a minor fraction of individuals received chemotherapy (permutation p > 0.05; Table 2). We complemented our statistical evaluations by modeling the relationship between DSS and variant burden using cubic splines. Interestingly, ABCC1 variability and DSS were significantly associated with the BRCA cohort (p = 0.042) with the relative hazard increasing specifically in individuals with more than 20 variants in ABCC1 (Fig. 3a). In contrast, no significant associations were observed for ccRCC (p = 0.46) and HCC patients (p = 0.61; Figs. 3b and 3c). Variability in ABCB1 or ABCG2 was not predictive for DSS in any cancer type analyzed (Supporting Information Fig. S1). Combined, these results corroborate the gene-specific association of ABCC1 variant burden with clinical outcomes in BRCA patients.
Notably, cubic spline modeling of the BRCA drug subgroups confirmed the significant correlations between variant burden in ABCC1 and ABCG2 with cyclophosphamide (p = 0.008) and doxorubicin treatment success (p = 0.001), respectively, observed with dichotomized burden (Supporting Information Fig. S2). However, these associations were nonmonotone, suggesting additional gene-gene or gene-drug interactions. In contrast to the analysis with categorized variant burden, the association between ABCC1 and DSS in the doxorubicin treatment group was not significant (p = 0.27).
Multiple variants with small effect sizes modulate diseasespecific survival Next, we explored whether specific genetic variants were associated with DSS and variant burden in ABCC1 and ABCG2 for BRCA patients. Notably, 34 ABCC1 variants were significantly associated with variant burden after Bonferroni correction (p < 8.6 × 10 −5 using linear-by-linear association tests), of which 21 were intronic, 10 were synonymous and three (G671V, C1047S and V1146I) resulted in amino acid exchanges (Table 3).   Linkage analysis revealed complex haplotype structures among significantly associated variants in Europeans (Supporting Information Fig. S3). In addition, four missense (V12M, R56Q, Q141K and F489L), one synonymous and 12 intronic variants in ABCG2 were significantly associated with variant burden in BRCA. Ten and two variants in ABCC1 and ABCG2, respectively, showed significant associations with DSS after Bonferroni correction, all of which were infrequent with minor allele frequencies ≤1.6%. The most significant associations were identified for the synonymous variant rs45607431 (H986H) in ABCC1 (log-rank p = 2.4 × 10 −14 ) and the stop-gain variant rs200190472 in ABCG2 (log-rank p = 3.6 × 10 −13 ). However, as statistical analyses of such rare variant associations are not robust, we want to emphasize that these results require further validation. In summary, our data show that multiple variants with individually small effect sizes, which are missed when using conventional variant interrogations, contribute to clinical outcomes and can explain a significant fraction of clinical outcomes in BRCA patients undergoing chemotherapy.

Structural mapping of the genetic ABC variability indicates mechanisms of drug resistance
To obtain further insights into the effects of genetic variability on transporter function, we used structural analyses in which we mapped variants of interest onto the tertiary structures of the corresponding human transporter proteins MRP1 (ABCC1) and BCRP (ABCG2). For BCRP, we used available highresolution crystal structures, 30,31 whereas we used homology modeling to derive the structure of human MRP1. Specifically, we mapped all variants found to be significantly associated with the high-risk group. Furthermore, as statistics for very rare variants, particularly singletons, are problematic, we also included all variations that were enriched in the high-risk group by more than threefold.
MRP1 is organized as two transmembrane domains (TMDs) that form the channel for substrate translocation and two cytoplasmic nucleotide-binding domains (NBDs; Fig. 4a), of which however only NBD2 has a functional ATP hydrolysis site. 32 In addition, MRP1 contains an N-terminal TMD0 domain with unclear functional relevance that is found exclusively in members of the ABCC subfamily. Based on available homologous crystal structures, we accurately modeled the structure of human MRP1, as indicated by an excellent alignment with the experimentally determined crystal structure of bovine MRP1, 33 which is 91% amino acid identity (Supporting Information Fig. S4).
R1066W and C1047S are located within the cytoplasmic loop CL6 in TMD2 or at its interface on the outside of the helix projecting into the lipid bilayer, respectively. CL6 interacts with the six-amino group of ATP thereby stabilizing the conformational change upon ATP binding to the TMDs. 34 Cysteine commonly stabilizes hydrophobic interactions and differs drastically in interaction profiles from the polar amino acid serine. 35 As such, the C1047S variant likely causes changes in CL6 positioning, which in turn might translate into differential ATP stabilization in the active ATP hydrolysis site in NBD2. R230Q and G231C is located in the Lasso motif that links TMD0 to TMD1, specifically within a section that has been shown to be essential for MRP1 folding and function. 36 MRP1 transports many substrates that are conjugated with large anionic groups and the amino acid exchange variant R433S is positioned close to the Ppocket at the interface between TMD1 and TMD2 that coordinates such moieties, as shown for glutathione, 33 and the exchange of basic arginine with uncharged serine likely destabilizes these interactions.
In contrast to MRP1, BCRP is organized as a homodimer of two half-transporters (Fig. 4b). The common Q141K polymorphism resides in an amphipathic α-helix of BCRP and leads to instability in the NBD, misfolding and ubiquitin-mediated proteasomal degradation, likely due to electrostatic repulsion. 37 In addition, we detected rare variants exclusively in patients with reduced DSS that result in amino acid exchanges of residues lining the substrate translocation pore, such as polar threonine at position 542 to nonpolar alanine (T542A) and of cysteine 438 to Only variants that were significantly associated with gene-specific variability burden after Bonferroni correction (p < 8.6 × 10 −5 ) are shown. 1 Variant carriers versus noncarriers; endpoint DSS. 2 Variant burden versus genotype groups, both coded as ordinal variables. 3 Alternative alleles were merged beforehand. 4 Global allele frequency in gnomAD (gnomad.broadinstitute.org).
charged arginine (C438R), as well as variants that are in <8 Å proximity to the substrate-binding pocket, such as the substitution of bulky aromatic phenylalanine at position 489 to hydrophobic, small leucine (F489L). 38

Discussion
Multidrug resistance transporters bind to >100 clinically relevant drugs and affect their disposition and pharmacokinetics, which is of particular importance for cancer drug resistance. 39 Our analyses of genomic data from 138,632 individuals indicate that the genetic complexity in these transporter genes is extensive and comparable to the variability in other highly variable pharmacogene families, such as CYPs, 40 SLCOs 41 and UGTs. 42 While a multitude of reports have implicated germline variants in multidrug resistance transporters in differences in chemotherapy response and dose-limiting toxicity, these studies interrogated only candidate polymorphisms and failed to incorporate the genetic complexity of these transporter genes into their associations, which might be insufficient for the functional interpretation of complex ABC genotypes. In agreement with this hypothesis, we found significant associations of ABCC1 (MRP1) variability with DSS in BRCA patients, whereas no significant correlations were identified when only individual candidate SNPs were considered. The majority of BRCA patients in our cohort received cyclophosphamide, doxorubicin or tamoxifen. Whereas the cytotoxic agents cyclophosphamide and doxorubicin exert direct cytotoxicity by DNA crosslinking or intercalation, respectively, tamoxifen acts cytostatic by modulation of estrogen receptor signaling. Importantly, when we stratified our cohort by treatment regimen, ABCC1 variability was significantly correlated with DSS in cyclophosphamide and doxorubicin subgroups, whereas no association was identified in an identically sized subgroup of patients treated with tamoxifen. These results align with the pharmacokinetics of these drugs, as cyclophosphamide and doxorubicin are high-affinity substrates of MRP1, 43 whereas tamoxifen and its active metabolite endoxifen are primarily transported by P-gp. 44 Moreover, a previous clinical study of 516 premenopausal hormone receptor-positive BRCA patients demonstrated that MRP1 expression correlates with resistance to cyclophosphamide, Figure 4. Structural analysis of genetic variants in important multidrug resistance genes. Schematic topology models as well as 3D protein structures are shown for MRP1 (a) and BCRP (b). Different domains in the topology models are shaded based on the identified number of deleterious variants per amino acid in the respective domain identified in NGS data from 138,632 individuals. MRP1 consists of two pseudosymmetrical transmembrane and nucleotide-binding domains (TMD and NBD) encoded in a single polypeptide, colored in orange and blue, respectively. In contrast, BCRP is a half-transporter that forms homodimers. Detailed 3D structure of key protein domains with functionally relevant variants (sticks in cyan or magenta) and substrates (sticks in yellow) are shown as insets under the respective topology models. In the 3D models, all variants significantly associated with the high-risk group of breast cancer patients, as well as all variations that were enriched in the high-risk group by more than threefold are shown. The corresponding residues in the secondary structure are colored in red for nonsynonymous variants and in blue for synonymous variants. methotrexate and fluorouracil chemotherapy, whereas no association with response to adjuvant endocrine treatment with tamoxifen was observed, 45 thus providing further clinical support for the reported associations. In addition to ABCC1 variability, drugspecific analyses revealed a significant correlation between ABCG2 variant load and DSS in patients treated with doxorubicin, an association that was not identified at nominal significance level in the overall BRCA cohort. BCRP encoded by ABCG2 has been implicated in the emergence of doxorubicin resistance in BRCA in vitro, which might lend mechanistic plausibility for this finding. 46 Importantly, ABC variability was not predictive of outcomes in patients with cancer types that are primarily treated by surgery, such as clear cell renal cell carcinoma and hepatocellular carcinoma (permutation p > 0.05), further corroborating the validity of our approach. In addition, our previous work supports the data for ccRCC since germline and somatic variants could not explain aberrant expression of ABCG2 and consequently patient survival. 47 Notably, while variant burden in ABCC1 and ABCG2 showed a monotone association with DSS in the overall BRCA cohort, subgroup analyses revealed nonmonotone relationships in the drug subgroups. As the genomic landscape of ABCC1 and ABCG2 is complex with a multitude of lowfrequency haplotypes, we speculate that specific haplotypes of low-effect size variants might drug-specifically modulate DSS. However, due to the relatively low numbers of individuals in the subgroups, analyses of combinatorial effects require replication in additional cohorts of BRCA patients treated with cyclophosphamide and doxorubicin.
Structural mapping of significant variants based on available high-resolution crystal structures of the respective transporters found multiple variants to reside in functionally important domains. The most common ABCC1 variant identified was G671V (global allele frequency 3.8%), which was previously associated with acute anthacycline induced cardiotoxicity. 48 G671 localizes near the Walker A motif in NBD1 and impairs extrusion of the reactive electrophile 4-hydroxy-2-trans-nonenal, a cardiotoxic product of doxorubicin redox cycling. 49 While the C1047S variant has to our knowledge not yet been identified in clinical studies as a modulator of chemotherapy outcomes, in vitro studies demonstrated the effects of this variant on chemotherapy resistance. 50 In addition to these missense variants, we identified a variety of synonymous variants in ABCC1 that associate with poor DSS in BRCA. Multiple of these variants have previously been associated with chemotherapy response or toxicity. N354N correlated with increased risk of hematological toxicity in colorectal cancer patients treated with FOLFOX4, suggesting increased function of this transporter variant. 51 Similarly, V275V, which is in linkage disequilibrium with N354N, is implicated in anthracycline-induced cardiotoxicity and febrile neutropenia in BRCA patients treated with fluorouracil, epirubicin and cyclophosphamide. 16,52 Synonymous variants can affect the structure of substrate and inhibitor interaction sites in ABCB1, likely due to the introduction of a rare codon, which affects the timing of cotranslational folding. 53,54 While such ribosomal stalling has to our knowledge not yet been reported for ABCC1, these studies raise the possibility that similar mechanisms might exist for this structurally related transporter.
Response to neoadjuvant chemotherapy is a strong prognostic indicator of survival in BRCA patients. 55 However, nonresponders can benefit from systemic chemotherapy or preoperative radiation. 56,57 It is thus of fundamental importance to identify nonresponders as early as possible to flag those individuals for salvage therapies. Importantly, current monitoring methods, such as conventional ultrasound or Xray mammography, can only detect major changes in morphological tumor features, whereas histological analyses of biopsies are more sensitive but are generally conducted only months after chemotherapy. Our results demonstrate that BRCA patients with high variant burden in ABCC1 are less prone to respond appropriately to pharmacological therapy with MRP1 substrates. Thus, such at-risk individuals are likely to benefit from earlier and more frequent monitoring, ideally using more sensitive methods, such as magnetic resonance imaging 58 or quantitative ultrasound spectroscopy, 59,60 which can detect early events (tumor cell death) that precede changes in tumor size.
In our study, we analyzed the impact of genetic variability of multidrug resistance transporters in the germline genome. While our data indicate important roles in chemotherapy resistance, we anticipate that somatic factors are further important contributors. Various genetic, epigenetic and posttranslational mechanisms have been reported to underlie the somatic multidrug resistance phenotype, including somatic gene rearrangements that result in the fusion of the ABCB1 gene with a more active promoter, 61 chromatin reorganization and hypomethylation of ABC promoters, 62 posttranscriptional regulation of expression through miRNAs, 63 transcriptional activation of other ABC family members 64 and posttranslational stabilization of transporter proteins. 65 We thus suggest that models to predict personalized chemotherapy response will benefit from an integration of both germline and somatic factors to maximize predictive accuracy.
In summary, we comprehensively mapped the genetic variability of the human ABC superfamily of transporters and revealed surprising genetic complexity with a plethora of rare single nucleotide and copy number variations. Importantly, we identified germline variant burden as a gene-and drug-specific predictor of DSS in BRCA patients undergoing chemotherapy, whereas no significant associations were observed for individual candidate markers. These results incentivize the consideration of sequencing-based genotyping for precision cancer medicine and suggest increased surveillance and monitoring for BRCA patient with high variant load in multidrug resistance transporters. Furthermore, we anticipate that a deeper mechanistic understanding of the functional consequences of ABC transporter variability can provide guidance for personalized pharmacogenetic treatment decisions and the optimization of precision public health programs.