Cancer Exacerbates Chemotherapy-Induced Sensory Neuropathy

Chemotherapy can achieve high rates of survival in patients with common cancers (1) but often causes severe side effects, including neuropathy, that can limit its use (2, 3). Debilitating sensory disorders, including pain, paraesthesia, and somatosensory loss, reduce quality of life for many patients and can persist for months or years after discontinuing chemotherapy (4, 5). These neurologic disorders occur for up to 80% of patients receiving commonly used antineoplastic agents, notably antitubulins and proteasome inhibitors, as well as platinum-based compounds (5, 6), which are the prescribed adjuvant treatment in 50% of cancer cases worldwide (7, 8). There are few options available for the prevention or treatment of sensory disorders and those that exist largely focus on symptomatic management not prevention (9). Thus, there is an urgent need to better understand the pathogenesis of sensory dysfunction to aid development of mechanism-based therapies.

Preclinical studies of chronic neuropathy in experimental models of cancer are unavailable, possibly because of presumption that chemotherapy alone is sufficient to explain the neuropathology. A recently published meta-analysis identified 341 preclinical studies of chemotherapy induced neuropathy, and none assessed interaction between chemotherapy and cancer (10). Omitting cancer from preclinical study produces a fundamental gap in knowledge that may explain why treatments for neuropathic side effects of chemotherapy have been unsuccessful in patients with cancer. Cancerous tumors located outside the nervous system can induce cognitive disability and peripheral nervous system dysfunction independent of chemotherapy, surgery, or comorbidities (11, 12). Moreover, cancer provokes dysregulation in immune, metabolic, oxidative, and neuronal excitability (13, 14), all of which are identified as targets through which chemotherapy produces neuropathy (2). Convergence of cancer and chemotherapy on the same biological processes seems likely to yield nonlinear interactions, leading us to hypothesize that clinically relevant neuropathy emerges from codependent actions of cancer and chemotherapy.

We tested our hypothesis by uncoupling the independent and combinatorial effects of cancer and chemotherapy on sensory neurons in studies of rats that are not feasible in human patients. Here we present original findings to show that sensory neuropathy depends on complex systemic interactions between cancer and chemotherapy. By combining cancer and chemotherapy, we reproduce the relevant clinical condition in preclinical study of rats to produce the first transcriptional profile of nervous tissue, to discover recurrent metabolic reprogramming and perturbed inflammatory and ion channel processes distinct from the independent effects of cancer or chemotherapy. Parallel codependency is found in novel neuropathic responses of mechanosensory neurons, both in their IHC for selected voltage-gated potassium ion channels and in electrophysiologic measures of encoding naturalistic stimuli in vivo. Collectively, our data show for the first time that sensory neuropathy emerges from codependencies between the systemic effects of cancer and chemotherapy. Apart from identifying new molecular targets for restoring function to mechanosensory neurons, our findings provide the fundamental insight that accounting for chemotherapy–cancer combinatorial effects on the somatosensory nervous system is prerequisite for developing meaningful treatment or prevention of neuropathy.

All procedures and experiments were approved by the Georgia Institute of Technology Institutional Animal Care and Use Committee. Adult male rats (250–350 g) were studied in terminal experiments only and were not subject to any other experimental procedures. All animals were housed in clean cages and provided food and water ad libitum in a temperature- and light-controlled environment in Georgia Institute of Technology's Animal facility.

We studied the rat model of colorectal cancer developed in Fisher 344 (F344) rats and known as polyposis in rat colon (Apc^Pirc/+; ref. 15). Pirc rats generated from a germline mutation in the Apc gene (Apcam^1137/+) approximate human colorectal cancer. We studied age-matched male rats (4 months). This age has been previously validated to result in 100% incidence of animals (16) developing colon cancer and was consistent with our population (Supplementary Fig. S1). Apc^Pirc/+ rats exhibit regional distribution of tumors that is best seen in human colorectal cancer (16), histopathology and morphology closely resembling human tumors (adenocarcinomas; ref. 17), upregulation of well-established progrowth factors, for example, β-catenin and EGFR, and proliferative markers (Ki-67) seen in human carcinogenesis (Supplementary Fig. S1; ref. 15), dramatic increases in Wnt signaling believed to be a key factor in human colorectal carcinogenesis, presence of apoptotic-resistant cells unique to carcinogenesis (15, 16), oxaliplatin induction of apoptosis in colon tumors, survivability (12–15 months) enabling long-term investigation of oxaliplatin's chronic effects (15).

Oxaliplatin was injected intraperitoneally once a week (10 mg/kg, 1 mL 5% dextrose in DMSO) in age-matched rats (4 months) to achieve a cumulative dose of 70 mg/kg over 7 weeks (Fig. 1A), which scales to a human dose of 420 mg/m² (conversion based on rat body surface area and Km = 6; ref. 18). This dose minimizes nerve degeneration in patients (19, 20). Throughout the entirety of the experimental timeline, rats were visually and physically (palpation) inspected for pain, distress, or signs of local infection. No individual rat reached set criteria established for early removal from the study: greater than 10% weight loss, pain, or distress assessed by vocalization, absence of grooming, and rats show normal exploratory and feeding behaviors (Supplementary Fig. S2). In contrast with previous reports (21), no animals exhibited signs of local inflammation, bloating, or infection, nor did we find any evidence of fibrino-purulent peritonitis during necropsy to recover tumors and other tissues.

Terminal in vivo experiments were performed 5 weeks after achieving clinically relevant chemotherapy doses. They were designed to measure the firing of individual intact sensory neurons in response to physiologically relevant muscle contraction and stretch with electrophysiologic techniques. All in vivo procedures are well established in our laboratory and have been extensively described in previous publications (22–24). Briefly, rats were deeply anesthetized by inhalation of isoflurane (5% in 100% O₂), intubated via a tracheal cannula, and then maintained for the remainder of the experiment (up to 12 hours by 1.5%–2.5% in 100% O₂). Respiratory rate, pCO₂, core temperature, pulse rate and pO₂ were continuously monitored and maintained by adjusting anesthesia and adjusting heat sources. Dorsal roots (lumbar L4–6), muscles, and nerves in the left hindlimb were prepared for stimulation and recording with the rat fixed in a rigid frame at the snout, vertebral bodies, distal tibia, and distal femur (knee angle 120°). Triceps surae muscles (lateral and medial gastrocnemii and soleus) were partially freed of surrounding connective tissue and marked for their resting length (L_o) at ankle angle 90° before their common Achilles tendon was severed at the calcaneus and tied directly to the lever arm of a force- and length-sensing servomotor (model 305B-LR; Aurora Scientific). Triceps surae nerves were loosely positioned in continuity on a unipolar silver stimulating electrode, and all other hindlimb nerves were crushed to reduce extraneous neuronal activity. Dorsal rootlets were carefully freed in continuity from overlying connective tissue and supported on bipolar hook electrodes positioned close to the rootlet's entry into the dorsal spinal cord. Exposed tissues were covered with warm mineral oil in pools formed by attaching the edges of severed skin to the recording frame.

