Angiogenic gene networks are dysregulated in opioid use disorder: evidence from multi-omics and imaging of postmortem human brain.

Opioid use disorder (OUD) is a public health crisis in the U.S. that causes over 50 thousand deaths annually due to overdose. Using next-generation RNA sequencing and proteomics techniques, we identified 394 differentially expressed (DE) coding and long noncoding (lnc) RNAs as well as 213 DE proteins in Brodmann Area 9 of OUD subjects. The RNA and protein changes converged on pro-angiogenic gene networks and cytokine signaling pathways. Four genes (LGALS3, SLC2A1, PCLD1, and VAMP1) were dysregulated in both RNA and protein. Dissecting these DE genes and networks, we found cell type-specific effects with enrichment in astrocyte, endothelial, and microglia correlated genes. Weighted-genome correlation network analysis (WGCNA) revealed cell-type correlated networks including an astrocytic/endothelial/microglia network involved in angiogenic cytokine signaling as well as a neuronal network involved in synaptic vesicle formation. In addition, using ex vivo magnetic resonance imaging, we identified increased vascularization in postmortem brains from a subset of subjects with OUD. This is the first study integrating dysregulation of angiogenic gene networks in OUD with qualitative imaging evidence of hypervascularization in postmortem brain. Understanding the neurovascular effects of OUD is critical in this time of widespread opioid use.

Introduction

In the United States, the prevalence of opioid use disorder (OUD) has rapidly increased during the last two decades and has now become a public health crisis (1). Despite measures taken to reduce opioid misuse, opioids remain the leading cause of overdose deaths, claiming the lives of over 50,000 individuals annually (2). This presents a strain on the U.S. health-care system with an annual economic burden estimated of $78.5 billion (1,3).

In addition to the detrimental effects of addiction and drug-seeking behavior, opioid use can cause short- and long-term neurotoxic effects that stem from the drugs and their metabolites. Long-term opioid exposure is known to cause epigenetic modifications, including changes in histone methylation and acetylation, DNA methylation, and miRNA expression (4). However, most studies examining opioid-induced gene-network alterations have been performed in animal models, and research in postmortem brain of human subjects with OUD is very limited. Therefore, the molecular mechanisms of gene deregulation in brains of subjects exposed to opioids are not well understood.

While 70 to 90% of the mammalian genome is transcribed into RNA, <2% of the genome represents protein-coding genes (5,6). One emerging area of research that can explain mechanisms of gene (de)regulation is non-coding RNAs, such as long non-coding (lnc)RNAs, which are largely unexplored in opioid exposure. lncRNAs are a type of regulatory RNA that are more than 200 nucleotides long, lack an open reading frame, and regulate many biological processes including embryonic stem cell differentiation and neurogenesis (7). Few studies to our knowledge have examined changes in both coding and lncRNAs following opioid exposure, and no studies have performed a comprehensive integration of transcriptome and proteomic alterations in the human brain. Importantly, most drugs currently available are designed to target proteins. Innovative targeting of specific lncRNAs greatly extends the number of druggable targets and facilitates new classes of mechanistic discoveries.

In this study, we performed comprehensive transcriptomic and proteomic profiling of dorsolateral prefrontal cortex (Brodmann area 9, BA9) in subjects with OUD and controls and found perturbations in genes and networks involved in endothelial cell function, cytokine signaling and angiogenesis, as well as differential expression of Immediate Early Genes (IEGs) which converge on angiogenic pathways. Further, using ex-vivo magnetic resonance imaging (MRI), we identified increased vascularization in postmortem brains from a subset of subjects with opioid exposure.

Materials and methods

Collection of Brain Samples

Postmortem brains were obtained from the University of Texas Health Science Center at Houston (UTHealth) Brain Collection in collaboration with the Harris County Institute of Forensic Science, with approval from the Institutional Review Board. Information obtained for each subject included demographics, autopsy and toxicology reports, and medical and psychiatric notes. A detailed psychological autopsy was obtained for each subject by interviewing the next of kin, where information regarding psychiatric clinical phenotypes (evidence of depression, mania, psychosis), age of onset of drug use, types of drugs used, smoking and alcohol use history, and any co-morbidities was obtained. A consensus diagnosis of OUD or non-psychiatric control was reached for each subject according to DSM-5 criteria by an independent panel of three trained clinicians after review of all available information.

