Differential methylation and expression of genes in the hypoxia-inducible factor 1 signaling pathway are associated with paclitaxel-induced peripheral neuropathy in breast cancer survivors and with preclinical models of chemotherapy-induced neuropathic pain

Survivors and settings

The methods for this analysis, which is part of a larger study of CIPN, are described in detail elsewhere.¹¹ In brief, survivors were recruited from throughout the San Francisco Bay area and met pre-specified inclusion and exclusion criteria. The National Coalition for Cancer Survivorship’s definition of cancer survivor was used in this study (i.e., a person is a cancer survivor from the moment of diagnosis through the balance of life).⁴³ Of the 1450 survivors who were screened, 754 enrolled, and 623 completed the self-report questionnaires and study visit. Data from a randomly selected sample of breast cancer survivors with (n = 25) or without (n = 25) chronic PIPN were used in this analysis.

Study procedures

Research nurses screened and consented the survivors over the phone; sent and asked them to complete the self-report questionnaires prior to their study visit; and scheduled the in-person assessment. At this assessment, written informed consent was obtained, responses to questionnaires were reviewed for completeness, and objective measurements were obtained. Blood samples were drawn, processed, and stored for subsequent molecular analyses in PAXgene® Blood RNA tubes (Qiagen, Venlo, the Netherlands). This study was approved by the institutional review board of the University of California, San Francisco.

Study measures

Biospecimen processing, quantification of methylation status, and quality control

