Integrative multiomic analyses of dorsal root ganglia in diabetic neuropathic pain using proteomics, phospho-proteomics, and metabolomics.

Diabetic peripheral neuropathy (DPN) is characterized by spontaneous pain in the extremities. Incidence of DPN continues to rise with the global diabetes epidemic. However, there remains a lack of safe, effective analgesics to control this chronic painful condition. Dorsal root ganglia (DRG) contain soma of sensory neurons and modulate sensory signal transduction into the central nervous system. In this study, we aimed to gain a deeper understanding of changes in molecular pathways in the DRG of DPN patients with chronic pain. We recently reported transcriptomic changes in the DRG with DPN. Here, we expand upon those results with integrated metabolomic, proteomic, and phospho-proteomic analyses to compare the molecular profiles of DRG from DPN donors and DRG from control donors without diabetes or chronic pain. Our analyses identified decreases of select amino acids and phospholipid metabolites in the DRG from DPN donors, which are important for cellular maintenance. Additionally, our analyses revealed changes suggestive of extracellular matrix (ECM) remodeling and altered mRNA processing. These results reveal new insights into changes in the molecular profiles associated with DPN.

Introduction

Chronic pain is a growing health burden, impacting an estimated 20% of adults in the United States[1]. Among chronic pain cases, an estimated 25% can be attributed to neuropathic pain[2]. A common painful neuropathic condition is Diabetic Peripheral Neuropathy (DPN), which is caused by Diabetes Mellitus (DM). Type II Diabetes Mellitus (T2D) is an epidemic that affects 462 million people globally, a number expected to grow markedly, with 50% of patients developing DPN[3-5]. DPN is defined as a length-dependent sensorimotor polyneuropathy. DPN is associated with hyperglycemic pathology as well as hyperlipidemia; however, a mechanism of DPN development remains unclear[5-7]. DPN typically manifests in a ‘stocking and glove’ pattern in which patients’ lower limbs are most affected, consistent with disease of long sensory neurons. Additionally, DPN is associated with axonal demyelination and degeneration leading to nerve dysfunction and possible cell death[8]. Exponential growth in the diabetes epidemic and a lack of an effective curative treatment for DPN make it increasingly imperative to gain insights into the mechanisms of disease progression[2,9,10].

Dorsal root ganglia (DRG) contain the cell bodies of peripheral sensory neurons that relay pain from the periphery to the central nervous system[11]. It is well established that transcriptomic changes occur in DRG in response to injury and in association with neuropathic pain development[11]. Our previously published mRNA sequencing (RNA-seq) study suggested an association between DPN, increased expression of immune-related genes, and decreased expression of neuronal genes in the human DRG[12]. However, due to translational and post-translational regulation, expression of a specific mRNA is not a perfect indicator of the abundance of the corresponding protein, especially in the context of neurological dysfunction[13-18]. Despite this, to the best of our knowledge, no study to date has explored translational and post-translational regulation in human DRG. Protein abundance and post-translational modifications (PTM) are central to cellular function, but easily deregulated in response to physiological changes[19]. Phosphorylation is an important PTM upon which many cellular processes are reliant. Additionally, metabolites can be influenced directly by endogenous and exogenous factors and are direct players in biochemical processes, thus offering a strong descriptor of molecular phenotype. A few studies have identified metabolomic and proteomic trends suggestive of disrupted energy metabolism and mitochondrial function in DRG of diabetic rodents[20-23]. Better understanding of molecular regulation with painful conditions in human systems is vital for development of novel therapeutics[24].

Therefore, we aimed to integrate previous transcriptomic findings[12] with metabolomic, proteomic, and phospho-proteomic data to offer a comprehensive understanding of the biochemical changes associated with DPN in human DRG. Presently, we profiled small molecular weight metabolites and lipids, as well as proteins and phospho-peptides in human DRG samples from DPN donors and non-diabetic controls. We identified decreased amino acids and phospholipid metabolites, as well as alterations in abundance of extracellular matrix (ECM) proteins and phosphorylation of RNA binding proteins. From this, we suggest a role of disrupted amino acid and protein metabolism in the development of DPN.

Results

Tissue donor cohort characteristics

For metabolomic analysis, DRG were obtained from 7 DPN donors and 10 non-diabetic control donors. For proteomic and phospho-proteomic analyses, DRG were obtained from 5 DPN donors and 5 non-diabetic control donors. Within both cohorts of donors, no significant differences were observed between donor groups in age, sex, or BMI (Table 1). A significant difference was observed in donors on ACE inhibitors in the metabolomic cohort, but no ACE inhibitor metabolites were detected in our analysis[25]. A non-significant trend was observed in differences in age in our metabolomic cohort, this was addressed in our statistical analysis of metabolite data.

Table 1

Donor characteristics for DRG tissues used in metabolomic and proteomic analyses.

	Metabolomic cohort			Proteomic cohort
	Control, N = 10	DPN, N = 7	p value	Control, N = 4	DPN, N = 5	p value