In total, brains were collected from 29 OUD and 18 non-psychiatric controls. Of these, samples of 20 OUD subjects and 12 controls were used to generate proteomics data, samples of 27 OUD subjects and 14 controls were used to generate ribonucleic acid (RNA) sequencing data (RNAseq), and tissue from 11 OUD subjects and 5 controls were used for MRI analysis (Supplementary Table 1). All OUD subjects except 1 were positive for opioids on toxicology at time of death. Cause of death of OUD subjects was opioid or polysubstance overdose (N=22), cardiovascular disease (N=4), blunt force trauma (N = 1), hepatic fibrosis (N = 1), and undetermined (N=1). Cause of death for controls was cardiovascular disease (N=13), pulmonary embolism (N=2), chest gunshot wound (N=2), and splenic artery aneurysm (N = 1).

Postmortem interval (PMI) was calculated from time of death until tissue preservation. Upon receipt, the right and left hemispheres were coronally sectioned in ~1.5 cm slices. Right hemisphere sections were immediately frozen and stored at -80ºC. Left hemisphere slices were fixed in 4% paraformaldehyde (PFA) and stored at 4ºC. BA9 was defined within the prefrontal cortex between the superior frontal gyrus and the cingulate sulcus. Dissections were obtained using a 4mm punch through the cortex, yielding approximately 100mg of tissue. Cerebellar pH was measured (8) (Table 1). As described below, all analyses were adjusted for the covariates of age, sex, PMI, cerebellar pH, and RNA integrity number (where applicable).

Table 1:

Summary of Patient Demographics

Demographics	Control	OUD
Total	18	29
Males (%)	16 (89%)	16 (58%)
Age mean (sd)	53 (14)	39(12)
PMI mean (sd)	28.3 (6.9)	26.2 (8.9)
Ethnicity (White/Black/ Hispanic/Asian)	10/5/2/1	24/5/0/0
Toxicology at Time of Death	Control	OUD
Negative	18 (100%)	0
Codeine/ morphine/ heroin	0	18 (62%)
Oxycodone/ hydrocodone/ hydromorphone/ oxymorphone	0	12 (41%)
Tramadol	0	4 (14%)
Fentanyl/ Methadone/ U-47700	0	6 (21%)
Cocaine	0	5 (17%)
Amphetamines	0	1 (3%)
Benzodiazepine	0	13 (45%)
Cannabis	0	5 (17%)

Proteomics quantification, differential expression analysis, and LGALS3 ELISA

We first performed proteomics analysis in BA9 from 20 subjects with OUD and 12 controls, based on sample availability at the time of this analysis. Fifty milligrams of BA9 tissue was lysed in radio-Immunoprecipitation assay (RIPA) buffer. Extracted proteins were reduced in 200mM tris(2-carboxyethyl)phosphine (TCEP) (Thermo Fisher Scientific) diluted with 50 mM tetraethylammonium bicarbonate (TEAB) (Sigma) at 55°C for 30min, and then alkylated with 375 mM iodoacetamide in 50 mM TEAB in the dark at room temperature for 45 min. Proteins were diluted with acetone and centrifuged at 15000g for 20 min at 4°C. Precipitated proteins were washed with 1 ml of ice-cold tri-n-butylphosphate/acetone/methanol (1:12:1 by volume) for 90min and centrifuged at 2800g for 15min at 4°C. 50ug of protein were solubilized with 5% SDS, 50mM TEAB, pH 7.55 and sample was centrifuged at 17,000g for 10 min. A solution of 12% phosphoric acid and 90% methanol in 100mM TEAB was then added, and the proteins were passed through a S-Trap spin column using a bench top centrifuge (30s spin at 4,000g). Trypsin was added to the protein mixture in a ratio of 1:25 in 50mM TEAB, pH=8, and incubated at 37○C for 4 hs. Peptides were eluted sequentially with 50mM TEAB, 0.2% formic acid, 50% acetonitrile, and 0.2% formic acid. The combined peptide solution was then dried in a speed vac and resuspended in 2% acetonitrile, 0.1% formic acid, 97.9% water and placed in an autosampler vial.

Nanoflow liquid chromatography tandem mass spectrometry (NanoLC MS/MS) was performed using a nano-LC chromatography system (UltiMate 3000 RSLCnano, Dionex), coupled on-line to a Thermo Orbitrap Fusion mass spectrometer (Thermo Fisher Scientific, San Jose, CA) through a nanospray ion source (Thermo Scientific). A trap and elute method was used with C18 PepMap100 (300um X 5mm, 5um particle size) and Acclaim PepMap 100 (75um X 25 cm) columns (Thermo Scientific). Samples were eluted (400 nL/min) by gradient elution onto the C18 column.