Dorsal rootlets positioned in continuity on bipolar recording electrodes were selected for sampling sensory neurons when they produced robust action potential activity in response to stretch of triceps surae muscles. Individual axons penetrated in these rootlets by glass micropipettes (∼30 MΩ filled with 2 mol/L K+ acetate) were selected for study when electrical stimulation of triceps surae nerves produced orthodromic action potentials that were readily resolvable and had conduction delay of <2 ms. Continuous intracellular recordings from sensory neurons were acquired with Spike2 software (version 8.02). Sensory neurons were classified as Ia by their perfect entrainment to 1-second bouts of high-frequency, small-amplitude vibration (100 Hz, 80 μm), pause in firing during rising twitch force response (muscle shortening), and by responding with an initial burst of high-frequency firing [>100 pulses per second (pps)] at the onset of muscle stretch (24). Two stretch paradigms were used to characterize the firing responses of afferents to physiologically relevant mechanical stimuli. In both paradigms, triceps surae muscles were stretched by 3 mm from L_o (7% strain). Ramp–hold–release stretches tested afferent encoding of both fast dynamic (20 mm/second, 47% strain rate) and static stimuli; successive triplets of triangular stretch tested slow (4 mm/second, 9% strain rate) dynamic and activity-dependent encoding in dynamic stretch known to be influenced by recent signaling history (24). Strains and strain rates fall within values expected for animals engaged in normal activities, for example, locomotion, and have been previously used in our laboratory (24).

At the conclusion of data collection, rats were overdosed with isoflurane inhalation (5%) then immediately transcardially perfused with cold vascular rinse (0.01 mol/L PBS pH 7.4) followed by room temperature fixative (2% paraformaldehyde in 0.1 mol/L PBS, pH 7.4). Muscles and dorsal root ganglia were quickly dissected and postfixed for 1 hour in the same fixative. After a brief wash in a 0.1 mol/L PBS, tissues were incubated in 0.1 mol/L PBS containing 20% sucrose at 4°C overnight. Three to four animals from each of the four groups were selected for whole-transcriptome analysis alone and did not undergo further analysis. Dorsal root ganglia (DRG) from the lumbosacral (L5-S2) spinal cord were surgically removed (4–6 per animal), washed with sterile saline and immediately flash frozen (fresh) in 2-methylbutane (isopentane) and prechilled in liquid nitrogen. DRG were then stored at −80°C for further analysis (see below). After DRG removal, colonic tumors from Apc^Pirc/+ rats were dissected and counted prior to postfixation for IHC analysis.

Details of our IHC procedures have been described previously (25). Briefly, 50-μm–thick sections of skeletal muscles, DRGs, and tumors were cut using a Cryostat (Leica). Tissues sections from all treatment groups were processed simultaneously. All tissue sections were incubated overnight in primary antibodies diluted in blocking buffer (5% normal goat serum, 0.3% Triton100 in PBS). The following primary antibodies used were: rabbit polyclonal anti-NaV1.6 (ASC-009, Alomone Laboratories, 1:200), rabbit polyclonal anti-NaV1.1 (ASC-001, Alomone Laboratories, 1:300), rabbit polyclonal anti-NaV1.7 (ASC-008, Alomone Laboratories,1:100), rabbit polyclonal anti-Kv3.3 (APC-102, Alomone Laboratories, 1:300), mouse monoclonal anti-Kv1.1 (K36/, NeuroMab, University of California Davis, Davis, CA; 1:300), guinea pig polyclonal anti-VGLUT1 (135–104 Synaptic Systems, 1:300), chicken polyclonal anti-neurofilament protein (NF-H, Aves Laboratories, 1:200), mouse monoclonal anti-ASIC2, (E-20, Santa Cruz Biotechnology,1:200), mouse monoclonal anti-EGRF (SC-120, Santa Cruz Biotechnology, 1:100), mouse monoclonal anti-KI67 (SC-23900, Santa Cruz Biotechnology, 1:100), mouse monoclonal antiIL6 (SC-57315, Santa Cruz Biotechnology, 1:100). After washing with PBS, tissue sections were incubated with appropriate fluorescent-conjugated secondary antibodies (Jackson ImmunoResearch Laboratories), diluted in blocking buffer, for 1 hour at room temperature. Following washes with PBS, slides were mounted using Vectashield containing DAPI (Vector Laboratories) to label cell nuclei.

RNA extraction and amplification were performed according to our previously described methods (26). Briefly, RNAs were isolated and purified from sensory neurons of DRG (L5-S2) using the mRNeasy Micro KIT (Qiagen). Total RNA concentration and integrity was assessed on the Bioanalyzer RNA Pico Chip (Agilent Technologies). Labeling of RNAs was performed with the FlashTag Biotin HSR RNA Labeling Kit (Affymetrix, Thermo Fisher Scientific). Gene expression was determined by microarray analysis on the Affymetrix platform (Rat Transcriptome Array 2.0).

In total, 13 [four groups with biologic triplicates per group (4 biologic replicates from Apc^Pirc/++ oxaliplatin group)] global transcriptional expression datasets were generated in this study. Raw mRNA expression data were processed using Affymetrix Expression Console (EC) software Version 1.4. Raw data probes were normalized using SST-RMA algorithm. Differentially expressed genes (DEG) were determined through a linear fixed effects model (limma; ref. 27) with a robust empirical Bayes (eBayes) framework to moderate the residual variances (27). This has the effect of increasing the effective degrees of freedom by which gene-wise variances are estimated, thus reducing the number of false positives for genes with small variances and improving power to detect DEG with larger variances (28).