DRG storage			0.9a
Snap frozen	5 (50%)	4 (57%)		4 (100%)	4 (100%)
RNAlater	5 (50%)	3 (43%)
Sex			0.9a			0.5a
Female	5 (50%)	3 (43%)		3 (75%)	2 (40%)
Male	5 (50%)	4 (57%)		1 (25%)	3 (60%)
Age	46.8 (7.6)	53.4 (9.8)	0.1b	51.8 (2.5)	53.4 (9.8)	0.7b
BMI	27.4 (6.1)	30.6 (8.6)	0.7b	35.0 (11.4)	30.6 (8.6)	0.5b
Race/ethnicity			0.5a			0.5a
African American	2 (20%)	3 (43%)		1 (25%)	2 (40%)
Hispanic/Latino	2 (20%)	3 (43%)		0 (0%)	1 (20%)
White	6 (60%)	2 (29%)		3 (75%)	2 (40%)
Cause of death			0.6a			0.6a
Anoxia/cardiovascular	2 (20%)	3 (43%)		1 (25%)	2 (40%)
CVA/ICH/stroke	4 (40%)	4 (43%)		1 (25%)	2 (40%)
MVA/head trauma/blunt injury	4 (40%)	1 (14%)		2 (50%)	2 (40%)
Medications
Vasodilators	3 (30%)	1 (14%)	0.6a	3 (60%)	2 (40%)	0.9a
Opioids	1 (10%)	3 (43%)	0.2a	0 (0%)	2 (40%)	0.4a
Gabapentin	0 (0%)	2 (29%)	0.2a	0 (0%)	2 (40%)	0.4a
Insulin	0 (0%)	2 (29%)	0.2a	0 (0%)	3 (60%)	0.2a
Metformin	0 (0%)	2 (29%)	0.2a	0 (0%)	2 (40%)	0.2a
Ca2+ channel blockers	1 (10%)	4 (57%)	0.1a	0 (0%)	2 (40%)	0.4a
β Blockers	1 (10%)	4 (57%)	0.1a	0 (0%)	3 (60%)	0.2a
ACE Inhibitors	1 (10%)	5 (71%)	0.035a,*	0 (0%)	3 (60%)	0.2a
Statins	1 (10%)	4 (57%)	0.1a	1 (20%)	2 (40%)	0.9a
Comorbidities
Cardiac disease	1 (10%)	4 (57%)	0.1a	1 (20%)	2 (40%)	0.9a
Hypertension	4 (40%)	6 (86%)	0.13a	3 (60%)	4 (80%)	0.9a
Neurological disease	2 (20%)	1 (14%)	0.058a	1 (20%)	3 (60%)	0.5a
Cancer	0 (0%)	1 (14%)	0.4a	0 (0%)	0 (0%)
Diabetes	0 (0%)	7 (100%)	5e−05a,*	0 (0%)	5 (100%)	0.008a,*

Results given as N (%) or mean (standard deviation).

aFisher’s Exact Test.

bT-test.

* indicates variables which are significantly different between control and DPN donor groups (p-value < 0.05).

Post-mortem tissue donation networks offer the best opportunity for analyses of human DRGs. L4, L5, and S1 DRGs were obtained as these ganglia contain the primary afferent neurons that innervate the distal extremities of the foot where patients with DPN often experience pain (Fig. 1a). Donor medical history confirms neuropathy for 10 years or more in the DPN cohort as observed from the available timeline. However, with current limitations in obtaining human DRG[26], we were unable to obtain detailed clinical measures of pain hypersensitivity in DRG donors. We have earlier reported pathological evaluation of DRG tissues from the DPN and control donors showing variation in pathology in both donor groups[12]. As reported earlier[12], within the non-diabetic control group, pathology ranges from apparently normal tissue to moderate ganglionic cell loss, whereas within the DPN group, pathology ranges from apparently normal tissue to moderate to severe ganglionic cell loss.

Figure 1

Study design and evaluation of DRG quality (a) Schematic representation of study design. Innervation pattern for L4, L5, and S1 DRG used in analyses overlap with painful region in DPN. Created with BioRender. (b) Heat map showing scaled and centered relative abundance data for select metabolites. Abnormally high profiles of all metabolites here would be indicative of tissue degradation. Two samples show higher abundance of some metabolites. However, these levels are reasonable for normal biological variation, and profiles of all noted metabolites are not observed to be consistently high. As such, DRG samples were determined to be of quality allowable for analysis. Row-wise annotations indicate degradative process which might produce indicated metabolites. CON control samples, DPN diabetic peripheral neuropathy.

For our present study, profiles of protein, DNA, and RNA degradation products as well as indicators of cellular damage were evaluated to characterize tissue quality (Fig. 1b). Abundance of these degradation products would be abnormally high if tissues were poorly handled or stored. We do note that one sample has exceptionally high levels of hypoxanthine, inosine, and certain amino acids (Fig. 1b). A second sample is noted to have high levels of certain amino acids as well (Fig. 1b). However, for both samples, these increases were not consistent across all amino acids. Additionally, diphosphates ADP and GDP were particularly low in these two samples. As such, we determined no DRG samples with abnormal degradation were included in our analysis.

Metabolomic profiling suggests disrupted amino acid and phospholipid metabolism

Metabolomic profiling was performed with 7 DRG from donors with DPN and 10 DRG from non-diabetic controls (Table 1). Capillary electrophoresis mass spectrometry (CE-MS) and liquid chromatography mass spectrometry (LC–MS) collectively identified 327 metabolites in donor DRG tissues. 165 metabolites were detected in more than 80% of samples and included in further analysis (Fig. 2a).

Figure 2

Decreased amino acids and phospholipid metabolites with DPN (a) Volcano plot showing metabolite abundance according to ratio of DPN/control and p-values before covariate control. Significantly regulated metabolites are considered to be those with p-value < 0.05. Metabolites with trends in regulation were identified with p-value < 0.1. Dashed lines represent p-values 0.1 and 0.05 before covariate control. (b) Heat map showing scaled and centered relative abundance of differentially regulated metabolites (p-value after covariate control < 0.05). Asterisks (*) indicates metabolites which were validated with quantitative panels. Row-wise annotations indicate metabolite class, metabolite classes represented by only one metabolite are grouped together as other. CON indicates control samples, DPN indicates diabetic peripheral neuropathy. (c) Metabolite enrichment analysis with SMPDB library shows enrichment in phospholipid biosynthesis and amino acid metabolism. Annotations shown are significant with p-value < 0.05. Enrichment ratio is observed number of metabolites divided by expected number metabolites in metabolite set. (d) Schematic representation of PC and PE biosynthetic? pathways as well as derivation of N-acylethanolamines, created with BioRender. (e) Box plot showing trends in decreased PE derived lipids, SEA and OEA. Indicated p-value was calculated after covariate control.

Dimensionality reduction with principal component analysis (PCA) suggested the influence of donor pain condition on metabolite profiles to be confounded by variables such as sample storage condition (Supplementary Fig. S1). Linear regression models were applied to candidate metabolites to account for variation due to sample storage condition, donor age, and donor sex. 21 metabolites were identified as significantly altered as a result of neuropathic pain (p-value adjusted for covariates < 0.05) after this covariate adjustment (Fig. 2b; Supplementary Table S1). All differentially regulated metabolites were decreased in the DRG from donors with DPN.