LC-MS/MS data were acquired using XCalibur version 2.1.0 (Thermo Fisher Scientific) in positive ion mode using a top speed data-dependent acquisition (DDA) method with a 3 sec cycle time. The survey scans (m/z 350–1500) were acquired at 120,000 resolution (at m/z = 400) in profile mode, with a maximum injection time of 50 msec and an AGC target of 400,000 ions. The raw mass spectrometry data files were processed using MaxQuant. Peptide identifications were accepted if they could be established at greater than 95.0% probability. MS/MS spectra were searched against the Swiss-Prot/UniProt human database (9). In total, 4,612 unique proteins were identified in all the samples. The mass spectrometry proteomics raw data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD025269.

Data was normalized and imputed, and DE analysis was performed using protein-wise linear models and empirical Bayes statistics implemented in the DEP (Differential Enrichment analysis of Proteomics data) (10) and limma (11) R packages, followed by correction for multiple testing using Benjamini-Hochberg. The linear model design formula consisted of control/OUD diagnosis, age, PMI, sex, and cerebellar pH (~0+Diagnosis+age+PMI+pH+sex). Differentially expressed proteins were identified based on the raw or the adjusted p-value of <0.05 combined with a fold change absolute value |FC| > 1.5 cutoff. Levels of the protein LGALS3 were validated on the same samples used for LC-MS at 625 nanograms of protein per sample using the Human LGALS3 enzyme-linked immunoassay (ELISA) kit (abcam ab269555), following manufacturer’s protocol.

Next-generation RNA Sequencing (RNAseq) and Differential Expression Analysis

We next performed RNAseq including additional available samples to those used in the proteomics analysis. RNA was extracted from 50mg of BA9 tissue using RNeasy Plus Mini kits (Qiagen, Hilden, Germany) and RNA integrity number (RIN) was measured for RNA quality (Agilent Bioanalyzer 2100 system, Agilent Technologies, Santa Clara, CA). 27 OUD and 14 controls passed RNA quality measures and were used for RNAseq. 1 μg RNA per sample was used for mRNA library construction using NEBNext® Ultra™ RNA Library Prep Kit (Illumina Inc, San Diego, CA). Paired-end sequencing reads (150bp) were generated on an Illumina HiSeq2000 platform (Novogene Bioinformatics Institute, Chula Vista, CA).

Sample sequencing quality was assessed using FastQC. Reads were mapped to the GRCh38 genome with TopHat v2.1.0 4 and Cufflinks v2.2.1 using default parameters. The annotation file was downloaded from NCBI (GCF_000001405.39_GRCh38.p13.gtf). Gene and transcript counts were obtained using htseq-count. From the raw counts, gene lengths and library size, counts per million (CPM), fragments per kilobase per million (FPKM), and transcripts per million (TPM) were calculated. To avoid inflation ratio, the expression of genes with FPKM and TPM values < 0.1 were set to 0.1.

Differential expression (DE) analysis was performed using R Bioconductor packages edgeR and limma (11,12). Genes with CPM >1 in more than 5 individuals were kept for analysis. Read counts were normalized by Trimmed Means of M (TMM) (13). The linear model design formula consisted of control/OUD diagnosis, age, PMI, RIN, sex, and cerebellar pH (~0+Diagnosis+age+PMI+RIN+pH+sex). Differentially expressed genes (DEGs) were determined using the following cutoffs: absolute fold change |FC|>1.5, Benjamini-Hochberg adjusted p-value with false discovery rate (FDR) <0.2 or FDR <0.05 and at least one samples’ FPKM >1.

Quantitative reverse transcription polymerase chain reaction

Quantitative reverse polymerase chain reaction (qPCR) for EGR1 and TGFB2 was performed. Complementary DNA (cDNAs) were synthesized from using the High Capacity cDNA Reverse Transcription kit (Thermo Fisher Scientific). Targets were amplified in triplicate using Taqman® primer/probe sets (EGR1 assay ID:Hs00152928_m1, TGFB2 assay ID: Hs00234244, GAPDH assay ID: Hs99999905_m1) and TaqMan® Gene Expression Master Mix (Thermo Fisher Scientific). Relative mRNA expression was calculated (14) with GAPDH as the reference gene. Pearson correlation was calculated between the targets’ logged qPCR expression (-ΔCT) and RNAseq logTPM.

