Environmental enrichment alters protein expression as well as the proteomic response to cocaine in rat nucleus accumbens

Prior research demonstrated that environmental enrichment creates individual differences in behavior leading to a protective addiction phenotype in rats. Understanding the mechanisms underlying this phenotype will guide selection of targets for much-needed novel pharmacotherapeutics. The current study investigates differences in proteome expression in the nucleus accumbens of enriched and isolated rats and the proteomic response to cocaine self-administration using a liquid chromatography mass spectrometry (LCMS) technique to quantify 1917 proteins. Results of complementary Ingenuity Pathways Analyses (IPA) and gene set enrichment analyses (GSEA), both performed using protein quantitative data, demonstrate that cocaine increases vesicular transporters for dopamine and glutamate as well as increasing proteins in the RhoA pathway. Further, cocaine regulates proteins related to ERK, CREB and AKT signaling. Environmental enrichment altered expression of a large number of proteins implicated in a diverse number of neuronal functions (e.g., energy production, mRNA splicing, and ubiquitination), molecular cascades (e.g., protein kinases), psychiatric disorders (e.g., mood disorders), and neurodegenerative diseases (e.g., Huntington's and Alzheimer's diseases). Upregulation of energy metabolism components in EC rats was verified using RNA sequencing. Most of the biological functions and pathways listed above were also identified in the Cocaine X Enrichment interaction analysis, providing clear evidence that enriched and isolated rats respond quite differently to cocaine exposure. The overall impression of the current results is that enriched saline-administering rats have a unique proteomic complement compared to enriched cocaine-administering rats as well as saline and cocaine-taking isolated rats. These results identify possible mechanisms of the protective phenotype and provide fertile soil for developing novel pharmacotherapeutics. Proteomics data are available via ProteomeXchange with identifier PXD000990.