Of the 21 identified metabolites, 8 are classified as carboxylic acids and derivatives (N-acetylaspartic acid (NAA), creatine (Cr), glycine (Gly), aspartic acid (Asp), glutamine (Gln), glutamate (Glu), guanidoacetic acid, and cysteine glutathione disulfide) (Fig. 2b). We noted additional non-significant decreases in carboxylic acids and derivatives including asparagine (ratio (DPN/control) = 0.46, p-value adjusted for covariates = 0.057) and ornithine (ratio (DPN/control) = 0.38, p-value adjusted for covariates = 0.058) (Supplementary Table S1). Similarly, our enrichment analysis revealed significantly enriched metabolite sets related to arginine and proline metabolism, aspartate metabolism, glycine and serine metabolism, and glutamate metabolism (Fig. 2c). Because these observations suggest alterations in amino acids, the concentrations of proteinogenic amino acids were then examined with quantitative panels (Supplementary Fig. S2). Three differentially regulated amino acids (Gln, Gly, and Glu) were observed to be among the highest abundance proteinogenic amino acids in the DRG (Supplementary Fig. S2). Of note, non-proteogenic amino acid NAA, was the most significantly decreased metabolite (ratio (DPN/control) = 0.37, p-value adjusted for covariates = 0.00995). Given that NAA is highly abundant in the nervous system and neuronally important[27], regulation of NAA was confirmed with a quantitative panel (Supplementary Table S1).

Among our significantly enriched metabolite sets, we also observed an enrichment for biosynthetic pathways of two phospholipids, phosphatidylcholine (PC) and phosphatidylethanolamine (PE) (Fig. 2c). Organonitrogen compounds choline (Cho, ratio (DPN/control) = 0.71, p-value adjusted for covariates = 0.049, Fig. 2b) and ethanolamine (Etn, ratio (DPN/control) = 0.48, p-value adjusted for covariates = 0.0068, Fig. 2b) are involved in both pathways (Fig. 2d). Ethanolamine phosphate (Etn-P, ratio (DPN/control) = 0.56, p-value adjusted for covariates = 0.027, Fig. 2b) is also involved in both pathways (Fig. 2d), whereas S-adenoylmethionine (SAM) (ratio (DPN/control) = 0.52, p-value adjusted for covariates = 0.0164, Fig. 2b) is involved in PC synthesis (Fig. 2d). To further support altered phospholipid dynamics, we observed trending losses in two N-acylethanolamines, oleoyl-ethanolamide (OEA, ratio (DPN/control) = 0.69, p-value adjusted for covariates = 0.082, Fig. 2e) and stearoyl-ethanolamide (SEA, ratio (DPN/control) = 0.66, p-value adjusted for covariates = 0.052, Fig. 2e). Given that PE is the sole endogenous source of OEA and SEA, these metabolites were selected for validation. Quantitative panels confirmed the identity and trend in differential abundance of OEA and SEA (Fig. 2e; Supplementary Table S1), as well as differential regulation of Cho and SAM (Supplementary Table S1).

Proteomic profile shows extracellular matrix remodeling

Proteomic profiling was performed with 5 DRG derived from donors with DPN and 5 DRG derived from control donors without diabetes or chronic pain (Table 1). After sample prep, total peptide abundance in 1 non-diabetic control sample was abnormally low, this sample was excluded from further analysis. 6186 proteins were identified in human DRG samples. 247 proteins were identified as significantly differentially regulated (Benjamini–Hochberg adjusted p-value < 0.05, Fig. 3a). Dimensionality reduction analyses with PCA and t-distributed Stochastic Neighbor Embedding (t-SNE) did not show clear sample grouping by neuropathic pain (Supplementary Fig. S3). We believe this was due to some influence of covariate variables and as such used linear modeling to adjust for age and sex. 41 proteins were identified as significantly altered after covariate control (p-value adjusted for covariates < 0.05, Supplementary Table S2), of these 5 are regulated by fold change greater than 1.5 (|log2(DPN/control)|> 0.58, Fig. 3b).

Figure 3

Extracellular matrix remodeling with DPN (a) Volcano plot showing protein abundance according to ratio of DPN/control and adjusted p-values calculated with ProteomeDiscoverer. Vertical lines at 0.58 and − 0.58 represent an abundance ratio of 1.5 and 0.67 respectively. Proteins are considered up-regulated or down-regulated with adjusted p-value < 0.05 and fold change > 1.5. Proteins are considered trending up or down with 0.05 < p-value < 0.1 or p-value < 0.05 and fold change < 1.5. (b) Heat map showing relative abundance of differentially regulated proteins, identified after covariate control (p-value after covariate control < 0.05, abundance ratio > 1.5 or < 0.67). CON indicates control samples, DPN indicates diabetic peripheral neuropathy. (c) Results of enrichment analysis with differentially regulated proteins (p-value after covariate control < 0.05). Top gene ontology annotations include biological adhesion and extracellular matrix. (d) Western blot shows full length TNR (150 kDa) and heavy TNR products (> 150 kDa) only for DPN samples. Displayed image is from one single membrane, with dividing line added to separate sample groups. (e) IHC results showing TNR (red) localization surrounding neurons (NeuN, green) in DRG sections. Identified neurons surrounded by TNR are indicated with a white arrow. Red staining inside the neuronal cell body is likely lipofuscin, indicated with grey arrowheads. DAPI is labeled in blue. A z-test compared proportion of neurons surrounded by TNR in control and DPN derived DRG sections (p-value = 7.6e−11).