DNA samples from buffy coats were archived using the PUREGene DNA isolation kit (Invitrogen, Carlsbad, CA) as previously described.^52,53 DNA samples were quantitated with a NanoDrop UV spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and normalized to a concentration of 50 ng/μL. The DNA preparation and microarray work were performed at the UC Berkeley Vincent J. Coates Genomics Sequencing Laboratory. DNA was bisulfite converted using the Zymo EZ-96 DNA Methylation Kit (Catalog #D5004) Deep-Well Format (Zymo Research, Irvine, CA) and used as input for the Illumina Infinium HD Methylation Assay (Illumina, San Diego, CA). Processed DNA was dispensed onto the Infinium MethylationEPIC BeadChip and scanned on the Illumina iScan (Illumina). Preliminary analysis and quality control (QC) of the data were performed using GenomeStudio (Illumina). Target success rates were determined. Samples that had <90% of their targets methylated with a detection p-value of ≤0.01 were flagged for review. Sample replicates and Jurkat control replicates were checked to ensure an r² value of >0.99.

Subsequent data analyses were done using well-established protocols in R (version 3.6.1).^54,55 Corrections for Infinium I and II probes, balance correction, background correction, and quantile normalization were performed using the minfi package in R (version 1.30.0).^56,57 Probes that contained a single nucleotide polymorphism at a CpG or flanking site and probes that aligned to multiple places on the genome were excluded.⁵⁸ Methylation scores were quantified as M-values.⁵⁹ Probes were annotated for genes based on mapping to the Genome Reference Consortium Human Build 38 (GRCh38) assembly.⁶⁰

Because DNA methylation differs among blood cell types,^61–63 cell types were estimated using the estimateCellCounts2() function in the FlowSorted.Blood.EPIC R package (version 1.2.0).⁶⁴ Cell type deconvolution was performed using the IDOL L-DMR library for CD8 and CD4 T-cells, natural killer (NK) cells, B cells, monocytes, and neutrophils. Differences between PIPN groups in estimates of cell type composition were evaluated using t-tests and assessed for significance at a p-value of <0.05. Any cell type composition estimates that were significantly associated with PIPN group membership were included as covariates in the final model.

Surrogate variable analysis (SVA, R package version 3.32.1)^65,66 was used to identify technical variations that contributed to heterogeneity in the sample (e.g., batch effects) that were not due to the variable of interest (i.e., PIPN), significant demographic (i.e., age, employment status) and clinical (i.e., AUDIT score, body mass index (BMI), KPS score) characteristics.⁶⁶ SVA can control for cell type heterogeneity in the methylation analysis.^62,63 The “be” method was used to identify surrogate variables.^65,66 Two surrogate variables associated with RIN, extraction date, and cell type composition estimates were included in the analysis.⁶⁷

Overlapping DM and gene expression

To understand the association of PIPN and biological variation in HIF-1 pathway genes, we performed a multi-staged integrated analysis using complementary layers of molecular data.^68,69 In the first stage, DM analysis of probes between cases and non-cases was performed using the limma R package (version 3.40.6).⁷⁰ A linear model was fit using the “ls” method that included significant demographic (i.e., age, employment status) and clinical (i.e., AUDIT score, BMI, KPS score) characteristics, significant cell type composition associated with PIPN (i.e., CD4T cells), and surrogate variables as covariates in the models. We estimates that a minimum of 25 replicates in each group would provide at least 80% power to detect a mean difference in M-values of ≥0.09 at a nominal single locus threshold (type I error rate of 0.01).⁷¹

In the second stage, differential gene expression (DGE) was evaluated for each candidate gene using a generalized linear model with edgeR⁷² as previously described.¹⁴ Briefly, these DGE analyses were adjusted for demographic (i.e., age, employment status) and clinical (i.e., AUDIT score, BMI, KPS score) characteristics that differed between the PIPN groups, as well as for technical variability (e.g., potential batch effects) using SVA. We estimated¹⁴ that, at a type I error rate of 0.01,⁷³ we were powered to detect 1.5-fold changes for 83% of genes across the whole transcriptome.

Overlapping gene expression and methylation probe pairs were identified at the gene level using the annotated Human Genome Organisation (HUGO) approved symbols using the HUGO Gene Nomenclature Committee (HGNC) database.⁷⁴ By using a systems genomics design integrating data from multiple biological sources, we would have increased power to identify and better interpret the omics-phenotype relationships relative to an analysis that used only a single source of omics data.^69,75 For this exploratory study, candidate genes were defined by a significant overlap assessed at a p-value <0.05 for both DM and expression tests. No minimal fold-change was utilized.

Functional analyses of candidate genes

To characterize potential functional roles and potential interactions among candidate genes with differentially methylated probes and DGE beyond the HIF-1 signaling pathway, we evaluated for protein–protein interaction (PPI) network connectivity using the Search Tool for the Retrieval of Interacting Genes (STRING).⁷⁶ The functional enrichment of the PPI was evaluated using the Kyoto Encyclopedia of Genes and Genomes (KEGG)⁷⁷ and Reactome⁷⁸ pathway databases. We assessed for significance of functional enrichment tests using a false discovery rate of 5 under the Benjamini–Hochberg (BH) procedure.^79,80 In addition, we evaluated for similar expression patterns of our candidate genes in various cell types in previously identified peripheral neuroimmune interactions⁸¹ from RNA-seq data sets using the “NIPPY—Neuro-Immune Interactions in the Periphery” database (https://rna-seq-browser.herokuapp.com/).

Differential expression of our candidate genes in DRG of mice with and without PIPN

Differential expression of our candidate genes was evaluated using a publicly available data set generated from mouse dorsal root ganglia (DRG) (NCBI GEO GSE113941). Data were collected from nine pools of mice 10 days after treatment with paclitaxel (n = 5) in a dose that was sufficient to induce PIPN (i.e. “CIPN model”) compared to not treated mice (n = 4) (i.e., “normal”). For our study, FASTQ files of the RNA-seq data from these experiments were downloaded from the NCBI Short Read Archive (GSM3124261, GSM3124262, GSM3124263, GSM3124264, GSM3124269, GSM3124270, GSM3124271, GSM3124272, and GSM3124273).

Following our previously described protocols,¹⁴ these sequences were trimmed with Trimmomatic and aligned to the mouse genome assembly (GRCm38.p6) annotated by GENCODE vM23 with the STAR aligner.⁸² Gene level counts were generated using the featureCounts tool in the subread R package (version 2.0.0).⁸³ QC and DGE analyses were performed as previously described.¹⁴ The “be” method was used to identify surrogate variables.^65,66 No outlier samples were identified and all samples were included in subsequent analyses. Genes with ≤1 read per million in at least four samples were excluded. Two surrogate variables were identified. Neither surrogate variable was associated with neuropathy and both were used as covariates in the final model. DGE was evaluated between the mouse DRG sample pools who did (n = 5) and did not (n = 4) received paclitaxel. These findings were compared to the differentially methylated and expressed candidate genes identified in our cancer survivors. Homologous mouse gene Ensemble ID was determined by the HUGO Symbol using the HGNC Database.⁷⁴ Significant DGE was assessed at a type I error rate of 0.05. No minimal fold-change was evaluated.

DM of our candidate genes in pre-frontal cortex and T-cells of rats with and without SNI

DM of our candidate genes was evaluated in two publicly available data sets generated from rat pre-frontal cortex (PFC) in sham (n = 8) vs. sciatic nerve injury (SNI)⁸⁴ (n = 7) animals (GEO GSE70006) and from T-cells in sham (n = 8) vs. SNI (n = 8) animals (GEO GSE70007) that were measured on a customized array (Aglient Technologies, Santa Clara, CA) nine months after the procedure. The PFC and T-cells were harvested from the same animals. Methylation levels were downloaded from GEO as the quantile normalized log2 ratio of the bound (Cy5) and input (Cy3) microarray channel intensities. DM in both experiments was evaluated using limma. Homologous rat genes were identified from their HUGO symbols using the HGNC Database⁷⁴ and Rat Genome Database.⁸⁵ Significant DM was assessed at a type I error rate of 0.05. No minimal fold-change was utilized.

Differences in demographic, clinical, and pain characteristics

Cancer survivor characteristics were reported previously.¹⁴ In brief, breast cancer survivors with chronic PIPN were significantly older (p = 0.006) and were more likely to be unemployed (p = 0.022) (Supplemental Table 1). In terms of clinical characteristics, survivors with PIPN had a lower AUDIT score (p = 0.012), a higher BMI (p = 0.011), and a lower KPS score (p < 0.001) (Supplemental Table 2). Of note, no between-group differences were found in the number of years since cancer diagnosis, the total dose of paclitaxel received or in the percentage of patients who had a dose reduction or delay due to PIPN. Worst pain severity was 6.3 (±2.1) and duration of PIPN was 3.8 (±3.9) years (Supplemental Table 3).

DM and DGE associated with PIPN in survivors

For the 100 genes that were identified in the KEGG HIF-1 signaling pathway, we mapped 732 probes that had a regulatory feature group classified as “promoter associated” and evaluated for DM. Only CD4+ T-cell composition estimates were associated with PIPN group membership (Supplemental Table 4). Eleven surrogate variables were identified. One was associated with PIPN group membership and four were identified as significantly associated with cell type composition. The final regression model included five significant demographic and clinical characteristics (i.e., age, employment status, AUDIT score, BMI, KPS score) and four surrogate variables. We had methylation and expression data for 81 genes in the HIF-1 signaling pathway that were candidates for evaluation for both DGE and DM. Twelve probes across eight genes were identified as both differentially methylated and expressed between survivors with and without PIPN (Table 1).

Functional analysis identified seven KEGG and three Reactome pathways that were enriched for the eight genes listed in Table 1 (excluding the HIF-1 signaling pathway itself; Table 2). These eight genes represent a significantly enriched PPI network (p = 0.002) and the network identified evidence of protein–protein interactions between six genes (Figure 1).

Differential expression of candidate genes in mouse DRG associated with PIPN

Of the eight candidate genes identified as both differentially methylated and expressed in breast cancer survivors, eight homologs were identified in the mouse DRG data set (Table 3). Of these eight candidates, one gene (i.e., Mknk1) was found to be differentially expressed between PIPN and “normal” pools of mice DRG.

DM of candidate genes in rat PFC and T-cells associated with SNI

Of the eight candidate genes with overlapping DM and DGE in breast cancer survivors, seven homologs were identified in the rat PFC and T-cell data set (i.e., Mknk1, Map2k2, Egln1, Rbx1, Pfkl, Cul2, Ldha). For these seven genes, we identified 72 probes (Egln1 n = 14, Map2k2 n = 10, Mknk1 n = 8, Pfkl n = 12, Rbx1 n = 5, Cul2 n = 12, Ldha n = 11) in the promoter regions and evaluated them for DM. Of these seven candidates, three genes (i.e., Mknk1, Ldha, Egln1) were differentially methylated between SNI and sham rats in both PFC and T-cells (Table 4). While no genes were differentially methylated only in T-cells, two genes (i.e., Cul2, Map2k2) were differentially methylated only in the PFC.

HIF-1 signaling pathway genes associated with DM and DGE

The mitogen-activated protein kinase (MAPK) 1 interacting serine/threonine kinase 1 (MKNK1) gene is involved in the control of the cellular proteome⁹⁰ and in inflammatory responses.⁴¹ MAPK interreacting protein kinases (Mnks) are broadly expressed across different tissues types. In terms of neuroinflammation, Mnks play critical roles in cytokine receptor signaling.⁴¹ MKNK1 mediates cytokine production in macrophages through post-transcriptional regulation and can mediate TNF-α translation.⁴¹

In terms of effecting translation, MKNK1 phosphorylates the eukaryotic initiation factor 4E (eIF4E) and is activated downstream of the p38 or MAPK pathways.⁴¹ The regulation of the eIF4E complex is a key signaling pathway that controls the cellular proteome.⁹¹ We found that the MKNK1 gene was differentially methylated and expressed in peripheral blood from breast cancer survivors with PIPN, differentially methylated in PFC and T-cells of rats with SNI, and differentially expressed in DRG of mice with PIPN. Recent work suggests that changes in Mknk1 expression and phosphorylation of elF4E in mice contribute to nociceptor plasticity^92–94 and that inhibition or elimination of Mkn1 (Mknk1) attenuates PIPN in a murine model.⁹³ In addition, inhibition of Mkn1 and Mnk2 reduces spontaneous pain in a SNI murine model.⁹⁵ Interestingly, the phosphorylation of elF4E by Mknk1 was identified in a neurological model animal, Aplysia californica, which suggests that this regulatory mechanism for neuroplasticity is highly conserved⁹⁶ across Bilateria.⁹⁷

In terms of neuroinflammation, Mknk1 is expressed in both DRG and peripheral blood and across peripheral and neuroimmune cell lines (Supplemental Figures 1 and 2) which is consistent with our findings of expression across tissue types. Compared to survivors without PIPN, MKNK1 had lower expression and higher methylation in survivors with PIPN. In contrast to our findings in survivors, Mknk1 had higher expression in mice with PIPN compared to normal mice. One possible explanation for this difference is the timing of the measures. For the mice with PIPN, gene expression was measured 10 days after CTX compared to approximately four years in our cancer survivors (Supplemental Table 2). Similarly, in a recent preclinical study⁸¹ of traumatic nerve injury, findings from macrophage and T-cells dissected from DRG suggests that Mknk1 expression levels decline 10 weeks after sciatic nerve injury (Supplemental Figure 4). Based on the fact that inflammation persisted unresolved for 3.5 months (the duration of the study) in these animals, the authors argued that the categorization of pain as either inflammatory or neuropathic may not be mechanistically appropriate. Given that DNA methylation reprogramming was observed in preclinical models of chronic neuropathic pain following nerve injury,³⁶ future research should evaluate for changes in the methylation and expression levels of Mknk1 over time.

MAP2K2 is a kinase which phosphorylates and activates MAPK1 and MAPK3.⁹⁸ The MAPK1 gene is involved in eukaryotic signal transduction⁹⁹ and is associated with the development of neuropathic pain.^100,101 We found that MAP2K2/Map2k2 was differentially methylated and differentially expressed in peripheral blood from breast cancer survivors with PIPN and differentially methylated in PFC of rats with SNI. ERK2 (MAPK1) activation by the three MAPK family pathways (i.e., phosphoinositide 3-kinases (PI3Ks), adenosine monophosphate-activated protein kinase (AMPK), mammalian target of rapamycin (mTOR)) promotes inflammatory and neuropathic pain conditions (reviewed in Wood et al.¹⁰²). Evaluations of the effects of paclitaxel on cancer cell lines found that paclitaxel activated these MAPK family pathways.^103–105 Thus, the regulation of MAPK1 was an early therapeutic target for the treatment of neuropathic pain induced by peripheral nerve injury.¹⁰⁶ More recently, inhibition of a number of MAPK family pathways (e.g., p38, MAPK, and JNK) was shown to attenuate neuropathic and inflammatory pain in preclinical models (reviewed in Ji et al.¹⁰¹). While no studies have examined associations between PIPN and changes in gene expression or methylation of MAP2K2, expression of Map2k2 increased in rat L4 and L5 DRGs cells with sciatic nerve lesions as compared with sham controls three days after injury.¹⁰⁷

The transferrin receptor (TFRC) gene encodes for a protein important for cellular iron uptake and homeostatis.¹⁰⁸ TFRC (TFR1) is a key iron metabolism gene that fine tunes cellular iron access, storage, and utilization.¹⁰⁹ Iron overload can lead to increased oxidative stress¹¹⁰ and may be regulated by HIF-1 TFs.¹¹¹ We found that TFRC was differentially methylated and expressed in peripheral blood of breast cancer survivors with PIPN. However, data were not available to test for DM in T-cells or PFC of rats with SNI or DGE in DRG of mice treated with paclitaxel. While no studies have examined associations between PIPN and changes in gene expression or methylation of TFRC, polymorphisms in TFRC decreased the risk for distal neuropathic pain in HIV-infected patients on combination antiretroviral therapy.¹¹² In addition, iron overload is a risk factor for diabetic neuropathy.¹¹³ In a preclinical central pain model, mRNA expression of TFR1 (an alias for TFRC) was upregulated in activated microglia after spinal cord injury.¹¹⁴ Finally, our previous finding of perturbations in iron homeostasis pathways (i.e., Ferroptosis) in this sample¹⁴ suggests that future studies should investigate the role of TFRC in the regulation of iron homeostasis in survivors with PIPN.

The phosphofructokinase, liver type (PFKL) gene codes for the liver subtype of an enzyme that catalyzes the phosphorylation of D-fructose 6-phosphate to fructose 1,6-bisphosphate. Therefore, PFKL is a key metabolic gene in glycolysis.¹¹⁵ We found that PFKL was differentially methylated and expressed in peripheral blood from breast cancer survivors with PIPN. However, it was not DM in T-cells or PFC of rats with SNI or differentially expressed in DRG of mice treated with paclitaxel. No studies were found that evaluated for associations between gene expression or methylation of PFKL and CIPN. However, in a preclinical model of diabetic neuropathy,¹¹⁶ PFKL was induced in IMS32 Schwann cells following chronic treatment with high levels of glucose (>8 weeks). Given that oncology patients with diabetes are at increased risk for developing CIPN,¹¹⁷ future studies should evaluate this mechanism in patients with both comorbidities.

The lactate dehydrogenase A (LDHA) gene catalyzes the forward and backward conversion of lactate to pyruvate in glycolysis. LDHA is activated by HIF^118,119 and the reduction of LDHA can induce oxidative stress.¹²⁰ We found that LDHA was differentially methylated and expressed in peripheral blood of breast cancer survivors with PIPN and differentially methylated in both T-cells and PFC of rats with SNI. Although no studies were found on associations between gene expression or methylation of LDHA and PIPN, LDHA-driven aerobic glycolysis was associated with pain in bortezomib-induced CIPN.¹²¹ Given that we found an enrichment in the glycolysis/gluconeogenesis KEGG pathway for these eight candidate genes, future research should investigate the role of LDHA in oxidative stress and cellular respiration in survivors with PIPN.

The Egl nine homolog 1 (EGLN1) gene produces a protein which acts as a cellular oxygen sensor that catalyzes the formation of 4-hydroxyproline in HIF alpha proteins and hydroxylates HIF1A and HIF2A. Then, hydroxylated HIFs are degraded through a ubiquination complex. The ring box 1 (RBX1) and cullin 2 (CUL2) genes encode for proteins that act as an E3 ubuquitin-protein ligase, which are involved in mediating the ubiquitination and degradation of proteins. We found that EGLN1/Egln1 was differentially methylated and expressed in peripheral blood of breast cancer survivors with PIPN and differentially methylated in both T-cells and PFC of rats with SNI. In addition, RBX1 and CUL2 were differentially methylated and expressed in peripheral blood of breast cancer survivors with PIPN. However, although Cul2 was differentially methylated in PFC of rats with SNI, we did not find DM of Rbx1 in T-cells or PFC of rats with SNI. In addition, we did not find Rbx1 or Cul2 to be DGE in DRG of mice treated with paclitaxel. In a pharmacological network-based analysis of drug-induced peripheral neuropathy (including Paclitaxel),¹²² Cul2 was identified as a highly connected significant intermediator between drugs and their pharmacological targets. Although no studies were found on associations between gene expression or methylation of RBX1 and neuropathy, in oxaliplatin-treated embryonic kidney cells, the overexpression of PHD2 (EGLN1) inhibits cold-induced hTRPA1 activation.¹²³ Future studies should evaluate this mechanism in patients having received taxanes, platinums, or both.

Epigenetic regulation of transcription

Translational control is an essential component in the regulation of gene expression¹²⁴ and may be a core mechanism for the development and maintenance of persistent pain.¹²⁵ Translational control of molecular processes associated with neuropathy allows for localized translation outside of the cell body in the axons^126,127 (i.e., activity-dependent regulation of mRNA translation¹²⁸) including the regulation of the eIF4E complex by MNKN1. However, the Mnks are effectors of other biological processes (i.e., receptor tyrosine kinase activity, TNF-α mRNA translation, arachidonate release).⁹⁰

Alternative splicing is an important step in post-transcriptional regulation of gene expression¹²⁹ and the diversity of various transcript isoforms of a gene are associated with environmental perturbations.¹³⁰ Alternative splicing is known to occur in the nucleus either during (co-transcriptional splicing) or after transcription (post transcriptional splicing)¹³¹ and DNA methylation can play a role in the regulation of alternative splicing.^132,133 Recent work describing variations in splicing across immune cells found that the distribution of isoforms was on/off, which suggests switch-like regulatory control.¹³⁴ Although not currently identified in CIPN, in a model of pre-diabetic polyneuropathy, errors in splicing factors were associated with neuronal dysfunction.¹³⁵ In addition, alternative splicing appears to play an important role in the development and function of the nervous system,¹³⁶ including axon guidance.¹³⁷ In terms of PIPN, two isoforms of the MKNK1 gene can occur through alternative splicing (i.e., MNK1a and MNK1b). However, the MKN1b isoform is not activated by either p38 MAPK or MAPK1.^138,139 This finding suggests that the regulation of MKN1b occurs through other mechanisms.⁴¹ Given that gene expression can be modulated at multiple levels, from chromatin folding to mRNA translation, future research should continue to evaluate the relative contributions of epigenetic regulation of translation to the development and maintenance of PIPN.