Quantification of spatial pharmacogene expression heterogeneity in breast tumors

Abstract Background Chemotherapeutic drug concentrations vary across different regions of tumors and this is thought to be involved in development of chemotherapy resistance. Insufficient drug delivery to some regions of the tumor may be due to spatial differences in expression of genes involved in the disposition, transport, and detoxification of drugs (pharmacogenes). Therefore, in this study, we analyzed the spatial expression of 286 pharmacogenes in six breast cancer tissues using the recently developed Visium spatial transcriptomics platform to (1) determine if these pharmacogenes are expressed heterogeneously across tumor tissue and (2) to determine which pharmacogenes have the most spatial expression heterogeneity. Methods and Results The spatial transcriptomics technology sequences the transcriptome of 55 um diameter barcoded sections (spots) across a tissue sample. We analyzed spatial gene expression profiles of four biobank‐sourced breast tumor samples in addition to two breast tumor sample datasets from 10× Genomics. We define heterogeneity as the interquartile range of read counts. Collectively, we identified 8887 spots in tumor regions, 3814 in stroma, 44 in lymphocytes, and 116 in normal regions based on pathologist annotation of the tissues. We showed statistically significant differences in expression of pharmacogenes in tumor regions compared to surrounding non‐tumor regions. We also observed that the most heterogeneously expressed genes within tumor regions were involved in reactive oxygen species (ROS) handling and detoxification mechanisms. GPX4, GSTP1, MGST3, SOD1, CYP4Z1, CYB5R3, GSTK1, and NAT1 showed the most heterogeneous expression within tumor regions. Conclusions The heterogeneous expression of these pharmacogenes may have important implications for cancer therapy due to their ability to impact drug distribution and efficacy throughout the tumor. Our results suggest that chemoresistance caused by expression of GPX4, GSTP1, MGST3, and SOD1 may be intrinsic, not acquired, since the heterogeneity is not specific to chemotherapy‐treated samples or cell type. Additionally, we identified candidate chemoresistance pharmacogenes that can be further tested through focused follow‐up studies.

cell type. Additionally, we identified candidate chemoresistance pharmacogenes that can be further tested through focused follow-up studies.

K E Y W O R D S
chemotherapy, pharmacogene, resistance, spatial, transcriptomics 1 | BACKGROUND Heterogeneity is one of the hallmarks of cancer 1 and is associated with resistance to treatment, poor prognosis, and treatment failure. [2][3][4] Intratumor heterogeneity exists at multiple levels including genetic, epigenetic, metabolic, and transcriptional. [5][6][7][8][9][10] Spatial transcriptomic heterogeneity, that is, a wide range in gene expression between neighboring regions of a tissue, can lead to phenotypic variability in traits such as cellular morphology, 11 growth, 12 metabolism, 13 immune function, 14 angiogenicity, 15 reactive oxygen handling, 16,17 and solute transport. 18 These phenotypic changes that result from transcriptionally distinct regions of the tumor (due to spatial heterogeneity) are significant because they are subsequently able to confer resistance to anti-cancer chemotherapies. 19 Success of systemic anti-cancer chemotherapy is dependent upon the ability to deliver cytotoxic concentrations of the drug to all cells within a tumor. Drug concentrations vary across different regions of tumors 20 and is thought to be involved in development of resistance. 21 Insufficient drug delivery to some regions of the tumor may be due to spatial differences in expression of genes involved in the disposition and transport of drugs (pharmacogenes). Drug transporters such as the solute-carrier protein (SLC) and the ATP-binding cassette (ABC) family are membrane proteins that transport a wide variety of molecules, including chemotherapeutic agents, and have been linked with multidrug resistance in cancers. 22,23 This mechanism of resistance is thought to involve survival of certain tumor cells with strong drug efflux mechanisms, [24][25][26][27] which seed resistant cancer cells that clonally expand and regenerate the tumor after chemotherapy treatment. 28 Transporters are not the only mechanism of chemotherapeutic resistance. For example, several chemotherapeutic agents alter the redox homeostasis of cancer cells by elevating levels of reactive oxygen species (ROS) and inducing ROS mediated cell injury. 29,30 Aberrant expression of ROS handling enzymes, including detoxifying enzymes such as superoxide dismutase (SOD) 31 and glutathione peroxidase (GPX) 32 can also promote development of regions across the tumors that are resistant to chemotherapeutic agents. Spatial expression of other pharmacogenes involved in the pharmacokinetics (PK) of chemotherapeutic agents, like the cytochrome P450 (CYP) enzymes, could also be involved in spatial differences in tumor drug concentrations. Therefore, in this study, we analyzed the spatial expression of 286 pharmacogenes in 6 breast cancer tissues using the recently developed Visium spatial transcriptomics platform 33 to determine if these pharmacogenes are expressed heterogeneously across tumors and to determine which pharmacogenes have the most spatial expression heterogeneity, as defined by interquartile range of read counts. The long-term goal of this research is to better understand whether and how spatial heterogeneity in pharmacogene expression impacts chemotherapy resistance. This study is the next step toward this goal and is designed to identify candidate chemoresistance pharmacogenes that have the most variation in expression from region to region (1) across combined samples and within tumor or tumor-adjacent breast tissue, and (2) between samples in tumor regions only. More specifically, the purpose of this study is to enable future studies to focus on a few candidate chemoresistance genes (narrowed down from 286 pharmacogenes) that can be perturbed in model systems of tumorigenesis and tested for their contribution to tumor sensitivity and resistance development to a single chemotherapeutic agent.