From statistical overrepresentation test using Panther, differentially regulated proteins (p-value adjusted for covariates < 0.05) showed enrichment in gene ontology annotations including ECM and biological adhesion (Fig. 3c). Changes to the ECM and cellular adhesion is highlighted by upregulation of the neural ECM protein, tenascin R (TNR, ratio (DPN/control) = 1.59, p-value adjusted for covariates = 0.048). Using an alternative method, we sought to further explore regulation of ECM proteins as identified by mass spec. TNR is 150 kDa protein and subject to PTMs such as glycosylation which increase its molecular weight. Western blot analysis of TNR revealed truncated TNR protein products (< 150 kDa) in all DRG samples (Fig. 3d). These truncated TNR products (< 150 kDa) were excluded from our analysis. Clear evidence of full-length TNR (150 kDa) and heavy TNR protein products (> 150 kDa) was exclusive to a subset of DPN samples (Fig. 3d). However, upon quantification of TNR (150 kDa) and heavy TNR protein products (> 150 kDa), this difference was not statistically significant (ratio (DPN/control) = 2.40, p-value = 0.1, Fig. 3d). Regardless, this reveals a TNR regulation pattern unique to samples within the DPN group.

Considering TNR is an ECM protein, we aimed to further explore regulation of this protein with visualization of its localization in the DRG. To do so, we performed immunohistochemistry with DRG tissue sections. NeuN was used to label neurons in the DRG. TNR fluorescence was observed surrounding select neurons in the DRG (Fig. 3e). We quantified the proportion of neurons in the DRG that are surround by TNR. In control DRG sections, of 453 neurons evaluated, 53 were surrounded by TNR. In DPN DRG sections, of 433 neurons evaluated, 128 neurons were surrounded by TNR. In total, 11.7% of neurons in DRG derived from control donors were surrounded by TNR, whereas 29.6% of neurons in DRG derived from the DPN donors were surrounded by TNR (p-value = 7.6e−11, Fig. 3e).

Phospho-proteomic results show changes to mRNA processing

Using the same samples as in proteomic profiling, samples were enriched for phospho-peptides and subjected to profiling by mass spectrometry. 7842 phosphorylated peptides mapping to 2533 master proteins were identified. Of these, 153 phospho-peptides were identified as differentially regulated between DPN and non-diabetic control groups (Benjamini–Hochberg adjusted p-value < 0.05, Fig. 4a). Dimensionality reduction analyses with PCA and t-SNE again did not show clear separation of samples by neuropathic pain (Supplementary Fig. S4). We again assumed an influence of covariates and we used linear regression to adjust for age and sex. 27 phospho-peptides were identified after covariate control (p-value adjusted for covariates < 0.05, Supplementary Table S3), 20 of these were regulated by a fold change greater than 1.5 (|log2(DPN/control)|> 0.58, Fig. 4b).

Figure 4

Disrupted protein translation in DPN (a) Volcano plot showing phospho-peptide abundance according to ratio of DPN/control and adjusted p-values calculated in ProteomeDiscoverer. Vertical lines at 0.58 and − 0.58 represent 1.5-fold change. Vertical lines at 1 and − 1 represent an abundance ratio of 2 or 0.5 respectively. Phospho-peptides are considered up-regulated or down-regulated with adjusted p-value < 0.05 and fold change > 1.5. Phospho-peptides are considered trending up or down with 0.05 < p-value < 0.1 or p-value < 0.05 and fold change < 1.5. (b) Heat map showing relative abundance of differentially regulated phospho-peptides (p-value after covariate control < 0.05, abundance ratio > 1.5 or < 0.67). CON indicates control samples, DPN indicates diabetic peripheral neuropathy. (c) 2D gel showing different migration patterns of SPP1 in DPN and control samples, isoelectric point (pI) is noted at the top. Due to the nature of 2D gels, each sample was run on a separate gel. Samples were run with loading controls, the positions of which were visualized with a Coomassie stain. Overlay of Coomassie stained membrane and SPP1 labeled membrane were used for pI estimation and alignment of separate membranes. Contrast was adjusted to improve visualization. See Supplementary Fig. S7 for images of full membranes. Within DPN samples, SPP1 is observed to have acidic, phosphorylated residues, whereas SPP1 in control samples has a more basic migration pattern. (d) Results of enrichment analysis with differentially regulated phospho-proteins (p-value after covariate control < 0.05). Top gene ontology annotations include regulation of mRNA stability and regulation of translation. (e) Western blot results confirmed increased phosphorylation of eEF2 with DPN. Displayed image is from one single membrane, with dividing line added to separate sample groups. p-eEF2 and eEF2 were visualized on the same membrane, with a stripping step in between.

Differential regulation of phospho-peptides could be explained by either differential phosphorylation or differential regulation of master proteins abundance. To explore this, we looked for overlap in proteomic and phospho-proteomic data. We found that not only are SPP1 phospho-peptides increased, but SPP1 protein is also slightly increased (ratio (DPN/control) = 1.32, p-value after covariate control = 0.012, Supplementary Table S3). However, abundance ratios of SPP1 phospho-peptides show larger fold changes, demonstrating phosphorylation differences. All other identified phospho-proteins were either not detected, or not significantly regulated in background proteome data (Supplementary Table S3). This confirms that identified differentially regulated phospho-proteins are in fact differentially phosphorylated.

3 phospho-peptides derived from SPP1 were upregulated by abundance ratio greater than 1.5 (pS219, ratio (DPN/control) = 1.6, p-value after covariate control = 0.017; pS263 and pS270, ratio (DPN/control) = 2.0, p-value after covariate control = 0.012; pS254 and pS263 and pS270, ratio (DPN/control) = 1.7, p-value after covariate control = 0.036, Fig. 4b). Additional phospho-peptides derived from SPP1 were also increased, however with a smaller fold change (pS375 and pS280, ratio (DPN/control) = 1.5, p-value after covariate control = 0.012; pS303, ratio (DPN/control) = 1.4, p-value after covariate control = 0.0218; pS254, ratio (DPN/control) = 1.478, p-value after covariate control = 0.02, Supplementary Table S3). Due to phosphorylation changes at multiple sites, we hypothesized that overall phosphorylation state and therefore isoelectric point of SPP1 is altered with DPN. As expected, 2D gel electrophoresis confirmed differential migration patterns of SPP1 in DPN and control groups, with more acidic migration of SPP1 in DPN samples (Fig. 4c). As phosphorylation is an acidic PTM, this supports increased phosphorylation of SPP1 in the DRG with DPN.