DEGs were subjected to paired multivariate analysis to discover high-level data structure. By serially linking two analytic techniques, principal components analysis (PCA)-driven Hierarchical clustering, we stabilized clustering performance and reduced dimensionality to tractably few continuous variables containing the most important information on which clustering algorithms are focused, which leads to a better clustering solution (29, 30). PCA and visualization was applied to DEGs with the factoextra (31) and FactoMineR (29) packages. PCs capturing high-level structure were then employed to drive hierarchical clustering. Hierarchical clustering and heatmaps were visualization with the gplots package (heatmap.2; ref. 32) in the R environment (3.5.0). Hierarchical clustering, using the ward.D2 method (33), was applied to the eBayes-filtered gene set to obtain unsupervised visualization of the gene clusters coordinately expressed among the different experimental groups. Prior to hierarchical clustering, genes were standardized across the samples to a mean of zero and an SD of 1. Subclusters were derived using the R function cutree (34).

Hypergeometric testing was conducted on subclusters of named genes using DAVID (Database for Annotation, Visualization, and Integrated Discovery) and g:Profiler to pinpoint significantly enriched pathways or gene ontology terms (35). We increase statistical power by hypergeometric testing on subclusters, which separates coregulated (up or down) probe sets (36). All biological processes with adjusted P < 0.05 were considered significantly enriched. g:Profiler 3–500 (min–max) Benjamini–Hochberg FDR correction intersection terms (2) only annotated genes were searched. Parallel identification of pathways enrichment was performed by gene-set enrichment analysis (GSEA). GSEA uncovers novel sets of functionally related genes by focusing on all detected (above-background) genes. GSEA is unbiased and more sensitive, thus allowing for the detection of even subtle enrichment signals. The Molecular Signature Database v 6.2 (C2-CP, canonical pathways; C3- MIR, miRNA targets; C3- TFT, transcription factor targets; C5-BP, GO biological process; C5- CC, GO cellular component; C5- MF, GO molecular function C7: immunologic signatures gene sets) with the permutation type set to “gene set,” 15–500 (min–max), to calculate statistical significance, as suggested for fewer than seven replicates; default settings were applied to all the other options. For GSEA, a false discovery rate (FDR) of <0.2 was considered statistically significant.

Complementary supervised analysis of DEGs leveraged knowledge of proteins known to be expressed at sensory neuron receptor endings that directly (25, 37, 38) or indirectly (39, 40) regulate excitability. In addition, due to incomplete knowledge of the full suite of proteins that are necessary and sufficient for normal function (40), we expand targeted transcriptional analysis to incorporate gene families, including ion channels, GPCRs, and channel regulators.

Z-axis stacks of images of muscle spindle receptor endings, DRGs, and tumors were constructed by sequentially imaging using a confocal microscope (LSM 700A, Zeiss). Muscle spindles terminals z-stacks (1-μm steps) were captured with a 40× oil immersion objective (N.A 0.6) at 0.6 digital zoom. DRG and tumors z-stacks (1-μm steps) were captured with a 20× objective (N.A 0.4) at 0.6 digital zoom. Stacks of images were processed and analyzed using Amaris (Bitplane) or ImageJ (NIH, Bethesda, MD) imaging software. Figures present images as flat maximal projections of the sum of the z-axis optical slices. Analysis of neuron receptor endings (spindles) and soma (DRG) were performed by an investigator blinded to treatment group. Immunoreactivity in nerve terminals and DRG was determined by measuring the mean pixel intensity of on each en-face portions of the terminal (visualized by the VGLUT1 staining) observed on every 1-μm single image or on the medium and large size cell somas observed on the flat maximal projections of the z-stack images, respectively.

Intra-axonal recordings of action potentials together with records of muscle length and force were digitized (20 kHz) and were monitored online and stored on computer for later analysis with Spike2 and custom-written MATLAB scripts. From raw intracellular data of both paradigms, we extracted 31 measured and derived features (Supplementary Fig. S3) that provide a comprehensive quantification of neuronal signaling characteristics. We then categorized features into four broad clusters that represent functional features encoded by these neurons, containing sensitivity (Thr), dynamic (Dyn), static (Stat), and history-dependent (Hx) signaling information. Supplementary Figure S3 identifies and describes the measured and computed features and functional clustering used for inference.

Linear discriminant analysis (LDA) provided supervised dimensionality reduction for the multiple features of neuron signaling, in an attempt to find a linear combination of features that separated and characterized independent and combinatorial treatment effects. Data were first demeaned, normalized to unit SD, and tested by Bayesian one-way ANOVA (stan_glm; ref. 41). The derived covariance matrix was then normalized by within-group, pooled covariance. The eigenvectors of that modified covariance matrix defined three canonical variables that characterized and separated the four treatment groups identified a priori as Apc^WT+control, Apc^WT+oxaliplatin, Apc^Pirc/++control and Apc^Pirc/++oxaliplatin. LDA and 10-fold cross validation of model performance (repeated holdout method) was performed with the MASS (7.3–51.1) library in the R environment (3.5.0). Individual features characterizing the firing responses of neurons sampled from multiple rats were tested for statistically significant differences with Bayesian one-way ANOVA (stan_glm). Bayesian parameters estimation derives the entire joint posterior distribution of all parameters simultaneously thus circumventing the need to correct for multiple tests on the data.

Bayesian hypothesis testing was performed by analyzing the probability mass of the parameter region in question. This leads to a direct probability measure that defines values inside the 95% of the highest posterior density (HDI) more credible than outside values. HDI was used to make unbiased inferences by directly comparing the posterior probability distributions between two (or more) contrasts of interests, for example, mean comparisons testing. For example, an HDI of a credible difference distribution that does not span zero indicates that the model predictions for the two conditions of interest are different from each other, resulting in evidential support that parameters for both populations are unequal.

All models were developed in fully hierarchal Bayesian framework (Supplementary Fig. S4) with the rstanarm package (2.18.1; ref. 41) in the R environment (3.5.0). Rstanarm implements regression models in stan, which are fit using Hamiltonian Markov Chain Monte Carlo sampling to compute credible parameter values (42, 43). For intercepts and predictors, we use Student t distribution with mean zero and four degrees of freedom as the prior distribution. The scale of the prior distribution is 10 for the intercept and 2.5 for the predictors. Each model was run with four independent chains for 400 warm-up, 4,000 sampling steps, and every second sample (thinning) was discarded. For all parameters, the number of effective (n_eff) samples was >500. Convergence was assessed and assumed to have reached the stationary distribution by ensuring that the Gelman–Rubin shrinkage statistic (rhat, R^ˆ) statistic for all reported parameters was <1.05 (and Monte Carlo SE of the parameter means was <0.002; ref. 44). We report the expected mean parameter values alongside 95% credible intervals using the HDI (Supplementary Fig. S5).