| Tissue collection
Tissues were obtained from the Indiana Biobank under an approved IRB protocol at Indiana University. Consent for research and publication of de-identified research information was provided by subjects of the Indiana Biobank. Surgically resected breast tumors were weighed and divided into $150 mg tissue pieces that were flash frozen in liquid nitrogen and placed in a cryovial. The tissues were bio-banked and stored in liquid phase of nitrogen until experimental use.

| Cryosectioning
A cryomold was placed in an ethanol-dry ice slurry. Tissue was transferred from the cryovial on dry ice to the cryomold using pre-cooled forceps. The pre-cooled Optimum Cutting Temperature (OCT) compound was added to the cryomold, completely covering the tissue, and was allowed to freeze for approximately 2 min. The frozen tissue block was removed and placed on a cryotome chuck using OCT and immediately transferred to a pre-chilled cryotome chamber (Leica CM3050 S). Once frozen, the cryotome chuck was placed on the stage and the tissues were cryosectioned at a thickness of 14-16 μm. The chamber temperature was maintained at À20 C and the specimen head was maintained at À28 to À32 C. Tissue sections were placed within the frames of a pre-chilled Visium Spatial slide.
The slide was transferred to a prechilled slide box and stored on dry ice.

| Visium spatial tissue optimization and library preparation
Visium Spatial Tissue Optimization was done as per manufacturer's instructions using a breast tumor tissue sample sectioned at 16 μm.
Based on the tissue optimization, the optimal permeabilization time was determined to be 12 min. Imaging was done using a Keyence BZ-X microscope. Visium Spatial library preparation was done as per manufacturer's instructions. H&E-stained tissue images were captured at 10Â magnification. Imaging time was 14 min per slide. cDNA and library preparation were done at the Center for Medical Genomics core at Indiana University School of Medicine. Visium-prepared sequencing runs and to convert barcode and read data to FASTQ files. Spaceranger count was then used to combine the FASTQ file and optical H&E stained tissue image, generate feature-spot matrices, determine clusters, and perform initial gene expression analysis. A Loupe Browser file was created for an interactive visualization functionality.

| Data analysis
Broad descriptions of the data analysis are given here, and the scripts used to analyze the data are provided in the supplement.
We used R (version 3.6.0) for all analyses. First, pathologist annotated cluster assignments (barcode annotation index) were downloaded from the Loupe Browser and read into R along with the full unique molecular identifier (UMI) read count matrix for each sample. Data was either normalized to (divided by) the total reads per spot or kept as is, depending on the analysis. Data that was normalized was done so similar to the transcripts-per-million (TPM) method, 34 without the unnecessary multiplication by 1 million. Log transformation was not performed on read counts for certain analyses (as indicated in the results section) to yield a more intuitive interpretation of the data. Log2 transformation was performed for fold-change analysis. For differential expression analysis, read counts that were zero in one group and nonzero in another were retained in the analysis as positive or negative "infinity," and shown above or below the solid horizontal line in the respective figure. While it is not inappropriate to discard read counts of zero, we kept them since these values are indicative of low expression and therefore are more informative to retain (e.g., a gene with four reads is likely more highly expressed than a gene with zero reads). Pharmacogenes were filtered based on a list of 298 pharmacogenes as described by www.pharmaadme.org, resulting in 286 pharmacogenes that were present in the read count matrix. The pharmacogene list is based on input from seven major pharmaceutical companies as to which genes perform or regulate drug metabolism or transport. Spots were categorized by sample or by annotation as needed for each analysis. Interquartile range of read counts (75th-25th percentile) was used to measure heterogeneity because it best characterizes the spread (variability) of data without being too heavily influenced by outlier values. Additionally, it allowed better visual representation than standard deviation or coefficient of variation. Spatial serial correlation (autocorrelation) was not used since our aims were focused on the spread of expression values, not the correlation of expression between neighboring spots.
To define which genes were related to reactive oxygen species handling (ROS), we used MGI's Batch Query (http://www.informatics. jax.org/batch/summary) to categorize pharmacogenes corresponding to ROS GO terms using key words "oxidative," "stress," "oxygen," and "reactive," followed by removal of CYP and ABC genes, to yield a final list of 35 "ROS genes" (corresponding to 27 ROS GO terms that can be found in Table S1).