Using a statistical overrepresentation test in Panther, enrichment analysis of differentially regulated phospho-peptides revealed enrichment in annotations involving post-transcriptional regulation of gene expression and regulation of mRNA stability (Fig. 4d). Central to this is eukaryotic translation elongation factor (eEF2), which shows increased phosphorylation at the T57 residue (ratio (DPN/control) = 1.77, p-value adjusted for covariates = 0.046, Fig. 4b). Because this phospho-site is well studied and known to inhibit translation activity of eEF2[28,29], we selected this phospho-protein for further study. We used western blot to validate mass spectrometry findings. Western blot analysis confirmed increased phosphorylation of eEF2 at the T57 residue in DRGs from DPN donors (ratio (DPN/control) = 1.78, p-value = 0.0016, Fig. 4e).

Multi-omics integration suggests contribution of disrupted amino acid and protein metabolism in neuronal dysfunction

Using an IPA, we looked for commonalities in disease and function annotations among our omics datasets. Synthesis of protein, metabolism of protein, and synthesis of amino acids were suggested to be decreased, while uptake of amino acids was possibly identified as increased (Table 2). As expected, these annotations are drawn from decreases in amino acids (Table 2). Phospho-proteins eEF2 and SPP1 are also related to decreased synthesis and metabolism of protein (Table 2). Additionally, SPP1, eEF2, TNR, certain amino acids (NAA, Glu), and certain phospholipid metabolites (SAM, Cho) are annotated as contributing to increased progressive neurological disorder (Table 2). Altogether, this links amino acid and protein metabolism with neurological dysfunction in DPN (Fig. 5a, b). This also links changes to ECM proteins and phospholipid metabolism to neurological dysfunction in DPN (Fig. 5b).

Table 2

Top results of integrative enrichment analysis with IPA.

Diseases or functions annotation	p-value	Activation z-score	Metabolites	Proteins	Phospho-proteins
Abnormality of cerebral cortex	8.3E−5	− 1.951	Glutamate*	CD81, CELSR3, CYGB, F12, TNR*	VIM*
Apoptosis	6.1E−5	− 1.458	O-acetyl carnitine*, alpha-tocopherol*, choline, ethanolamine*, glycine*, glutamate*, glutamine*, ethanolamine phosphate*, S-adenosylmethionine*	B2M, CYGB, GHRH, HAPLN1*, IFIT3, ISG15, MTSS1*, PCP4*, SERPINB9, SPP1, STAR, VCAN	ADAM10*, EIF4B*, HNRNPC*, MAP2K4, PPP1R2*, SPP1*, UBE2V2*, VIM*, ZFP36L1
Synthesis of protein	1.9E−7	− 1.315	Creatine*, glycine*, aspartate*, glutamate*	IGFBP2, ISG15, SPP1	EEF2*, EIF4B*, LARP1, MAP2K4, PPP1R2*, SPP1*, VIM*, ZFP36L1
Metabolism of protein	7.9E−08	− 0.781	Creatine*, glycine*, aspartate*, glutamate*	B2M, CD81, F12, F9, IGFBP2, ISG15, SPP1, VCAN	ADAM10*, EEF2*, EIF4B*, LARP1, LONP1*, MAP2K4, PPP1R2*, SPP1*, VIM*, ZFP36L1
Synthesis of amino acids	4.4E−5	− 0.555	Creatine*, glycine*, aspartate*, glutamate*	VCAN
Transmembrane potential of mitochondria	7.4E−5	0.958	Glutamate*, S-adenosylmethionine*	B2M, STAR	LONP1*, MAP2K4, VIM*
Conversion of lipid	2.0E−7	1.476	O-acetyl carnitine*, alpha-tocopherol*, creatine*, glycine*, guanidoacetic acid*, glutamate*, glutamine*, S-adenosylmethionine*	CYGB, STAR
Peroxidation of lipid	2.6E−08	1.82	O-acetyl carnitine*, alpha-tocopherol*, creatine*, glycine*, guanidoacetic acid*, glutamate*, glutamine*, S-adenosylmethionine*	CYGB
Uptake of amino acids	6.5E−6	2.173	Choline, glycine*, guanidoacetic acid*, aspartate*, glutamate*	TNR*
Progressive neurological disorder	6.0E−6	2.186	O-acetyl carnitine*, choline, creatine*, glutamate*, N-acetylaspartic acid*, S-adenosylmethionine*, glycerophosphocholine*	CHCHD2, F12, IGFBP2, ISG15, SPP1, STAR, TNR*	ADAM10*, EEF2*, SPP1*, VIM*, ZFP36L1

*Indicates molecules with |log2(ratio)|> 1.5.

Figure 5

Suggested association between amino acids and DPN development. (a) Schematic showing direct relation between decreased metabolite NAA and decreased gene transcript NAT8L, created with BioRender. NAA is synthesized in the mitochondria from Asp by NAT8L. NAA is then transported to Schwann cells, where it is a precursor for myelin lipids. (b) Schematic showing suggested interplay between identified molecular changes, created with BioRender. Briefly, amino acid starvation can impair mitochondrial function as well as inhibit protein translation, resulting in ER stress. Impaired function of the mitochondria and ER can also lead to disrupted crosstalk between these organelles and thus disrupted phospholipid synthesis. Disrupted mitochondrial function can also impair NAA synthesis, impacting glial-neuronal structures. Changes to the extracellular matrix can be related to ER stress and amino acid uptake. Finally, both ER stress and amino acid depletion can induce hyperphosphorylation of eEF2, thus inhibiting protein translation. Known molecular relationships are marked with solid lines, while relationships which are suggested by literature, but not well studied are indicated with dashed arrows.

Discussion