Models performance was validated (Supplementary Fig. S5) by computing out-of-sample predictive accuracy using Pareto-smoothed importance sampling (PSIS) to perform leave-one-out cross (LOO) validation. Although alternative validation strategies (WAIC) are asymptotically equal to LOO, PSIS-LOO provides a more accurate and stable estimate of model performance and is more robust in the finite cases with diffuse priors (weakly informative). We obtained SEs for predictions and for comparison of predictive errors between models to definitively test and rank models based on predictive accuracy (elpd = expected log pointwise predictive density).

We used the ladder rung walking task as a validated outcome to detect and describe sensorimotor deficits (45). The ladder rung apparatus was built as described previously (45). Secure hind foot placement was assessed using the seven-category scale scoring system introduced by Metz and Whishaw (45) with modifications that did not distinguish whether or not an error was associated with a deep fall. The total number of errors in forelimb and hindlimb placement was calculated and the mean error/step ratio was calculated per animal and represented as percent of error where an increase in the number of errors/steps corresponds to a decrease in secure foot placement of the hind and/or fore foot.

Microarray datasets supporting the conclusions of this article are available in the Gene Expression Omnibus (GEO) repository under accession code GSE126773. Data supporting the conclusions of this article are included within the article and its Supplementary Files. Other datasets, code, and models that support the findings of this study are available from the corresponding author upon reasonable request.

We first profiled the transcriptomes of lumbosacral (L4-S1) DRGs, where sensory neurons (Fig. 1B) are vulnerable to chemotherapy (2). We studied wild-type (Apc^WT) rats and rats carrying an Apc gene mutation (Apc^Pirc/+; ref. 15) associated with development of colon cancer (Supplementary Fig. S1; ref. 46), not only in these rats, but also in 80% of patients with colorectal cancer (46). Presence of malignant colonic tumor burden was validated in all Apc^Pirc/+ rats by IHC evaluation of key indicators of malignancy, for example, Ki67 and EGFR (Supplementary Fig. S1) from tumors collected during terminal experiment dissection. Rats in each group were randomized to receive injections of oxaliplatin, a platinum-based compound known to induce neuropathy in patients with cancer and shown in our previous studies to induce chronic movement disorders (23). Oxaliplatin was administered by intraperitoneal injection on a human-scaled dose schedule (Fig. 1A; ref. 23). The oxaliplatin treatment schedule resulted in four experimental groups used for all studies: Apc^WT+control, Apc^WT+oxaliplatin, Apc^Pirc/++control, and Apc^Pirc/++ oxaliplatin (Fig. 1A).

We employed an empirical Bayesian (eBayes) linear model (27) to modify residual variances on a total of 31,042 (22,875 named) genes and identified 3,426 DEGs across all groupwise contrasts (P_eBayesAdj < 0.005). Subsequent post hoc analysis identified DEG from all possible contrasts [fold change (FC) ≥1.2 and ≤−1.2]. Overlap between DEGs retrieved from individual comparisons is illustrated in Fig. 1C. In keeping with the notion that chemotherapy and cancer act through similar molecular mechanisms, we found overlap of 112 DEG affected independently by cancer and chemotherapy (Fig. 1C). Of note, we found that 73.8% (n = 1,704) of all DEGs (n = 2,307) were uniquely affected when cancer and chemotherapy interacted (Apc^Pirc/++oxaliplatin; Fig. 1C). Figure 1D lists the top DEGs (FC ≥ 4 and ≤ −4) identified in Apc^Pirc/++oxaliplatin in comparison with the effects observed in Apc^WT+oxaliplatin or Apc^Pirc/++control. Collectively, these findings demonstrate the first evidence that cancer exacerbates genetic dysregulation induced by chemotherapy, unmasking a more than 4-fold increase in DEGs.

Next, we performed two independent downstream analyses (Supplementary Fig. S6) to test for interaction and to determine which biologic processes might be involved in pathogenesis of sensory disorders. First, we subjected the eBayes-filtered database to unsupervised PC analysis. DEGs were visualized in the new 3D composite space created by PC1–3 (Fig. 1E) that explained 83.98% of the variance and clearly segregated the four experimental groups. Inspection of individual PCs revealed significant evidence of nonlinear interactions between cancer and chemotherapy captured in PC1 (65.6% of explained variance), whereas PC2 (11.14%) and PC3 (7.24%) separated the effects of cancer and chemotherapy, respectively (Fig. 1F). Second, we drove unsupervised hierarchical clustering onto the PC1 solution to identify groups of coregulated genes that express interaction (Fig. 2A). Analysis identified two clusters of genes displaying distinct, but internally consistent, expression profiles clearly representing patterns of interaction (Fig. 2A, C1 and C2; Supplementary Table S1).

We then subjected clusters of named genes (Supplementary Table S1) to pathway enrichment analysis (35) and identified unique subclusters (Fig. 2A) expressing dysregulation in distinct biological processes (Fig. 2B). In keeping with the notion that combinatorial effects of cancer and chemotherapy (oxaliplatin) reflect convergence onto shared signaling pathways, genes mediating inflammatory response, for example, Il6 (P_eBayesAdj = 7.82 × 10⁻⁶), Ptgs2 (P_eBayesAdj = 0.0028), and Cxcr4 (P_eBayesAdj = 3.74 × 10⁻⁷), and reactive oxygen species processes, for example, Hmbox1 (P_eBayesAdj = 0.001), were found to be among the most upregulated in cluster 1a (C1a; Fig. 2B). Notably, we observed substantial induction of genes encoding glycolytic and carbohydrate metabolic processes (FDR ≤ 1.83 × 10⁻³, C1a; Fig. 2B) in the Apc^Pirc/++ oxaliplatin neurons. Coincident with these changes was a statistically significant downregulation of genes mediating lipid metabolism (C2a), suggestive of pervasive metabolic dysfunction beyond the small subset of dysregulated mitochondrial-related processes that we found in response to Apc^WT+ oxaliplatin or Apc^Pirc/++control alone (Fig. 2B; Supplementary Fig. S7A and S7B) and similar to dysfunctional effects previously suggested (2). Among additional unique transcripts downregulated in Apc^Pirc/++ oxaliplatin neurons, Fig. 2B highlights multiple deficits in processes related to extracellular matrix control, mechanical stimulus transduction, (C2a), DNA-binding, RNA polymerization, (C2b), transcription factor activity, hormone signal transduction (C2c; Supplementary Table S1). Furthermore, we observed genes in C1b and 2a that represent specific dysregulation of neuronal and sensorimotor function, for example, ion-channel genes determining neuronal excitability along with coincident downregulation of locomotion. These data provide the first evidence that codependent neuropathy induces broad metabolic reprogramming, fro example, increased glycolytic and decreased lipid metabolism, consistent with cancers capacity to corrupt metabolic processes (47, 48), which extends knowledge beyond the limited impairment observed in mitochondria (Supplementary Fig. S8).

