Age-related Effects of Heroin on Gene Expression in the Hippocampus and Striatum of Cynomolgus Monkeys.

OBJECTIVE: The aim of this study was to investigate differentially expressed genes and their functions in the hippocampus and striatum after heroin administration in cynomolgus macaques of different ages.
METHODS: Cynomolgus monkeys were divided by age as follows: 1 year (A1, n = 2); 3 to 4 years (A2, n = 2); 6 to 8 years (A3, n = 2); and older than 11 years (A4, n = 2). After heroin was injected intramuscularly into the monkeys (0.6 mg/kg), we performed large-scale transcriptome profiling in the hippocampus (H) and striatum (S) using RNA sequencing technology. Some genes were validated with real-time quantitative PCR.
RESULTS: In the hippocampus, the gene expression of A1H was similar to that of A4H, while the gene expression of A2H was similar to that of A3H. Genes associated with the mitogen-activated protein kinase signaling pathway (STMN1, FGF14, and MAPT) and γ-aminobutyric acid-ergic synapses (GABBR2 and GAD1) were differentially expressed among control and heroin-treated animals. Differential gene expression between A1S and A4S was the least significant, while differential gene expression between A3S and A2S was the most significant. Genes associated with the neurotrophin signaling pathway (NTRK1 and NGFR), autophagy (ATG5), and dopaminergic synapses (AKT1) in the striatum were differentially expressed among control and heroin-treated animals.
CONCLUSION: These results suggest that even a single heroin exposure can cause differential gene expression in the hippocampus and striatum of nonhuman primates at different ages.

INTRODUCTION

Heroin is an opioid, interacting with receptors on neuronal cells in the brain to allow euphoric effects and pain relief [1]. Heroin is an illicit drug, unlike other opioids that can be purchased by prescription. Heroin use has increased across most age groups in the United States, especially ages 18 to 25 years. This increased use has elevated the possibility of heroin abuse leading to death [1]. Asia is currently the biggest market for heroin; two-thirds of people who abuse opioids (including heroin) live in Asia [2]. Previous study has shown that heroin abuse causes deleterious effects on cognitive function and impulse control [3].

The hippocampus, where neurogenesis occurs, is responsible for learning, memory, and spatial cognition [4]. Therefore, chronic exposure to heroin can destroy learning, memory, and cognitive function [5]. Researchers have demonstrated that acute heroin exposure negatively impacts attention, motivation, and impulse control in mice, suggesting that even short-term exposure affects various functions in the brain [6,7]. In addition, mice that experience prenatal heroin exposure show impairment of learning and memory. These mice display increased neuronal apoptosis and alterations of protein kinase C activity in the hippocampus, indicating that indirect exposure to heroin negatively affects the brain [8]. In a study by Benoilid et al. [9], a single heroin inhalation caused cortical laminar necrosis in the hippocampus, as well as amnesia. However, the underlying mechanisms remain uncertain.

The striatum is composed of the caudate nucleus, putamen, and nucleus accumbens (NAc). As a core motor control center, the striatum undergoes cellular and molecular adaptations on exposure to addictive drugs [10]. Acute administration of heroin induces a rapid increase of dopamine in the striatum of animal models [11]. In contrast, the chronic use of heroin reduces the level of dopamine transporter in the striatum of humans [12]. Chronic use also impairs transcriptional events associated with glutamatergic neurotransmission in the striatum of heroin addicts [13]. The size of the NAc in heroin-dependent patients is decreased compared to that of healthy controls, indicating that alterations in apoptosis or neurotoxicity may affect volumetric changes [14]. In humans, it is difficult to analyze molecular mechanisms in areas of the brain exposed to heroin, aside from postmortem analyses. For this reason, the effects of addictive drugs are often studied by analyzing drug-induced cellular and molecular changes in the brains of animals similar to humans.

The genome of the long-tailed cynomolgus macaque (Macaca fascicularis) shows 92.8% similarity to humans [15]. In particular, their mesocortical dopaminergic system, neuroanatomical structures, neural circuits, and aging processes have many similarities to those of humans [16]. Based on these similarities, macaques have been used to study the molecular mechanisms associated with neurodegenerative diseases and drug abuse in the human brain [17,18]. To investigate the age-related effects of a single exposure to heroin on gene expression in the hippocampus and striatum, we injected heroin (0.6 mg/kg) into cynomolgus monkeys, performed large-scale transcriptome profiling using RNA sequencing (RNA-Seq), and analyzed certain functional regulatory networks. Based on our findings, we validated the messenger RNA (mRNA) expression levels of genes related to the mitogen-activated protein kinase (MAPK) signaling pathway, γ-aminobutyric acid (GABA)ergic synapses, and retrograde endocannabinoid signaling in the hippocampus. We further validated the mRNA expression levels of genes related to the neurotrophin signaling pathway, autophagy, oxytocin signaling pathway, and dopaminergic synapses in the striatum.

METHODS

Animals

Ten female cynomolgus monkeys with no history of previous participation in drug studies were used. On the basis of lifespan data for caged monkeys, we used a ratio of 1:3 when comparing ages of cynomolgus monkeys to humans [19]. The monkeys originated from Suzhou Xishan Zhongke Laboratory Animal Co. (Suzhou, China). They were maintained in individual indoor cages at the National Primate Research Center at the Korea Research Institute of Bioscience and Biotechnology (KRIBB), as described previously [20]. All procedures were approved by the KRIBB Institutional Animal Care and Use Committee (approval No., KRIBB-AEC-17079).

Heroin Treatment

Diamorphine (Johnson Matthey Macfarlan Smith, Edinburgh, Scotland) was freshly dissolved in 0.9% saline immediately before administration. Cynomolgus monkeys for heroin treatment were grouped by age as follows: 1 year (A1), 3 to 4 years (A2), 6 to 8 years (A3), and older than 11 years (A4). Ten monkeys were divided into a control group (n = 2) and four heroin-treated groups (A1, n = 2; A2, n = 2; A3, n = 2; A4, n = 2). The control group received an intramuscular injection of 0.1 ml of 0.9% saline. The heroin groups received a single-dose intramuscular injection of heroin (0.6 mg/kg) at 10:00 am. The dosage for heroin administration was determined based on previous research [21]. Control and heroin-injected groups were sacrificed after they were maintained for 4 weeks to investigate whether a single injection of heroin had an effect on brain impairments even after a certain period of time.