DPN is a common form of chronic pain for which there is no effective curative treatment. Demyelination, axonal degeneration, and neuronal dysfunction and death are linked to progression of pain in DPN[8]. Here, we used a multi-omics approach to gain a comprehensive understanding of molecular alterations associated with this pathology. Given the difficulties in obtaining human DRG samples, this descriptive omics-based study reveals unique information necessary for further elucidating the etiology and progression of DPN. Our metabolomic study revealed depletion of select amino acids in the DRG in our DPN cohort (Fig. 2b, c) as well as decreased phospholipid-related metabolites (Fig. 2c, d). In proteomic and phospho-proteomic data we find changes to structural ECM proteins (Fig. 3b, c) and differential phosphorylation of RNA binding proteins (Fig. 4b, d) in the DRG with DPN. A TNR regulation pattern unique to DPN samples was also observed in western blot analysis (Fig. 3d) and IHC (Fig. 3e). Hyperphosphorylation of eEF2 (Fig. 4e) and SPP1 (Fig. 4c) in DPN were confirmed using western blot. Integration of this data using IPA suggests decreased amino acids and altered phosphorylation of RNA binding proteins contribute to altered protein synthesis and neurological disorder in the DRG in DPN (Table 2; Fig. 5a, b).

In our metabolomics analysis, we record decreases in 8 carboxylic acid and derivatives, including proteogenic amino acids Glu, Gln, Gly, and Asp (Fig. 2b). Of note, IPA reveals a link between amino acids Glu and NAA, and progressive neurological disorders (Table 2; Fig. 5a, b). Although homeostasis for most amino acids is maintained by transport processes, intracellular Glu and Asp levels are influenced most by energy metabolism[30] (Fig. 5b). Observations of decreased Gly and Gln are consistent with alterations in plasma levels of diabetic patients[31,32]. Such states of amino acid depletion have been linked to cell death and neurodegeneration[33-35]. However, whether amino acid depletion in the DRG is sufficient to produce degeneration observed in DPN is unclear. Additionally, Gly can function as an inhibitory neurotransmitter; however, a causal relation between decreased Gly in diabetes and neuronal hyperexcitability in DPN has not been studied[31].

Moreover, the non-proteinogenic amino acid NAA is decreased in our data (Fig. 2b). NAA is synthesized in neuronal mitochondria by the enzyme NAT8L[36] (Fig. 5a). Interestingly, the gene transcript for the enzyme was decreased in RNA-seq data[12]; however, this same change was not detected in our proteomic data. This discrepancy could be related to post-transcriptional regulation or sampling effects. Nonetheless, decreased NAA is important in the peripheral nervous system as it is utilized by Schwann cells to synthesize myelin[37] (Fig. 5a). Neuronal injury and degeneration decrease NAA, suggesting the differential regulation of NAA in DRG is specific to the neuropathic phenotype in DPN[27,38]. NAA’s role in myelin synthesis also provides a possible mechanism for demyelination in DPN (Fig. 5a, b).

We also report decreases in metabolites linked to phospholipid biosynthesis (Fig. 2c–e). Specifically, differentially regulated metabolites are metabolically related to PC and PE (Fig. 2c). We also see a trend with DPN towards decreased SEA, a PE derived compound that has been reported to have an analgesic effect[39]. Phospholipids are also important in mitochondrial function and myelin maintenance[40,41] (Fig. 5b). PC is a high abundance phospholipid in peripheral myelin[41]. The short half-life of this molecule[41] suggests a link between PC biosynthesis and demyelination in DPN (Fig. 5b).

Among significantly altered proteins, we observe enrichment of proteins in the ECM (Fig. 3c). The ECM has been implicated in painful conditions by regulating the neuronal microenvironment, influencing synaptic plasticity, and modulating cell signaling[42-44]. We observe a unique TNR regulation pattern with DPN, especially a greater proportion of DRG neurons surrounded by TNR with DPN (Fig. 3d, e). Western blot results show heavy TNR products (> 150 kDa) in the DRG with DPN, suggestive of glycosylation (Fig. 3d). Glycosylation of TNR is reported to vary with development, possibly impacting TNR’s adhesive properties and interactions in the ECM[45]. TNR is also noted to be related to amino acid uptake and progressive neurological disorder in IPA (Table 2; Fig. 5b). Specifically, TNR has been reported to indirectly regulate Glu uptake[46]. A direct relation between amino acid depletion and TNR upregulation in the DRG is not clear. However, increased TNR does inhibit axonal regeneration[47], thus, implicating this protein in neuronal degeneration observed in the DRG with DPN.

Among differentially regulated phospho-proteins, we report enrichment in annotations relating to RNA binding and translation (Fig. 4c). Changes in phosphorylation of RNA binding proteins is highlighted by hyperphosphorylation of eEF2 at the T57 residue (Fig. 4b, e; Supplementary Table S3). Phosphorylation at this site is well known to inhibit protein translation[29] (Fig. 5b). Interestingly, eEF2 hyperphosphorylation can occur in response to ER stress[48], a phenomenon suggested by our previous transcriptomic analysis[12]. Further evidence of ER stress can be noted by the ER’s role in phospholipid biosynthesis[49] and the evidence of decreased phospholipid biosynthesis in our data (Figs. 2c, d; 5b). ER stress is well known to contribute to chronic pain and neurodegeneration[50-52]. With the ER being a single continuous organelle, spanning the entity of the axon, long sensory neurons are particularly vulnerable to ER stress[53]. However, the myelinating activities of glial cells also make them susceptible to ER stress[54]. Further work is needed to explore cell-type specific ER-stress in DPN.

We also report hyperphosphorylation of SPP1 (Fig. 4b, c). We note a slight increase in total SPP1; however, this result was not validated, and the fold change was minimal (ratio (DPN/control) = 1.32); Supplementary Table S2). SPP1 is marker for proprioceptors in the DRG and is linked to nerve injury and mechanical pain[55,56]. SPP1 has also been reported to be linked to ER stress[57,58], inflammation in diabetes and neurodegeneration[59,60], and adhesion in the extracellular matrix[61] (Fig. 5b). However, these reports looked at total SPP1, less is known about the function of phosphorylated SPP1.