Markers for nerve degeneration, inconsistently observed in patients treated with antineoplastic drugs (49), exhibited mixed results. While we observed few significantly downregulated genes associated with myelination processes (C2b; P ≤ 7.12 × 10⁻³, Fig. 2B), we also found significant induction of processes involved in the preservation of neurons (e.g., enrichment of negative regulation of apoptosis P ≤ 5.96 × 10⁻³) in C1a. Mixed results from these data suggests heterogeneous neuronal populations in the DRG may contain discrete subpopulations with selective vulnerability, whereas others are resistant to degeneration, which is consistent with mixed results seen in human studies (49) and highlights the need for future cell-type specific investigations that are currently unavailable.

Because inflammatory pathways play multiple roles in regulating neuron function (50), we focused attention on their significant response to codependent neuropathy. We selected the most dysregulated gene in the inflammatory pathways, Il6, for molecular validation by assessing gene—protein correspondence. Sensory neurons in dorsal root ganglia were immunolabeled with mAbs targeting IL6 (Fig. 2C). As predicted from the transcriptional profiling (Fig. 2D), the IL6 protein was constitutively expressed in Apc^WT+control rats at lower levels than in Apc^Pirc/+ control animals or Apc^WT+ oxaliplatin animals, and at significantly higher levels in Apc^Pirc/++ oxaliplatin neurons (Fig. 2E), demonstrating that codependent neuropathy was conserved in protein expression (Fig. 2E; Supplementary Fig. S9). Previous studies implicate IL6's capacity to mediate hyperexcitability and morphologic changes associated with neuron death (51). However, studies restricted investigation to acute neurotoxicity, which limits generalizability to symptoms that persist long after treatment. Our data present the first evidence that Il6 (IL6) expression levels remain elevated long after treatment cessation without observable morphologic changes. We show that codependent neuropathy significantly increases Il6 (IL6) when compared with chemotherapy or cancer alone. In light of IL6’s capacity to decrease excitability following chronic exposure (51), our data suggest that neuronal dysfunction (see below) as a result of codependent neuropathy may be expressed as hypoexcitability, in contrast with previous work (5).

We then used GSEA as a second independent tool to investigate the cellular processes underlying codependent neuropathy. We initially focused on the comparison of Apc^WT+ oxaliplatin and Apc^Pirc/++ oxaliplatin (Fig. 3A) and Apc^Pirc/++control and Apc^Pirc/++ oxaliplatin (Fig. 3B) and interrogated entire datasets against the Molecular Signatures Database (Supplementary Table S2). Two novel conclusions can be drawn from the GSEA findings. First, GSEA independently validated codependencies of cancer and chemotherapy in the DEGs identified by transcriptional profiling using paired unsupervised analysis (Fig. 3A and B). GSEA corroborated downregulation of lipid metabolic pathways [normalized enrichment score (NES) = 1.73] and upregulation of glycolytic pathways (NES = −1.94) in Apc^Pirc/++ oxaliplatin sensory neurons when compared with either chemotherapy (Fig. 3A) or cancer alone (Fig. 3B). Focusing analysis on pathways directly involved in neuronal signaling confirmed dysfunction in peripheral neuron ensheathment in the absence of evidence supporting neuron-specific apoptotic response (FDR > 0.83; Fig. 3A). Second, GSEA exposed dysregulation of genes relevant to neuronal excitability particularly potassium ion channels and transporters, and others participating in synaptic communication, for example, SNARE and glutamatergic transmission (Fig. 3A and B). Dysregulation of these and other gene sets (e.g., proinflammatory chemokines) was not predicted from the independent effects of cancer or chemotherapy alone (Supplementary Table S2). Unlike previous studies that consistently report dysregulated voltage-gated sodium channels (52), we find little evidence of differentially expressed sodium channels or disrupted regulatory pathways (52). Instead, we find targeted DEGs related to potassium ion channels unlike those previously identified. We draw two conclusions from these data. First, ion channel clusters are differentially vulnerable to codependent neuropathy. Second, chronic neuropathy may be mechanistically linked to dysregulation of ion channels distinct from those identified in acute preparations.

Taken together, our analyses indicate significant dysregulation of genes and proteins mediating core cellular processes of Apc^Pirc/++ oxaliplatin sensory neurons. Our integrated transcriptomics and protein-level findings confirmed that cancer and chemotherapy affect common processes (13, 53, 54). Moreover, our findings support our central hypothesis that clinically relevant neuropathy depends on systemic interaction between cancer and chemotherapy by showing for the first time that codependence exacerbates neuropathy and revealing it operates through distinct mechanistic pathways, unpredicted from cancer or chemotherapy alone.

Sensory deficits, for example, somatosensory and hearing loss, are commonly reported by patients long after cessation of chemotherapy, but knowledge of the combined effects of cancer and chemotherapy is missing. Mechanosensory systems, for example, cochlear hair cells, and various somatosensory neurons in skin and muscle are responsible for encoding this information and exhibit deficits following chemotherapy (6). We therefore next focused on identifying the combined effects of cancer and chemotherapy on a representative member of mechanosensory neurons, called muscle spindle neurons. Muscle spindles share many of the molecular mechanisms (40) underlying mechanotransduction and encoding found in this broad class of sensory neurons [e.g., hair cells (55), Merkel (56)]. Furthermore, signaling by these neurons encodes sensory features of muscle mechanics necessary for perceiving (proprioception) and coordinating body position and movement, functions, which, when impaired, have the potential to explain persistent disorders in patients (57). We performed targeted analyses of our transcriptome database by querying DEGs identified in Apc^Pirc/++ oxaliplatin rats compared with all other groups against genes encoding proteins known to be constitutively expressed in mechanosensory neurons that mediate unique contributions to neuronal signaling (Supplementary Fig. S10; refs. 25, 37, 40). We found no evidence for independent nor combinatorial effects of cancer or chemotherapy at the genetic level (Fig. 4A) for channels mediating mechanotransduction (e.g., Asic2, Piezo2, ENaCs), signal amplification, and spike encoding [e.g., Scn1a (Nav1.1), Scn8a (Nav1.6), Scn9a (Nav1.7), Cacna1s,-c,-d, -f (Cav1.1-4), Kcnn2 (SK2), Kcna1 (Kv1.1)].

