Evaluating the Stability of RNA-Seq Transcriptome Profiles and Drug-Induced Immune-Related Expression Changes in Whole Blood

Methods were developed to evaluate the stability of rat whole blood expression obtained from RNA sequencing (RNA-seq) and assess changes in whole blood transcriptome profiles in experiments replicated over time. Expression was measured in globin-depleted RNA extracted from the whole blood of Sprague-Dawley rats, given either saline (control) or neurotoxic doses of amphetamine (AMPH). The experiment was repeated four times (paired control and AMPH groups) over a 2-year span. The transcriptome of the control and AMPH-treated groups was evaluated on: 1) transcript levels for ribosomal protein subunits; 2) relative expression of immune-related genes; 3) stability of the control transcriptome over 2 years; and 4) stability of the effects of AMPH on immune-related genes over 2 years. All, except one, of the 70 genes that encode the 80s ribosome had levels that ranked in the top 5% of all mean expression levels. Deviations in sequencing performance led to significant changes in the ribosomal transcripts. The overall expression profile of immune-related genes and genes specific to monocytes, T-cells or B-cells were well represented and consistent within treatment groups. There were no differences between the levels of ribosomal transcripts in time-matched control and AMPH groups but significant differences in the expression of immune-related genes between control and AMPH groups. AMPH significantly increased expression of some genes related to monocytes but down-regulated those specific to T-cells. These changes were partially due to changes in the two types of leukocytes present in blood, which indicate an activation of the innate immune system by AMPH. Thus, the stability of RNA-seq whole blood transcriptome can be verified by assessing ribosomal protein subunits and immune-related gene expression. Such stability enables the pooling of samples from replicate experiments to carry out differential expression analysis with acceptable power.

Introduction time in the control group. In the current study, stability was assessed across four time points (ranging over a 2 year period) for a control group and a treatment group. Stability was evaluated primarily through comparisons of the expression profile of ribosomal subunit proteins and secondarily through comparisons profiling the entire transcriptome.