We previously reported a pro-inflammatory signature in the DRG with DPN[12]. TNR, of which we observe a unique regulation pattern with DPN (Fig. 3b, d, e), is tightly linked to neuroinflammation[62]. SPP1 is involved in inflammatory processes as well[59-61,63] although the role of SPP1 phosphorylation in these mechanisms is not known. We also find minor changes in interferon induced protein with tetratricopeptide repeats 3 (IFIT3, ratio = 1.25, p-value after covariate control = 0.024; Supplementary Table S2) and interferon stimulated gene 15 (ISG15, ratio = 1.44, p-value after covariate control = 0.017; Supplementary Table S2), both of which play a role in immune signaling. Inflammation is known to contribute to progression of both T2D and DPN[64,65]. Lack of an identification of a stronger inflammatory signature in proteomic data is likely a result of detection limits of our methods.

This study is limited by a small sample size brought on by the lack of available human DRG tissue in both disease and healthy conditions. Additionally, a lack of available clinical pain assessment of DRG donors limits the power of this study. A study with the resources to overcome these limitations has yet to be powered and would allow valuable information to build off these preliminary results[26]. Out of an abundance of caution, we used strict covariate controls, however this may have added to limitations in our statistical power. Although determined to be minimal, postmortem changes likely limited the power of this study, as did detection limits of mass spectrometry. Despite this, our study offers an important preliminary look into metabolite and protein dynamics in a disease state in the human DRG. Further work is needed to explore this regulation in larger cohort of donors. Additionally, future studies are needed to explore mechanistic contributions of identified molecules in the development of pain in DPN.

In this pilot study, we explored for the first time, metabolite, protein, and phospho-protein regulation in the human DRG. Furthermore, comparison of non-diabetic and DPN derived DRG suggests a role of disrupted amino acid and protein metabolism in neuronal dysfunction in DPN. We propose use of findings here to advise further studies with DPN. In depth understanding of molecular regulation identified here could offer novel insights into the etiology and development of DPN.

Methods

Sample acquisition

DRG used in this study were acquired from the cadaveric donors with informed consent of the next of kin (Anabios, San Diego, CA). We obtained approval for carrying out these studies from the National Institutes of Health (NIH) Office of Human Subjects Research Protection (OHSRP) and Biosafety Committee, Bethesda, MD, USA. All methods were performed in accordance with the guidelines and regulations approved by NIH Biosafety Committee, Bethesda, MD, USA. L4, L5 and S1 DRGs were collected under cold ischemic conditions, within 3 h of aorta cross-clamp. Donor medical history was obtained by trained interviewers from donor family members to the best of their knowledge. Despite lack of clinical data, organ donation offers unique and valuable insight molecular regulation in the soma of sensory neurons impacted in DPN.

Metabolomic analysis

Data acquisition

S1 DRG samples (30–50 mg) were used for metabolite analysis (Human Metabolome Technologies-America, Boston, MA). For polar metabolites, frozen tissues were homogenized in acetonitrile in water (50%) with internal standards (20 μM). The supernatant was filtered through a 5 kDa filter, centrifugally concentrated, and resuspended in ultrapure water (50 μl). CE-TOFMS was carried out using Agilent CE-TOFMS system (Agilent Technologies Inc, Waldbronn, Germany) and fused silica capillary (50 μm × 80 cm) with HMT electrophoresis buffer and HMT sheath liquid. Capillary electrophoresis was run at 30 kV. The spectrometer scanned from mass/charge (m/z) 50–1000.

For non-polar metabolites, frozen tissue was homogenized by vortexing with zirconium beads in 1% formic acid in acetonitrile containing internal standards (10 μM) and centrifuged (2300×g, 4 °C, 5 min). The supernatant was collected. The pellet was homogenized in 1% formic acid in acetonitrile and MilliQ-water (167 ul) followed by centrifugation (2300×g, 4 °C, 5 min). Supernatants were combined and filtered through a 3 kDa filter and a far filtered through phospholipid affinity column (Hybrid SPE phospholipid 55261-U, Supelco, Bellefonte, PA, USA). Filtrate was desiccated and resuspended in isopropanol (50%) in Mili-Q water. LC-TOFMS was carried out with Agilent 1200 series RRLC system SL (Agilent Technologies Inc, Wadbronn, Germany) with ODS column (2 × 50 mm, 2 μm) coupled to Agilent LC/MSD TOF MS system (Agilent Technologies Inc, Wadbronn, Germany). For chromatographic separation, mobile phase A was 0.1% HCOOH In H2O. Mobile Phase B was 0.1% HCOOH and 2 mM HCOONH4 in 65:30:5 isopropanol: acetonitrile: H2O. The gradient condition was 1% mobile phase B for 0.5 min, followed by a 13-min ramp from 1% mobile phase B to 100% mobile phase B, followed by 100% mobile phase B for 6.5 min. Flow rate was 0.3 mL/min, column temperature was 40 °C, injection volume was 1 μL, MS capillary voltage was 4 kV and 3.5 kV in ESI positive and ESI negative mode respectively, nebulizer pressure was 40 psi, gas flow was 10 L/min, and gas temperature was 350 °C. The spectrometer scanned from mass/charge (m/z) 100 to 1700.

CE-TOFMS and LC-TOFMS data was processed with MasterHands (v2.17.1.11) for identification of peaks and quantification based on internal standards. Metabolites detected in fewer than 80% of samples were excluded from downstream analysis. Remaining missing values were imputed using half minimum method[66]. Identity and concentration of select metabolites was validated using external unlabeled standards. Human Metabolome Database was used to identify metabolite class and directly related enzymes.

Proteomics

Sample prep

Tissues were prepared for proteomic analysis using EasyPep Mini MS Sample Prep Kit according to manufacturer’s instructions (Thermo Fisher Scientific). Protein concentration was measured using Pierce BCA Protein Assay (Thermo Fisher Scientific). 100 μg aliquots of protein from each sample were reduced, alkylated, and digested. Remaining protein extracts were saved for electrophoresis and western blot. Samples were labeled with TMTpro 16-plex (Thermo Fisher Scientific) according to the manufacturer’s instructions. 10 μg of each sample was combined, fractionated using high pH fractionation kit (Thermo Fisher Scientific), and used for full proteome analysis. The remaining 90 μg of each sample was pooled for phospho-peptide enrichment with PTMScan Phospho-Enrichment IMAC FE-NTA Magnetic Beads (Cell Signaling Technology) according to the manufacturer’s instructions. The flow through was then enriched according to the High-Select SMOAC protocol (Thermo Scientific Scientific). In the SMOAC protocol, samples were sequentially enriched with High-Select TiO2 Phosphopeptide Enrichment Kit (Thermo Fisher Scientific) and High-Select Fe-NTA Phosphopeptide Enrichment Kit (Thermo Fisher Scientific).