Next, we performed cell-specific IHC analyses to test gene–protein correspondence. We found no evidence to suggest dysregulated protein expression nor did evidence emerge to suggest altered protein distributions (Supplementary Fig. S11–S13) for most proteins, including those mediating mechanotransduction (ASIC2), signal amplification (Nav1.1, Nav1.6, Nav1.7), spike encoding (Kv1.1), and mechanosensory sensitivity (VGLUT1). These findings demonstrate that many of the molecular mechanisms responsible for mechanotransduction and determining excitability of mechanosensory neurons were preserved in the Apc^Pirc/++ oxaliplatin rats. These findings narrow the field of candidate mechanisms for which treatment might reasonably restore normal function.

Having exhausted examination of proteins known to be expressed in mechanosensory neurons, we extended our search to include ion channels yet to be identified, but known to regulate neuronal signaling in other classes of neurons (58–60). As one of the most dysregulated voltage-gated ion channels, Kcnc3 (encoding protein Kv3.3) emerged as a novel candidate (Fig. 4A and B; Supplementary Fig. S14). Pharmacologic and genetic perturbation of Kv3.3 (Kcnc3) impairs neuronal signaling (60) and is causally linked to ataxias (61) that are consistent with functional deficits observed in patients with cancer long after cessation of chemotherapy. We then immunolabeled sensory cell bodies in DRG and their receptor endings taken from Apc^WT+control with mAbs targeting Kv3.3, which revealed the first evidence of this ion channel in all large diameter cell bodies in the DRG that supply mechanosensors (Fig. 4C; Supplementary Fig. S15) and endings (Fig. 4E). Results of comparison across experimental groups indicate that Kv3.3 expression was significantly lower in Apc^Pirc/++ oxaliplatin compared with Apc^WT+control (Fig. 4D), Apc^Pirc/++control, or Apc^WT+ oxaliplatin (Supplementary Fig. S15), demonstrating that codependent dysfunction is conserved in cell-type–specific protein expression. Moreover, conserved Kv3.3 expression in neuromuscular junctions suggests that Kv3.3 dysfunction was constrained to sensory neurons (Supplementary Fig. S16). This is the first evidence implicating Kv3.3 in the development of neuropathy and the first data, to our knowledge, demonstrating a channelopathy persists long after treatment cessation. The capacity of Kv3.3 to drive fast spiking and to enhance transmitter release–processes that are required for mechanotransduction and spike encoding–suggests that this ion channel might provide both a potential novel mechanism for the dysfunction of mechanosensory neurons and a potential target for therapy.

Next, we examined whether the codependent neuropathy we discovered at gene and protein levels had negative consequences on mechanosensory function in living animals (Supplementary Fig. S17). Applying electrophysiologic methods to rats in vivo, we recorded spiking activity from single mechanosensory neurons responding to naturalistic mechanical stimuli (Fig. 5A). From these spiking responses, we collected 31 measured and derived parameters (average of four trials; n = 11 Apc^WT+control, n = 19 Apc^WT+OX, n = 20 Apc^Pirc/++control, n = 10 Apc^Pirc/++ oxaliplatin; Supplementary Fig. S3) from which we extract four functional features encoded by these neurons—sensitivity, dynamic, static, and history-dependent signaling information (Supplementary Fig. S3). In Apc^WT+control rats, we found that muscle stretch elicited the spiking expected from sensory neurons in normal animals (24) and humans (62), specifically high frequency initial bursting at stimulus onset, increasing dynamic firing during increasing stimulus, sustained static firing with slow accommodation during the hold phase of stretch (Fig. 5B), and history-dependent reduction in spike number (Supplementary Fig. S18).

We next tested whether Apc^Pirc/+ or oxaliplatin treatment alone induce signaling dysfunction. Although we find transcriptional dysregulation in Apc^Pirc/++control (Fig. 1C; Supplementary Fig. S7A), we found a remarkable degree of concordance with signaling in Apc^WT+control (Fig. 5B). Apc^Pirc/++control neuronal signaling displayed comparable sensitivity to stimuli and showed similar static behavior across all neurons (n = 20; Fig. 5B and C) as the Apc^WT+control. We then found that oxaliplatin treatment of Apc^WT rats induced mild signaling deficits that were less pronounced, occurred in a small proportion of neurons (5/19), and were primarily restricted to sensitivity functional clusters (Fig. 5B and C; Supplementary Fig. S19D), validating our previous findings of restricted deficits observed in a different strain of healthy rats treated with oxaliplatin alone (23).

In contrast to the observations in Apc^WT+control, Apc^WT+ oxaliplatin, or Apc^Pirc/++control rats, we found drastically impaired neuronal signaling in Apc^Pirc/++ oxaliplatin rats (Fig. 5B). High frequency initial bursting was attenuated 4-fold and sensitivity to stimuli was reduced by 3-fold at both high (Supplementary Fig. S19D) and low velocity stretch (Supplementary Fig. S18 and S19D). During dynamic stimuli, we found marked reduction in the number of spikes over both stretch conditions (Fig. 5C). Neurons in Apc^Pirc/++ oxaliplatin rats also failed to sustain firing (Fig. 5B). In addition, we observed fewer spikes (Fig. 5C) during rapid accommodation immediately preceding signal deletion, despite constant stimulus (Fig. 5B). We then simulated naturalistic compensation by subjecting Apc^Pirc/++ oxaliplatin neurons to repeated trials at 1, 2, and 3× background stimulus (L_o strain) intensity, which failed to completely rescue signaling back to Apc^WT+control levels (Supplementary Fig. S20B–S20D). These data are the first to directly demonstrate a damaging functional interaction between the systemic effects of cancer and chemotherapy in sensory neurons in a living animal. While validating our previous findings of signaling deficits in chemotherapy alone (23), these data indicate that restrictive signaling deficits induced by chemotherapy or cancer alone were insufficient to reproduce the magnitude, direction, and number of dysfunctional parameters induced by codependent neuropathy.