For TGFB2-OT1, cDNAs were synthesized from using iScript gDNA Clear cDNA Synthesis Kit (Bio-Rad, Cat #1725035) and random primers. PowerUp SYBR Green Master Mix (Applied Biosystems, Cat# A25780) was used for qPCR detection. Primer pair sequences and thermocycling conditions previously reported for TGFB2-OT1 were used (15).

Cell-type deconvolution analysis and Gene-Cell type correlation

Cell type ratio scores were computed using CIBERSORTx (16,17), with the following settings: relative quantification, 100 permutations, quantile normalization disabled. We modified a published signature matrix reference of cell types (18) adapted from single-cell RNAseq of human brain (19), which included constituent cell types of the adult prefrontal cortex: astrocytes, neurons, oligodendrocytes, microglial cells, and endothelial cells (Supplementary Table 2). The R package BRETIGEA (20) was used to compute relative cell type proportion surrogate variables for the same cell types, using the top 50 cell-type specific genes for each cell type (Supplementary Table 2). Gene-cell type linear correlation was calculated using the Pearson correlation coefficient of each surrogate variable and the TPM-normalized counts for each gene, considered significant if correlation was positive and BH FDR-corrected p-value < 0.05.

Gene Ontology, Pathway Analysis, and WGCNA

Gene ontology (GO) overrepresentation analysis was performed on the DE RNAs and proteins using R package TopGO (21) to query GO annotation in 3 categories: biological processes (BP), molecular function (MF), and cellular component (CC). Significant GO terms in each category were determined by Fisher’s exact tests using the TopGO weight correction algorithm and filtered for GO terms with > 2 genes represented. Pathway analysis and corresponding network diagrams was conducted using Ingenuity® Pathway Analysis (IPA) software (QIAGEN, Redwood City, CA) on the combined list of differentially expressed (DE) RNAs and proteins.

Weighted correlation network analysis (WGCNA) analysis was performed using the R WGCNA package (22) on normalized transcripts per million (TPM) values after covariates (age, sex, PMI, cerebellar pH, RIN) were regressed out using a linear mixed model of covariates and diagnosis created by the jaffelab cleaningY function (23). Briefly, after filtering for low variance genes (variance >3) we derived a soft thresholding power of 8 with a scale-free model fit cutoff of 0.9. A total of 9,558 genes were used in the WGCNA package with parameters power=8, minModuleSize=30, reassignThreshold = 0, mergeCutHeight=0.25. Modules were correlated with BRETIGEA cell type surrogate variables and diagnosis, and heatmaps were created using WGCNA package and pheatmap (24). Hub genes were defined as the top 10% most connected genes in each module that also were correlated with a diagnosis of OUD (abs(correlation) >0.4) and had high module membership (MM > 0.85). Network diagrams of hub genes were constructed in Cytoscape (25) for visualization. GO overrepresentation of biological processes (GO-BP) of module hub genes was performed using the same methods as described above for the DE RNAs and proteins.

Ex vivo Magnetic Resonance Imaging

Fixed (> 1 month) postmortem brain slices from 11 OUD and 5 controls (four coronal slices, covering from the prefrontal to occipital lobes) were imaged. Time in fixative (PFA) was calculated from date of tissue preservation to date of imaging. The 4 brain slices of each subject were grouped together with spacers in between each slice. Brain slices from up to three subjects were securely caged in a cylindrical acrylic imaging container with acrylic plates separating the group of slices from different subjects. The container was then filled with Fluorinert (FC-3283, 3M, Saint Paul, MN) and put under a vacuum of ~0 bar for 15 min to minimize air bubbles. Imaging was performed on a 3T human scanner (Prisma, Siemens, Erlangen, Germany) using a 64-channel head/neck coil. The protocol included multi-echo gradient echo (ME-GRE) acquisitions (0.47 mm isotropic resolution, TR / TE1 / TE2 / TE3 / TE4 = 22.0 / 2.40 / 7.48 / 12.56 / 17.64 ms, flip-angle = 30° / 20° / 10°, 55 min. acquisition time for each flip-angle). ME-GRE images were combined by weighted root-sum-of-squares, and T1, T2*, and apparent proton density (S0) maps were jointly fitted based on the GRE signal equation using Python (Python 3.6.9, SciPy 1.4.1 module, curve_fit function).

T2* weighted MRI is highly sensitive to the magnetic susceptibility changes introduced by iron in the blood clots, which can be leveraged to visualize blood vessels in ex vivo MRI (26). Blood vessels were visually identified as contiguous linear hypointensity in the combined weighted ME-GRE image or the T2* map by consensus of two radiologists who were blinded to subject diagnosis, led by an experienced (>10 years) neuroradiologist.

Results

Differential expression of genes and proteins identified enrichment of angiogenesis pathways

Table 1 shows demographic information of subjects. Groups were matched for PMI, pH, and RIN and both consisted of primarily White/Caucasian race/ethnicity. OUD subjects were significantly younger than controls and had a higher female/male ratio. Variance partition analysis of normalized (TPM) RNAseq and protein expression (Supplementary Fig. 1b and c, respectively) revealed that the covariates age, sex, pH, PMI, and RIN, all adjusted for in the differential expression analyses, each contributed 1 to 5 percent of the variance, while OUD diagnosis contributed to less than 1 percent of gene variance.

We used two sets of criteria for differential expression analysis. Setting the threshold as adjusted p-value < 0.2 and |FC| > 1.5, we identified 345 coding genes (230 upregulated, 115 downregulated) and 49 long-noncoding RNAs (lnc-RNAs) (35 upregulated, 14 downregulated) (Supplementary Table 3). As depicted in the volcano plot (Fig. 1a), top downregulated RNAs included the immediate early genes (IEGs) EGR1, EGR2, EGR4, NR4A1, NR4A2, ARC, and NPAS4, previously identified in transcriptome studies of long-term opioid exposure, and the lncRNA CNIH3-AS2 (4). Top upregulated RNAs included growth factor genes FGF2, TGFB2, and TGFBI, and lncRNAs TGFB2-OT1 and SOX2-OT (Fig. 1a).

Fig. 1

Differential expression of RNAs and proteins in BA9 of OUD subjects.

(A) Fold change versus adjusted p-value for upregulated (red) and downregulated (blue) DEGs with selected genes from referenced networks and pathways labeled. Significance threshold for pathway analysis (padj < 0.2 and abs(Log2FC) > 0.58) indicated by dashed line. (B) Top: Select enriched gene ontology (GO) terms in biological processes (BP), cellular component (CC) and molecular function (MF). Bottom: Upregulated (red) and downregulated (blue) DE RNAs attributed to angiogenesis, cytokine signaling, and extracellular matrix GO terms. (C) Top: Venn diagram showing number of unique genes represented by RNAseq pipeline, proteomics pipeline, or both. Bottom: Four genes differentially expressed in both RNAseq and proteomics data, with columns indicating upregulation (red) or downregulation (blue) in respective dataset. (D-H) Correlation of RNAseq log2TPM-normalized counts plotted against logged qPCR expression (-ΔCT) values for EGR1, TGFB2, TGFB2-OT1, VAMP1 and LGALS3. (I) Correlation of normalized LC-MS/MS protein expression with galectin-3 (LGALS3) protein concentration as determined by ELISA. (J) Fold change versus unadjusted p-value for upregulated (red) and downregulated (blue) DE proteins, selected genes from referenced networks and pathways labeled. Threshold used for pathway analysis (p-value < 0.2 and abs (Log2FC) > 0.58) indicated by dashed line. (K) Top: Select enriched gene ontology (GO) terms in biological processes (BP), cellular component (CC) and molecular function (MF). Bottom: DE proteins attributed to endothelial function and stress fiber GO terms, all upregulated (red).

Functional enrichment analysis of the DE RNAs revealed enrichment of pathways involving cytokine signaling, angiogenesis, and extracellular matrix organization (Fig. 1b, Supplementary Table 4). Genes involved in angiogenesis included TGFBI, TGFB2, and FGF2. Of interest, TGFB2 and the lncRNA TGFB2 overlapping transcript 1 (TGFB2-OT1) were both upregulated in OUD (Fig. 1a). Analysis with a more stringent threshold of adjusted p < 0.05 and |FC| > 1.5 reduced the number of differentially expressed genes to 10, including EGR1 (Supplementary Table 3, bold).

We next explored expression patterns in the proteome. 90.7% of proteins identified by LC-MS had their corresponding coding RNAs detected by RNAseq (Fig. 1c). We identified 160 upregulated and 53 downregulated proteins (Fig. 1j, Supplementary Table 5). Of these, 4 were also identified at the RNA level: LGALS3, SLC2A1 and PLCD1 exhibited the same direction of change in protein and RNA expression; VAMP1 was downregulated in RNA and upregulated in protein (Fig. 1c). Enrichment analysis revealed DE proteins were related to processes of endothelial cell migration and proliferation and inflammatory response (Fig. 1k, Supplementary Table 6). Analysis with adjusted p-value < 0.05 and |FC| > 1.5 revealed differential expression of 3 proteins: ATP5J2, LGALS3, and TAF15 (Supplementary Table 5, bold).

The highly differentially expressed early response gene EGR1, the angiogenic TGFB2 and TGFB2-OT1, and the RNA/Protein convergent genes LGALS3 and VAMP1 were selected for validation of significant findings. qPCR expression levels of the 5 genes correlated with normalized transcript counts (logTPM) (Fig. 1d-h). Additionally, LGALS3 protein levels as determined by ELISA correlated with normalized and imputed LC-MS protein expression data (Fig. 1i).

Using Ingenuity Pathway Analysis (IPA) we identified networks and pathways enriched in the combined lists of DE RNAs and proteins (Supplementary Table 7). The top network produced by this analysis was enriched in cellular functions of apoptosis (p = 4.5E-13), migration of endothelial cells (p = 9.8E-10), and angiogenesis (p = 1.3E-7) (Fig. 2, Supplementary Table 7 bold). Several important mediators in this network, including DE RNA TGFB2, were also identified by GO to be related to angiogenesis and endothelial cell function (Fig. 1b). IPA also showed enrichment of DE RNAs and proteins involved in the canonical opioid receptor signaling pathway (Supplementary Fig. 2). Other top enriched canonical IPA pathways of the combined list of DE RNAs and proteins included Protein Kinase A signaling and ERK/MAPK signaling. (Supplementary Table 7)

Fig. 2

Convergent pathways of differentially expressed genes and proteins.

Left: Network identified Ingenuity Pathway Analysis (IPA). Node shapes represent the functional class of the molecule: kinase (inverted triangle), enzymes (vertical rhombus), transcription regulators (horizontal ellipse), cytokines/growth factor (square), transmembrane receptor (vertical ellipse), and complex/group/others (circle). Solid and dotted lines denote direct or indirect interaction of two genes, respectively. Arrows indicate action of a gene on a target. Right: Differentially expressed genes/proteins are labeled as upregulated (red) or downregulated (blue).

Cell type deconvolution reveals correlation of genes with specific cell types

Cell type composition can influence RNAseq data obtained from bulk brain tissue (27). Using CIBERSORTx (16) and BRETIGEA (20), we obtained absolute cell type ratios and relative cell type proportions for astrocytes, neurons, oligodendrocytes, endothelial cells, and microglia (Supplementary Fig. 3, Supplementary Table 2). Although no significant microglia population was detected using CIBERSORT, we calculated relative microglial proportions using the BRETIGEA algorithm. We found no significant differences in the absolute or relative cell type proportions of the major detected cell types between groups (Fig. 3b, Supplementary Fig 3b, Welch’s T-test, p > 0.05). However, we found that cell type contributed to a large portion of the RNA expression variance (Fig. 3a, Supplementary Fig. 1b). We therefore explored the linear correlation between each BRETIGEA relative cell type proportion and the DE RNAs. Hierarchical clustering revealed that most DEGs correlated with either neurons, astrocytes, microglia or endothelial cells (Fig. 3c, Supplementary Table 2), and many of these formed part of the angiogenesis pathway identified by GO (TGFB2, FGF2, SEMA3E, IFITM1, IFITM2, IL33, COL4A1 and COL4A2, NR4A1 and TBX1) (Fig. 3d).

Fig. 3

Cell type deconvolution and WGCNA.

(A) Principal component analysis (PCA) of RNAseq data reveals that subjects align along the first principal component (PC1, 15.05% of data variance) based on neuronal RNA composition (Neuron Ratio) identified by CIBERSORT. (B) Left: CIBERSORT absolute cell type ratios for astrocytes, endothelial cells, neurons, and oligodendrocytes. Right: Boxplots of OUD versus control for each cell type. Horizontal line indicates median, box represents interquartile range (25th to 75th). Whiskers represent 1.5 times IQR above the 75th or below the 25th quartile. (C) Heatmap showing strength of positive (red) or negative (blue) Pearson correlations between TPM of each DE RNA and BRETIGEA relative cell type proportions or OUD diagnosis. Hierarchical clustering of the genes by correlation reveals clusters positively correlated with astrocytes or endothelial cells (cluster 2), and neurons (cluster 3). Cluster 1 contains genes correlated with oligodendrocytes or no particular cell type. (D) Select DE RNAs that are positively correlated with individual cell types, from pathways of interest. (E) Heatmap of Module eigengenes (ME) identified from WGCNA, colored by strength of positive (red) or negative (blue) Pearson correlation between MEs and BRETIGEA cell type proportions or diagnosis. Labels indicates Pearson correlation with p values in parentheses. (F) Networks of WGCNA hub genes identified in 3 modules highly correlated with OUD diagnosis and endothelial/astrocyte/microglia cells (Blue), neurons (Green), or endothelial cells alone (Brown). DE RNAs outlined in black. Tables indicate top 3 enriched GO BP terms and p-values of corresponding network gene hubs for Blue and Green hub genes. Brown module had no significantly enriched GO-BP terms.

WGCNA reveals hub genes related to OUD pathology and cell-type specific networks

WGCNA identified 5 modules (blue, green, brown, tan, and cyan) significantly correlated with OUD and with endothelial, neuronal, microglial, or astrocyte cell type (p-value < 0.05; Fig. 3e, Supplementary Table 8). Hub genes of the three modules exhibiting the strongest correlation with the OUD phenotype (blue, green, brown) were found to incorporate many of the DE RNAs that were involved in signaling pathways identified by GO (Fig. 3f). The blue endothelial/astrocyte/microglia module and the green neuron modules were enriched in cytokine signaling and synaptic vesicle GO biological pathway (BP) genes, respectively (Supplementary Table 9). The brown endothelial-correlated module was not significantly enriched in any GO-BP pathway. The tan and cyan modules were not enriched in any GO-BP pathway and there were no correlated DEGs as hub genes in these modules (data not shown).

Ex vivo MRI indicates hypervascularization in OUD subjects

To explore structural changes in brains of opioid users, anatomical and quantitative (T1 and T2*) MRI was performed on fixed brain slices. Specifically, blood vessels can be visualized in ex vivo MRI (26). Based on qualitative analysis of images, we identified a hypervascularization phenotype, shown as distinct and prominent linear hypointensities, in 5 out of 11 OUD subjects (Fig. 4, Supplementary Fig. 5), particularly in the basal ganglia (e.g., globus pallidus, nucleus accumbens, and putamen). Most of the observed vascularization resided within the white matter, with some projection to the cortex, including the frontal cortex. This phenotype was not observed in any control brains imaged. One control showed clustered periventricular vessels radiating from the ventricle, appearing as a large block of vessels in coronal view. Proliferating periventricular blood vessels were also a prominent feature in OUD, although they were not existent in all cases. Supplementary Table 1 provides a description of the radiological observations that were used to define the vascularization phenotype. There were no differences in the time of fixation of imaged tissues between OUD and controls or between OUD cases with or without the observed vascularization phenotype (p > 0.1, two-tailed Wilcoxon rank sum, Supplementary Fig. 4).

Fig. 4

Ex Vivo MRI of Postmortem Brain.

Representative ex vivo MRI T2* map of postmortem brains slices (midbrain to prefrontal lobe from left to right in each case) show increased vasculature in OUD (A, 40yr and B, 70yr) as compared to age matched control (C, 53yr and D, 69yr) cases. Brown arrows point to most discernable blood vessels shown as contiguous hypointensity in T2* map.

Discussion

Preclinical studies show substantial neurobiological alterations caused by opioid exposure, including changes in gene expression, likely mediated by epigenetic modifications, in various regions of the brain (4). However, the long-term consequences of opioid exposure in the human brain remain poorly understood. In this study, we utilized RNAseq and proteomic technologies, coupled with a sophisticated bioinformatics pipeline, to comprehensively investigate the PFC transcriptome and proteome of OUD subjects. Going beyond profiling, we used ex vivo MRI to examine potential structural alterations in a subset of the same postmortem brains.

Differential RNA and protein expression followed by pathway analysis revealed enrichment of pathways related to angiogenesis, cytokine signaling and endothelial cell migration and proliferation in OUD. Common to the enriched networks are the IEGs, including EGR1, EGR2, EGR4, NR4A1, NR4A2, ARC, and NPAS4, which have been implicated in OUD pathology (28,29). In line with our finding, a recent study found NPAS4 and EGR4 downregulated in dorsolateral PFC of opioid-exposed subjects (30). IEGs are known to modulate endothelial cell function and angiogenesis downstream of inflammatory and growth factor signaling pathways (31,32). In particular, EGR1, downregulated in our study, regulates angiogenesis via modulation of pro-angiogenic FGF2 signaling (33–36). Importantly, we found angiogenesis-related differential expression at both the RNA and protein levels, including LGALS3, SLC2A1, PCLD1 and VAMP1. The galectin genes, including LGALS3, increase endothelial cell proliferation and induce angiogenesis (37–40). SLC2A1 is a major glucose transporter in the brain, which is upregulated in response to oxycodone in animal studies, and is likely involved in regulating angiogenesis and endothelial cell metabolism (41,42). Altogether, these findings suggest angiogenesis dysregulation in OUD.

In addition to coding genes, we evaluated differential expression of long-noncoding RNAs, which have been largely unexplored in addiction disorders. We identified several differentially expressed lncRNAs, including TGFB2-OT1, SOX2-OT, and CNIH3-AS2. TGFB2-OT1 overlaps with the differentially expressed coding RNA TGFB2. It regulates endothelial cell functions, including inflammation, autophagy, and survival via micro-RNA interaction (15,43). While the functions of SOX2-OT and CNIH3-AS2 are not currently known, CNIH3 is known to be involved in glutamatergic AMPA receptor signaling, which is strongly implicated in OUD (44); SOX2, an important developmental gene that modulates cell stemness, is upregulated in the presence of morphine (45,46). Further studies are now warranted to investigate the mechanisms by which the lncRNAs identified in this study regulate transcriptional activity in OUD.

Opioids are known to modulate the immune system and cytokine signaling in CNS and peripheral tissues (47). Recent studies (48,49) have also identified transcriptional and epigenetic dysregulation of cytokine signaling in postmortem dorsolateral prefrontal cortex of opioid users, corroborating our findings. Cytokine signaling is well-known to mediate angiogenesis, particularly by targeting endothelial cells (50). We found that many of the DEGs involved in angiogenesis and cytokine signaling correlated with endothelial cells, microglia, and astrocytes. The IEGs NR4A1 and NR4A2 strongly correlated with the endothelial cell type. These genes regulate vascular endothelial cell proliferation and angiogenesis (51–54), likely through VEGF signaling mechanisms (55). The DEGs TGFB2 and FGF2 function via signaling of brain endothelial cells by astrocytes, are involved in cytokine signaling, and may be important for modulating cerebral angiogenesis (56–58). TGFBI, a microglial-correlated gene, interacts with the endothelial and astrocyte-correlated DEGs COL4A1 and COL4A2, and this interaction is important in signaling endothelial cell migration and proliferation (59). In addition, we found the neuroinflammatory microglial correlated gene KCNK13, known to be crucial for modulating microglial surveillance and cytokine signaling in the brain (60).

WGCNA further revealed correlation of DEGs and networks with specific cell types. In particular, many of the DE RNAs were found to be hub genes within an astrocyte/endothelial/microglia correlated module. We also identified an endothelial-correlated module containing the hub gene NPAS3, which has been previously implicated in psychiatric disorders (61).

The qualitative ex vivo MRI evidence of increased vascularization in a subset of OUD brains, shown as distinct and prominent linear hypointensity in ME-GRE images or T2* maps, could be interpreted as a visual corroboration of the molecular evidence of angiogenesis dysregulation in some OUD subjects. Angiogenesis is considered an important process in the recovery from ischemic stroke (62). In vivo MRI studies have previously reported neurovascular complications characterized by ischemia in opioid abusers (63,64). Interestingly, the globus pallidus is the most commonly involved location of the ischemic events (65) and this is one of the main locations of our findings in ex vivo tissue analysis. Ischemia is the most commonly observed acute neurovascular complication caused by heroin (66) and opioid use is associated with increased risk of stroke and cardiovascular events (67). Opioid antagonists have even been suggested as potential therapeutics for ischemic stroke (68). In addition to acute ischemia, some studies have demonstrated an increase in white matter T2-weighted hyperintensities in opioid abusers, which has been attributed to the chronic vasoconstrictor effects of the mu-opioid receptor on the vascular smooth muscle (69,70). Finally, the vasculitis-type changes observed in some studies are thought to be the result of immune-mediated responses (62,70), which could be mediated by genes involved in cytokine signaling, such as those identified in this study. In this regard, studies in animal models suggest that opioids induce angiogenesis via MAPK/ERK signaling (71–73). The angiogenic effect may depend on interaction with specific receptors, such as mu versus kappa opioid receptors (70, 74–76).

While the ex vivo MRI findings likely represent long-term ischemic changes in brain vasculature, we cannot definitely determine whether they are ischemic in nature. Further, it cannot be clearly determined whether in opioid abusers they are caused by opioids, or by other disorders or factors not considered here. In this regard, it is important to highlight that evidence of vascularization was only observed in a subset of OUD brains imaged. Therefore, the imaging findings could perhaps represent an extreme phenotype of OUD, and are meant as a qualitative corroboration to the molecular results of our study. Blood vessel characterization with larger sample sets are needed to validate the imaging evidence of brain pro-angiogenic effects of opioids observed in this study.

Limitations of our study include a relatively small sample size and the use of bulk brain tissue for analysis. Studies with larger sample sizes are necessary to investigate the effects of covariates controlled for in our analysis, including sex and age, as well as the potential effects of the different opioids or other substances used by some subjects. Future studies using single-cell RNA sequencing technology are needed to validate the gene correlation findings in specific cell types. In addition, it is not possible in postmortem brain to determine whether the gene expression and network alterations we observed are the result of an acute overdose or of chronic exposure to opioids. Additional studies are needed to clarify the causal mechanisms of angiogenic dysregulation and neurovascular damage in OUD. Finally, a more detailed investigation of the observed vascularization using ex vivo MRI may reveal clinically relevant information for radiological practice.

In summary, we have provided convergent transcriptomic and proteomic evidence for disruption of genes and pathways related to CNS angiogenesis in opioid-exposed brains. To our knowledge, this is the first study relating dysregulation in angiogenic gene networks with qualitative evidence of hypervascularization in human postmortem brains. Our results pointing towards vasculature damage as a consequence of opioid abuse provide fundamental new insights into the molecular mechanisms underlying neurobiological alterations in OUD.