| Summary of tumor tissues
We analyzed spatial gene expression profiles of four biobank-sourced breast tumor samples in addition to two breast tumor sample datasets from 10Â Genomics. The patient demographics and disease characteristics of the biobank samples are listed in Table S2.
Spatial gene expression patterns were generated from a total of 13 600 spots across the six tissues with 27 542 genes detected. Figure S1 shows that pharmacogenes were evenly distributed across the relative expression range ( Figure S1A) and that many of the CYP and SLC genes were generally expressed at lower levels ( Figure S1B).
Pharmacogenes involved in ROS handling (ROS genes) and the ABC transporters are generally seen toward the top half of relative expression among the pharmacogenes.   In addition to heterogeneous expression, we observed statistically significant differences in the overall expression of pharmacogenes in F I G U R E 1 Boxplot of all pharmacogenes, across all samples, in tumor regions only, with an interquartile range greater than zero. Each data point on the y-axis is the number of unique-UMI reads from a single barcoded spot. The genes on the X-axis are sorted by interquartile range in descending order tumor regions compared to adjacent regions. Figure 4A illustrates differential expression of pharmacogenes between tumor regions and adjacent non-tumor regions. We noticed pharmacogenes tend to be downregulated in tumor regions as can be seen in Figure 4B, however there were still many significantly upregulated genes. Figure 4C shows the top 35 most differentially expressed genes, with ADH1B and GPX3 being significantly downregulated in all six tissue samples. Figure 4D,E,F show subsets of differential expression for ROS genes, ABC transporter genes, and SLC transporter genes, respectively.   Table S3. We found a predominance of genes involved in reactive oxygen species handling in addition to drug and metabolite transport.