All fractions were analyzed with nano LC–MS/MS with Thermo Scientific Fusion Lumos Tribrid mass spectrometer interfaced to a UltiMate3000 RSLCnano HPLC system (Thermo Fisher Scientific, San Jose, CA). For each analysis, 1 micro of corresponding fraction was loaded and desalted in an Acclaim PepMap 100 trap column (75 µm × 2 cm) at 4 μl/min for 5 min. Peptides were then eluted into a 75 μm × 250 mm Accalaim PepMap 100 column (3 μm, 100 Å) and chromatographically separated using a binary solvent system consisting of A: 0.1% formic acid and B: 0.1% formic acid and 80% acetonitrile, at a flow rate of 300 nl/min. A gradient was run from 1% B to 42% B over 150 min, followed by a 5-min wash step with 80% B and a 10-min equilibration at 1% B before the next sample was injected. Precursor masses were detected in the Orbitrap at R = 120,000 (m/z 200). HCD fragment masses were detected in the orbitrap at R = 50,000 (m/z 200). Data-dependent MS/MS was carried out with top of speed setting, cycle time 2 s with dynamic exclusion of 20 s.

Proteome Discoverer Software (v2.5, Thermo Fisher Scientific, San Jose, CA) processed mass spectrometry data. MS spectra were searched against Homo sapiens and contaminant databases using the SEQUEST HT with PhosphoRS node for verification of phosphorylation sites. The following search parameters were used: enzyme: trypsin; maximum missed cleavage sites: 2; precursor mass tolerance: 10 ppm; fragment mass tolerance: 0.02 Da; dynamic modifications: oxidation (M), phosphorylation (S, T, Y), acetylation (protein N-terminus); static modifications: TMTpro (peptide N-terminus), carbamidomethyl (C), TMTpro (K); percolator strict FDR: < 0.01, percolator relaxed FDR < 0.05.

For quantification, samples were normalized by total peptide abundance. Spectra with > 50% isolation interference were excluded. Protein quantification was performed using unique and razor peptides, protein ratio calculation was protein abundance based, missing data was imputed using low abundance resampling. Proteome Discoverer reported background-based t-tests, with the Benjamini Hochberg correction.

2-D gel electrophoresis

Two-dimensional electrophoresis was performed with 500 ug of protein by Kendrick Labs, Inc. (Madison, WI) according to the carrier ampholine method of isoelectric focusing[67]. Isoelectric focusing was carried out in a glass tube of inner diameter 3.3 mm using 2.0% pH 3–10 isodalt Servalytes (Serva, Heidelberg, Germany) for 20,000 V-h. Tube gels were then sealed to the top of a stacking gel that overlaid a 10% acrylamide slab gel (1.0 mm) and SDS slab gel electrophoresis was carried out. The gel was then transblotted onto PVDF membrane.

Western blotting

Proteins were separated with 4 to 12% bis–tris gels for eEF2 or 7% tris–acetate gel for TNR and transferred to a PVDF membrane. Membranes were blocked and incubated overnight with primary antibodies against pT57-eEF2 (1/1000, 2331, Cell Signaling Technologies), eEF2 (1/1000, 2332, Cell Signaling Technologies), TNR (1/2000, AF3865, R&D systems), alpha-tubulin (1/10,000, ab7291, Abcam), or SPP1 (1/1000, ab8448, Abcam). Membranes were washed and incubated with HRP conjugated secondary antibodies against rabbit IgG (1/10,000, 711-035-152, Jackson Immunoresearch), goat IgG (1/10,000, 705-035-003, Jackson Immunoresearch), or mouse IgG (1/10,000, ab6789, Abcam). Chemiluminescence was detected using ECL substrate and imaged with FlourChem M (Biotechne) or Amersham Imager 600 (GE Healthcare Life Sciences) for membranes from 2D gel. For detection of total eEF2 after detection of p-eEF2, the membrane was stripped with Restore PLUS Western Blot Stripping Buffer (Thermo Fisher Scientific). Quantification was performed using ImageJ software. Contrast was adjusted using adobe illustrator to improve visualization.

Immunohistochemistry

Formalin fixed paraffin embedded DRG tissue sections were deparaffinized and subject to citrate buffered antigen retrieval. Sections were stained using a primary antibody raised against TNR (1/50, AF3865, R&D systems) and secondary antibody conjugated to Rhodamine Red-X raised against goat IgG (1/100, 705-296-147, Jackson Immunoresearch). Sections were also stained with antibody raised against NeuN and conjugated to Alexa Flour488 (1/100, MAB377X, Sigma Aldrich). Sections were mounted using DAPI fluoromount-G (17984-24, Electron Microscopy Sciences). Fluorescence was imaged using a Nikon A1R HD25 Spectral microscope. Representative images at 10× magnification were used for analysis, images for figures were acquired at 20× magnification. The proportion of NeuN expressing cells which are also surrounded by TNR was compared between groups.

Statistical analysis

Statistical analyses were performed in R, unless otherwise specified. A T-test was used to assess statistical differences in metabolomic and western blot data, a z-test was used to assess statistical differences in immunohistochemistry data. For omics datasets, the R package stats (v3.6.2) was used to apply linear regression models to candidate molecules for covariate control. Enrichment analyses were performed with molecules with p-value after covariate control < 0.05. Metabolite enrichment was performed using metaboanalyst with Small Molecule Pathway Database (SMPDB). Proteome and phospho-proteome enrichment analysis was performed separately with Panther with p-values resulting from fisher’s exact test. Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems, Mountain View, CA) was used for integrated enrichment of candidate metabolites proteins, and phospho-proteins.

Supplementary Tables.

Supplementary Information.