We next tested whether the impaired neuronal signaling of Apc^Pirc/++ oxaliplatin neurons might be explained by dying-back degeneration of sensory nerve terminals. While the underlying mechanisms for chemotherapy induced neuropathy are not understood, current opinion identifies “dying-back” axon degeneration as a major pathology in this disorder (19). Blinded reviewers assessed the structural integrity of sensory afferents and receptor endings immunolabeled against neurofilament protein (NEFH; ref. 25). We then tested for functional evidence of degeneration by comparing differences in conduction delays, a physical measure of nerve degeneration that is used clinically (49). Neither histologic observation of sensory afferents and receptors (Fig. 5D) nor axon conduction tests (Supplementary Fig. S21) revealed evidence of dying-back degeneration of sensory nerve terminals (Apc^Pirc/++ oxaliplatin: 1.41 ± 0.13 ms, n = 10. Apc^WT+control: 1.52 ± 0.14 ms, n = 11; Supplementary Fig. S21). Our findings demonstrate that transcriptional changes we observed for some markers of nerve degeneration (Fig. 3) were insufficient to yield nerve degeneration. While our results indicate selective resistance of mechanosensory neurons in Apc^Pirc/++ oxaliplatin rats, they are in line with the inconsistency of physical evidence for nerve degeneration reported in clinical studies (49). Furthermore, the ability of muscle vibration to entrain firing of these neurons in all treatment groups as they do in Apc^WT+control rats further refutes the existence of degeneration's necessity in signaling disorders (Supplementary Fig. S21). Overall, our findings suggest that the structure and core ability of mechanosensory neurons to produce action potentials remain unimpaired by codependent neuropathy and does not depend on dying-back degeneration.

Data presented in Fig. 5 suggests damaging cancer–chemotherapy interactions are conserved in vivo and have consequences in sensory function. To assess the extent of interaction quantitatively, we took an unbiased statistical approach by subjecting all neurons (n = 60, 240 total trials) from all experimental groups to a machine learning algorithm [linear discriminant (LD) analysis]. Our approach reduced complex feature space into canonical variables giving us a high-level understanding of where interaction emerges, without biased feature selection a priori. Our analysis yielded three canonical variables (Fig. 6A) that achieved overall 94.7% classification accuracy (Supplementary Fig. S19A). We then visualized neuronal signaling in the new 3D composite space created by LD1-3 (Fig. 6A). By projecting high-dimensional parameter and feature data onto a simplified 3D canonical space, statistically significant nonlinear interaction between cancer and chemotherapy clearly emerged in the first dimension (LD1; 54.23% proportion of variance; Fig. 6A; Supplementary Fig. S19B). Notably, dysfunction in LD1 induced by the cancer—chemotherapy interaction occurred in the opposite direction to that predicted by independent effects, and its magnitude was amplified greater than their sum (Fig. 6A; Supplementary Fig. S19B and S19C). To test the statistical significance of the cancer—chemotherapy interaction, we conducted Bayesian model comparison with full factorial and all restricted models using leave-one-out cross validation. We then quantified and validated each model's predictive performance by computing the expected log predictive densities (ELPD; measure of a model's out-of-sample predictive accuracy in Fig. 6B–E). We found decisive evidence in favor of a model including a cancer–chemotherapy interaction predictor (ELPD diff ≥ 48 SE <8.1; Fig. 6E). Moreover, our data and generative modeling conclude that codependent interaction (Fig. 6B–D) is necessary to accurately and reliably reproduce clinically relevant neuronal signaling deficits observed in Apc^Pirc/++ oxaliplatin rats (Fig. 6E)

Having determined the significance of an interactive effect, we then asked which functional clusters are highly enriched in LD1. By examining the discriminant function coefficients, we found that static signaling and sensitivity functional clusters represent a large portion of the explained variance in LD1. These findings suggest the sensitivity and static signaling clusters are most susceptible to codependent neuropathy (Supplementary Fig. S19D).

While dynamic signaling provided relatively modest contributions to LD1 (lower coefficient weights in Supplementary Fig. S19D), examining the parameter-level data allowed us to discover the expression of distinct classes of interactions. Independently, oxaliplatin and Apc^Pirc/+ induce highly conserved (87.5%) opposing effects exclusively within the dynamic functional cluster. We consistently found that Apc^Pirc/+ mutation increased and oxaliplatin treatment decreased dynamic firing (Fig. 6G). This led us to predict that their combination would nullify dysfunction and approximate Apc^WT+control signaling. Instead, we observed opposing interactions that resulted in drastic reduction in dynamic firing properties. In contrast, for history-dependent functional clusters, we found that cancer–chemotherapy interaction emerged exclusively in Apc^Pirc/++ oxaliplatin neurons, since neither Apc^WT+ oxaliplatin nor Apc^Pirc/++control neurons were disrupted (Fig. 6H).

Taken together, we found that 58% of the neuronal signaling parameters in Apc^Pirc/++ oxaliplatin rats showed high-confidence interaction effects, in that those deficits were greater than those observed in either Apc^WT+ oxaliplatin nor Apc^Pirc/++control rats. Our findings demonstrate that complex systemic perturbations of chemotherapy and cancer result in synergistic interactions for sensitivity signaling characteristics (Fig. 6F), opposing interactions for peak firing rate during dynamic signaling (Fig. 6G), and emergent interactions for history-dependent and static signaling characteristics (Fig. 6H). In summary, Apc^Pirc/++ oxaliplatin rats express widespread neuronal signaling dysfunction that targets all functional clusters, was independent of structural or physiologically detected degeneration and were not fully accounted for by a global decrease in sensitivity.