Animal Dosing and Sacrifice
This study was carried out in accordance with the declaration of Helsinki and the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the National Institutes of Health. The use of animal testing in this study was done under protocols E7295 and E7519 (issued to John Bowyer) that were approved by the NCTR institutional animal care and use committee (IACUC) which is fully accredited (Food and Drug Administration-National Center for Toxicological Research Accreditation #A4310-01) by NIH-OLAW. Nine-to ten-weekold ( 65 days) male Sprague-Dawley rats, ninety five total, were obtained from the Charles River Laboratories [Crl:CD(SD)]. Upon arrival at NCTR they received tail tattoos for identification. Prior to testing, rats were housed 2 per cage with food and water available ad libitum. Rats were housed on a daily 12 hr light cycle with lights on at 6:00 am and off at 6:00 pm. During housing, the temperature (23°C) and humidity (53%) were controlled. The rats were tested between 12 and 13 weeks of age ( 87 days old).
The data used in this study (a cohort of a larger study contained in GSE62368 and GSE64778) is comprised of data from n = 22 saline dosed "control" rats and n = 21 rats given a neurotoxic exposure to AMPH. The larger study was carried out to identify mRNA biomarkers and to determine transcriptional immune-related responses present in circulating blood after neurotoxic exposures to hyperthermia, neurotoxic AMPH exposure, and non-neurotoxic AMPH exposure. Herein, we focus solely on the control and neurotoxic AMPH treatment groups. Multiple replicates of the same experiment were conducted over a 2 year period. This was done to ensure that any significant effects seen in transcript expression would be robust and reproducible over time and hopefully subsequently observable in other laboratories. The sample size of the control group at each time period was: Con1: n = 3, Con2: n = 7, Con3: n = 6, and Con4: n = 6. The size of the AMPH groups were: AMPH1: n = 5, AMPH2: n = 3, AMPH3: n = 6, and AMPH4: n = 7. The saline control and AMPH groups were sacrificed on: Con1 and AMPH1, September 2012; Con2 and AMPH2, June 2013; Con3 and AMPH3, Aug 2013; Con4 and AMPH4, July 2014.
Dosing commenced at 8:00 a.m. and ended at 2:00 p.m. AMPH-exposed animals were given 4 doses of amphetamine comprised of a sequential exposure to 5, 7.5, 10, and 10 mg/kg AMPH subcutaneously with 2 hr between each dose. The d-amphetamine-sulfate (Sigma-Aldrich, St. Louis, MO) dose was dissolved in normal saline (1 ml/ kg injected). This AMPH dosing paradigm was modified slightly from an original acute methamphetamine dosing paradigm that was shown to produce neurotoxicity [20,21]. The saline animals received 4 injections of 1 ml/kg s.c. In all groups, the behavior and body temperature of each rat were monitored hourly during saline or AMPH exposure and until at least 3 hr after the last dose (time of sacrifice). The lethal effects of hyperthermia and hyperpyrexia in the AMPH groups when body temperatures exceeded 41.4°C were prevented by placing the animals unrestrained on crushed ice for 15 to 30 min in a clean, wood chip free cage to allow their temperatures to drop below 40.0°C. withdrawn from the heart using an 18 gauge needle attached to a 5 ml syringe, and the rats were then euthanatized with decapitation. Approximately 1 ml of the blood was immediately separated into approximately three 300 μl aliquots, frozen on dry ice in cylindrical foil capsules, and then stored at -70°C for later RNA-seq analysis. The remainder of the blood was injected into a 7.0 ml BD (Franklin Lanes, NJ) Vacutainer containing 12 mg EDTA (data from these aliquots are to be reported at later date). One ml was used to count the total RBC and WBC numbers as well as the different types of WBC present. Within 12 h of collection of the cardiac blood in EDTA coated vacutainers, one ml was used to determine total RBC, WBC counts, and the numbers of lymphocytes (T-and B-cells), monocytes, neutrophils, basophils and eosinophils present. Complete blood counts were determined on an ABX Pentra 60 C+ analyzer (ABX, Irvine CA). Maintenance and calibration was done according to the manufacturer's recommendations. Three levels of assayed controls were included in daily analyses as internal controls.
RNA isolation and RNA-seq processing RNA isolation was performed using RNAzol BD (Molecular Research Center, Inc., www. mrcgene.com) [9] with modified procedures for frozen blood. Just prior to RNA isolation, the frozen 300-400 μl RBC or whole blood aliquots were placed on dry ice and cut into approximately 100 to 150 μl sections. Two or three of these smaller sections were immediately (were not allowed to thaw) homogenized in 1.6 ml of RNAzol BD containing 45μl of 5 N acetic acid using a 2 ml Teflon pestle homogenizer. Approximately 300 to 400 μl of frozen whole blood was used for RNA isolation. Separation of the upper phase containing the extracted RNA with lower phases containing protein and DNA was accomplished by adding 350 μl of 4-bromoanisole and centrifugations at 12,000 x g for 10 min. To precipitate the RNA, 600 μl 2-propanol was added to 600 to 700 μl the extracted upper phase. After 15 min, centrifugation at 12,000 x g was performed. The RNA ppt was washed once with 70% ethanol and water and a second wash of 100% ethanol. After 1 min air drying the RNA was dissolved in 25 to 35 μl of RNase-free water. The final total cellular RNA recovered ranged from 5 to 25 μg and was stored at -70°C until shipment for data processing. Two to four μl of each aliquot was set aside for purity analysis at NCTR using an Agilent 2200 TapeStation System (Agilent Technologies, Palo Alto, CA).
The isolated total RNA was then shipped overnight on dry ice to Expression Analysis Inc. [EA; Durham, NC] for globin transcript removal and sequencing. Alpha and Beta Globin mRNA were substantially depleted from total RNA samples using the GlobinClear-Mouse/Rat Kit (Life Technologies # AM1981), essentially as described by the vendor. Briefly, 1.25 μg of total RNA isolated from whole blood was combined with biotinylated capture oligonucleotides complementary to globin mRNAs and the mixture incubated at 50°C for 15 minutes to allow duplex formation. Streptavidin magnetic beads were added to each specimen, and the resulting mixture was incubated for an additional 30 minutes at 50°C to allow binding of the biotin moieties by Streptavidin. These complexes, comprising Streptavidin magnetic beads bound to biotinylated oligonucleotides that are specifically hybridized to globin mRNAs, were then captured using a magnet. The globin-depleted supernatant is transferred to a new container and further purified using RNA binding beads. The final globin mRNA-depleted RNA samples are quantitated by spectrophotometry using a NanoDrop ND-8000 spectrophotometer. synthesis methods, as implemented via Illumina technology, were used to generate the RNAseq data. Globin mRNA-depleted RNA samples were converted into cDNA libraries using the TruSeq Stranded mRNA Sample Prep Kit (Illumina, #RS-122-2103). Starting with 300 ng of globin mRNA-depleted RNA, polyadenylated RNA (primarily mRNA) was selected and purified using oligo-dT conjugated magnetic beads. This mRNA was chemically fragmented and converted into single-stranded cDNA using reverse transcriptase and random hexamer primers, with the addition of Actinomycin D to suppress DNA-dependent synthesis of the second strand. Double-stranded cDNA was created by removing the RNA template and synthesizing the second strand in the presence of dUTP in place of dTTP. A single A base was added to the 3' end to facilitate ligation of sequencing adapters, which contain a single T base overhang. Adapter-ligated cDNA was amplified by polymerase chain reaction to increase the amount of sequence-ready library. During this amplification the polymerase stalls when it encounters a U base, rendering the second strand a poor template. Accordingly, amplified material used the first strand as a template, thereby preserving the strand information. Final cDNA libraries were analyzed for size distribution using an Agilent 2200 TapeStation (D1000 Screentape, Agilent # 5067-5582), quantitated by qPCR (KAPA Library Quant Kit, KAPA Biosystems # KK4824), then normalized to 2 nM in preparation for sequencing.
The 2nM normalized samples were pooled. cDNA templates were denatured using fresh 0.1N NaOH, and then diluted to a final loading concentration of 13pM. Cluster generation was performed on an Illumina cBOT (v1.5.12.0) using an Illumina TruSeq Paired-End Cluster Kit v3 (Illumina # PE-401-3001) to cluster an Illumina Paired End Flow Cell. Templates were attached to the flowcell via a dense lawn of oligonucleotides that bind to the sequencing adapters added during sample preparation, which are extended and then denatured. The flowcell was then sequenced through 51 bases, paired end, with an 8 base index cycle on an Illumina HiSeq 2000 (HiSeq Control Software v 1.5.15.1). During sequencing cycles, fluorescent reversible terminator dNTPs were added to the clusters, with only a single base per target being incorporated. Following imaging of the clusters, the terminator and fluorescent tag were cleaved so that the next base could be incorporated. Quality control files contained read length and depth results (before and after clipping), presence of artifact/ duplicate sequences, distribution of base quality and base frequency by sample. Also, flow cell total yield, PF reads, barcode quality, alignment summaries and number of genes detected were determined. The basecall files were converted to fastq files using CASAVA 1.8.2. The fastq files were clipped using fastqmcf with the parameters "-max-ns 4-qual-mean 25-H-p 5-q 7-l 25" [22]. The fastq files were aligned to the rat Ensembl release 70 transcriptome (rn5) using bowtie 0.12.9 with the parameters "-e 500 -m 100 -chunkmbs 256". The alignments were quantified using RSEM v 1.2.0 with no special parameters.

RNA-seq Data Analysis
RNA-seq by Expectation-Maximization (RSEM) [23] data generated by EA was rounded to produce digital expression values. In total, 16,881 transcripts were mapped to the reference rat genome in all samples. Transcript-level outliers within each experimental group were identified by the iLOO method [24] and subsequently replaced with the trimmed mean. All reports of average transcript expression are summarized using DESeq2 normalized counts [25]. DESeq2 was also used to carry out differential expression analysis. Significance was assessed using pvalues adjusted by the Benjamini-Hochberg method for multiple comparison testing. Transcripts with adjusted p-values lower than 0.05 and absolute log 2 (FC) values exceeding specific cutoffs were declared significant. All reported values of correlation reflect Pearson's correlation coefficient and are significant at α = 0.05.
A number of criteria were used to identify the most stable transcripts for RNA-seq data. Ideally, a stable transcript would show minimal variability, be consistently expressed in the high region of expression, and demonstrate consistent expression across treatment groups. Thus, candidates for stable expression were transcripts with no detected outliers, no measurements in the low expressed region [26], and no significant difference between the control and treatment group.
All analyses were carried out in R (http://www.r-project.org).

Results
The RNA integrity number (RIN) for the control whole blood samples (n = 22) ranged from 7.0 to 9.4, averaging 8.3 ± 0.15 (mean ± SEM). The RIN for the AMPH group ranged from 6.9 to 9.3, averaging 8.1 ± 0.15 (mean ± SEM), and did not significantly differ from control. Although globin depletion is reported to reduce expression [2], hemoglobin alpha (Hba) remained the highest expressed transcript. Also, as would be expected, the expression of olfactory receptor genes (approximately 1,200 unique transcripts [10]) was virtually absent in blood with only 14 genes having counts greater than 1. In addition, we assessed whether high levels of myoglobin transcript (Mb) were present in the blood samples. Since cardiac muscle contains very high levels of myoglobin [27], high expression of the myoglobin transcript could indicate that cardiac muscle was "contaminating" the samples as a result of the cardiac punctures to obtain whole blood. Two of the control animals from the Con4 group had very high levels (>5000 counts) for Mb. They were removed from further analysis resulting in n = 20 for the control group. None of the animals in the AMPH group had to be removed from the study due to high Mb levels. Mb expression for the remaining control and AMPH animals was relatively low with mean expressions of 29.8 ± 3.4 and 35.8 ± 1.8, respectively. Average expression for the 70+ genes necessary to encode the 80s ribosome are shown in Table 1. The large 70s subunit transcripts are labeled with an NCBI gene symbol beginning with "Rpl," while the small 30s subunit transcripts begin with "Rps." Average expression for most of the transcripts was ranked in the top 5 th percentile of expression levels (i.e. rank of 844 or higher). All of the ribosomal subunit genes, except Rpl22, had at least one transcript (splice variant) in the top 5 th percentile; Rpl22 ranked 1010 in whole blood. RNA-seq expression of the ribosomal proteins for control animals were highly correlated (r0.98 for log2 expression data) across the first 3 time points (Fig 1).
Initial sequencing of the Con4 group (denoted by Con4A) yielded poorer correlation with the other groups (r0.74 for log2 expression data) (Fig 1). To some extent, diminished correlation between Con4A and Con1 -Con3 was also evident when assessing profiles of the entire transcriptome (Fig 2). Interestingly, this was the case despite the fact that total clusters, genes detected, % genes detected, % transcriptome mapped, skewness, median coverage, and % of all genes mapped did not differ between Con1, Con2, Con3, and Con4A. One noticeable fluctuation was the slight decrease in %(G+C), which ranged from 43% to 48% for Con4A, compared to the usual 48% to 53% range for the other groups. The less than optimal sequencing results shown with Con4A were due to an Illumina library preparation kit that contained a less than optimal PCR Master Mix reagent.
Subsequently, aliquots of the same RNA from the 4 animals in the Con4A group were resequenced (denoted by Con4B), yielding a more accurate expression of the ribosomal proteins (Fig 1) and comparable %(G+C) values. Differential expression (DE) analysis of pairwise comparisons of the four control groups indicated that error rates (i.e. DE transcripts) were well below the nominal level. For example, DE analysis comparing Con1 and Con4B, the most distant pairwise time points, returned only 3 significant transcripts. Thus, although the groups were dosed and sacrificed two years apart, stability of transcript expression in control animals was maintained. When profiles of the ribosomal subunit expression for the four AMPH groups (AMPH1, AMPH2, AMPH3 and AMPH4) were compared to each other, pairwise correlation values were comparable to those observed for the control groups (correlation values among the four pairs ranged between 0.97 and 0.98). Somewhat surprisingly, correlation between each control group and the time-matched AMPH group was at least 0.97 (Fig 3).
To further characterize repeatability of the control expression profiles, 29 genes specifically related to the immune system and found predominately on leukocytes were evaluated. The mean (and standard error) of each gene is reported in Tables 2 and 3 for both groups (all four subgroups were included in the calculations). Most of the selected genes are classified as cluster of differentiation antigens (i.e. classification determinant) and are located on the outer membranes of leukocytes. The seven transcripts primarily expressed on monocytes/macrophages (Ccr1, Cd14, Cd163, Cd68, Itgax, Nos2 and Arg1) had a 1000-fold range if indeed interleukin 1 β is exclusively produced in monocytes (Il1b counts reached 18,000). Ten genes preferentially expressed in either T-cells (Cd2, Cd3e, Cd8a, Cd28, Cd247) or B-cells (Cd19, Cd22, Cd79b, Cd180 and Ly86) had a median expression of 1921 with only slightly more than a 10-fold range. The median expression of 6 selected transcripts found in almost all leukocytes (Cd300a, Cd44, Cd48, Cd84, Cd97 and Sell) was 6342 but ranged over 70-fold. The neutrophil related genes (Camp, Elane, Mpo, Mme, Prtn3 and Ctsg) had expression levels below 100 except for Camp.
The expression of these immune-related genes had near-1 correlation (r0.97) among the four control groups (Con1-Con3 and Con4B) (Fig 4). However, as expected, the correlation between each control group and its paired time-matched AMPH group was not as strong because of the known effect of the drug on the immune system (Fig 5). Fig 6 presents the log 2 fold change of AMPH to control expression for the selected immune-related genes, except for the three neutrophil-related genes with very low expression. Relative expression is presented for all four time-matched groups and also for all of the data (denoted by a purple diamond). AMPH has a significant effect on 15 of the 26 genes. Relative to control, all of the genes primarily expressed in T-cells decreased approximately 1.5-to 3.0-fold, while 4 of the 7 monocytespecific transcripts increased between 2.0 and 6-fold. Expression of 3 of the 5 genes relatively specific for B-cells decreased around 1.5-fold, and genes that are normally expressed in most leukocytes were not significantly different (less 1.5-fold), with only two of the six genes (Cd300a and Cd84) reaching significance between AMPH and control. The levels of the major types of leukocytes in blood for control and AMPH are presented in Table 4. The levels of the lymphocytes (72.3%), neutrophils (23.6%), monocytes (3.5%), eosinophils (0.4%) and basophils (0.2%) were as would be expected based on standard reported levels Monocytes were significantly upregulated by AMPH; the levels of neutrophils were also upregulated by AMPH and marginally significant. Lymphocytes were significantly down-regulated by AMPH.

Discussion
The results of our study demonstrate the stability of whole blood RNA-seq transcriptome profiles after globin transcript depletion. Globin depletion was implemented since it has been observed to improve the transcriptome profile in whole blood in previous studies using RNAseq methods [2]. Any loss of RNA due to globin removal was reasonably consistent since comparisons of transcriptome profiles over a period of two years were highly correlated for both the control and treatment group. The level of agreement is also attributed to consistent expression of each experimental group across time. Our study provides evidence that the consistency and repeatability of RNA-seq of whole blood can be validated by profiling expression of the ribosomal subunits. However, further studies will be necessary to determine whether this approach is more sensitive than methods such as the detection of aberrant G-C profiles [28]. The utility of RNA-seq data from whole blood was further demonstrated by the completeness  Circulating Whole Blood Transcriptome Stability and Drug Response of the expression profile of immune-related genes as well as the consistent effect that drug treatment had on the expression of these genes, at least in the case of a toxic exposure to AMPH [29].
The 80s ribosome is a group of more than 70 transcripts, of which are all expressed at relatively very high levels [14]. Note that there are 10 transcripts with lower expression; however, the gene for each of these transcripts produced another transcript (splice variant) that measured at much higher levels. As seen in this study, the one anomalous group of RNA-seq data (Con4A) that was less than optimal showed large variability in these transcript levels (Fig 1). This was solely due to a very rare instance of amplification kit variability. Although the observed variability in Con4A expression was somewhat apparent in the entire transcriptome (Fig 2), the deviant variation was readily apparent when profiling the ribosomal subunits. As an example, Pearson correlation values for Con4A were at least 0.90 for the whole genome but were as low as 0.74 when comparing expression levels of the ribosomal transcripts. An additional benefit of using ribosomal transcripts rather than the entire genome is that the majority of ribosomal transcripts do not change much with AMPH treatment. Only two of the ribosomal transcripts were significantly differentially expressed between control and AMPH (p<0.05 and FC>2). However, it is not yet known how other types of drug/toxicant treatments will affect subunit expression relative to control. Several ribosomal proteins have been used as "housekeeping" genes to normalize RT-PCR expression across individual RNA samples, although they may not be optimal in all tissues [11][12][13][14]. As an example, in our study, ribosomal protein L7, Rpl7, was identified as the highest consistently expressed transcript across all 41 samples. Finally, a relatively high expression for at least one transcript of each gene in the ribosomal subunit is almost as important as consistent expression. This is because expression of all of the genes is necessary to form one of the most numerous protein complexes in an organism. Our data showed that at least one transcript for each ribosomal gene ranked in the upper 5 th percentile of expression. One could argue that the use of a set of 10 or so genes commonly used to normalize RT-PCR data would provide a better way to determine the potential for significant variation in whole blood RNA-seq. However, such analysis would be limited in some respects because the specific genes in a set of housekeeping genes for RT-PCR are not linked by a common process or protein complex. Thus, unlike the set of ribosomal transcripts, Table 3. Citations for Gene Expression within Leukocyte Cell Types.
In this study, we profiled the expression of immune-related genes (mostly cell determinant genes that are expressed on leukocytes). The stability and consistency of the expression of these immune-related genes are also important as these are the genes in whole blood, along with cytokines and chemokines protein levels in blood, that are used to evaluate how the immune system is responding to a given disease, treatment or toxic insult [30][31][32][33][34][35][36][37]. These genes are major effectors in the chemokine and cytokine signaling system of leukocytes. There have been previous reports on how gene expression, as determined by microarray or RNA-seq, is affected in human blood in anorexic patients or after acute ethanol exposure [38,39]. However, these studies focused primarily on genes not necessary directly related to cytokines, chemokines, or receptors that are exclusively or predominately localized to leukocytes. One very recent RNA-  The pooled expression of AMPH to control over all subgroups (AMPH/Con) is represented by a purple diamond. Differential expression was assessed for the latter comparison. Log2FC expression for genes with adjusted p-values smaller than 0.05 and 0.58 < |log2FC| < 1 is marked by one asterisk. Likewise, expression for genes with adjusted p-values smaller than 0.05 and |log2FC| 1 is marked with two asterisks.  seq study investigated the effects of immunosuppression on the expression of genes related to cytokine/chemokine signaling in mononuclear cells isolated from peripheral blood in patients involved in [40]. A pronounced 10-fold or more down regulation of these genes was observed due to the immunosuppression. Our data clearly indicates that RNA-seq methods are also a sensitive method for detecting changes in gene expression related to cytokine and chemokine signaling resulting from toxic exposures, in this case AMPH. We observed that the levels of transcripts relatively specific for monocytes, T-cells, B-cells, neutrophils or present in all leukocytes were roughly as expected in the circulating blood of rat. Transcript levels relatively specific for monocytes, the least abundant cell type (3%), were readily detectable in the control and amphetamine groups, and some such as Ccr1 had very high expression. All but one of the genes specific for neutrophils, which are 8 to 15% of the all WBCs in rat, were at low levels. This would be expected due to the nature of mature neutrophils present in circulating blood. However, some genes, mostly those also present in other types of WBCs, can be induced [19]. Therefore, some of the expression that we observed in blood of interleukin 1b (Il1b) and neutrophil cytosolic factor 1 (Ncf1), which were at very high levels (data not shown), could have been due to expression in neutrophils as well as monocytes. Thus, the overall profile of immune expression was well represented.
There was a significant difference between the control and amphetamine groups for several of the selected immune-related genes. Specifically, 10 of the 26 immune-related genes had 2-fold or greater change in expression and 5 additional genes demonstrated at least 1.5-fold change in expression. Two genes, Cd14 and Arg1, increased expression more than 4-fold as a result of amphetamine exposure. Importantly, the response (change in transcript levels) of all these genes, except Cd163, was very consistent across all four treatment groups. In general, many of the changes observed in these immune related genes might be expected to be due to AMPH exposure, including the increase in β-adrenergic stimulation produced by AMPH releasing norepinephrine [41][42][43][44]. Finding treatment-induced changes in immune-related transcripts is of the utmost importance since previous, and likely subsequent, animal studies involving gene expression in blood or serum have focused on changes related to immune system function.
Our data on cell counts of monocytes and lymphocytes (T-cells and B-cells) in circulating blood indicate that the transcript expression changes may be partially explained by increased or decreased numbers of specific types of leukocytes present in circulating blood after AMPH. The monocyte count increased slightly over two fold while the lymphocyte levels decreased approximately 30%. In regards to the monocytes, the increase in their numbers could partially explain the increase in expression of genes specifically found in them. Nonetheless, the changes in numbers/counts cannot explain all the gene expression changes seen. The expression increases in monocyte-specific genes range from a less than 1.5-fold increase (Ccr1 and Itgax) to almost a 6-fold increase for Cd14, which would strongly argue that there is likely a relative change in the relative expression of the genes present within the circulating monocytes. Likewise, the 30% decreases in the number of circulating T-cells and B-cells may not explain the much greater fold change decreases, particularly in T-cell-related genes. Overall, the increases in the circulating blood expression profile and monocyte and neutrophil cells numbers indicate an activation of the innate immune system by AMPH [45][46][47][48].
The stability of transcript expression in four control groups of animals (Sprague-Dawley rats in this study) over a two-year timeframe was strong. However, more importantly, the changes in gene expression of the animals given a neurotoxic dosage of AMPH relative to the controls were very consistent across all four time points. Thus, the work presented here demonstrates that under these conditions control and treatment-specific animals can be pooled for analysis. The increase in sample size will subsequently increase power to detect differentially expressed transcripts. The ability to show reproducibility of treatment effect on the transcriptome over time lends to the validation of the expression changes that are reported.

Author Contributions
Conceived and designed the experiments: RPS JFB NIG NMC JPH KMT. Performed the experiments: JFB KMT NMC. Analyzed the data: JFB NIG RPS NMC JPH KMT. Contributed reagents/materials/analysis tools: JFB KMT NIG NMC. Wrote the paper: JFB NIG RPS JPH NMC.