Because corticostriatal circuit dysfunction contributes to the impairments that are a hallmark of substance abuse (21), we determined the impact of OUD on transcriptional differences in the human NAc. In line with our previous work, we found that many more NAc DE transcripts were downregulated than upregulated in OUD when collapsing all human subjects across biologic sex (Fig. 1A; 1085 down vs. 221 up). However, sex differences have been described across several domains in OUD including comorbidity with psychiatric symptoms (22), cue-induced craving/neural activity (23), and DNA methylation of OUD-related genes (24). To identify potential sex-specific transcriptional signatures of OUD, we examined the sexes separately. When disentangling the cohorts by sex, the pattern of downregulation was preserved. More transcripts were downregulated than upregulated in females (134 down vs. 14 up) and males (110 down vs. 18 up) (Fig. 1A). Most DE transcripts were protein coding (89.2%) with a much smaller fraction representing long non-coding RNAs (lncRNAs; 7.5%) and pseudogenes (3.3%).
We then investigated the effect of opioid administration in two rodent models (mice and rats) at three stages of opioid use and withdrawal. For the mouse study, we examined DE in three groups compared to naïve controls: 1) chronic self-administration and intoxicated at time of sacrifice (Intoxication); 2) chronic self-administration and in withdrawal for 24 hours (Withdrawal); 3) chronic self-administration and abstinent for two weeks (Abstinence). In the mouse NAc, opioid exposure was associated with more down than upregulated genes across all three opioid exposed groups (Fig. 1B). When separating by sex, both sexes recapitulated the pattern except for Withdrawal males, where the pattern was reversed (88 down vs 144 up). The majority of DE transcripts were protein coding genes followed by lncRNAs across the three stages of opioid exposure.
In our rat model (25), we examined the NAc core and shell separately. We examined DE in three groups compared to naïve controls: 1) chronic self-administration and intoxicated at time of sacrifice (Intoxication); 2) chronic self-administrations and in withdrawal for 24 hours (Withdrawal); 3) chronic self-administration and abstinent for 4–5 weeks (Abstinence). Unlike humans or mice who died with opioids on board (i.e., intoxicated), more DE transcripts were upregulated than downregulated in the NAc core of rats across the three opioid exposed groups (Fig. 1C). When separating by sex, the pattern of upregulation remained in females across the three stages of exposure. However, the ratio between DE transcripts was more balanced in males which could be a product of the relatively small number of DE transcripts in the NAc core of male rats. Similarly, in the rat NAc shell, there were more upregulated DE transcripts than downregulated across the three groups (Fig. 1D). However, when disentangling groups by sex, this pattern was reversed in Intoxicated and Abstinent male rats (Fig. 1D). For both rat NAc core and shell, greater than 97% of DE transcripts were protein coding across all three stages of opioid administration.
Next, we interrogated whether the pattern of gene expression in human OUD mapped onto the signature from our mouse opioid exposure groups using threshold-free RRHO analysis. These analyses have the potential to identify which transcriptional changes in human OUD are a product of the acute effects of opioid intoxication versus the long-term effects of chronic use. When the sexes in our human cohorts were pooled, there was substantial overlap in both upregulated and downregulated transcripts between humans and Intoxicated mice. However, when splitting by sex, OUD females but not males showed strong concordance with sex-matched Intoxicated mice (Fig. 2A). The same pattern was true for mice in Withdrawal (Fig. 2C) and Abstinence (Fig. 2E). Pathway enrichment was performed on shared transcripts between OUD females and female mice across the three mouse groups (Fig. 2B, D, F). Top shared pathways included Wnt/β-catenin, myelination, CLEAR (Coordinated Lysosomal Expression and Regulation), Platelet-derived growth factor (PDGF)-mediated signaling, and synaptogenesis signaling. We identified genes which act as potential upstream regulators in our mouse groups. In Intoxicated and Withdrawal mice, we identified levodopa (L-DOPA), the central precursor in the synthesis of the catecholamine neurotransmitters dopamine, norepinephrine, and epinephrine as an upstream pathway regulator (Fig. 2B, D). We also identified the gene for Estrogen receptor alpha (ERα), ESR1, as a potential upstream regulator in Intoxicated and Withdrawal mice (Fig. 2B, D). When activated by its ligands, ERα acts a transcription factor with roles across the organism (26). The HNF4A gene which also encodes a transcription factor, emerged as an upstream regulator across all three opioid exposed groups.
We then turned to opioid exposed rats and used the same approach as used for mice to interrogate the relationship of gene expression to human OUD and probe how that relationship differs from that of opioid exposed mice. Unlike in humans and mice, we looked at the NAc shell and core separately in our rat groups. In the NAc shell, the pattern of gene expression was strongly discordant between humans and Intoxicated rats when the sexes in our human cohorts were pooled. When the sexes were separated, gene expression was still discordant but more strongly in women than in men (Fig. 3A). Gene expression was strongly concordant between humans and Withdrawal rats. However, this seems to be driven by OUD men given that when the sexes were examined separately, only OUD men showed concordance with male Withdrawal rats. OUD women were weakly discordant with female Withdrawal rats (Fig. 3B). The same overall patterns held consistent between humans and abstinence rats. Gene expression was weakly concordant between all humans and Withdrawal rats but when separated by sex, males were weakly concordant, and females were weakly discordant (Fig. 3D). Pathway enrichment was performed on shared transcripts between OUD males and male Withdrawal rats (Fig. 3C). Top pathways included protein ubiquitination and mTOR signaling. We identified upstream regulators for the Withdrawal group. Surprisingly, two genes, ESR1 and HNF4A, which we identified as upstream regulators in our concordant mouse pathways, popped up here as well.
In the NAc core, identical to the NAc shell, the pattern of gene expression was strongly discordant between humans and Intoxicated rats when the sexes in our human cohorts were pooled. However, when the sexes were separated, gene expression was strongly discordant in females and weakly discordant in males (Fig. 3E). Also, like the NAc shell, gene expression was strongly concordant in OUD males and male Withdrawal rats, weakly discordant between OUD females and female Withdrawal rats, and weakly concordant between all human OUD subjects and Withdrawal rats (Fig. 3F). The same was true between human OUD subjects and Abstinence rats (Fig. 3H). Pathway enrichment was performed on shared transcripts between OUD males and male Withdrawal and Abstinence rats (Fig. 3G, I). While there were no shared pathways between Withdrawal and Abstinence groups, two of the top pathways in the Abstinence comparison, synaptogenesis and CLEAR signaling, emerged in the Withdrawal comparison in the rat NAc core, and across the three comparisons in the mouse groups, suggesting the importance of both pathways in the development of tolerance across species (Fig. 3I). We identified upstream regulators for the Withdrawal and Abstinence groups. HNF4A, emerged again across both groups which indicates the importance of this gene’s regulatory role across sex and species. We identified ESR1 as an upstream regulator in the NAc core of Withdrawal rats but not Abstinence rats (Fig. 3G, I). Upstream regulators were highly consistent across the NAc shell and core. TP53, a tumor suppressing gene which acts as a transcriptional activator, and bestradiol (estradiol), a ligand for the estrogen receptor, were identified as upstream regulators across both regions (Fig. 3C, G, I).
To investigate gene correlations in human OUD and our two rodent models, we interrogated our gene networks using WGCNA. We first built the WGCNA modules using the human OUD data and then looked for enrichment in OUD-associated DE transcripts as well as cross-species concordance patterns within each module. There was strong enrichment of DE transcripts across all human OUD subject groups in the darkorange, yellowgreen (Supplementary Fig. 4), and lightgreen modules (Fig. 4A). In addition to being enriched for DE transcripts in all our human OUD subject groups, the lightgreen module showed a high prevalence of DE transcripts that were concordant between human female OUD and mouse as well as human male OUD and rat (Fig. 4A). Because this module showed the strongest enrichment for DE transcripts, we focused our attention there. We visualized the pathway network associated with the lightgreen module and identified clusters with significantly enriched terms based on pathways (Fig. 4B). The two largest clusters were categorized as RNA splicing/spliceosome and acetylation lysine/acetyltransferase. Other network clusters included ‘DNA transcription preinitiation’, ‘regulation gene heterochromatin’, ‘cytosolic ribosome structural’, and ‘transcription coactivator activity’. Together, this network analysis underscores the profound effect opioid exposure has on the regulatory mechanisms of DNA transcription from epigenetic modifiers through post-transcriptional regulation of RNA splicing. Other notable network clusters involved synaptic activity: ‘synaptic vesicle localization’ and ‘signal release synapse.’
To distinguish potential regulators of this co-expression network, we identified highly connected hub genes within the lightgreen module that were predicted to drive the expression of other genes within the module (Fig. 4B). Of the 11 identified hub genes, 5 were specific to OUD: ANKRD11, CAMTA2, IQSEC1, RIPOR1, and SKI. In addition, we also identified sex-specific hub genes as part of our network analysis. Three of those genes overlapped OUD hub genes: CAMTA2, IQSEC1, and ANKRD11 while the other two, SEC16A and WNK2, were uniquely sex specific. The central hub gene in our network is KMT2D, a histone methyltransferase with a direct relationship with the estrogen receptor (27). This is particularly interesting given that one of the largest clusters in the lightgreen module is acetylation lysine/acetyltransferase, further suggesting opioids lead to heterochromatin modifications (15, 28).
The lightgreen module represents a potential interacting, correlative gene network that is integrated across mouse, rat, and human. To further deconvolve the cell types of represented by the genes in the lightgreen module, we assessed whether module genes were enriched for specific cell types in human striatum (29, 30). Overall, the lightgreen module was enriched for genes in neuron subtypes compared to glial cells (Fig. 5A), suggesting the module is composed mainly of genes active in neurons across species in response to opioids. In OUD, module genes were significantly reduced in enrichment within medium spiny neurons relative to interneurons and glial cells (Fig. 5B). Enrichment changes were independent of sex (Fig. 5B), occurring in both females and males with OUD. Thus, module genes were decreased in expression in medium spiny neurons in OUD, reflecting similar genes and consistent correlative expression across the gene network in medium spiny neurons across rodent models of opioids and individuals with OUD.