Apc^WT+control and Apc^Pirc/++control rats achieved a high incidence of rear and fore foot placement, when walking on the rungs of an uneven horizontal ladder, with a foot-drop error rate of only 2.4% (±2.7% SD, 6 animals; Fig. 7A; Supplementary Movie S1) and 2.5%, respectively (±2.9% SD, 4 animals; Fig. 7A; Supplementary Movie S2). In contrast, Apc^Pirc/++ oxaliplatin exhibited a 19.2% foot-drop error rate [±5.6% SD, 7 animals; hierarchical Bayesian model (stan_glm) vs. Apc^WT+control and Apc^Pirc/++control; Fig. 7A; Supplementary Movie S3]. Apc^Pirc/++ oxaliplatin error rates significantly exceeded our previous findings following oxaliplatin treatment alone [Wistar+ oxaliplatin) of 8.4% (±3.1% SD, 7 animals; hierarchical Bayesian model (stan_glm) vs. Apc^Pirc/++ oxaliplatin; Fig. 7]. Independent analysis of fore- (Fig. 7B and D) and hind-foot (Fig. 7C and E) error rates revealed that skilled motor behavior degradation was significantly exacerbated by the emergence of fore-limb errors (10.4% ± 3.9 SD) not present in Apc^WT+control (0%), Apc^Pirc/++control (2.3%±2.8SD), or Wistar+ oxaliplatin (0%; hierarchical Bayesian model (stan_glm)).

Here we present original evidence that cancer transforms the nature and magnitude of neuropathy induced by chemotherapy alone. Our preclinical study of rats is the first to compare the neuropathic effects of chemotherapy and cancer, both independently and in combination, impractical in human study. These comparisons for global transcriptional analysis of sensory neurons in DRG revealed dysregulation of genes uniquely induced, amplified, or suppressed by the combination of cancer and chemotherapy. Codependence was conserved as impaired spike encoding of mechanosensory stimuli and novel ion channelopathy. While the current report is limited to rats, extensive conservation of core molecular and cellular processes across mammalian species (63) leads reasonable expectation that cancer–chemotherapy codependency, whether different in detail, extends to humans. All considered, we conclude that inattention to codependencies necessarily prevents the development of mechanism-based treatments for sensory neuropathy, which remains unexplained and unabated in patients receiving chemotherapy for cancer.

Our transcription analyses expose two potential targets for treating sensorimotor disorders among the debilitating patient symptoms that persist following treatment of various cancers with platinum-based compounds and other antineoplastic agents, for example, taxanes. Patients display deficits in the spatiotemporal parameters (speed and stride) of walking gait, in control of posture, and in balance relying on proprioception (64). Acknowledging that these disabilities may arise from a variety of lesions in the nervous system, the deficits we observe here at the receptor origin of detection and encoding of muscle mechanics would necessarily impair movements and postures and cannot be fully compensated by other senses, for example, vision. Furthermore, signaling by mechanosensory neurons in rats closely resembles that in human (62). For these reasons, we assign special attention to depressed signaling by muscle spindles and to the associated decrease in expression of the voltage-gated ion channel Kv3.3 and its gene Kcnc3. Reports that Kcnc3 gene knockout impairs firing responses of neurons and results in ataxia (61) promote this gene and its ion channel as potential contributors to signaling deficits observed in this study. At this time, we are uncertain about the sufficiency of Kv3.3 channelopathy to explain deficient mechanosensory signaling in cancer treated by chemotherapy. Uncertainty arises, in part, because comprehensive understanding of the molecular mechanisms underlying the function of these sensory neurons is lacking, although advanced by our discovery of Kv3.3 in mechanoreceptors of healthy animals. Another candidate target for treating sensorimotor disorders is the inflammatory signaling molecule IL6, shown here to express large increases in both gene and protein expression. As a result of its effect in suppressing neuronal firing behavior (51), IL6 has the potential to explain decreased signaling by mechanosensory neurons. We identify both Kv3.3 and IL6, therefore, as targets worthy of testing for their potential value in treating movement disorders. Additional points of interaction likely exist, which remain either undiscussed (e.g., FosB) or potentially masked due to transcriptional profiling of DRG that may result in variable differential expression across heterogenous cell types (Fig. 1D).

We observed abnormal skilled sensorimotor behavior in Apc^Pirc/++ oxaliplatin rats. Deficits were significantly greater than those previously observed in our lab from animals treated with oxaliplatin alone (23). The emergence of novel fore-limb deficits typically not observed in this challenge to hind-limb capacities provides further evidence that exacerbation effects are unpredicted from studying chemotherapy alone. Furthermore, our data indicates that intact proprioceptive sensory input from muscle spindles is necessary to ensure limb placement accuracy in a skilled behavioral task and may be sufficient to explain sensorimotor deficits experienced by cancer survivors.

Our study, being restricted to a single time point following chemotherapy, does not assess the stability of potential treatment targets that emerge from codependent neuropathy. It is reasonable to expect, however, that the codependence undergoes dynamic change. Cancer, presumably also its associated systemic effects, undergoes complex progression as subpopulations of cancer cells differentially resist, adapt, or succumb to chemotherapy. Similarly, surgical resection of tumors or surgical associated analgesics might impose additional points of interaction that have the potential to impact DRG neurobiology. Moreover, biological systems themselves initiate dynamic responses to perturbations. In the nervous system, homeostatic regulation initiates compensatory mechanisms to offset perturbations in neuronal excitability (59). It is thought provoking in this regard to consider the possibility that the chronic hypoexcitability, which follows the acute hyperexcitability (2, 19), that develops during chemotherapy reflects the actions of a dysregulated compensatory mechanism. Understanding these nonlinear interactions and their effects on potential treatment are both challenging and necessary for developing effective cancer treatment.

No potential conflicts of interest were disclosed.

Conception and design: S.N. Housley, P. Nardelli, T.C. Cope

Development of methodology: S.N. Housley, P. Nardelli, J.F. McDonald, T.C. Cope

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.N. Housley, P. Nardelli, D.I. Carrasco, T.M. Rotterman, E. Pfahl, L.V. Matyunina, J.F. McDonald, T.C. Cope

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.N. Housley, P. Nardelli, J.F. McDonald, T.C. Cope

Writing, review, and/or revision of the manuscript: S.N. Housley, P. Nardelli, J.F. McDonald, T.C. Cope

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S.N. Housley, D.I. Carrasco, T.M. Rotterman, T.C. Cope

Study supervision: S.N. Housley, T.C. Cope

This work is supported by NIH grant R01CA221363 and R01HD090642. We thank the staff of the Physiological Research Laboratory at Georgia Institute of Technology. We are grateful to Dr. Thomas Burkholder for technical assistance, discussions, and comments. We thank Dr. Mengnan Zhang and Dr. Evan Clayton for discussions and comments.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.