INTRODUCTION
It is clear that humans display vast individual differences in susceptibility to drug addiction. Some people become addicted after a single exposure yet others can be resistant to addiction even after many exposures to high doses of drug. Understanding the mechanisms leading to resistance to addiction will provide new targets for the treatment or even the prevention of addiction. To identify novel targets, addiction science needs good animal models of resistance to addiction.
The environmental enrichment paradigm is an animal model of resistance to addiction. In this paradigm, rats are assigned either to an enriched condition (EC) with daily exposure to novelty (children's toys), exercise, and social contact with conspecifics (i.e., group housing), or an isolated condition (IC; single housed with no novelty). Our work and others' has shown that enrichment increases sensitivity to the locomotor activating and conditioned place preference effects of stimulants yet produces a decrease in intravenous self-administration of cocaine or amphetamine (Bowling and Bardo, 1994;Bardo et al., 1995;Green et al., 2002Green et al., , 2010Thiel et al., 2009Thiel et al., , 2010. This decrease in self-administration has been described for acquisition, maintenance, extinction and reinstatement phases, using fixed and progressive ratios and in both prevention (i.e., before drug exposure) and treatment (during abstinence) models.
Currently there are no FDA approved pharmacotherapeutics for cocaine addiction in the United States despite decades of study of currently known targets. Thus, identifying completely novel targets for pharmacotherapeutic and genetic intervention is paramount to the success of developing a viable treatment for addiction. Accordingly, the first step in the process is determining intrinsic differences between the EC and IC brain and then examining how each group responds to cocaine. Our prior EC/IC cocaine research has focused on gene transcription and transcription factors (Green et al., 2010;Pavlovsky et al., 2013); the current project utilizes an unbiased discovery-based proteomic method looking for differential expression of proteins in the nucleus accumbens of EC and IC rats self-administering cocaine or saline with the goal of discovering proteins and/or biochemical pathways that could serve as targets in the treatment of cocaine addiction.
The analyses employed for this study exploit known functional relationships to identify over-representation of related protein groups in two complementary analyses. First, the Ingenuity Pathways Analysis (IPA) interrogates a curated database of published protein functional relationships (Kramer et al., 2014) to identify coordinated regulation of protein lists derived using a traditional p-value cutoff. Second, a Gene Set Enrichment Analysis (GSEA; in this case using proteins) uses a running-sum statistic of all rank-ordered proteins (i.e., no p-value threshold) (Subramanian et al., 2005).

SUBJECTS
Male Sprague-Dawley rats (40) arrived at 21 days of age and were randomly assigned to either an EC or an isolated condition to remain throughout the experiment. Enriched rats were housed 10 per cage in a large cage (77 × 78× 60 cm) with 14 hard plastic children's toys replaced daily. Isolated rats were single housed in standard polycarbonate cages. Isolated rats were chosen as the control group rather than pair-housed rats because pair-housed rats represent an intermediate enrichment condition. Rats had ad lib access to food and water and were kept on a 12 h light/dark cycle. All experiments were conducted during the light phase of the cycle. All animal procedures were approved by the UTMB Institutional Animal Care and Use Committee and conform to the NIH Guide for the Care and Use of Laboratory Animals.

SELF-ADMINISTRATION
To minimize differences in self-administration among EC and IC rats, all rats were first trained to bar-press for sucrose pellets (Noyes 45 mg). Prior to sucrose training (beginning PND51) rats were brought down to 85% of their free-feed body weight over 5 days. Rats were first trained on a continuous schedule of reinforcement in standard two-lever operant chambers (Med-Associates, St. Albans, VT); the schedule was incremented daily until the rats were on an FR5 schedule. All sessions lasted 15 min. After rats learned the FR5 schedule, they were allowed to regain their 100% free feed body weight prior to catheter surgery. For surgery, rats were anesthetized with ketamine (100 mg/kg IP) and xylazine (10 mg/kg IP) and a Silastic catheter (Fisher Scientific, Pittsburgh, PA) was fed into the jugular vein down to the heart. The catheter exited from the back from a stainless steel threaded pedestal. After 1 week of recovery, rats were placed back into the operant chamber and allowed to self-administer 0.5 mg/kg cocaine (National Institute on Drug Abuse, Bethesda, MD) or saline given in 0.1 ml over 5.8 s under a continuous schedule of reinforcement. Sessions lasted 2 h daily for 14 days, and the session terminated when the rat received 30 infusions to eliminate EC/IC differences in cocaine intake that could confound proteomic results. Control rats self-administered saline. Thus, the study was a 2 (EC/IC) X 2 (cocaine/saline) design. The number of sessions (14) was chosen because sensitization, tolerance and transcriptional changes are all stable at 14 days (Miller et al., 1998;Alibhai et al., 2007;Green et al., 2008). For the final session (PND88), nucleus accumbens tissue was harvested 3 h after the start of the self-administration session.

PROTEIN EXTRACTION
Nucleus accumbens (NAc) tissue samples from the right side of two rats in same group were combined for protein extraction, washed with ice cold tris-buffered saline (TBS), homogenized on ice in a buffer containing TBS pH 7.4, 1% Igepal-CA630 (NP-40), 1× protease inhibitor cocktail, 20 mM NaF, 1 mM Na 3 VO 4 , 10 mM DTT, and 5 mM EDTA, and then centrifuged at 750 g for 20 min at 4 • C. The pellet, containing the nuclear fraction, was set aside. The top fraction was transferred to new tubes and centrifuged at 20,000 g for 20 min at 4 • C. Methanol/chloroform (v/v 1:4:1) was added to the resulting supernatant, and it was kept at room temperature 15-30 min, vortexing every 5 min, and centrifuged at 16,000 g for 20 min at 4 • C. Finally, 500 uL of acetone was added to remove the methanol and chloroform. The 20,000 g pellet fraction was treated similarly to remove lipids (membrane fraction). All fractions were then dissolved in a buffer containing 6 M urea, 1% NP-40, 20 mM Tris-HCl pH 7.4, 1× protease inhibitor cocktail, 10 mM DTT and 5 mM EDTA, and stored at −80 • C until further analysis.

PREPARATION OF SAMPLES FOR LC-MS
Equal amount of protein extracts from each stock (100 µg) was digested with trypsin (sequence grade modified trypsin, Promega, trypsin:protein = 1:333 molar ratio), and incubated at 37 • C overnight. The peptide mixture was extracted with a C18 tip (Pierce), and eluted with 90% acetonitrile in 0.1% TFA. The extracted peptides were dried and resuspended in 5 mM Tris, pH 7.4 (25 µL) for LC-MS analysis.
All LC-MS/MS data were acquired using XCalibur, version 2.1.0 (Thermo Fisher Scientific). For high-resolution data-driven analyses (DDA), the survey scans (m/z 300-2000) (MS) were acquired in the Orbitrap at 60,000 resolution (at m/z = 400) in profile mode, followed by six CID fragmentation MS/MS spectra, acquired in centroid mode. The automatic gain control targets for the Orbitrap were 5 × 10 5 for the MS scans and 1 × 10 4 for MS/MS scans. The maximum injection times for the MS1 and MS/MS scans were 300 and 500 ms, respectively. For MS/MS acquisition, the following settings were used: parent threshold = 50,000; isolation width = 2.0 Da; normalized collision energy = 35%; activation time = 30 ms. Monoisotopic precursor selection, charge state screening, and charge state rejection were enabled, with rejection of singly charged and unassigned charge states. Dynamic exclusion was used to remove selected precursor ions (±10 ppm) for 60 s after MS/MS acquisition. A repeat count of 1 and a maximum exclusion list size of 500 was used. The following ion source parameters were used: capillary temperature 275 • C, source voltage 2.0 kV, source current 100 uA, and S-lens RF level 63%.

PROTEOMIC DATA ANALYSIS
MS files (.raw) were imported into Progenesis LC-MS (version 4.1, Nonlinear Dynamics) for m/z and retention time alignment. The top 5 spectra for each feature were exported (charge deconvolution, top 1000 peaks) as a combined.mgf file for database searching in PEAKS (Ma et al., 2003) (version 6, Bioinformatics Solutions Inc., Waterloo, ON) against a combined UniProtKB/SwissProt rat-mouse database (September 2013 version, 20,264 proteins), appended with the common Repository of Adventitious Proteins (cRAP) contaminant database (The Global Proteome Machine, http://www.thegpm.org/cRAP/index.html). PEAKS DB (Zhang et al., 2012) and Mascot (version 2.3.02, Matrix Science) searches were performed with a parent ion tolerance of 20 ppm, fragment ion tolerance of 0.8 Da, and variable modifications of oxidation (M) and phosphorylation (STY). Trypsin was specified as the enzyme, allowing for 2 missed cleavages and a maximum of 3 PTMs per peptide. An additional search for unexpected modifications was performed with the entire Unimod database (Zhang et al., 2012). Finally, homology searching was performed using the SPIDER algorithm (Han et al., 2005) to identify peptides resulting from non-specific cleavages or amino acid substitutions. Mascot and PEAKS SPIDER searches were combined (inChorus), using a 1% false discovery rate cutoff for both search engines. The resulting peptide-spectrum matches (95% peptide probability) were imported into Progenesis LC-MS. Conflict resolution was performed manually to ensure that a single peptide sequence was assigned to each feature by removing lower scoring peptides. The resulting normalized peptide intensity data were exported, and the peptide list was filtered to remove non-unique peptides and all modified peptides except methionine oxidation. For quantification, the filtered list of peptide intensities was imported into Dante-R (Polpitiya et al., 2008;Karpievitch et al., 2009) (version 0.1.1), log 2 transformed, and combined to protein abundances (RRollup) using the default settings, including one-hit wonders (50% minimum presence of at least one peptide, minimum dataset presence 3, p-value cutoff of 0.05 for Grubbs' test, minimum of 5 peptides for Grubbs' test). The resulting proteins were quantified by Two-Way ANOVA, and p-value adjustment was performed according to the Benjamini and Hochberg method (Benjamini and Hochberg, 1995).

GENE SET ENRICHMENT ANALYSIS (GSEA) ANALYSES OF PROTEIN
The IPA analyses involve deciding on a P-value cutoff for deciding if proteins are significantly differentially expressed. Highly stringent cutoffs introduce Type II error (false negatives) and less stringent cutoffs introduce Type I errors (false positives). To mitigate this problem, normalized protein intensities were analyzed via GSEA using the Broad Institute's GSEA software (http://www. broadinstitute.org/gsea). Protein names were converted to gene names by IPA. The advantage of the GSEA is that all proteins are ranked (not just the statistically significant) and significance is determined using a running-sum statistic for the whole gene set (Subramanian et al., 2005). Analyses included proteins from all 3 fractions. Gene sets were taken from the Broad Institute's (v.4.0) set of curated gene sets [C2.all; including Biocarta, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome], Gene Ontology gene sets (C5.all) and transcription factor target genes (C3.tft). Analyses were run with 1000 permutations of gene sets (size 15-500) using the Signal2noise ranking metric.

RNA SEQUENCING
The left nucleus accumbens from each rat was purified using the RNeasy kit using the manufacturer's directions. cDNA libraries were created by reverse transcribing the RNA and creating the second strand. Blunt ends were phosphorylated and "a-tailed" so that adapters could be ligated to both ends. RNA was sequenced with a HiSeq 1000 system from Illumina. Four samples were added to each lane. cDNA was amplified using "bridge" amplification. Base calls were made using fluorescently labeled nucleotides. Over 100 million reads were taken at 50 bp (paired-end reads). FastQC (v0.9.1) (Andrews, 2014) was used to check the quality of the reads. Reads were mapped to the rat reference genome (RN4) using Tophat2 (v2.0.4) (Kim et al., 2013) and Bowtie2 (v2.0.0.6) (Langmead and Salzberg, 2012) software packages. The R package EdgeR (v.3.0.8) (Anders and Huber, 2010;Robinson et al., 2010) was then used for analysis using the "trimmed mean for M-values" (TMM) method for normalization and tag-wise dispersion using "count" data. A likelihood ratio F-test was used for generating P-values to compare EC vs. IC rats. Transcripts significantly regulated (p < 0.05) were overlaid onto the IPA energy metabolism pathway for comparison to protein values.

RESULTS
The sessions were limited to 30 infusions to eliminate differences in self-administration among groups. The two rats with the slowest acquisition from each group were dropped from the experiment, leaving eight rats in each group. The remaining cocaine rats all acquired self-administration within 3 days and all rats received the maximum 30 infusions for each of the last five sessions (data not shown). There was also no EC/IC difference in time to complete 30 infusions. In several instances, data acquisition problems in the LC-MS/MS runs prevented proper alignment of samples for label-free quantification. In those cases, three samples per condition (corresponding to protein from six mice) were used for quantification and statistical analysis. Table 1 contains a selection of significantly regulated proteins for the cocaine main effect, environmental enrichment main effect and interaction. A complete list of all identified proteins can be found in Supplementary Materials. Overall, 52 proteins were significantly regulated (p < 0.05; fold change >1.5) as a function of cocaine self-administration, with 11 proteins decreased and 41 proteins increased. In the case of environmental enrichment, 117 proteins were significantly regulated, with 50 decreased and 67 increased. The interaction between environmental enrichment and cocaine produced 128 significant proteins. For both the cocaine and environmental enrichment main effects, the protein exhibiting the greatest decrease was cysteine sulfinic acid decarboxylase (CSAD). Moreover, this protein had the largest positive effect size of the significant proteins in the Enrichment X Cocaine interaction term.

Confirmatory evidence
In discovery-based science it is important to have independent (i.e., orthogonal) validation of a subset of the results to provide greater confidence in the experiment as a whole. Luckily, the effect of cocaine on protein expression has been studied extensively. As a result, the published literature provides many validating measures. For cocaine main effect, the top-ranked biological functions ( Figure 1A) were Amino Acid Metabolism and Molecular Transport [−log(p) = 5.18 for each]. Amino acid metabolic pathways lead to production of the neurotransmitters glutamate and dopamine, known to be important in cocaine addiction (see Figure 1C). Another highly ranked function, Cellto-Cell Signaling and Interaction [−log(p) = 4.84], represents another important facet of the effects of cocaine on the brain and is also manifested in the significant canonical pathways as determined by IPA ( Figure 1B). The top-ranked canonical pathway in IPA was RhoA Signaling [−log(p) = 5.05; Figure 1B], a pathway which has been demonstrated to show decreased activity in the NAcc by treatment with cocaine (Kim et al., 2009;Gourley et al., 2011). This result was also supported by the GSEA protein analysis showing significant regulation of the RhoA gene set curated by the Protein Interaction Database [Normalized Enrichment Score (NES) = −1.8, p < 0.005; Figure 1D; Table 2]. The 4thand 7th-ranked canonical pathways were CREB Signaling in Neurons [−log(p) = 3.1] and ERK/MAPK Signaling [−log(p) = 2.9], respectively ( Figure 1B; also see network in Figure 1F). These two signaling pathways are well known addiction-related   signaling pathways (Berhow et al., 1996;Green et al., 2006;Larson et al., 2011). Further, the significant biological function named Behavior was comprised of Conditioning (p < 0.001), Emotional behavior (p < 0.005), Cocaine-seeking behavior (p < 0.01) and Place preference (p < 0.05) (see Figure 1E). Lastly, BDNF was identified as a significant upstream regulator (p < 0.01), confirming what is known about this protein in cocaine self-administration (Graham et al., 2007).

Glutamate/dopamine release
Our results show a significant increase for several proteins involved in glutamate and dopamine signaling-notably, we have significant increases in vesicular glutamate transporters (VGLUTs) 1 and 2 (SLC17A7 and SLC17A6, respectively) as well as the vesicular monoamine transporter VMAT2 (SLC18A2; Figure 1C).

Protein ubiquitination
Ubiquitin C (UBC) is a major hub in the network shown in Figure 1G. Note that expression for 6 of the 8 ubiquitin target proteins was increased by cocaine.

Confirmatory evidence
Although not at the very top of the list of regulated pathways, the enrichment main effect IPA analysis also produced significant canonical pathway regulation of CREB signaling in neurons [−log(p) = 2.3] and ERK/MAPK signaling [−log(p) = 2.1; Figure 2B], supporting our own prior research (Green et al., 2010;Fan et al., 2013a,b). Also of interest is the regulation of multiple proteins associated with mood disorders ( Figure 2C) as a subset of Psychological disorders (Figure 2A). In addition, the energy metabolism results described below confirm and greatly extend upon our earlier work (Fan et al., 2013a,b).

Ubiquitination and heat-shock proteins
EC rats had increased expression of several proteins in the Protein ubiquitination canonical pathway [−log(p) = 3.0; Figure 2D], particularly heat-shock proteins. These results agree with our prior research showing increases in heat-shock proteins in EC rats (Fan et al., 2013a).

Splicing and mRNA processing
EC rats exhibit decreased expression of proteins related to splicing, as measured by the GSEA protein analysis using the KEGG Spliceosome gene set (NES = −2.13, p < 0.001; Figure 2E; Table 3). This effect was driven, to a large degree, by U2-interacting small nuclear ribonucleoproteins and splicingfactor 3 proteins. In addition to the KEGG spliceosome gene set, the largely overlapping Gene Ontology (GO) gene set for Ribonucleoprotein complex was similarly regulated (NES = −2.28, p < 0.001; Figure 2F). In fact, the top seven scoring GO gene sets downregulated in EC rats were all related to splicing (NESs range from −2.28 to −1.59, p-values all less than 0.006).

Transporters
One gene set robustly upregulated in EC rats was the GO gene set Substrate-specific transporter activity (NES = 1.74, p < 0.001; Figure 2G; Table 3). Proteins contributing to this effect were largely solute carriers (SLCs) and subunits of the Na+/K+ exchange pump. In contrast, calcium channels (CACNs) tended to be downregulated in EC rats.

Kinase activity
Another gene set robustly upregulated in EC rats is the GO gene set Protein kinase cascade (NES = 1.68, p < 0.005; Figure 2H; Table 3). Two of the top three proteins in the set (MINK1 and DBNL) are involved in receptor trafficking and endocytosis.

Energy metabolism
One of the most regulated biological functions in EC vs. IC rats was Energy production (p = 2.12E-10; Figure 2A). Likewise, 4 of the top 5 canonical pathways were related to energy production (see Figure 2B) Figure 3A depicts the consistent upregulation of these proteins in EC rats (particularly EC control rats) for the 4 canonical pathways. In addition, the significantly-regulated biological function Energy production was driven to a large degree by proteins related to ATP synthesis and degradation ( Figure 3B). In support of the IPA results, the gene set with the highest normalized gene set enrichment score (NES) in the GSEA protein analysis was the Reactome TCA cycle and respiratory electron transport gene set (NES = 1.96, p < 0.001; Figure 3F; Table 3). The IPA upstream regulator analysis identified PPARGC1A as being activated by environmental enrichment (activation Z-score = 3.26, p = 6.8E-7). In addition, ESRRA (activation Z-score = 2.21, p = 2.96E-08) was also identified as an upstream regulator, with similar energy metabolism proteins in both results ( Figure 3C). A subsequent literature search revealed that PPARGC1A interacts with ESRRA to increase expression of energy metabolism proteins (Ichida et al., 2002;Luo et al., 2003), offering good evidence for this mechanism driving the robust increases in energy metabolism proteins shown in Figure 3A. In support of the ESRRA results, a GSEA protein analysis identified a significantly upregulated gene set of proteins induced by ESRRA in a paper by Stein et al. (2008) (Figure 3G). It is important to note that the Stein et al. data are not part of the Ingenuity Knowledgebase, so the complementary IPA and GSEA results are completely orthogonal. In addition to PPARGC1A and ESRRA, the insulin receptor (INSR) was identified as a significant upstream modulator, with a clear majority of these proteins being related to energy metabolism (activation Z-score = 3.57, p = 9.4E-8; Figure 3D). The INSR, when activated, increases glucose uptake into cells (Weinstein et al., 2009) and signals through the PI3K/and ERK signaling cascades, both of which were identified by IPA as significant hubs [−log(p) = 1.8 and 2.13, respectively]. It is likely through transcriptional effects of these two mechanisms that INSR also increases expression of energy metabolism proteins.

ENRICHMENT BY COCAINE INTERACTION
The interaction of enrichment X cocaine is potentially important because it shows how EC and IC rats react differently to cocaine.

Rho signaling
The second-and third-ranking canonical pathways were RhoA signaling [−log(p) = 7.0] and Signaling by Rho family GTPases [−log(p) = 4.6], respectively ( Figure 4B). Inspection of protein intensities revealed that IC control rats have less Rho-related proteins than the other groups and that IC rats showed induction after cocaine ( Figure 4C).

Kinase signaling cascades
Several protein interaction networks were identified for the interaction term, including those in Figures 4D,E scoring 52 and 33, respectively. Both of these networks have hubs of known kinases (ERK, PI3K, p38, etc.) and downstream transcription factors (CREB, NFKB) all known to be important in the response to cocaine (Berhow et al., 1996;Russo et al., 2009;Larson et al., 2011). Further, our prior research has shown some of these to be important in environmental enrichment (Green et al., 2010;Fan et al., 2013a,b). The interpretation of these results is that EC and IC rats show differential kinase responses to cocaine, a result intriguing given that EC rats have a protective addiction phenotype.

Ubiquitination
Prior research from this laboratory demonstrated that EC and IC rats had differential expression of ubiquitin target proteins and these rats exhibited differential regulation of ubiquitin targets after restraint stress. The IPA analysis of the current data also suggests ubiquitination may be important for the differential expression of proteins in EC and IC rats after cocaine. Two different high-scoring networks (scoring 29 and 24) were identified with UBC as the major hub (the higher-scoring network is shown in Figure 4F).

Energy metabolism
Closer inspection of the heat maps from Figures 3F,G reveals that the EC main effect for energy metabolism proteins is mostly driven by EC saline-administering rats (left 3 columns). Not surprisingly, the IPA analysis for the interaction of enrichment and cocaine produced a significant result for Energy production (p = 6.14E-4; Figure 4A). For canonical pathways, TCA cycle II [−log(p) = 3.3] and Glycolysis I [−log(p) = 3.2] ranked 11th and 12th among regulated pathways, respectively ( Figure 4B). When the GSEA data are reanalyzed comparing the EC saline vs. all other groups, 20 gene sets related to energy production and mitochondria were statistically significant, all with NES scores greater than 1.75 and p-values < 0.005. Thus, EC control rats had drastically higher levels of energy metabolism proteins than all other groups. In addition, a set of genes comprising genes with the ESRRA response element in the promoters also produced a significant result (NES = 1.49, p < 0.005). It is notable that this gene set is completely orthogonal to the Stein et al. and the IPA knowledgebase, meaning that 3 independent analyses suggest EC saline rats have higher expression of ESRRA target proteins. Figure 5 shows the same IPA pathway as Figure 3A, only this time overlaid with environmental enrichment main effect RNAseq data. As with the protein data, one of the top-regulated canonical pathways was Mitochondrial dysfunction [−log(p) = 9.22]. Further, ESRRA was identified as a significant (p < 0.05) upstream regulator with an activation Z-score of 2.29, nearly identical to the protein data (Z = 2.21).

DISCUSSION
The results of this experiment provide compelling evidence that EC and IC rats have significant detectable differences in protein expression in the nucleus accumbens and that EC and IC rats show differential protein regulation after 14 days of cocaine self-administration. Our confidence in the novel EC/IC results is extremely high because the cocaine analysis provided more than a dozen results confirming what is already known about cocaine exposure. Further, in most cases the GSEA results support the IPA results even though the data sets comprising the IPA and GSEA functional databases are often completely orthogonal. The most clear and striking difference between EC and IC rats is in the expression of proteins involved in energy production. Increased energy demand might be expected given that the brains of EC rats are more active than IC rats, but the clear Enrichment X Cocaine interaction for these proteins suggests that EC rats respond to cocaine differently than IC rats. In looking at the heat maps it is clear that EC control rats have the highest amounts of energy-producing proteins and that cocaine self-administration decreases these proteins to levels similar to those of IC rats. IC rats, by contrast had lower levels already and no further decrease after cocaine. Energy metabolism proteins were also differentially regulated by psychological stress in EC and IC rats in our prior experiments (Fan et al., 2013a,b). Thus, it is becoming increasingly clear that EC and IC brains respond differently to perturbations that affect energy metabolism. The protein results are backed by RNA-seq data confirming upregulation of the energy metabolism pathway at the mRNA and protein level. Thus, protein and mRNA both suggest the energy metabolism pathway as a whole is upregulated in EC rats. Given that energy metabolism proteins are expressed at extremely high levels and have a rapid turnover (a result of the oxidative stress of the mitochondria), one might expect a better correlation between individual protein and mRNA expression. The lack of an exact match between protein and mRNA illustrates the main weakness of -omics studies: high Type II error (i.e., false negatives). Regardless, the analysis strategy employed in this study is designed to be somewhat resistant to both Type I and Type II error. The fact that the separate protein and mRNA analyses both identified Mitochondrial Dysfunction as a major factor attests to the resilience of the analysis. The fact that IPA and GSEA concur despite differential statistical strategies and functional annotation databases adds further confidence to the overall analysis strategy.
One recurring theme in this study as well as our upcoming transcriptomics study using mRNA from these same rats and our previous EC/IC stress proteomic study (Fan et al., 2013a,b) is proteins involved in neurodegeneration (Alzheimer's, Parkinson's, and Huntington's diseases). In the current study, MAPT was upregulated in EC rats, and MAPT, PSEN1, APP and HTT were among the most significant upstream regulators in multiple analyses. In trying to make sense of the results, it is interesting to note that mitochondrial dysfunction and energy metabolism are key players in the pathology of these neurodegenerative diseases (Quintanilla et al., 2013;Wang et al., 2013;Camilleri and Vassallo, 2014). Given that the mitochondrial proteins of EC rats are upregulated, one might hypothesize that environmental enrichment might alter the course of these neurodegenerative diseases. Published literature confirms that enrichment is protective in animal models of all three neurodegenerative diseases (Hu et al., 2010;Wood et al., 2010).
Much study has focused on the separate roles of VMAT2 and VGLUTs in addiction, but recent research suggests glutamate and dopamine are co-released in the NAcc of mice by mesolimbic dopaminergic neurons (Tecuapetla et al., 2010). Although it has not been shown that these neurotransmitters are packaged into the same vesicles, the current results showed concordant increases of VMAT2, VGLUT1, and VGLUT2 protein levels after cocaine. A literature search reveals that knockout of either VMAT2 or VGLUT2 can increase sensitivity to cocaine (Wang et al., 1997;Alsio et al., 2011). It has been shown previously that short-term cocaine exposure increases VMAT2 (Schwartz et al., 2007) yet long-term cocaine users have decreased VMAT2 (Little et al., 2003). Intriguing opportunities for regulation arise when one imagines fast ionotropic glutamate neurotransmission overlaid with slower metabotropic dopaminergic neuromodulation. Thus, it is possible that the coordinated increases in expression in the current study are relevant to the addictive nature of cocaine. Increased expression of these vesicular transporters might potentiate dopamine and glutamate release in the NAcc. Indeed, it has been shown that repeated cocaine administration potentiates glutamate release in addition to potentiating dopamine release (Reid and Berger, 1996). The authors propose that this pathway represents a fast mechanism for reward signaling.
For both main effects and the interaction, the largest effect size of regulated protein was cysteine sulfinic acid decarboxylase (CSAD). CSAD is the rate-limiting enzyme in the biosynthesis of taurine, an antioxidant and agonist of glycine receptors (Wang et al., 2005). Acute cocaine exposure has been shown to decrease taurine tissue concentration in the rat NAcc at 2 h post injection (Li et al., 2012). Additionally, repeated cocaine was shown to increase extracellular release of taurine at 30 min post injection (Yablonsky-Alter et al., 2009). The significance of these findings is unclear at this point, but the evidence at hand suggests taurine may play some role in addiction.
Taking a step back to get a broad view of the results, this study further supports the idea that the neurobiology of enriched rats differs significantly from that of isolated rats, particularly expression of proteins related to energy metabolism, RNA processing, membrane transporters and kinases. Further, EC rats show fundamentally distinct responses to cocaine compared to IC rats, particularly proteins related to RhoA signaling, energy metabolism, a variety of kinase cascades and ubiquitination. Interestingly, the EC-saline group was most distinct from the other three groups. Current and future experiments will begin to uncover the importance of these proteins in mediating this protective phenotype in an attempt to illuminate why many people are resistant to addiction despite exposure to drugs of abuse, thus identifying new targets for treatment of those who do become addicted. Future experiments teasing apart individual aspects of the enrichment condition (e.g., novelty, social contact, and exercise alone) will shed additional light on the exact nature of environmental enrichment.

AUTHOR CONTRIBUTIONS
Cheryl F. Lichti assisted in experimental design for proteomics experiments, performed proteomics data analysis and assisted in writing of manuscript. Xiuzhen Fan performed experiments and assisted in experimental design. Robert D. English assisted in proteomics experimental design, performed LC-MS analysis and assisted in writing materials and methods. Yafang Zhang and Dingge Li were involved with behavioral, surgical, dissection and design aspects of the project. Fanping Kong, Heidi Spratt, Mala Sinha and Bruce A. Luxon were involved with design and execution of the primary and secondary RNA-seq analyses. Thomas A. Green participated in all aspects of experimental design, performed GSEA analysis, generated relevant figures, and participated in writing of the manuscript. All authors read, provided feedback and approved the final version of the paper.