| DISCUSSION
We analyzed the spatial expression patterns of pharmacogenes in six human breast tumor samples (including tumor and adjacent tissue) Our study evaluated spatial heterogeneity in tumor pharmacogene expression but did not evaluate the downstream impact of this heterogeneity. Therefore, our results should be regarded as hypothesis-generating and not as direct evidence behind chemoresistance mechanisms. However, the literature supports the hypothesis that spatial heterogeneity in tumor pharmacogene expression contributes to chemotherapeutic resistance mechanisms. For example, certain drug and metabolite transporters, such as those included in the multidrug resistance protein (MRP) group of the large ABC family, are known to play a role in chemoresistance when overexpressed as they contribute to the efflux of chemotherapeutic compounds from the cell. 35 One of these transporters, ABCC5, was among our most heterogeneously expressed pharmacogenes. ABCC5 transports cyclic nucleotides, including the metabolites of 5-fluorouracil (5-FU), a common anticancer agent used in both breast cancer and colon cancer treatment. ABCC5-transfected cells have been shown to have a nearly 9-fold increase in 5-FU resistance. 35 Another transporter in this same family also found to be heterogeneously expressed is ABCB8; it transports compounds from the mitochondria to the cytosol. One study showed inhibition of ABCB8 with shRNA resulted in increased doxorubicin-induced mitochondrial DNA damage. 36 Combined with literature knowledge that these genes can affect 5-FU (or capecitabine) and doxorubicin resistance, respectively, we are now able to hypothesize that large degrees of expression heterogeneity (e.g., in ABCC5 and GPX4 has been shown to be critical for survival of lapatinib-resistant cancer cells, but not lapatinib-naïve cancer cells. 32 In colorectal cancer, increased GPX3 expression resulted in increased resistance to oxaliplatin and cisplatin. 37 This evidence, combined with our detection of spatial heterogeneity in GPX3 and GPX4 expression within tumors, could explain why certain cells survive platinum-agent or lapatinib treatment, leading to resistant tumors. GSTP1 was also highly heterogeneously expressed in our data. Enzymes in the glutathione s-transferase (GST) family catalyze the conjugation of polar glutathione groups that enhance systemic elimination of chemotherapeutic agents and toxic metabolites. GST activity has been shown to be inducible by treatment with vincristine, doxorubicin, or topotecan. 38 When GSTP1 enzymatic activity is impaired, as is the case with the rs1695 missense variant, 39 platinum-based chemotherapy-induced granulocytopenia was shown to be more common in a meta-analysis of 12 case control trials. 40 Additionally, GSTP1 expression was found to be higher in adriamycin-resistant cells, and higher GSTP1 expression was also found in breast cancer tissues from subjects with progressive/stable disease versus those with partial/ complete response. Interestingly, this finding was true for tissue collected before and after anthracycline/taxane treatment, indicating intrinsic mechanisms of resistance. These data indicate that spatial differences in GSTP1 expression could be involved in chemoresistance and seeding of resistant tumor cells following chemotherapy treatment.
MGST3 was another highly heterogeneously expressed gene in our data. This enzyme is involved in immune function by catalyzing the conjugation of reduced glutathione and leukotriene A4, producing leukotriene C4. Overexpression of MGST3 was found in cisplatin resistant lung adenocarcinoma cells compared to non-resistant progenitor cells, and when MGST3 expression was increased (via antagonism with its microRNA regulator, mir-432-5p) the progenitor cells demonstrated increased survival to cisplatin treatment. 41 Although it is unclear if this relates to immune modulation, these data and other similar studies 42,43 indicate that spatial differences in MGST3 expression could be involved in chemoresistance and seeding of resistant tumor cells following chemotherapy treatment.
Another highly heterogeneously expressed gene in our data was superoxide dismutase 1 (SOD1). SOD1 eliminates damaging superoxide radicals by converting them to less toxic molecular oxygen and hydrogen peroxide. SOD activity counteracts superoxide-induced autophagy 44 and inhibition of SOD1 with siRNA or small molecule inhibitors results in increased cisplatin sensitivity in ovarian cancer cells. 45,46 These studies combined with our results suggest spatial differences in SOD1 expression could be involved in chemoresistance and seeding of resistant tumor cells following chemotherapy treatment.
We found that reactive oxygen species handling (by gene families like GPX, GSTP, and SOD), is a pathway that is overrepresented among the pharmacogenes that were heterogeneously expressed in our data. The over-representation analysis (ORA) approach identified the top pathways associated with our dataset's most heterogeneously expressed genes.
Our study also assessed differential pharmacogene expression between tumor versus tumor-adjacent regions. This analysis was secondary to the goal of characterizing pharmacogene expression heterogeneity, but does provide some insight into the effect of pharmacogene regulation in tumor regions. Interestingly, ADH1B and GPX3 were downregulated in all six samples. As discussed earlier, GPX3 may play a role in resistance to platinum agents, and while it is not clear how, tumor-induced downregulation could somehow be involved. A potentially more relevant conclusion is that tumor regions do not appear to consistently induce expression of pharmacogenes that could predispose to chemotherapeutic resistance.
We acknowledge some additional limitations to our study. First, we have a small sample size and more information will likely be gained from bigger studies. Second, the lower expressed genes may contain additional variability that we did not observe due to the low read numbers. Third, the resolution between spots is not fine enough to capture single cells, so there may be some overlap in cell types. Lastly, we were not able to identify somatic mutations that may account for variability in pharmacogene expression.

| CONCLUSIONS
Our data show substantial heterogeneity in the expression of many pharmacogenes across areas of the breast tumors. These results provide more quantitative measurements of expression heterogeneity across spatial regions of tumors. Our results provide more detail on the current understanding of chemoresistance development by providing evidence that there is heterogeneity in the expression of these chemoresistance genes across tumor sub-regions. Our results suggest that chemoresistance caused by GPX4, GSTP1, MGST3, and SOD1 may be intrinsic, not acquired, since the heterogeneity is not specific to chemotherapy-treated samples or cell type. Additionally, we demonstrate the utility of spatial transcriptomics to identify candidate chemoresistance pharmacogenes that can now be further tested through focused follow-up studies.

CONFLICT OF INTEREST
Joseph Ipe began working at 10Â Genomics after this study was completed, and while no conflict of interest existed during the study, this statement is provided for full disclosure. None of the other authors have potential conflicts of interest related to this work.

DATA AVAILABILITY STATEMENT
Data files are available upon reasonable request from the corresponding author.

ETHICS STATEMENT
Tissues were obtained from the Indiana Biobank under the approved IRB protocol (number 1803676586) at Indiana University.

CONSENT FOR PUBLICATION
Consent for research and publication of de-identified research information was provided by subjects of the Indiana Biobank.