(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .
Anthracyclines induce cardiotoxicity through a shared gene expression response signature [1]
['E. Renee Matthews', 'Department Of Biochemistry', 'Molecular Biology', 'University Of Texas Medical Branch', 'Galveston', 'Texas', 'United States Of America', 'Omar D. Johnson', 'Biochemistry', 'Cellular']
Date: 2024-03
TOP2 inhibitors (TOP2i) are effective drugs for breast cancer treatment. However, they can cause cardiotoxicity in some women. The most widely used TOP2i include anthracyclines (AC) Doxorubicin (DOX), Daunorubicin (DNR), Epirubicin (EPI), and the anthraquinone Mitoxantrone (MTX). It is unclear whether women would experience the same adverse effects from all drugs in this class, or if specific drugs would be preferable for certain individuals based on their cardiotoxicity risk profile. To investigate this, we studied the effects of treatment of DOX, DNR, EPI, MTX, and an unrelated monoclonal antibody Trastuzumab (TRZ) on iPSC-derived cardiomyocytes (iPSC-CMs) from six healthy females. All TOP2i induce cell death at concentrations observed in cancer patient serum, while TRZ does not. A sub-lethal dose of all TOP2i induces limited cellular stress but affects calcium handling, a function critical for cardiomyocyte contraction. TOP2i induce thousands of gene expression changes over time, giving rise to four distinct gene expression response signatures, denoted as TOP2i early-acute, early-sustained, and late response genes, and non-response genes. There is no drug- or AC-specific signature. TOP2i early response genes are enriched in chromatin regulators, which mediate AC sensitivity across breast cancer patients. However, there is increased transcriptional variability between individuals following AC treatments. To investigate potential genetic effects on response variability, we first identified a reported set of expression quantitative trait loci (eQTLs) uncovered following DOX treatment in iPSC-CMs. Indeed, DOX response eQTLs are enriched in genes that respond to all TOP2i. Next, we identified 38 genes in loci associated with AC toxicity by GWAS or TWAS. Two thirds of the genes that respond to at least one TOP2i, respond to all ACs with the same direction of effect. Our data demonstrate that TOP2i induce thousands of shared gene expression changes in cardiomyocytes, including genes near SNPs associated with inter-individual variation in response to DOX treatment and AC-induced cardiotoxicity.
Anthracycline drugs such as Doxorubicin are effective treatments for breast cancer; however, they can cause cardiotoxicity in some women. It is unclear whether women would experience the same toxicity for all drugs in this class, or whether specific drugs would be better tolerated in specific individuals. We used an in vitro system of induced pluripotent stem cell-derived cardiomyocytes from six healthy females to test the effects of five breast cancer drugs on cell heath and global gene expression. We identified a strong shared cellular and gene expression response to drugs from the same class. However, there is more variation in gene expression levels between individuals following treatment with each anthracycline compared to untreated cells. We found that many genes in regions previously associated with Doxorubicin-induced cardiotoxicity in cancer patients, respond to at least two drugs in the class. This suggests that drugs in the same class induce similar effects on an individual’s heart. This work contributes to our understanding of how drug response, in the context of off-target effects, varies across individuals.
Funding: This work was funded by a Cancer Prevention Research Institute of Texas (
https://www.cprit.state.tx.us/ ) Recruitment of First-Time Faculty Award (RR190110) to M.C.W. O.D.J was supported by a Jeane B. Kempner Predoctoral Fellowship administered through UTMB (
http://www.kempnerfund.org/ ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. S.R.P was supported by an NSF (
http://www.nsf.gov ) grant (NSF1933321) to Dr. Stefan Bossman (PI).
We thus designed a study to compare the cardiomyocyte response to five different breast cancer drugs in the same set of individuals. To do so, we established an in vitro model of cardiotoxicity using iPSC-CMs from a panel of six healthy female individuals and drugs that inhibit TOP2 (TOP2i) including the ACs DOX, EPI and DNR, the anthracenedione MTX, and TRZ. We collected data for cell viability, calcium handling, and global gene expression to understand the gene regulatory response to these drugs over time. This framework allowed us to identify a gene expression signature associated with the TOP2i response and gain insight into expression variability across individuals in response to different TOP2i.
There is conflicting evidence about whether analogs such as EPI are equally or less likely to cause cardiotoxic effects than DOX [ 27 , 28 ]. Cardiomyopathy incidence based on pediatric cancer patient follow-up has shown that for every case of cardiomyopathy in patients treated with DOX, there are 0.8 cases for patients treated with EPI, 0.6 for DNR, and 10.5 for MTX [ 29 ]. However, clinical trials in breast cancer patients suggest that MTX is less cytotoxic than DOX [ 30 ]. Unrelated breast cancer drugs such as the HER2 receptor agonist Trastuzumab (TRZ) have also been shown to induce cardiotoxicity in as many as 6% of patients [ 31 ]; however these effects appear to be reversible [ 32 ]. It is challenging to compare drug-induced cardiotoxicity across different populations of individuals where each individual is administered a different treatment.
The National Cancer Institute has approved tens of drugs for use in the treatment of breast cancer. This list includes ACs that are structural analogs of DOX, such as Epirubicin (EPI), which is an epimer of DOX, Daunorubicin (DNR), which differs from DOX by a hydroxyl group, and Mitoxantrone (MTX), which is an anthracenedione that is structurally related to ACs. All of these drugs are considered to be intercalating TOP2 poisons and show evidence of cardiotoxicity [ 24 – 26 ].
The precise mechanistic basis of DOX-induced cardiotoxicity is unclear. DOX intercalates into DNA forming a complex with DNA topoisomerase II (TOP2), which results in double-stranded breaks [ 16 ]. This leads to the aberrant activation of the p53 stress response pathway, mitochondrial dysfunction and apoptosis [ 17 ]. There are two TOP2 isoforms in humans–TOP2A and TOP2B. TOP2A expression is regulated in a cell cycle-dependent manner and is essential for cell division, while TOP2B contributes to transcriptional regulation in mitotic and post-mitotic cells [ 18 ]. TOP2B is essential for the cardiotoxicity observed in mice [ 19 ] and disruption of TOP2B in iPSC-CMs decreases the sensitivity to DOX [ 20 ]. In addition to inducing DNA damage, DOX has been shown to evict histones and initiate chromatin damage [ 21 ]. Indeed, in breast cancer patients, sensitivity to ACs is mediated through chromatin regulators [ 22 ]. It has been proposed that ACs, which have only chromatin-damaging activity and not DNA damage activity, do not induce cardiotoxicity [ 23 ].
Breast cancer survivors are more likely to suffer from heart failure and arrhythmia than ischemic heart disease, suggesting that the heart muscle itself is affected by the treatment [ 11 ]. Cardiomyocytes make up 70–85% of the heart volume [ 12 ], and are the target of DOX-induced toxicity given their high mitochondrial content and metabolic activity [ 13 ]; however these cells are challenging to obtain from humans. With the advent of iPSC technology, we are now able to acquire easily accessible cell types from blood from humans, and reprogram these cells into iPSCs, which can subsequently be differentiated into cardiomyocytes (iPSC-CMs). This in vitro iPSC-CM system has been shown to recapitulate the clinical effects of DOX-induced cardiotoxicity including apoptosis, DNA damage, and oxidative stress in cells from breast cancer patients treated with DOX [ 14 ]. iPSC-CMs generated from 45 healthy individuals and treated with DOX revealed hundreds of genetic variants that associate with the gene expression response to DOX treatment (eQTLs) [ 15 ]. These studies suggest that there is a genetic basis to DOX-induced cardiotoxicity, and highlight the need, and potential, for personalized medicine to reduce side effects of chemotherapeutic agents.
Anthracyclines (ACs) such as Doxorubicin (DOX), which is prescribed in ~32% of breast cancer patients, can cause left ventricular dysfunction and heart failure both during treatment, and years following treatment [ 5 ]. The risk of adverse cardiovascular events increases with higher doses of DOX [ 4 ]. DOX-induced congestive heart failure has been observed in 5% of patients treated with 400 mg/m 2 of the drug, 26% of patients treated with 550 mg/m 2 , and 48% of patients treated with 700 mg/m 2 [ 6 ]. However, patients with CVD risk factors treated at low doses are also at risk for cardiotoxicity [ 6 ], suggesting inter-individual variation in response to drug treatment. Indeed, genome-wide association studies (GWAS) have implicated a handful of genetic variants in AC-induced cardiotoxicity that are close to genes including RARG, SLC28A3, UGT1A6 and GOLG6A2/MKRN3 [ 7 – 10 ].
Globally, breast cancer is the most common cancer in women [ 1 ]. 13% of women in the United States will be diagnosed with breast cancer during their lifetime [ 2 ]. While the number of deaths attributed to the disease is decreasing, and the 5-year survival rate is 90.8%, there are an estimated 3.8 million women living with breast cancer [ 2 ]. Breast cancer survivors are now more likely to suffer from secondary conditions such as cardiovascular disease (CVD) than tumor recurrence [ 3 ]. This is likely due to shared risk factors for breast cancer and CVD, and the cardiotoxic side effects induced by chemotherapeutic agents [ 4 ].
Results
We obtained iPSCs from six healthy female individuals in their twenties and thirties with no known cardiac disease or cancer diagnoses and differentiated these cells into cardiomyocytes. iPSC-CMs were metabolically selected and matured in culture for ~28 days post initiation of differentiation (See Methods). The purity of the iPSC-CM cultures was determined as the proportion of cells expressing cardiac Troponin T by flow cytometry. The median iPSC-CM purity was 97% (range 63–100%) across individuals (S1 Fig).
TOP2-inhibiting breast cancer drugs decrease iPSC-CM viability We studied the response of iPSC-CMs to a panel of drugs used in the treatment of breast cancer (Fig 1A). Specifically, we chose drugs belonging to AC and non-AC classes that inhibit TOP2: DOX, EPI, DNR and MTX, and TRZ as an unrelated monoclonal antibody. Given the reported cytotoxic effects of several of these drugs, we first sought to measure iPSC-CM viability across six individuals following drug and vehicle (VEH) exposure at a range of drug concentrations (0.01–50 μM) for 48 hours. We observe a dose-dependent decrease in viability following DOX, EPI, DNR, and MTX treatments and no effects on viability for TRZ and VEH (Fig 1B and S1 Table). These effects are consistent across independent differentiations from the same individuals (S2 Fig), as well as across individuals (Fig 1B). Following analysis of the dose-response curves, we extracted the concentration which resulted in 50% cardiomyocyte cell death (LD 50 ) for each drug, across each replicate, in each individual, and calculated the median value across all six individuals. We find that the median DNR LD 50 and MTX LD 50 values are significantly lower than those for DOX (T-test; P < 0.05; median LD 50 in μM: DOX = 14.02, DNR = 0.98, EPI = 3.79, MTX = 0.98; Fig 1C and S2 Table). PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 1. TOP2i drugs affect cardiomyocyte viability in a dose-dependent manner. (A) Experimental design of the study. iPSCs derived from six healthy women aged 21 to 32 were differentiated into cardiomyocytes (iPSC-CMs) and exposed to a panel of drugs used in the treatment of breast cancer. The response to TOP2i Doxorubicin (DOX), Epirubicin (EPI), Daunorubicin (DNR) and Mitoxantrone (MTX) is compared to the non-TOP2i Trastuzumab (TRZ) and a water vehicle (VEH). Effects on cell viability, cellular stress, cell function, and gene expression are measured. (B) Proportion of viable cardiomyocytes following exposure to increasing concentrations of each drug. Cell viability in each individual (colored line) was assessed following 48 hours of drug treatment. Dose-response curves were generated using a four-point log-logistic regression with the upper asymptote set to one. Each data point reflects the mean viability from two independent differentiations per individual where each is measured in quadruplicate. (C) LD 50 values for each drug treatment for each individual. Each value is calculated as the mean across two independent experiments. LD 50 values could not be calculated for TRZ and VEH treatments. (D) EC 50 values for each drug treatment in ten breast cancer cell lines were obtained from the Depmap portal (
https://depmap.org) using the PRISM [33], CTRP CTD2 [34], and GDSC2 [35] databases. Asterisk represents a statistically significant change in LD 50 (C) or EC 50 (D) values between each drug and DOX treatment (*p < 0.05, **p < 0.01, ***p < 0.001).
https://doi.org/10.1371/journal.pgen.1011164.g001 In order to understand the cell-type specificity of the drug responses, we extracted high confidence EC 50 values following treatment with the same drugs in a panel of ten breast cancer cell lines from the PRISM, GDSC2 and CTRP databases [33–35]. Interestingly, we find that DNR is less cytotoxic in breast cancer cell lines than DOX and MTX (T-test; P < 0.05), and that MTX is the most cytotoxic drug in the panel (Fig 1D). This indicates that cancer drug doses tested on breast cancer cell lines may not be predictive of cardiotoxicity. To determine the physiological relevance of our model, we collated observed plasma concentrations of these drugs in patients being treated for cancer. We find the blood serum levels range from 0.002–1.73 μM, indicating that the concentrations we identified in vitro could be observed in vivo and therefore warrant further study (S3 Table). To determine whether the effects on viability are associated with cellular stress we measured activity of secreted lactate dehydrogenase. We observe a significant inverse correlation between cell viability and cellular stress across TOP2i concentrations (DNR rho = -0.54; DOX rho = -0.35; EPI rho = -0.34; MTX rho = -0.64; P < 0.05) and no correlation in TRZ- and VEH-treated cells (S3 Fig and S4 Table). This suggests that samples with lower viability have a higher level of cellular stress as expected given that cardiomyocytes are post-mitotic. Moving forward, we were interested in understanding the direct effects of these drugs on cardiomyocytes prior to the initiation of secondary effects leading to cell death. We used our dose response curve data at 48 hours (S4 Fig), together with the collated clinical serum drug concentration data, and published dose response data for DOX treatment in iPSC-CMs over time [36], to select a dose for deep characterization within the first 24 hours of treatment. We chose a treatment concentration of 0.5 μM, which is below the LD 50 calculated 48 hours post drug treatment for all drugs.
TOP2-inhibiting breast cancer drugs affect iPSC-CM calcium handling To elucidate the impact of cancer drugs on the calcium handling mechanism of cardiomyocytes, a fundamental determinant of cardiomyocyte contraction, we used the Fluo-4 AM fluorescent calcium indicator for real-time imaging and quantification of calcium transients. We randomly selected three individuals and treated these iPSC-CMs for 24 hours with each drug at 0.5 μM and quantified the calcium-associated fluorescence over time using spinning disc confocal microscopy (Fig 2A). We observe a significant decrease in peak amplitude for DNR, EPI and MTX compared to VEH (T-test; P < 0.05; Fig 2B and S5 Table), indicating decreased cytoplasmic calcium entry and therefore decreased contractility. PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 2. Calcium dysregulation occurs in iPSC-CMs post exposure to sub-lethal concentrations of TOP2i. (A) Schematic representation of the features measured to assess intracellular calcium flux in beating iPSC-CMs from three individuals treated with DOX, EPI, DNR, MTX, TRZ and VEH. Fluorescent intensity of the calcium-sensitive fluorescent dye Fluo-4 AM was measured over time by spinning disc confocal microscopy. (B) Mean amplitude of calcium peaks in each individual (orange dot: Individual 2, blue dot: Individual three, green dot: Individual 5) in response to DOX (mauve), EPI (pink), DNR (yellow), MTX (blue), TRZ (dark green), VEH (light green). (C) Mean rising slope of calcium peaks. (D) Mean decay slope of calcium peaks. (E) Mean peak width at half maximum peak height. (F) Mean contraction rate over ten seconds. (G) PCA representing five calcium flux features. TOP2i drugs are represented as triangles and non-TOP2i drugs as circles. Colors represent the specific drug treatment. Asterisk represents a statistically significant change in calcium flux measurements between drug treatments and vehicle (*p < 0.05, **p < 0.01, ***p < 0.001).
https://doi.org/10.1371/journal.pgen.1011164.g002 To estimate the rate of calcium influx and efflux from the cytosol, we examined the rising and decay slopes of the calcium transients. A reduced rate of cytosolic calcium influx during contraction was observed in DOX, EPI, and MTX samples compared to VEH (T-test; P < 0.05; Fig 2C). We observed a more gradual decay slope during the relaxation phase of contraction for DNR, DOX, EPI, and MTX samples compared to VEH (T-test; P < 0.05; Fig 2D). These data therefore suggest potential dysfunction in both cytosolic calcium entry and clearance. We did not observe effects of drug treatment on the duration of calcium transients, as determined by measurements of the full width of the calcium transient at half-maximum fluorescence intensity (Fig 2E). However, we did observe a significant increase in the rate of cardiomyocyte contraction in DNR-treated cells compared to VEH (T-test; P = 0.02; Fig 2F). Notably, contraction rate did not change in response to other drug treatments suggesting potential alterations in cardiac rhythm associated with DNR specifically. In order to determine the overall similarity in drug effects on calcium handling, we performed principal component analysis on the data obtained for the five calcium transient features described above. PC1 which accounts for 56% of the variation, separates the VEH- and TRZ-treated samples from the TOP2i-treated samples (Fig 2G). Taken together, these findings reveal that TOP2-inhibiting drugs at sub-lethal doses can significantly impact calcium handling in cardiomyocytes, which may contribute to cardiomyocyte dysfunction.
AC treatments converge on a shared gene expression response over time We were next interested in determining the gene expression changes that may lead to drug-induced effects on cardiomyocyte contraction and cell viability. To do so, we collected global gene expression measurements by RNA-seq following treatment of 80–99% pure iPSC-CMs with the five drugs at 0.5 μM, and a vehicle control, in six individuals. We assessed two timepoints to capture early (3 hours) and late (24 hours) effects on gene expression. We processed the 72 high quality RNA samples in treatment and time-balanced batches (S5 Fig). Sample processing and sequencing metrics are described in S6 Table. Following sequencing, we aligned reads to the genome (S5 Fig), counted the number of reads mapping to genes and removed genes with low read counts (See Methods), leaving a final set of 14,084 expressed genes. Correlation of read counts followed by hierarchical clustering across samples identifies two distinct clusters, which primarily separates the 24-hour AC treatments from the other samples (S6 Fig). Principal component analysis reveals that PC1 accounts for 29.06% of the variation in the data, and associates with treatment time and type, while PC2, accounting for 15.67% of the variation, associates with the individual from which the samples came (Figs 3A and S7). Together, these results indicate that treatment with TOP2 inhibitors affects global gene expression in cardiomyocytes. PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 3. TOP2i induce global gene expression changes in iPSC-CMs over 24 hours. (A) PCA of RNA-seq-derived expression measurements (log 2 cpm) across 72 samples representing six individuals (1,2,3,4,5,6), the length of exposure (three hours: circles; 24 hours: triangles), and six treatments (DOX (mauve), EPI (pink), DNR (yellow), MTX (blue), TRZ (dark green), VEH (light green)). Data is representative of 14,084 expressed genes. (B) Correlation between drug responses across drugs. The log 2 fold change of expression between each drug treatment and VEH was calculated for all 14,084 expressed genes at each time point. These values were compared across 12 treatment-time groups by hierarchical clustering of the pairwise Pearson correlation values. Color intensity indicates the strength of the correlation. Drug responses are colored by whether the drug is a TOP2i (TOP2i: dark green; non-TOP2i: light green), whether the drug is an AC (AC: dark orange; non-AC: light orange), and the exposure time (3 hours: pink; 24 hours: brown). (C) Comparison of differentially expressed genes in response to each drug treatment after three hours of exposure. n = the total number of differentially expressed genes per drug treatment. N = the total number of differentially expressed genes across drugs. (D) Comparison of differentially expressed genes in response to each drug treatment after 24 hours of exposure. (E) All biological pathways enriched amongst differentially expressed genes in response to drug treatments at three hours. Individual drug treatments as well as genes that respond to all ACs and all TOP2i are shown. KEGG pathways that are significantly enriched in at least one treatment group compared to all expressed genes are displayed (HSV-1 infection and p53 signaling pathways), where color represents the significance of the enrichment across groups and treatment groups where the pathways that are significantly enriched are represented with an asterisk. (F) Top biological pathways enriched amongst differentially expressed genes in response to drug treatments at 24 hours. BER: base excision repair; DNA rep.: DNA replication.
https://doi.org/10.1371/journal.pgen.1011164.g003 We next sought to identify gene expression responses to each drug treatment compared to vehicle at each timepoint using pairwise differential expression analysis. Following treatment for three hours we find tens to hundreds of differentially expressed (DE) genes at 5% FDR for the TOP2i drugs and no response to TRZ (DOX vs VEH = 19; EPI vs VEH = 210; DNR vs VEH = 532; MTX vs VEH = 75; TRZ vs VEH = 0; S8 Fig and S7–S11 Tables). The number of genes differentially expressed in response to TOP2i increases to thousands at the 24-hour timepoint, while there are still no gene expression changes in response to TRZ treatment (DOX vs VEH = 6,645; EPI vs VEH = 6,328; DNR vs VEH = 7,017; MTX vs VEH = 1,115; TRZ vs VEH = 0). Analyzing the overall magnitude of the effect of response to treatment, regardless of a p-value threshold, similarly reveals greatest responses amongst AC-treated cells at 24 hours (S9 Fig). When comparing the magnitude of the effect across all drug treatments and time by hierarchical clustering, we observe three predominant treatment clusters corresponding to the two treatment times of the TOP2i, and a TRZ treatment cluster including both timepoints (Fig 3B). We find a generally low correlation in the individual drug responses over time (DOX 3 vs 24 hour r2 = 0.27; EPI r2 = 0.35; DNR r2 = 0.28; MTX r2 = 0.23). Within the 3-hour treatment cluster, the response to DOX is moderately correlated with EPI (r2 = 0.68), DNR (r2 = 0.61) and MTX (r2 = 0.54). The similarity between drug responses increases at the 24-hour timepoint where the DOX response is highly correlated with EPI (r2 = 0.97) and DNR (r2 = 0.98), and moderately so with MTX (r2 = 0.70). We next compared drug responses within a timepoint by overlapping the sets of DE genes for each drug. After three hours of treatment, we find that out of 566 genes that respond to at least one TOP2i, 1% of genes respond to all four drugs and 2% respond to all three ACs (Fig 3C). After 24 hours of treatment, of the 8,188 genes that respond in at least one treatment, 11% are shared across TOP2i and 54% are shared across ACs suggesting a convergence in the response over time (Fig 3D). While most genes respond similarly across TOP2i after 24 hours of exposure, we were interested in investigating the subset of genes that respond specifically to one drug. There are 356 drug-specific response genes at 3 hours (DOX = 0; EPI = 18; DNR = 322; MTX = 16) and 1,536 at 24 hours (DOX = 355; EPI = 444; DNR = 615; MTX = 122) at 24 hours. However, this approach for identifying drug-specific responses with a significance threshold in each drug treatment is affected by incomplete power to detect responses in a single drug treatment. Indeed, when we investigate the drug-specific response genes we find evidence of signal in the other drug treatments, particularly in the ACs, suggesting an underestimate of the degree of sharing in response to drug treatments (S10A Fig). To identify a high confidence set of drug-specific response genes we selected genes based on the distribution of all p-values for each drug. This resulted in selection of a two-step threshold based on an adjusted p-value threshold of less than 0.01 for the drug of interest, and greater than 0.05 in all other drugs (S10B Fig). This analysis identified 104 drug-specific response genes at 3 hours (DOX = 0; EPI = 0; DNR = 100; MTX = 4; 29% of all response genes) and 305 drug-specific response genes at 24 hours (DOX = 68; EPI = 84; DNR = 112; MTX = 41; 4% of all response genes; S12 Table). Drug-specific response genes at 3 hours include the MAFK transcription factor for DNR, and the transcriptional regulator ZNF547 for MTX, while 24-hour drug-specific response genes include the lncRNA gene ZNF793-AS1 for DOX, the SLC16A9 plasma membrane protein for EPI, the epigenetic regulator gene SMYD4 for DNR, and the rRNA processing gene WDR55 for MTX (S10C Fig). Together, these results suggest that drugs within the TOP2i class induce minimal drug-specific effects. To determine the gene pathways responding to drug treatments, we performed KEGG pathway enrichment analysis on the set of genes DE in response to DOX, EPI, DNR and MTX at both timepoints and found cell cycle, p53 signaling, DNA replication, and base excision repair pathways to be amongst the most enriched pathways across drugs (adjusted P < 0.001 relative to all expressed genes; Fig 3E). p53 signaling and base excision repair pathways are similarly enriched amongst all drugs, while cell cycle and DNA replication genes are most enriched amongst MTX and TOP2i-shared response genes at 24 hours (Fig 3F). This set of genes includes CDKN1A, a p53 response gene and cell cycle regulator. These results corroborate work that investigated the response to 0.05–0.45 μM DOX over time, and identified DNA damage and cell cycle genes to be affected following treatment [37]. The only pathways enriched at three hours are p53 signaling for MTX and Herpes simplex infection for EPI, DNR and MTX. The Herpes simplex infection pathway consists of many genes related to p53 signaling and apoptosis [38], specifically C2H2-type zinc binding domain proteins suggested to be involved in the DNA damage response [39]. To determine if the stringently-identified drug-specific response genes are enriched for biological processes, we performed gene ontology analysis rather than KEGG pathway analysis, given the relatively low number of genes in this set. At three hours we find transcription-related terms to be most enriched amongst DNR- and MTX-specific response genes as well as metabolism-related terms for DNR specifically (the only two drugs that initiate a drug-specific response at this timepoint), while at 24 hours we find enrichment for terms related to cation and calcium channel activity amongst the 68 DOX-specific response genes, and terms related to DNA replication for the 41 MTX-specific response genes (adjusted P < 0.001 relative to all expressed genes; S11 Fig). We further investigated the lack of response to TRZ by using an alternative multiple testing correction approach, and again did not identify any gene expression changes. We identified only 36 genes that pass a nominal p-value threshold of 0.05 following either three or 24 hours of treatment including PHLDA1 and ANKRD2 (S12A and S12B Fig), Nominal TRZ DE genes at three hours are enriched in pathways related to transcriptional regulation in cancer and p53 signaling, while there are no pathways enriched amongst 24 hour response genes (S12C Fig). Pairwise comparisons of multiple drugs and treatment times makes it challenging to determine overall trends in the data, therefore we next sought to jointly model the data to identify gene expression response signatures that may be shared across drugs or specific to a drug.
AC-sensitive chromatin regulators are enriched amongst early response genes To determine the appropriate model for gene expression signature selection, we used Bayesian information criterion and Akaike information criterion analysis. Using this approach we determined that there are four gene clusters that explain the major gene expression patterns across timepoints and treatments (S13A Fig). TRZ treatment did not contribute to any of the gene expression response patterns and had a probability of differential expression of less than 0.1 at both timepoints. We assigned genes to each of the four clusters based on them having a probability > 0.5 of belonging to that cluster (Figs 4A and S13B). Using this approach, we were able to uniquely classify 99.6% of expressed genes into one of the four clusters (S13 Table). We categorized the 7,409 genes, which do not respond to any drug at either timepoint as non-response genes (NR), the 487 genes that respond to the TOP2i drugs only after three hours of treatment as early-acute TOP2i response genes (EAR), the 5,596 genes that respond to the TOP2i drugs only after 24 hours of treatment as late TOP2i response genes (LR), and the 589 genes that respond to the TOP2i drugs at three and 24 hours as early-sustained TOP2i response genes (ESR; S14 Table). There are no clusters driven by individual drugs or the AC drug class at either timepoint. In line with the pairwise differential expression analysis, most genes that respond to drug treatment, respond specifically at the late timepoint (83%; Fig 4B). Early-sustained response genes show a heightened response in ACs at 24 hours (Fig 4C). Late response and early-sustained response clusters show a lower probability of differential expression in the MTX-treated samples at 24 hours (p = 0.3) compared to the AC drugs (p = 1), suggesting a divergence between AC and non-AC treatments over time (Fig 4C). PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 4. Chromatin regulators that mediate breast cancer sensitivity to ACs are enriched amongst early TOP2i response genes. (A) Gene expression motifs identified following joint modeling of test pairs. Shades of black represent posterior probabilities of genes being differentially expressed in response to each drug treatment compared to VEH at each time point. Genes are categorized based on their posterior probabilities: genes with a p > 0.5 in only 24 hour TOP2i drug treatments are designated as ‘Late response genes’ (LR: blue), genes with a p > 0.5 in only three hour TOP2i drug treatments are designated as ‘Early-acute response genes’ (EAR: red), genes with a p > 0.5 in the three and 24 hour TOP2i drug treatments are designated as ‘Early-sustained response genes’ (ESR: green), and genes with a p < 0.5 across all tests are designated as ‘Non-response genes’ (NR: purple). (B) The total number of genes that are assigned to each TOP2i response category. (C) The mean absolute log 2 fold change of each gene in response to each drug at each timepoint within each TOP2i response category. (D) Enrichment of chromatin regulators amongst TOP2i response gene categories compared to the non-response gene category. A curated list of chromatin regulators (n = 408) and chromatin regulators that are sensitive to ACs in breast cancer patients (n = 54) was obtained from Seoane et al. [22]. Enrichment amongst response categories was calculated by a Chi-square test of proportions. Color represents the -log 10 P value for all tests. Significant tests are denoted with an asterisk.
https://doi.org/10.1371/journal.pgen.1011164.g004 Using gene ontology enrichment analysis, we find that early-acute response genes and early-sustained response genes are enriched in biological processes related to transcription and the regulation of transcription compared to all expressed genes, while late response genes are enriched in terms related to mitosis and the cell cycle indicating time-dependent effects of drug treatment on cellular processes (adjusted P < 0.001; S13C Fig). The non-response cluster is enriched in terms related to oxidative phosphorylation and metabolism. ACs inhibit TOP2 proteins and induce damage to both DNA and chromatin. Modified ACs, such as aclarubicin, that induce only DNA damage, are effective anti-cancer agents that show limited cardiotoxicity in mice [23], suggesting that it is the damage to chromatin that affects the heart. In large cohorts of breast cancer patients, and in breast cancer cell lines, chromatin regulator expression predicts response to AC treatment [22]. The authors of this study identified 54 chromatin regulators amongst a curated set of 404 chromatin regulators whose expression level associated with survival in breast cancer patients treated with ACs. We therefore investigated whether these regulators respond to TOP2i in cardiomyocytes by testing for overlap with our TOP2i gene expression response clusters. All three TOP2i response gene clusters are enriched in the full set of chromatin regulators compared to genes that do not respond to any treatment (Chi-square P < 0.05; Fig 4D). Genes that respond early in an acute or sustained manner are also enriched for AC-sensitive chromatin regulators (including KAT6B for example; P < 0.05) compared to genes that do not respond to any treatment. This suggests that chromatin regulators involved in AC sensitivity in breast cancer patients are also involved in mediating the response in cardiomyocytes. Interestingly, the early-acute gene expression signature is enriched in terms for histone modification, histone lysine methylation, histone H3K36 methylation, histone H3K36 dimethylation, regulation of chromatin organization, and heterochromatin organization from gene ontology analysis (adjusted P < 0.05). These terms are not enriched amongst the early-sustained and late response genes.
AC treatment induces transcriptional variation over time To gain insight into inter-individual variation in drug response, we next investigated the variation in gene expression levels of all 14,084 expressed genes across individuals in response to different drug treatments at each timepoint. Mean expression levels are largely stable across treatment groups (Fig 5A and S15 Table). There is a significant reduction in the variation of gene expression levels across individuals in response to all TOP2i drugs compared to VEH at three hours, while TRZ-treated cells show no change in variation (P < 0.05; Fig 5B). After 24 hours of treatment, the AC drugs show increased variation in gene expression levels, while MTX and TRZ show reduced variation (P < 0.005; Fig 5B). Given that all samples were collected and processed in treatment-balanced batches, this suggests that there is a robust response to drug treatments shortly after exposure, but that expression diverges across individuals over time. EPI shows the greatest increase in variance followed by DOX and DNR. To determine whether the increase in AC-induced variability following treatment is evident across a larger panel of individuals, we investigated variance in response to DOX treatment across iPSC-CMs using published data from 45 individuals [15]. We find that 24 hours of treatment with 0.625 μM DOX similarly increases variability in gene expression across individuals compared to vehicle treatment (S14 Fig). We next asked how the variation across individuals changes between each drug and vehicle treatment for each gene by calculating the F statistic. We tested for correlation between the F statistic across drugs and time and observed two distinct sample clusters based on time (Fig 5C). At three hours the variability is highly correlated across all drugs (r2 > 0.47), while at 24 hours the samples cluster based on whether the drugs are ACs or not. EPI-induced variance is distinct from DNR and DOX but is still highly correlated (r2 > 0.52 for both) compared to MTX (r2 = 0.39) and TRZ (r2 = 0.28). These results imply that AC treatment induces variability in expression across individuals after 24 hours of treatment, and these effects are replicated in a DOX-focused study. PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 5. AC treatments induce transcriptional variation across individuals over time. (A) Mean expression levels of all 14,084 expressed genes across individuals for each drug treatment at each timepoint. (B) Variance of gene expression levels of all 14,084 expressed genes across individuals for each drug treatment at each timepoint. (C) The F-statistic was used to test for differences in the variance between each drug treatment and vehicle for all 14,084 expressed genes at each timepoint. These values were compared across 12 treatment-time groups by hierarchical clustering of the Spearman correlation values. Color intensity indicates the strength of the correlation. Drug responses are colored by whether the drug is a TOP2i (TOP2i: dark green; non-TOP2i: light green), whether the drug is an AC (AC: dark orange; non-AC: light orange), and the exposure time (3 hours: pink; 24 hours: brown). Asterisk represents a statistically significant change between drug treatment and vehicle (*p < 0.05, **p < 0.005, ***p < 0.0005).
https://doi.org/10.1371/journal.pgen.1011164.g005
DOX response eQTLs are enriched amongst AC response genes We were next interested in understanding how the drug response genes relate to genes whose expression is known to vary across individuals i.e., eQTLs or eGenes. We collated the set of eGenes in left ventricle heart tissue (q-value < 0.05) identified by the GTEx consortium [40] and selected those that are expressed in our iPSC-CMs (n = 6,261). We first asked whether genes that respond to different drug treatments are likely to vary in their expression level across individuals in the absence of treatment. We find that following 24 hours of drug treatment DOX and DNR response genes are no more likely to be an eGene in heart tissue than not, and EPI and MTX response genes are in fact depleted amongst heart eGenes (Fig 6A). This supports previous studies that suggest that in order to identify genetic effects on gene expression in disease contexts, one must study the relevant context [41]. To do so, we took advantage of a study that investigated the association between genetic variants and the response to DOX in iPSC-CMs treated with a range of DOX concentrations (0.625–5 μM) after 24 hours of treatment [15]. The authors identified 417 baseline eGenes in these cells, as well as 273 response eGenes i.e. eQTLs where the variant associates with the gene expression response to DOX. Reassuringly, we find that genes that respond to 0.5 μM DOX at 24 hours in our study are enriched in response eGenes compared to baseline eGenes (142 of the 273 response eGenes; Chi-square test; P < 0.05; Fig 6B). We also find that EPI, DNR and MTX response genes are enriched in response eGenes compared to baseline eGenes (P < 0.05). Of the 142 DOX response eGenes, 93% respond to DNR, 83% to EPI and 18% to MTX (14% of all DOX response genes respond to MTX), corresponding to 96% of DOX response eGenes responding to at least one TOP2i (Fig 6C). This indicates that most DOX response eGenes respond to the other AC drugs and have the potential to be DNR and EPI response eGenes. Only one DOX response eGene is categorized as one of the 68 high-stringency set of DOX-specific response genes. JPH3, involved in intracellular ion signaling, shows a significant response to DOX but not to any other drug (Fig 6D). PPT PowerPoint slide
PNG larger image
TIFF original image Download: Fig 6. DOX response eGenes are enriched amongst TOP2i response genes. (A) Proportion of drug response genes amongst eGenes in left ventricle heart tissue. Heart eGenes (q-value < 0.05) were obtained from the GTEx database [40] and filtered to retain only those genes expressed in our data. All expressed genes that are not eGenes in heart tissue are denoted as ‘not eGene’. These gene sets were compared amongst genes classified as differentially expressed in response to a particular treatment (DE: light blue) at 24 hours or not differentially expressed (not DE: dark blue). Asterisk represents drug treatments where there is a significant difference in the proportion of DE genes amongst eGenes and non-eGenes. (B) Proportion of drug response genes amongst DOX response eGenes in DOX-treated iPSC-CMs. eGenes identified in iPSC-CMs (base eGenes) and in response to DOX (response eGene) were obtained from Knowles et al. [15]. (C) Overlap between DOX eGenes identified by Knowles et al. that are responsive to DOX in our data, and genes that are differentially expressed in response to at least one other TOP2i. (D) Example of one of the five DOX response eGenes that responds to DOX in our cells but not to any other drug treatment. JPH3 expression levels following 24 hours of treatment are shown. Asterisk represents drug treatments for which this gene is categorized as DE.
https://doi.org/10.1371/journal.pgen.1011164.g006
[END]
---
[1] Url:
https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1011164
Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.
via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/