RNA-Seq Library Preparation and Sequencing

Total RNA was extracted from the hippocampus and striatum in control (n = 1), A1 (n = 1), A2 (n = 1), A3 (n = 1), and A4 (n = 1) animals using TRIzol (Life Technologies, Carlsbad, CA, USA), according to the manufacturer’s instructions. Total RNA integrity was assessed using a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and an RNA integrity number value greater than 8. For RNA-Seq, RNA libraries were prepared with the TruSeq RNA library preparation kit (Illumina, San Diego, CA, USA) using a modified protocol. Briefly, mRNA was purified using oligo-dT-attached magnetic beads. The purified mRNA was cleaved and reverse-transcribed into first-strand complementary DNA (cDNA) using reverse transcriptase and random primers. The RNA template was removed, followed by second-strand cDNA synthesis. End repair, 3′ adenylation, and adapter ligation were performed. The purified cDNA templates were then enriched using polymerase chain reaction (PCR). The constructed libraries were 101-bp paired-end sequenced using an Illumina HiSeq 2500 sequencer.

Differential Gene Expression Analysis

Raw reads obtained from RNA-Seq underwent quality control analyses using FastQC (version 0.10.1; http://www.bioinformatics,babraham.ac.uk/projects/fastqc/). To remove low-quality data and artifacts (including adaptor sequences, contaminant DNA, and PCR duplicates), preprocessing of reads was performed using Trimmomatic version 0.32 [22]. The preprocessed reads were mapped into a reference genome using TopHat software (version 2.1.0), and aligned reads were produced. The transcripts of each sample were assembled by Cufflinks [23] based on the fragments per kilobase of transcript per million mapped reads (FPKM) method. If there was a transcript with 0 FPKM in more than one of five samples per brain region (hippocampus or striatum), the transcript was excluded from the following analysis. To facilitate the log2 transformation, 1 was added to each FPKM value of the filtered genes. Filtered data were log2-transformed and subjected to quantile normalization. Statistical significance of the differential expression data was determined using an independent t test and fold change (FC), in which the null hypothesis stated there are no differences among groups. The false discovery rate was controlled by adjusting the p value using the Benjamini–Hochberg algorithm. For the differentially expressed gene (DEG) set, hierarchical clustering analysis was performed using complete linkage and Euclidean distance as a measure of similarity for the expression patterns of differentially expressed transcripts, while satisfying |fc| ≥ 2 and independent t test raw p < 0.05. All data and visualization of DEGs were conducted using R 3.1.2 (www.r-project.org).

Functional Annotation and Pathway Analysis

The DAVID 6.8 tool (https://david.ncifcrf.gov) was used for functional annotation and gene enrichment of transcripts that were differentially expressed in response to heroin. Statistically over-represented gene ontology (GO) categories at p < 0.05 were considered significant. GO categories were classified into three subcategories: biological process (BP), cellular component, and molecular function. After subsequently subdividing the categories into hyperlinked GO categories, we selected GO terms that showed a value of p < 0.05. To analyze the biological responses and various canonical pathways associated with DEGs, the KEGG pathway database (http://www.kegg.jp/kegg/pathway.html) was employed. The KEGG pathway analysis produced molecular and cellular networks that might involve DEGs based on previously known associations between genes or proteins, if they are otherwise independent of established canonical pathways.

Real-time Quantitative PCR (RT-qPCR)

To validate the genes that were differentially expressed in RNA-Seq, we performed RT-qPCR. Total RNA extracted from the hippocampus and striatum of controls (n = 2) and heroin-treated animals (A1, n = 2; A2, n = 2; A3, n = 2; A4, n = 2) was reverse transcribed to cDNA using the Superscript RT III system (Thermo Fisher Scientific, San Jose, CA, USA). Details of the methods for qPCR have been described previously [20]. Experiments for qPCR were performed three times independently for all samples in each group (control, A1, A2, A3, and A4) to guarantee reliable results. The differential expression of genes in each sample was normalized to GAPDH expression. The relative expression differences among control and heroin-treated groups were calculated using the 2−ΔΔCT method [24]. The primers used for the amplification of candidate genes are presented in Supplementary Table 1.

Statistical Analysis

All data obtained from qPCR are expressed as mean ± standard error of the mean. The statistical significance of differences among groups was analyzed by one-way analysis of variance (ANOVA) using SPSS software (ver. 18.0; SPSS Inc., Chicago, IL, USA). ANOVA results (followed by Tukey’s honestly significant difference post-hoc test) with p < 0.05 were considered statistically significant.

RESULTS

Identification of Genes Expressed Differentially according to Age on Exposure to Heroin

We investigated DEGs in the hippocampus (H) and striatum (S) using RNA-Seq technology after exposing the monkeys to heroin. Control animals (ConH and ConS) were not treated with heroin. We extracted 33,886 and 30,513 transcripts in the hippocampus and striatum, respectively. We mined the data from the control and heroin-treated groups for each brain region (hippocampus and striatum) using a log2 fc cutoff of 1.5 (|fc| ≥ 2). We constructed heat maps of 6,540 transcripts in the hippocampus and 8,708 transcripts in the striatum. These transcripts satisfied the |fc| ≥ 2 conditions in at least one of the total comparison pairs. Heat maps were constructed by two-way hierarchical clustering based on similarities of expression (Fig. 1A, 1B). In the hippocampus, the gene expression pattern of A1H was similar to that in A4H, while the gene expression pattern of A2H was similar to that in A3H. Of the 6,540 transcripts in the hippocampus, 5,534, 2,879, 3,544, and 4,990 transcripts were up-/down-regulated in A1H, A2H, A3H, and A4H, respectively, compared to ConH (Fig. 1C). Among the heroin-treated groups, the number of DEGs was smallest between A1H and A4H (1,222 transcripts), while the number of DEGs was largest between A2H and A4H (2,229 transcripts). Of the 8,708 transcripts in the striatum, 7,382, 6,253, 8,978, and 7,242 transcripts were up-/down-regulated in A1S, A2S, A3S, and A4S, respectively, compared to ConS (Fig. 1C). When comparing genes differentially expressed in the striatum among the heroin-treated groups, the number between A1S and A4S was the smallest (1,102 transcripts), while that between A3S and A2S was the largest (2,915 transcripts), showing a pattern consistent with the hippocampus results.

Fig. 1

Age-related alterations of gene expression in the hippocampus and striatum after exposure of monkeys to heroin. Heroin was injected into monkeys of different ages (A1, A2, A3, and A4), and differentially expressed genes (DEGs) were analyzed using RNA sequencing technology. Genes that satisfied |fold change| ≥ 2 conditions in at least one of the total comparison pairs were selected. (A) Heat map analysis of DEGs among the hippocampus (H) of control (ConH) and heroin-treated (A1H, A2H, A3H, and A4H) animals. (B) Heat map analysis of DEGs among the striatum (S) of control (ConS) and heroin-treated (A1S, A2S, A3S, and A4S) animals. (C) Comparison of DEGs among groups.

A1, monkey aged 1 year; A2, monkey aged 3–4 years; A3, monkey aged 6–8 years; A4, monkey aged over 11 years.

Gene Enrichment and Functional Annotation

We analyzed the GO annotation and KEGG pathways of genes differentially expressed among control and heroin-treated groups. In the BP category, we selected the top five networks (including the KEGG pathway). As presented in Table 1, genes associated with actin cytoskeleton organization (GO:0030036) were differentially expressed in the hippocampus of all heroin-treated groups (A1H, A2H, A3H, and A4H) compared to ConH. On the basis of our KEGG pathway analysis, genes related to retrograde endocannabinoid signaling (mcf04723) and GABAergic synapse (mcf04727) were differentially expressed in the A2H, A3H, and A4H groups compared to ConH. On the other hand, when we analyzed the GO annotation and KEGG pathways of genes differentially expressed among A1H, A2H, A3H, and A4H groups, genes related to nervous system development (GO:0007399) and neuropeptide signaling pathways (GO:0007218) were differentially expressed in A2H compared to A1H (Table 2). In addition, we found that 19 genes (mcf05010) related to Alzheimer’s disease (AD) were differentially expressed between A1H and A4H.

Table 1

Gene ontology (GO) annotation of genes differentially expressed in the hippocampus of cynomolgus monkeys exposed to heroin relative to control (unexposed) animals

Category	Term	Count	p value
A1H vs. ConH
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	10	0.007
GOTERM_BP_DIRECT	GO:0008380–RNA splicing	11	0.019
GOTERM_BP_DIRECT	GO:0015031–protein transport	16	0.024
GOTERM_BP_DIRECT	GO:0006357–regulation of transcription from RNA polymerase II promoter	10	0.037
GOTERM_BP_DIRECT	GO:0007156–homophilic cell adhesion via plasma membrane adhesion molecules	6	0.046
KEGG_PATHWAY	mcf04010: MAPK signaling pathway	78	< 0.001
KEGG_PATHWAY	mcf00970:Aminoacyl-tRNA biosynthesis	21	< 0.001
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	72	< 0.001
KEGG_PATHWAY	mcf04141:Protein processing in endoplasmic reticulum	57	< 0.001
KEGG_PATHWAY	mcf01100:Metabolic pathways	313	< 0.001
A2H vs. ConH
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	7	0.018
GOTERM_BP_DIRECT	GO:0016192–vesicle-mediated transport	11	0.020
GOTERM_BP_DIRECT	GO:0051260–protein homooligomerization	8	0.035
KEGG_PATHWAY	mcf04723:Retrograde endocannabinoid signaling	28	< 0.001
KEGG_PATHWAY	mcf04727:GABAergic synapse	22	< 0.001
KEGG_PATHWAY	mcf05033:Nicotine addiction	14	< 0.001
KEGG_PATHWAY	mcf04720:Long-term potentiation	19	< 0.001
KEGG_PATHWAY	mcf04724:Glutamatergic synapse	26	< 0.001
A3H vs. ConH
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	9	0.002
GOTERM_BP_DIRECT	GO:0006357–regulation of transcription from RNA polymerase II promoter	10	0.003
GOTERM_BP_DIRECT	GO:0051260–protein homooligomerization	9	0.029
KEGG_PATHWAY	mcf04727:GABAergic synapse	27	< 0.001
KEGG_PATHWAY	mcf04723:Retrograde endocannabinoid signaling	30	< 0.001
KEGG_PATHWAY	mcf04724:Glutamatergic synapse	31	< 0.001
KEGG_PATHWAY	mcf05033:Nicotine addiction	16	< 0.001
KEGG_PATHWAY	mcf04720:Long-term potentiation	22	< 0.001
A4H vs. ConH
GOTERM_BP_DIRECT	GO:0008380–RNA splicing	11	0.009
GOTERM_BP_DIRECT	GO:0015031–protein transport	16	0.010
GOTERM_BP_DIRECT	GO:0046907–intracellular transport	4	0.038
GOTERM_BP_DIRECT	GO:0006397–mRNA processing	12	0.043
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	8	0.047
KEGG_PATHWAY	mcf04727:GABAergic synapse	31	< 0.001
KEGG_PATHWAY	mcf04010:MAPK signaling pathway	72	< 0.001
KEGG_PATHWAY	mcf04723:Retrograde endocannabinoid signaling	33	< 0.001
KEGG_PATHWAY	mcf04141:Protein processing in endoplasmic reticulum	50	0.001
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	62	0.002

A1H, hippocampus of the monkey aged 1 year; A2H, hippocampus of the monkey aged 3–4 years; A3H, hippocampus of the monkey aged 6–8 years; A4H, hippocampus of the monkey aged over 11 years; ConH, hippocampus of control; MAPK, mitogen-activated protein kinase; tRNA, transfer RNA; GABA, γ-aminobutyric acid; mRNA, messenger RNA.

Table 2

Gene ontology (GO) annotation of genes differentially expressed in the hippocampus of cynomolgus monkeys of different ages exposed to heroin

Category	Term	Count	p value
A2H vs. A1H
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	9	< 0.001
GOTERM_BP_DIRECT	GO:0007399–nervous system development	10	0.001
GOTERM_BP_DIRECT	GO:0007156–homophilic cell adhesion via plasma membrane adhesion molecules	5	0.030
GOTERM_BP_DIRECT	GO:0015031–protein transport	11	0.032
GOTERM_BP_DIRECT	GO:0007218–neuropeptide signaling pathway	4	0.043
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	46	< 0.001
KEGG_PATHWAY	mcf01100:Metabolic pathways	187	< 0.001
KEGG_PATHWAY	mcf00970:Aminoacyl-tRNA biosynthesis	14	< 0.001
KEGG_PATHWAY	mcf04144:Endocytosis	47	0.002
KEGG_PATHWAY	mcf04010:MAPK signaling pathway	45	0.003
A3H vs. A1H
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	9	< 0.001
GOTERM_BP_DIRECT	GO:0046907–intracellular transport	3	0.046
KEGG_PATHWAY	mcf04144:Endocytosis	39	< 0.001
KEGG_PATHWAY	mcf04010:MAPK signaling pathway	34	0.004
KEGG_PATHWAY	mcf00310:Lysine degradation	11	0.008
KEGG_PATHWAY	mcf04261:Adrenergic signaling in cardiomyocytes	22	0.008
KEGG_PATHWAY	mcf00511:Other glycan degradation	6	0.014
A4H vs. A1H
GOTERM_BP_DIRECT	GO:0006511–ubiquitin-dependent protein catabolic process	6	0.024
GOTERM_BP_DIRECT	GO:0000166–nucleotide binding	12	0.017
KEGG_PATHWAY	mcf00970:Aminoacyl-tRNA biosynthesis	8	0.005
KEGG_PATHWAY	mcf04024:cAMP signaling pathway	19	0.009
KEGG_PATHWAY	mcf05010:Alzheimer’s disease	19	0.018
KEGG_PATHWAY	mcf04144:Endocytosis	21	0.033
KEGG_PATHWAY	mcf04720:Long-term potentiation	8	0.042
A3H vs. A2H
KEGG_PATHWAY	mcf04010:MAPK signaling pathway	29	0.002
KEGG_PATHWAY	mcf04144:Endocytosis	29	0.004
KEGG_PATHWAY	mcf04666:Fc gamma R-mediated phagocytosis	12	0.024
KEGG_PATHWAY	mcf00511:Other glycan degradation	5	0.026
KEGG_PATHWAY	mcf05204:Chemical carcinogenesis	10	0.026
A4H vs. A2H
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	9	< 0.001
GOTERM_BP_DIRECT	GO:0007399–nervous system development	7	0.036
GOTERM_BP_DIRECT	GO:0006897–endocytosis	5	0.036
KEGG_PATHWAY	mcf00970:Aminoacyl-tRNA biosynthesis	14	< 0.001
KEGG_PATHWAY	mcf04010:MAPK signaling pathway	42	< 0.001
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	39	< 0.001
KEGG_PATHWAY	mcf01200:Carbon metabolism	24	0.002
KEGG_PATHWAY	mcf04066:HIF-1 signaling pathway	21	0.002
A4H vs. A3H
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	6	0.007
GOTERM_BP_DIRECT	GO:0015031–protein transport	9	0.009
GOTERM_BP_DIRECT	GO:0006412–translation	16	0.024
GOTERM_BP_DIRECT	GO:0046907–intracellular transport	3	0.037
KEGG_PATHWAY	mcf00310:Lysine degradation	9	0.025
KEGG_PATHWAY	mcf00970:Aminoacyl-tRNA biosynthesis	8	0.027
KEGG_PATHWAY	mcf00760:Nicotinate and nicotinamide metabolism	6	0.043

A1H, hippocampus of the monkey aged 1 year; A2H, hippocampus of the monkey aged 3–4 years; A3H, hippocampus of the monkey aged 6–8 years; A4H, hippocampus of the monkey aged over 11 years; tRNA, transfer RNA; MAPK, mitogen-activated protein kinase; cAMP, cyclic adenosine monophosphate; HIF-1, hypoxia-inducible factor 1.

In a GO analysis of the striatum, genes associated with endocytosis (GO:0006897) were differentially expressed in A1S, A3S, and A4S compared to ConS (Table 3). Furthermore, genes associated with RNA splicing (GO:0008380) were differentially expressed in A1S, A2S, and A4S compared to ConS. Autophagy-related genes (GO:0006914) were differentially expressed in A3S and A4S compared to ConS. In our KEGG pathway analysis of the striatum, genes related to long-term potentiation (LTP) (mcf04720) were differentially expressed in A1S, A2S, A3S, and A4S compared to ConS (Table 3). On the other hand, when we analyzed the GO annotation and KEGG pathways of genes differentially expressed among A1S, A2S, A3S, and A4S, circadian entrainment-related genes (mcf04713) were differentially expressed in A2S vs. A1S in response to heroin (Table 4). Besides, we found that 25 genes (mcf05010) related to AD were differentially expressed between A2S and A4S.

Table 3

Gene ontology (GO) annotation of genes differentially expressed in the striatum of cynomolgus monkeys exposed to heroin relative to control (unexposed) monkeys

Category	Term	Count	p value
A1S vs. ConS
GOTERM_BP_DIRECT	GO:0006897–endocytosis	10	< 0.001
GOTERM_BP_DIRECT	GO:0008380–RNA splicing	14	0.003
GOTERM_BP_DIRECT	GO:0007264–small GTPase mediated signal transduction	37	0.019
GOTERM_BP_DIRECT	GO:0015031–protein transport	19	0.020
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	10	0.033
KEGG_PATHWAY	mcf04720:Long-term potentiation	34	< 0.001
KEGG_PATHWAY	mcf04722:Neurotrophin signaling pathway	54	< 0.001
KEGG_PATHWAY	mcf04141:Protein processing in endoplasmic reticulum	71	< 0.001
KEGG_PATHWAY	mcf04713:Circadian entrainment	43	< 0.001
KEGG_PATHWAY	mcf04921:Oxytocin signaling pathway	66	< 0.001
A2S vs. ConS
GOTERM_BP_DIRECT	GO:0015031–protein transport	20	0.001
GOTERM_BP_DIRECT	GO:0008380–RNA splicing	13	0.004
GOTERM_BP_DIRECT	GO:0006810–transport	18	0.024
GOTERM_BP_DIRECT	GO:0006412–translation	38	0.048
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	9	0.049
KEGG_PATHWAY	mcf01200:Carbon metabolism	52	< 0.001
KEGG_PATHWAY	mcf04141:Protein processing in endoplasmic reticulum	67	< 0.001
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	84	< 0.001
KEGG_PATHWAY	mcf04921:Oxytocin signaling pathway	62	< 0.001
KEGG_PATHWAY	mcf04720:Long-term potentiation	30	< 0.001
A3S vs. ConS
GOTERM_BP_DIRECT	GO:0006897–endocytosis	9	0.018
GOTERM_BP_DIRECT	GO:0007399–nervous system development	14	0.018
GOTERM_BP_DIRECT	GO:0006914–autophagy	12	0.019
GOTERM_BP_DIRECT	GO:0015031–protein transport	21	0.026
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	11	0.035
KEGG_PATHWAY	mcf04261:Adrenergic signaling in cardiomyocytes	73	< 0.001
KEGG_PATHWAY	mcf04720:Long-term potentiation	38	< 0.001
KEGG_PATHWAY	mcf04728:Dopaminergic synapse	62	< 0.001
KEGG_PATHWAY	mcf04722:Neurotrophin signaling pathway	61	< 0.001
KEGG_PATHWAY	mcf04015:Rap1 signaling pathway	98	< 0.001
A4S vs. ConS
GOTERM_BP_DIRECT	GO:0008380–RNA splicing	14	0.003
GOTERM_BP_DIRECT	GO:0006897–endocytosis	9	0.005
GOTERM_BP_DIRECT	GO:0015031–protein transport	19	0.016
GOTERM_BP_DIRECT	GO:0006397–mRNA processing	15	0.036
GOTERM_BP_DIRECT	GO:0006914–autophagy	10	0.047
KEGG_PATHWAY	mcf04720:Long-term potentiation	33	< 0.001
KEGG_PATHWAY	mcf01200:Carbon metabolism	54	< 0.001
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	88	< 0.001
KEGG_PATHWAY	mcf05010:Alzheimer’s disease	82	< 0.001
KEGG_PATHWAY	mcf04261:Adrenergic signaling in cardiomyocytes	62	< 0.001

A1S, striatum of the monkey aged 1 year; A2S, striatum of the monkey aged 3–4 years; A3S, striatum of the monkey aged 6–8 years; A4S, striatum of the monkey aged over 11 years; ConS, striatum of control; GTPase, guanosine triphosphatase; mRNA, messenger RNA.

Table 4

Gene ontology (GO) annotation of genes differentially expressed in the striatum of cynomolgus monkeys of different ages exposed to heroin

Category	Term	Count	p value
A2S vs. A1S
KEGG_PATHWAY	mcf04144:Endocytosis	29	0.007
KEGG_PATHWAY	mcf04713:Circadian entrainment	14	0.007
KEGG_PATHWAY	mcf00330:Arginine and proline metabolism	9	0.017
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	25	0.017
KEGG_PATHWAY	mcf04530:Tight junction	17	0.021
A3S vs. A1S
GOTERM_BP_DIRECT	GO:0006897–endocytosis	6	0.003
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	7	0.003
KEGG_PATHWAY	mcf04071:Sphingolipid signaling pathway	22	< 0.001
KEGG_PATHWAY	mcf04728:Dopaminergic synapse	20	0.004
KEGG_PATHWAY	mcf05142:Chagas disease (American trypanosomiasis)	17	0.005
KEGG_PATHWAY	mcf00230:Purine metabolism	26	0.005
KEGG_PATHWAY	mcf05203:Viral carcinogenesis	29	0.006
A4S vs. A1S
KEGG_PATHWAY	mcf05203:Viral carcinogenesis	21	0.002
KEGG_PATHWAY	mcf04612:Antigen processing and presentation	9	0.019
KEGG_PATHWAY	mcf04144:Endocytosis	21	0.024
KEGG_PATHWAY	mcf00330:Arginine and proline metabolism	7	0.033
KEGG_PATHWAY	mcf00380:Tryptophan metabolism	6	0.044
A3S vs. A2S
GOTERM_BP_DIRECT	GO:0030036–actin cytoskeleton organization	7	0.026
KEGG_PATHWAY	mcf04142:Lysosome	31	< 0.001
KEGG_PATHWAY	mcf01100:Metabolic pathways	196	< 0.001
KEGG_PATHWAY	mcf04144:Endocytosis	49	0.002
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	44	0.003
KEGG_PATHWAY	mcf01230:Biosynthesis of amino acids	20	0.004
A4S vs. A2S
GOTERM_BP_DIRECT	GO:0006914–autophagy	6	0.006
GOTERM_BP_DIRECT	GO:0055114–oxidation-reduction process	5	0.038
KEGG_PATHWAY	mcf05203:Viral carcinogenesis	25	0.002
KEGG_PATHWAY	mcf05010:Alzheimer’s disease	25	0.002
KEGG_PATHWAY	mcf04961:Endocrine and other factor-regulated calcium reabsorption	9	0.004
KEGG_PATHWAY	mcf04720:Long-term potentiation	11	0.005
KEGG_PATHWAY	mcf04728:Dopaminergic synapse	15	0.017
A4S vs. A3S
GOTERM_BP_DIRECT	GO:0006914–autophagy	7	0.001
KEGG_PATHWAY	mcf00480:Glutathione metabolism	10	0.006
KEGG_PATHWAY	mcf01100:Metabolic pathways	106	0.008
KEGG_PATHWAY	mcf01230:Biosynthesis of amino acids	12	0.016
KEGG_PATHWAY	mcf01130:Biosynthesis of antibiotics	25	0.016
KEGG_PATHWAY	mcf00563:Glycosylphosphatidylinositol-anchor biosynthesis	6	0.019

A1S, striatum of the monkey aged 1 year; A2S, striatum of the monkey aged 3–4 years; A3S, striatum of the monkey aged 6–8 years; A4S, striatum of the monkey aged over 11 years.

Validation of DEGs

Of genes differentially expressed in the hippocampus, genes related to the MAPK signaling pathway (STMN1, FGF14, and MAPT ), GABAergic synapses (GABBR2 and GAD1), and retrograde endocannabinoid signaling (GRIA2) were validated using RT-qPCR. STMN1 was significantly downregulated in the A2H group compared to the control group, while it was significantly upregulated in the A1H group compared to the control group (Supplementary Table 2 and Fig. 2). FGF14 and GABBR2 were significantly decreased in three heroin-treated (A1H, A3H, and A4H) compared to the control group, while MAPT, GAD1, and GRIA2 were significantly downregulated in all heroin-treated groups. However, expression level of GABBR2 also decreased in all heroin-treated groups compared to the control group, showing a significant decrease in A1H, A3H, and A4H (Supplementary Table 2 and Fig. 2).

Fig. 2

Validation of differentially-expressed genes identified in the hippocampus by RNA sequencing. After a single injection of heroin in monkeys of different ages, differentially expressed genes in the hippocampus (H) were analyzed using RT-qPCR. Two monkeys per group (ConH, A1H, A2H, A3H, and A4H) were used for the experiment. Genes related to MAPK signaling (STMN1, FGF14, and MAPT), GABAergic synapses (GABBR2 and GAD1), and retrograde endocannabinoid signaling (GRIA2) were analyzed. Data are represented as the mean ± standard error of the mean (n = 3). Data were analyzed using one-way ANOVA followed by Tukey’s honestly significant difference post-hoc test.

mRNA, messenger RNA; RT-qPCR, real-time quantitative PCR; MAPK, mitogen-activated protein kinase; GABA, γ-aminobutyric acid.

*Significantly different from the control group (p < 0.05).

Of genes differentially expressed in the striatum, genes related to the neurotrophin signaling pathway (NTRK1 and NGFR), oxytocin signaling pathway (CACNG3 and NFATC4), autophagy (ATG5), and dopaminergic synapses (AKT1) were validated using RT-qPCR. Expression levels of NTRK1 and NGFR were significantly downregulated in heroin-treated groups (A1S, A2S, A3S, and A4S) compared to the control (ConS) by both RNA-Seq (Supplementary Table 3) and RT-qPCR (Fig. 3). On the other hand, the expression level of CACNG3 was increased in all heroin-treated groups compared to the control, while that of NFATC4 was decreased in all heroin-treated groups compared to the control by both RNA-Seq (Supplementary Table 3) and RT-qPCR (Fig. 3). ATG5 was increased in all heroin-treated groups compared to the control group, showing a tendency to decrease with age in heroin-treated groups. Conversely, AKT1 was significantly decreased in all heroin-treated groups compared to the control (Fig. 3). Taken together, expression patterns of all the validated genes in RT-qPCR were the same as the result obtained from RNA-Seq.

Fig. 3

Validation of differentially-expressed genes identified in the striatum by RNA sequencing. After a single injection of heroin to monkeys of different ages, differentially expressed genes in the striatum (S) were analyzed using RT-qPCR. Two monkeys per group (ConS, A1S, A2S, A3S, and A4S) were used for the experiment. Genes related to neurotrophin signaling (NTRK1 and NGFR), oxytocin signaling (CACNG3 and NFATC4), autophagy (ATG5), and the dopaminergic synapse (AKT1) were analyzed. Data are represented as the mean ± standard error of the mean (n = 3). Data were analyzed using one-way ANOVA followed by Tukey’s honestly significant difference post-hoc test.

mRNA, messenger RNA; RT-qPCR, real-time quantitative PCR.

*Significantly different from the control group (p < 0.05).

DISCUSSION

The short-term effects of heroin include dry mouth, skin flushing, vomiting, nausea, and severe itching. In addition, heroin induces slowness of heart function and breathing as well as somnolence. These changes can cause coma and permanent brain damage. Despite these detrimental effects, people of all ages are illegally using heroin to induce euphoria. To identify the effects of heroin on the brain at the molecular and pathophysiological levels depending on age, we investigated global alterations in gene expression and functional networks in the hippocampus and striatum of cynomolgus monkeys, which have a brain structure and function very similar to humans.

In the present study, we found that gene expression of A1H and A4H was similar, as was that of A2H and A3H. These findings are consistent with our previous study, in which cynomolgus monkeys exposed to the psychostimulant methamphetamine showed similar gene expression patterns according to age [20]. These findings imply that age affects the way in which addictive drugs alter gene expression, regardless of the type of drug (opiate or psychostimulant). Based on our results that expression patterns of genes in A1H were closer to those of A4H than in other animals, it is supposed that hippocampi of animals at these ages show a similar response to heroin. The ages of animals in A1H and A4H are equivalent to the ages of children and middle-aged or older humans, respectively. Human children are especially vulnerable to drugs [25]. Therefore, it is believed that the hippocampus of A1H and A4H animals is vulnerable to heroin, exposure to which causes similar gene expression in A1H and A4H at a molecular level. The ages of animals in groups A2H and A3H are comparable to human adolescence, when brain and body systems experience crucial maturational rearrangements [26]. Among the heroin-treated groups, the number of DEGs between A2H and A3H was somewhat smaller than that among other age groups. This indicates that the hippocampus tissues of adolescents and young adults exhibit common molecular changes on exposure to heroin. Considering that the use of addictive drugs (including heroin) is associated with dopaminergic activation in the striatum [27], the alteration of gene expression by activated dopaminergic synapses may occur. Unlike the hippocampus, the number of DEGs between A2S and A3S in the striatum was the largest among the heroin-treated groups. Because animals in the A2S and A3S age groups experience critical maturational rearrangements, it is likely that even a single exposure to heroin induces molecular changes in the striatum.

In the present study, heroin treatment led to up- or downregulation of genes related to actin cytoskeleton organization in the hippocampus of all heroin-treated groups. The actin cytoskeleton is responsible for regulating junctional integrity and cellular remodeling under pathological and physiological states [28]. Therefore, the heroin-induced reorganization of the actin cytoskeleton can affect a variety of cellular functions in the brain, such as motility, cellular signaling, cell adhesion, and intracellular trafficking. On the other hand, A2H, A3H, and A4H groups exhibited up- or downregulation of genes related to retrograde endocannabinoid signaling (mcf04723) compared to ConH. Excitatory and inhibitory transmission in the brain leads to feedback inhibition through the activation of presynaptic CB1 receptors by endocannabinoids (eCBs) [29]. A previous report has shown that acute exposure to cocaine and alcohol results in the disruption of eCB-mediated synaptic plasticity in animal models [30]. Our findings suggest that heroin, an addictive drug, may also induce an imbalanced eCB-mediated synaptic plasticity, promoting drug-seeking behaviors and the stress response. Taken together, these results suggest that it is likely that even a single exposure of heroin disrupts synaptic plasticity.

In this study, genes associated with nervous system development were differentially expressed between A1H and A2H groups. Monkeys in groups A1 and A2 are similar in age to humans in their early childhood and teenage years. These are ages in which structural brain development and neurophysiological mechanisms (such as increase of myelination and axon density) actively occur [31]. Therefore, exposure to heroin at these ages can cause alterations in genes related to nervous system development. On the other hand, we found that 19 genes associated with AD were differentially expressed between A1H and A4H. Animals in group A4 are equivalent to humans in their mid-thirties, whose brains have already stopped growing. Even a single exposure to heroin may lead to alterations in genes related to cognitive function in the hippocampus, possibly affecting the early onset of diseases like dementia.

In the present study, interestingly, autophagy-related genes (GO:0006914) were up- or downregulated in A3S and A4S compared to ConS. Autophagy is provoked by hypoxia, metabolic stress, and starvation and is responsible for removing intracellular aggregates and cytoplasmic organelles and promoting cellular senescence [32]. Previous researchers have demonstrated that chronic exposure to opiates (including heroin) promotes autophagy in rats [33] and in human addicts ranging in age from 19 to 40 years [34]. The present study also found that a single exposure to heroin induced the alteration of autophagy-related genes in the striatum tissues of monkeys of specific ages. Some research has demonstrated that autophagy-related genes are increased in brain-derived cells (including SH-SY5Y and C6 cells) and in the hippocampus tissues of rats in response to morphine, a major psychoactive opiate [33,35]. These results indicate that opiates (including heroin) encourage autophagy, promoting neurotoxicity via autophagic cell death in the striatum. In our study, the expression of LTP-related genes was also changed in the striatum of all heroin-treated groups (A1S, A2S, A3S, and A4S). Heroin addiction is associated with upregulated LTP in the NAc [36] and hippocampus [37], followed by impairments of synaptic plasticity. Based on previous studies and our current data, we speculate that heroin promotes LTP in the brain, regardless of exposure time. Further investigations are needed to determine how LTP-related genes affect synaptic plasticity in response to heroin.

In this study, circadian entrainment-related genes were differentially expressed between A1S and A2S. A previous study has reported that circadian preferences are different between children and teenagers [38]. Considering that animals in groups A1S and A2S are age-matched with human children and teenagers, respectively, it is possible that the circadian rhythms of A1S and A2S are inherently different. Further studies are needed to assess whether the differential expression of these genes between A1S and A2S is due to an alteration of circadian features by heroin or just the result of age-specific circadian differences. In addition, as a result of comparing DEGs between A2S and A4S, there were changes of genes associated with AD. This implies that heroin affects the expression of genes associated with neurodegenerative disease differently in the growing brain versus the mature brain. In particular, some genes known to be upregulated during the development of AD [39,40] were upregulated in A4S compared to A2S (data not shown). These included Fe65/amyloid β A4 precursor protein-binding family B member 1 and microtubule-associated protein tau (MAPT). These results suggest that adults may be more susceptible than teenagers to heroin-induced neurodegenerative disease.

In the present study, STMN1 was significantly downregulated in only A2H group. STMN1 controls microtubule dynamics, in which the mitotic spindle assembles or disassembles in the last stage of cell division [41]. It has been reported that neurotoxicants (such as copper) cause a decrease of STMN1 in the hippocampus of the mouse, resulting in neurotoxicity [42]. Similarly, heroin may cause age-specific neurotoxicity by inhibiting the expression of STMN1 in the hippocampus of the monkey. On the contrary, FGF14 in our study decreased in all heroin-treated groups except A2H. FGF14 is abundant in various regions of the brain, including the hippocampus, cerebellum, and neocortex. Previous studies have reported that mutations or genetic deletions of FGF14 in animal models not only induced impairments in hippocampal synaptic plasticity and cognitive function [43], but also caused reduction of food intake and aggressivity [44]. Based on previous results and our data, the expression patterns of STMN1 and FGF14 in the A2H group are opposite those of other heroin-treated groups, indicating that heroin differentially affects neurotoxicity and hippocampal synaptic plasticity depending on age. Considering that the A2H group is matched with human early teenagers, a single exposure to heroin in this age group may cause neurotoxicity through different mechanisms compared to other age groups. However, as a result, even a single injection of heroin may indirectly induce impairments in hippocampal synaptic plasticity.

In this study, genes related to GABAergic synapses (GABBR2 and GAD1) were decreased in the hippocampus of all heroin-treated groups. This suggests the possibility that pre- and post-synaptic GABAergic receptors are affected by even a single injection of heroin. The GABBR2 gene encodes GABA B receptor 2 (GABABR2), which forms heterodimers with GABABR1 [45]. The GABAergic system is inhibited presynaptically by GABABRs, and alteration of the system is related to addictive behavior [46]. In our study, the A1H and A2H groups exhibited particularly significant decrease of GAD1 mRNA. GAD1 encodes the glutamate decarboxylase 67 protein that regulates GABA level and signaling [47]. Expression level of GAD1 is decreased in the hippocampus of patients with major depression and bipolar disorder and also in cocaine addicts [48,49]. Taken together, these findings imply that even a single injection of heroin negatively affects the GABAergic system, possibly promoting drug addictive behavior and psychiatric disorders.

GRIA2, a subtype of the glutamate AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid) receptor, is the richest glutamate receptor in the central nervous system. It regulates excitatory neurotransmission in retrograde endocannabinoid signaling [50] and inhibits calcium influx through AMPA-receptor complexes [51]. In the present study, heroin treatment induced a decrease of GRIA2 mRNA expression in the hippocampus of all monkeys. In a study by Mead and Stephens [52], gria2 knock-out mice showed deficits in learning stimulus-reward associations. This implied that a GRIA2 deficit may negatively affect behavior, contributing to drug addiction. Therefore, even a single exposure to heroin is believed to decrease GRIA2 level in all ages, leading to impairments in retrograde endocannabinoid signaling.

In the present study, a single exposure of heroin caused decrease of NTRK1 and NGFR in the striatum of all animals. NTRK1 and NGFR are genes that encode receptor proteins responsible for the development and survival of neurons [53]. In particular, NGFR was drastically decreased in all groups exposed to heroin. Nerve growth factor receptor (also known as p75NTR) is a low-affinity receptor that binds neurotrophins, leading to the development, survival, and function of neurons [54]. In a study by Busch et al. [54] that examined physiological and morphological changes in p75NTR knockout mice, the deficiency of p75NTR led to an increase of LTP in the basolateral amygdala and promoted anxiety-like behavior. In our GO analysis, we observed that LTP-related genes were changed in the striatum by heroin exposure. We can speculate that heroin may regulate LTP and negatively affect the survival of neurons in the striatum by inhibiting p75NTR.

We identified upregulation of CACNG3 and downregulation of NFATC4 in all heroin-treated animals. These two genes are associated with the oxytocin signaling pathway. Oxytocin signaling controls the hypothalamo-pituitary-adrenal axis through behavioral responses in the brain [55]. CACNG3 encodes the voltage-dependent calcium gamma-3 subunit that mediates trafficking to the somatodendritic compartment [56]. Little has been studied about the features or functions of CACNG3. NFATC4 increases neurogenesis and inhibits anxiety [57]. A deletion of NFATC4 causes partial defects in the encoding of hippocampal-dependent spatial memories [58]. Based on studies by Quadrato et al. [57,58] and on our results, the downregulation of NFATC4 in response to heroin may have negative effects on neurogenesis and memory function.

Autophagy is a degradation mechanism mediated by lysosomes and induced by drug addiction [59]. It is an important pathway for cells to protect against pathological factors. A previous study found that neurons exposed to morphine promote ATG5-related autophagy [60]. In this study, ATG5 gene expression increased in all heroin-treated groups, particularly the youngest group (A1S). ATG5-related autophagy plays a protective role in cells undergoing apoptosis [61]. Apoptotic cells exhibited an accumulation of cytoplasmic inclusion bodies in the brains of atg5-deficient mice [62]. The increase of ATG5 in the striatum may be a protective mechanism to inhibit cell death or to block the accumulation of cytoplasmic debris due to heroin. In addition, it is believed that exposure to heroin at a young age is more toxic to nerve cells, highly activating a protective mechanism.

In our study, AKT1 was significantly decreased in animals of all ages exposed to heroin. The AKT1 gene has a dopaminergic synapse-associated function. On the other hand, the AKT1/PKB protein kinase interaction plays an important role in cell survival [63]. Phosphorylation of AKT1 leads to its activation, which facilitates antiapoptotic action [64]. In a study by Ramos-Miguel et al. [65], chronic opiate addicts exhibited a decrease of phosphorylated AKT1 in the brain, while total AKT1 was not changed. This result suggests that chronic exposure to opiates inactivates AKT1, but does not block its expression. In contrast, our research found that acute exposure to heroin decreased AKT1 expression. More studies are needed to examine whether acute exposure to heroin also affects the phosphorylation of AKT1.

In conclusion, we found that a single injection of heroin induced alterations in gene expression in the hippocampus and striatum of all tested ages of cynomolgus monkeys. When analyzing DEGs in the hippocampus, genes involved in actin cytoskeleton organization were differentially expressed in all ages of heroin-treated monkeys (A1H, A2H, A3H, and A4H) compared to the ConH group. This result suggests that heroin may affect a variety of cellular functions in the brain (such as motility, cellular signaling, cell adhesion, and intracellular trafficking) by reorganizing the actin cytoskeleton. When analyzing DEGs in the striatum, genes related to LTP were differentially expressed in heroin-treated groups of all ages (A1S, A2S, A3S, and A4S) compared to the ConS group. Interestingly, heroin decreased STMN1 (involved in microtubule dynamics) expression in the only hippocampus of the A2H group but did not induce expression change of FGF14 only in the group. On the other hand, heroin decreased the expression of NTRK1 and NGFR (genes associated with neurotrophin signaling) in the striatum of all age groups. In particular, the amount of increase in ATG5 (involved in autophagy) expression by heroin tended to decrease. However, there were a few limitations to this study. For the RNA-Seq analyses, only one heroin-treated animal was used in each age group. To partially overcome this limitation, two heroin-treated animals in each age group were used to validate the results obtained from RNA-Seq and genes presenting a greater than 2-fc were selected for further analysis. To the best of our knowledge, this is the first study to analyze heroin exposure in nonhuman primates of different ages. DEG in the hippocampus and striatum were analyzed simultaneously by RNA-Seq technology.