Research | Open | Published:
Epigenetic impacts of stress priming of the neuroinflammatory response to sarin surrogate in mice: a model of Gulf War illness
Journal of Neuroinflammationvolume 15, Article number: 86 (2018)
Gulf War illness (GWI) is an archetypal, medically unexplained, chronic condition characterised by persistent sickness behaviour and neuroimmune and neuroinflammatory components. An estimated 25–32% of the over 900,000 veterans of the 1991 Gulf War fulfil the requirements of a GWI diagnosis. It has been hypothesised that the high physical and psychological stress of combat may have increased vulnerability to irreversible acetylcholinesterase (AChE) inhibitors leading to a priming of the neuroimmune system. A number of studies have linked high levels of psychophysiological stress and toxicant exposures to epigenetic modifications that regulate gene expression. Recent research in a mouse model of GWI has shown that pre-exposure with the stress hormone corticosterone (CORT) causes an increase in expression of specific chemokines and cytokines in response to diisopropyl fluorophosphate (DFP), a sarin surrogate and irreversible AChE inhibitor.
C57BL/6J mice were exposed to CORT for 4 days, and exposed to DFP on day 5, before sacrifice 6 h later. The transcriptome was examined using RNA-seq, and the epigenome was examined using reduced representation bisulfite sequencing and H3K27ac ChIP-seq.
We show transcriptional, histone modification (H3K27ac) and DNA methylation changes in genes related to the immune and neuronal system, potentially relevant to neuroinflammatory and cognitive symptoms of GWI. Further evidence suggests altered proportions of myelinating oligodendrocytes in the frontal cortex, perhaps connected to white matter deficits seen in GWI sufferers.
Our findings may reflect the early changes which occurred in GWI veterans, and we observe alterations in several pathways altered in GWI sufferers. These close links to changes seen in veterans with GWI indicates that this model reflects the environmental exposures related to GWI and may provide a model for biomarker development and testing future treatments.
A coalition of 34 countries deployed approximately 956,600 troops during the 1990–1991 Gulf War , with the majority, ~ 700,000, from the USA. An estimated 25–32% of these veterans fulfil the requirements of a Gulf War illness (GWI) diagnosis . GWI is an archetypal, medically unexplained, chronic condition characterised by persistent sickness behaviour, with neuroimmune and neuroinflammatory components.
Symptoms of GWI include fatigue, musculoskeletal pain, cognitive dysfunction, chemical sensitivities, loss of memory and sleep disruption, which can be characterised as ‘sickness behaviour’ [3, 4]. ‘Sickness behaviour’ is normally a result of inflammatory response to illness or injury, which usually resolves itself over time after the initial insult is removed. Symptoms were reported within 6 months of the conflict [5,6,7].
Although the exact cause of GWI is still unknown, there is a consensus that exposure to environmental toxins is the most likely cause . GWI symptoms are highly heterogeneous, and specific symptoms may be related to specific experiences: for example, Gulf War veterans exposed to nerve agents or oil well fires are at increased risk of brain cancer compared to other Gulf War veterans .
A leading hypothesis for the cause of GWI is that the high physical and psychological stress of combat interacted with exposure to acetylcholinesterase (AChE) inhibitors [1, 4, 10,11,12,13,14,15,16,17,18,19]. Military personnel were exposed to a number of AChE inhibitors [1, 20], including pyridostigmine bromide (PB), a reversible AChE inhibitor, as a prophylactic against nerve agents; sarin, soman, and related nerve agents, irreversible AChE inhibitors, which combatants were inadvertently exposed to after demolition of Iraqi supply depots, such as at Khamisiyah; organophosphate pesticides, irreversible AChE inhibitors, which were widely used to prevent pest-borne diseases and irritation ; permethrin, an insecticide which may inhibit AChE [10, 21]; and DEET, an insect repellent and a weak AChE inhibitor which may enhance the activity of other AChE inhibitors . For example, an estimated 95,000 deployed personnel were exposed to the plume from the Khamisiyah demolition, and approximately 250,000 may have been exposed to low levels of nerve agents during aerial bombardments earlier in the conflict . Further, the number of nerve agent alarms heard is correlated with risk for GWI . Accumulating research has indicated that deleterious health effects of exposures to psychophysiological stress [25, 26] and environmental toxicants [27, 28] involve epigenetic modifications that affect transcriptional regulation.
The overall objective of this study was to examine genome-wide epigenetic transcriptional modifications in the brain using an established mouse model of GWI [4, 12, 15, 16]. Our previous research demonstrates that effects on neuroinflammatory pathways occur shortly after initial exposures. For example, pre-exposure with the stress hormone corticosterone (CORT) causes an increase in expression of specific chemokines and cytokines in response to diisopropyl fluorophosphate (DFP) , an irreversible acetylcholinesterase inhibitor  used here as a sarin surrogate. This corresponds well with the work in GWI study participants [30,31,32,33,34,35,36,37], which have shown immunological abnormalities, and a recent paper , which has shown specific immune-related biomarkers for GWI veterans. We hypothesized that epigenetic and transcriptomic changes upon initial exposures would identify gene pathways linked to poor health outcomes in GWI.
Adult male C57Bl/6J mice were purchased from Jackson Laboratory (Bar Harbor, ME, USA). A total of 79 animals were used for the analyses presented here. All procedures were performed under protocols approved by the Institutional Animal Care and Use Committee of the Centers for Disease Control and Prevention, National Institute for Occupational Safety and Health and the US Army Medical Research and Materiel Command Animal Care and Use Review Office. The animal facility was certified by AAALAC International. Upon receipt, the mice were housed individually in a temperature (21 ± 1 °C) and humidity-controlled (50 ± 10%) colony room maintained under filtered positive-pressure ventilation on a 12-h light/12-h dark cycle beginning at 06:00 h. The plastic tub cages were 46 × 25 × 15 cm; cage bedding consisted of heat-treated pine shavings spread at a depth of approximately 4 cm. Teklad 7913 irradiated NIH-31 modified 6% rodent chow, and water were available ad libitum.
The dosing paradigm is presented in Fig. 1. Mice were given CORT in the drinking water (200 mg/L in 0.6% EtOH) for 4 days. This CORT regimen is known to be immunosuppressive as evidenced by decreased thymus weight ; thymus and spleen weights were confirmed to be decreased (> 20%) in similarly exposed animals [4, 15]. On day 5, mice were given a single intraperitoneal injection of either DFP (4 mg/kg) or saline (0.9%).
Thus, there were four exposure groups: (1) saline: vehicle for 4 days, then saline injection on day 5; (2) CORT: CORT for 4 days with a saline injection on day 5;( 3) DFP: vehicle for 4 days with DFP injection on day 5; and (4) CORT + DFP: CORT for 4 days with a DFP injection on day 5.
Brain dissection and tissue preparation
Mice were killed by decapitation and the brains rapidly removed. The frontal cortex, consisting of the anterior portion of the cortex , and total hippocampus were dissected free-hand on a thermoelectric cold platform (Model TCP-2; Aldrich Chemical Co., Milwaukee, WI, USA) and immediately frozen at − 80 °C.
Differential gene expression
Frontal cortex mRNA-seq data was generated on the Illumina HiSeq 2000 by Q Squared Solutions Expression Analysis LLC (Morrisville, NC, USA), paired-end, with a read length of 100 bp (four samples of saline, CORT and DFP, five samples of CORT + DFP). Hippocampus mRNA-seq data was generated by Sickkids (Toronto, Ontario). Sequencing was carried out on the HiSeq 2500, paired-end, with a read length of 125 bp (n = 4 for all groups).
Fastq files were trimmed to remove adaptors and low-quality reads (q < 30) using TrimGalore version 0.4.1 , around Cutadapt  (version 1.9.1). The pre-processed reads were examined using FastQC .
Trimmed files were aligned with the STAR aligner  (version 2.5.2a) in a two-pass mode. The GENCODE GRCm38.p4 assembly (mm10) and annotations were obtained from the GENCODE website [44, 45] and used throughout. For the frontal cortex, there was an average of 33,547,161 reads per sample and 98.2% average mapped reads. For the hippocampus, there was an average of 35,429,647 reads per sample and 98.2% average mapped reads.
The resultant Bam files were analysed using the GenomicAlignments  and DEseq2  (version 1.10.1) R packages, using the DEseq2 standard pipeline. A recent comparison study has identified this as an appropriate tool to use with replicates and when relatively large biological effects are expected . Briefly, DESeq2 first fits a generalized linear model for each gene, with read counts modelled as a negative binomial distribution. An empirical Bayes approach is used for shrinkage of dispersion estimation, and the Wald test is used for significance testing, which is then adjusted for multiple corrections using the Benjamini and Hochberg method .
Samples were checked for similarity using Poisson dissimilarity matrix  with the R package PoiClaClu 1.0.2 and visualised with pheatmap 1.0.8.
The R Bioconductor package DeconRNASeq  was used to estimate the proportion of different cell types within the sample from the RNA-seq data. Data enriched for specific CNS cell types were downloaded from the Gene Expression Omnibus (GEO) [52, 53], Series GSE52564, which contains data from the Mus musculus cerebral cortex , to use as a reference of cell-type-specific gene expression. This RNA-seq data was trimmed and aligned and gene expression quantified as above. An expression signature for each of six cell types (astrocytes, neurons, oligodendrocyte precursor cells (OPC), myelinating oligodendrocytes (MO), microglia and endothelial cells) was obtained by finding those genes with a five-fold difference in expression in one cell type, compared to each of the others.
Using these expression signatures, the proportion of astrocytes, neurons, oligodendrocyte precursor cells (OPC), myelinating oligodendrocytes (MO), microglia and endothelial cells were estimated for each of the 17 cortex samples and 16 hippocampus samples using DeconRNASeq . DeconRNASeq is based on a linear model of a sum of pure tissue or cell-type-specific reads of all cell types, weighted by the respective cell-type proportions. To estimate the proportions of known tissue types in a sample, DeconRNASeq solves a non-negative least-squares constraint problem with quadratic programming to obtain the globally optimal solution for estimated fractions. It is accurate down to cell types making up only ~ 2% of the total cell populations .
DNA methylation modifications
DNA was extracted using the E.Z.N.A Tissue DNA Kit (VWR-Omega Bio-Tek). Bisulfite conversion was carried out using the Qiagen Epitect Fast Bisulfite Conversion Kit, and library preparation was performed using the Ovation NuGen RRBS Kit. Reduced representation bisulfite sequencing (RRBS) was carried out by the Princess Margaret Genomics Centre, part of the University Health Network, Toronto, on a NextSeq 500, using a single end, 70 base read length and multiplexed at 9–10 samples per flowcell. Samples for the cortex were 6 saline, 6 CORT, 6 DFP and 8 CORT + DFP, and for the hippocampus, 4 saline, 5 CORT, 3 DFP and 8 CORT + DFP.
RRBS fastq files were trimmed to remove adaptors and low-quality reads (q < 30) using TrimGalore version 0.4.1  around Cutadapt  (version 1.9.1). Trimmed files were then aligned to the GENCODE GRCm38.p4 (mm10) assembly, using Bismark (v0.16.0)  wrapped around bowtie2 (version 2.2.6) . For the frontal cortex, the average reads per sample was 34,962,881 with 57% mapping efficiency, and for the hippocampus, the average reads per sample was 37,512,807, 63.6% mapping efficiency.
The resultant bam files were analysed with MethPipe (3.4.2) [57, 58], using the suggested methods. The bisulfite conversion rate was > 98.9 for all samples. The methylation level for every cytosine site at single-base resolution is estimated as a probability based on the ratio of methylated to total reads mapped to that loci. Differential methylation was calculated by beta-binomial regression, with all batches and exposures included as factors, and CORT + DFP exposure set as the test factor. We examined both differentially methylated cytosines (DMCs) and differentially methylated regions (DMRs).
H3K27ac is a mark of active enhancers, strongly suggesting that genes with differential enrichment of H3K27ac will be differentially expressed [59, 60]. Native chromatin immunoprecipitation sequencing (ChIP-seq) using an MNase digestion and alignment to the mm10 genome using Burrows–Wheeler Alignment was carried out by the Genome Sciences Centre, BC Cancer Agency. Samples were sequenced single end, 75 base read length on a Hiseq 2500 platform. The average read per sample was 131,727,189 with 98.9% average mapped reads. There were four samples per group, with each sample having immunoprecipitated and input DNA sequenced.
PePr (Python) and diffReps (Perl) packages were used for ChIP-seq analysis, as a recent paper by Steinhauser et al.  suggested that both are good tools when biological replicates are available. The results from each were compared and analysed to provide a conservative list of sites showing differential enrichment.
PePr (1.1.14)  was run according to the authors’ suggested pipeline. PePr used a sliding window approach, shifting all reads toward their 3′ direction by half of the empirically estimated DNA fragment length and estimating window width based on the average width of the top pre-candidate peaks. The genome was then divided into consecutive widths that overlap by 50%, and the number of reads within each window was counted. This read count was then normalised based on total read count among ChIP and control sample and the relative average peak heights among ChIP samples. Read counts were modeled across replicates and between groups with a local negative binomial model. Genomic regions with less variable read counts across replicates were ranked more favourably than regions with greater variability, thus prioritizing consistently enriched regions . Narrow peaks were assumed, in line with previous literature on H3K27ac (e.g. ).
diffReps (1.55.6)  was run according to the authors’ suggested pipeline. Bam files first had to be converted to bed files, using Bedtools (v2.26.0) . Unlike PePr, diffReps used a set window size of 1000 bp for narrow peaks with a step size of 100 bp. The genome was pre-screened to remove regions with low read count to improve power and decrease computational time. Normalisation was carried out using the read count for a particular window over read count across all samples. An exact negative binomial test was used for differential analysis, which used biological replicates. p values were adjusted by the Benjamini-Hochberg method . Peaks were annotated to genes using region_analysis .
Identifying genes unique to the CORT + DFP exposure
For several computational tools, e.g. DESeq2, PePr and diffReps, only direct comparisons between any two exposure groups (1v1) could be made, rather than the multifactorial comparisons available for differential DNA methylation modifications with MethPipe (RRBS). Therefore, in these cases, a series of 1v1 comparisons were made to conservatively estimate which genes were differentially expressed (RNA-seq) or enriched for H3K27ac (ChIP-seq). To use the RNA-seq data as an example, the 1v1 comparisons were carried out as:
Genes differentially expressed between CORT and CORT + DFP and then include only the subset of genes which were not differentially expressed between saline and DFP.
Changes between CORT and CORT + DFP may be due to DFP or the combination of CORT + DFP; removing those differentially expressed between saline and DFP removes those which are due to DFP alone
Genes differentially expressed between DFP and CORT + DFP and then include only the subset of genes which were not differentially expressed between saline and CORT.
Changes between DFP and CORT + DFP may be due to CORT or the combination of CORT + DFP; removing those differentially expressed between saline and CORT removes those which are due to CORT alone
Intersect of the genes which appear in both list 1 and list 2
Both should be genes differentially expressed due to the combined CORT + DFP exposure
Those genes which appear in list 3 and are differentially expressed between saline and CORT + DFP
Final check, as they should be different in the CORT + DFP exposure vs saline
This provides a list of changes unique to the combined CORT + DFP exposure, which are not seen in either CORT or DFP exposure alone. Note that this is conservative, as genes with low but significant differential expression in either CORT or DFP alone, but a larger change in the combined CORT + DFP exposure, will be lost.
For the sets of significant genes identified by RNA-seq, RRBS and ChIP-seq, as well as those genes that were significant in two or more of these analyses, gene set enrichment analyses were carried out. The R package clusterProfiler 3.4.4  was used for Gene Ontology Biological Process (GO BP) [67, 68] and KEGG pathway [69, 70] enrichment analysis, with p and q value cutoffs of ≤ 0.05. Reactome pathway analyses were carried out using the ReactomePA 1.20.2 R package , with a p value cutoff of ≤ 0.05. All three packages reference the latest versions of their respective databases. These significantly enriched annotations were then visualised with the ‘enrichMap’ function of DOSE 3.2.0 R package , with parameters altered to aid legibility with different numbers of enriched annotations.
Overlap between gene sets
The overlap between gene sets was visualised using the UpSetR R package . This provides similar information to a Venn diagram, but in a way which makes proportions clear. Significance of overlap was determined by permutation analysis: random gene sets of the same size as our observed gene sets were taken from the same annotations, and the overlap between the two random gene sets was recorded. This was repeated 1,000,000 times, and the number of occurrences of an overlap equal to or larger than our observed overlap was divided by 1,000,000 to give an empirical p value.
Differential gene expression
Samples were clustered using a Poisson dissimilarity matrix to determine if samples from the same exposure group showed similar expression profiles. As can be seen in Additional file 1: Figure S1 and Additional file 2: Figure S2, the samples largely clustered by exposure group. The only sample that appeared to be out of place was a CORT + DFP sample in the frontal cortex that appeared intermediate between CORT and DFP alone.
In the frontal cortex, the RNA-seq analysis identified 206 GENCODE genes (204 with unique entrez IDs) that were uniquely differentially expressed in the CORT + DFP exposure group compared to all other groups (Additional file 3: Table S1). Enrichment analysis showed 12 enriched KEGG pathways (Additional file 4: Figure S3; Additional file 5: Table S2) and 24 enriched GO BP annotations (Fig. 2; Additional file 6: Table S3). These annotations formed several broad groups related to immune response, including chemokine production, oxidative stress and steroid biosynthesis.
In the hippocampus, 667 GENCODE genes (637 with unique entrez IDs) were uniquely differentially expressed in the CORT + DFP exposure group (Additional file 7: Table S4) compared to all other groups. Enrichment analysis showed 19 enriched KEGG pathways (Additional file 8: Figure S4, Additional file 9: Table S5) and 294 enriched GO BP annotations (Fig. 3, Additional file 10: Table S6). Similar to the frontal cortex, these annotations were grouped into several clusters (Fig. 3), including immune-related annotations (e.g. I-kappaB and NF-kappaB signalling), annotations related to nervous system differentiation, and development.
The two analyses revealed two overlapping KEGG annotations: cytokine–cytokine receptor interaction and rheumatoid arthritis and nine GO BP annotations (Additional file 11: Table S7). There are 32 genes differentially expressed under CORT + DFP priming and exposure found in both the cortex and hippocampus RNA-seq data (Additional file 12: Figure S5; Additional file 13: Table S8).
Enrichment analysis of these 32 genes showed 31 enriched KEGG pathways, including annotations such as rheumatoid arthritis and cytokine–cytokine receptor interaction (Additional file 14: Table S9) and 333 GO BP annotations, including positive regulation of steroid biosynthetic process, positive regulation of chemokine production and regulation of I-kappaB kinase/NF-kappaB signaling (Additional file 15: Table S10).
DNA methylation modifications
We next examined DNA methylation in the frontal cortex and hippocampus using RRBS to identify DNA methylation modifications associated with the exposures. The frontal cortex RRBS data showed 297 differentially methylated cytosines corresponding to 60 differentially methylated regions. Once these regions were annotated to genes (53 entrez IDs; 60 GENCODE; Additional file 16: Table S11), there was no significant enrichment for any KEGG or GO BP annotations. The hippocampus RRBS data showed 926 differentially methylated cytosines corresponding to 192 differentially methylated regions and annotated to 98 GENCODE genes (95 unique entrez IDs; Additional file 17: Table S12). Enrichment analysis was carried out for KEGG pathways or GO BP annotations, showing three significant GO BP enrichments: norepinephrine metabolic process (n = 3, p = 0.048, q = 0.045), cilium morphogenesis (n = 7, p = .048, q = 0.045) and cilium organization (n = 7, p = 0.049, q = 0.046). It is interesting to note in relation to the acetylcholinesterase action of DFP that a CpG site within the acetylcholinesterase gene (Ache) was significantly differentially methylated in the hippocampus (chr5: 137291317; adjp = 0.0296).
Given our hypothesis that DNA methylation modifications contribute to long-term changes in gene expression as a function of GWI exposures, this apparent lack of large, coordinated changes in DNA methylation was unexpected but could have at least two explanations. First, DNA methylation is thought to be relatively stable, and therefore, there may not have been an opportunity for substantial methylation changes to have occurred only 6 h after DFP exposure. A second possibility is that any changes were confounded by the number of different cell types within the brain. Consequently, methylation changes from any single cell type, especially cells making a small proportion of the tissue, may be lost in the ‘noise’. To investigate this second possibility, we used RNA-seq data to estimate the proportions of cells in our two tissues.
The estimated average proportion of each of our five cell types of interest for each exposure group is shown in Fig. 4. In the rat cortex, neurons make up ~ 40% of cells  and ~ 44% of whole mouse brain . In the human and mouse cortex, microglia make up ~ 5% of cells [76, 77]. These reports are in line with our estimates of ~ 40–50% of cells being neurons and ~ 4–6% of cells being microglia. This, therefore, may indicate that enriching for specific cell types, such as microglia, may enhance our ability to detect cell-type-specific methylation modifications due to these exposures. For example, currently only ~ 1:25 RRBS counts will come from microglia.
An interesting incidental finding in the cortex was that CORT exposure, with or without co-exposure with DFP, was associated with an increase in the proportion of neurons and a decrease in the proportion of myelinating oligodendrocytes (MOs) in the frontal cortex (Fig. 4b). As we would not expect neurogenesis to occur in the frontal cortex, this suggests that the increase in the proportion of neurons is driven by a decrease in the absolute number of myelinating oligodendrocytes. A reduced number of oligodendrocytes would be in line with previous work in rats where CORT was shown to reduce the proliferation of oligodendrocytes [78, 79]. We emphasize that these are estimated cell proportions; however, the results indicate that stereology to confirm this will be important in future studies.
H3K27ac is a mark of regions of the genome that are being actively transcribed [59, 60]. Our RNA-seq data showed enrichment for genes involved in histone modification, suggesting that changes in chromatin accessibility may play a role in the response to the exposures. H3K27ac ChIP-seq provides an additional layer of epigenetic regulation, which may respond more quickly than methylation. It also allows an indirect examination of current transcription in the largest cell population, as H3K27ac is found at actively transcribed regions. In ChIP-seq (and RRBS), every locus gives a single signal: either it is enriched for H3K27ac or it is not (or is methylation or is not). However, in RNA-seq, every locus could produce none, one or hundreds of RNA molecules, meaning that a small cell population with large changes in gene expression could mask the signal from a large population with small changes in gene expression. Therefore, using ChIP-seq will allow indirect examination of potential gene expression changes in neurons.
PePr identified 3294 GENCODE genes (3023 entrez IDs; Additional file 18: Table S13) with differential enrichment of H3K27ac, whereas diffReps identified 1518 GENCODE genes (1465 entrez IDs; Additional file 19: Table S14). The overlap between these two analyses was 563 GENCODE genes (557 entrez IDs; Additional file 20: Table S15) which were used for further analysis. However, gene annotation enrichment for each of the two gene sets (PePr and diffReps) demonstrated a large overlap in enriched annotations, suggesting that they are both finding changes in similar pathways but that the individual pathway members they find are different (85% of diffRep and 69% of PePr KEGG pathways (74) are found in both; 68% of diffReps and 54% of PePr GO BP annotations (521) are found in both).
The enrichment data indicated a clear bias towards neuronal-linked annotations, including neuronal morphology and synapse-related annotations (Additional file 21: Figure S6 and Fig. 5; Additional file 22: Table S16 and Additional file 23: Table S17). Of particular interest are the observed enrichment of GO BP annotations ‘cognition’ and ‘learning or memory’, both of which are observed to be disrupted in GWI study participants, ‘Circadian entrainment’ which may relate to observed sleep disruption, and ‘response to steroid hormone’, which likely relates to CORT (a steroid hormone involved in the response to stress).
These findings demonstrate that there are potential changes in neuronal-related gene expression in the frontal cortex, as was also seen in the hippocampus RNA-seq, highlighted by the fact that 33 genes were found in both the ChIP-seq frontal cortex analysis and the hippocampus RNA-seq analysis.
Overlap between genes found in different analyses
As shown in Fig. 6, there is not a large overlap in genes found between any of our analyses. However, this disparity may be partly explained from the aforementioned difference between mRNA and DNA, whereby one locus can produce many mRNA molecules, but DNA either has a modification or does not. This is reflected by the fact that the largest percentage overlap is between those genes found with RRBS and ChIP-seq, as these are both examining DNA modifications: 12% of genes found in frontal cortex RRBS, and 16% in hippocampus RRBS, are also found in the frontal cortex ChIP-seq, whereas this is only 1% and 5% for frontal cortex and hippocampus RNA-seq respectively. Similarly, 15.5% of genes found in the frontal cortex RNA-seq are also found in the hippocampus RNA-seq.
Because of the very few annotations enriched in genes differentially methylated in either of our tissues, we see very little overlap between annotations in methylation and ChIP-seq. However, we do see annotations and pathways enriched in both our ChIP-seq and RNA-seq data (Additional file 24: Figure S7 and Additional file 25: Figure S8).
To our knowledge, this is the first transcriptome- and epigenome-wide study to examine evidence for transcriptional, chromatin and DNA methylation modifications in a model of GWI. A previous study  examined 84 microRNAs (MiRNAs) and global but not gene-specific changes in DNA methylation and hydroxymethylation in rats subjected to restraint stress and a protocol of PB, DEET and permethrin to simulate troops’ chemical exposures. This study found differential expression of two miRNAs and in increase in global methylation in the hippocampus. Although this was an important step forward, our study was able to examine the whole transcriptome, a histone modification associated with chromatin accessibility and DNA methylation at a genome-wide, single cytosine level.
Overall, our results represent several interesting findings. First, as expected, there was a large change in the expression of immune-related genes in both the frontal cortex and hippocampus, building upon previous findings in this model . Second, many genes associated with synaptic function are changed in their activity, as shown by our frontal cortex H3K27ac ChIP-seq data and hippocampus RNA-seq. These changes in gene expression seem to be subtler than those found for immune-related genes (lower expression), but differential expression of genes related to synaptic function in mice is associated with impaired memory and cognition, consistent with impairments reported by GWI suffers. Interestingly, long-term potentiation- and depression-related genes are enriched in the ChIP-seq data. Finally, we see evidence of not just a change in gene activity but of a suggested change in cell proportions. This is in line with previous work with CORT [78, 79].
It is possible that microglia are responsible for this large change in immune-related gene expression, but we acknowledge that other cells, such as astrocytes, also express cytokines and chemokines (e.g. Kim et al. ). However, this is usually at a much lower level than that in microglia: in genes significantly differentially expressed uniquely with CORT + DFP in both the frontal cortex and hippocampus the six genes with the largest fold change are expressed in microglia (Additional file 26: Table S18). Therefore, although other cell types may be contributing to cytokine and chemokine expression, it is highly unlikely that microglia are not the cell type driving this change. Similarly, we attribute many of the transmembrane transporter-related annotations we see in the frontal cortex ChIP-seq data and hippocampus RNA-seq data to neurons; however, many of these transporters are also expressed in glial cells . For this reason, future work should be carried out to isolate or enrich specific cell populations, allowing these predictions to be tested.
Our findings of altered cholinergic neurotransmitter expression after DFP exposure are in line with previous studies of low dose sarin exposure in rats [81, 82]. Specifically, changes in M1 and M3 acetylcholine receptor expression were seen in the frontal cortex in a dose-dependent manner when the animals were maintained under hyper-thermic conditions (i.e. a stressor; Henderson et al. [81, 82]). In relation to this, we see choline-related annotations at all levels we analysed: our H3K27ac ChIP-seq analyses show changes in cholinergic synapse-related genes and differential binding related to Chrm2, the gene coding for M2; the KEGG annotation, ‘choline metabolism in cancer’, is enriched in our hippocampus RNA-seq genes and our H3K27ac ChIP-seq genes; finally, in our hippocampus RRBS, we see altered methylation of a CpG site within Ache, the gene coding for acetylcholinesterase.
These findings not only link to rat models but also to GWI subjects, where differences in prefrontal cortex working memory have been previously found and have been attributed to the cholinergic system . This again connects with our findings from ChIP-seq where both ‘learning and memory’ and ‘cholinergic synapse’ annotations are significantly enriched. Therefore, pre-exposure with CORT may not only be potentiating the immune response to AchE inhibitors but also causing longer term changes to the cholinergic signalling system in line with previous animal models. Human studies have been more varied, with evidence both for [84, 85] and against  long-term changes in the cholinergic system, perhaps related to the heterogeneity of symptoms or to differences in tissues being examined. Recent work by Locker et al.  has shown that increased neuroinflammatory response in this model is not directly related to AChE inhibition but instead may result from changes in other aspects of signaling (e.g. epigenetic alterations in gene expression).
In relation to the potential reduction in myelinating oligodendrocytes in the cortex, this may have an effect on some of the phenotypes seen in GWI: reduced oligodendrocytes have been linked to major depressive disorder (MDD), functional consequences in neurons and mood-related symptoms in rats . As this change in cell proportion would affect myelinating cells, it could also contribute to the reported alterations in white matter in GWI veterans [20, 88,89,90].
Our RNA-seq data suggests both immune dysfunction and oxidative stress. Previous reports have shown immune dysfunction and oxidative stress in other disorders with similar phenotypes to GWI such as myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) . However, whether the immune response is causing oxidative stress or oxidative stress is causing an immune response is less clear. This is directly relevant, as oxidative stress contributes to the toxicity of AchE inhibitors  and has been suggested as a cause of GWI [1, 92].
One of the six genes (Additional file 26: Table S18) we find strongly differentially expressed in both the frontal cortex and hippocampus, Tlr2, has previously been linked to sickness behaviour when expressed in hypothalamic microglia  and was reported to show increased expression in a model of GWI . Therefore, it is of great interest in terms of the sickness-like behaviour that defines GWI. TLR activation by immune insult appears to increase Tlr2 expression . We can speculate that the change in Tlr2 expression we see is upstream of the change in Il1b and TNF, as specific activation of TLR2 increases their expression [93, 96]. Further, TLR stimulation increases SLC15A3 expression in dendritic cells, which in turn regulates TNF and IL-1β expression . At the very least, these six, consistently, strongly upregulated genes make up a robust core of immune- and microglia-related genes for further investigation (Additional file 26: Table S18).
Our results correspond very well with those of Broderick et al.  in blood samples from GWI study participants, who found changes to NF-KappaB-related genes (which are enriched in our hippocampus RNA-seq data and in the genes significant in both the cortex and hippocampus RNA-seq), and in pathways under the broad theme of neuronal development and migration, which we see in our hippocampus RNA-seq and our cortex ChIP-seq. They also saw ligand–receptor interactions supporting neurotransmission, which again we see strongly in our cortex ChIP-seq data. Further, our KEGG pathway analysis of the transcriptomics data highlighted rheumatoid arthritis-related gene enrichment. Interestingly, a study in veterans with GWI also highlighted the possibility of medication used to treat rheumatoid arthritis also being used to treat GWI . The same study identified the TNF-alpha pathway (which we see in our transcriptomics) and the estrogen pathway (which we saw in our ChIP-seq) as potential drug targets . Therefore, data from our model corresponds to data identified from veterans with GWI and could be used as a model to test these compounds as potential therapeutics.
A recent paper described eight potential blood biomarkers for GWI, the strongest being a 9.27-fold increase in CaMKII protein in the blood of these veterans . In our study, Camk2b was found to be differentially methylated in the hippocampus and have differential H3K27ac in the cortex.
Further linking our study with human data was the finding that plasma CRP levels are increased in the blood of GWI veterans . CRP is often used as a biomarker of IL-6-mediated inflammation [38, 101], and IL-6 related annotations were enriched in the genes differentially expressed in both the cortex and hippocampus (Additional file 15: Table S10). This annotation was due to four genes, Il1b, Tlr2, Il1a and Tnf, four of our six most strongly and consistently differentially expressed genes.
We see minimal overlap in genes with significant differential methylation between the hippocampus and frontal cortex (two genes, Bcar3 and Tmem242). This is likely due, at least partly, to the two factors outlined above: cellular heterozygosity and a short time point after treatment. Of the two consistent genes, Bcar3 is involved in cell proliferation, which when overexpressed in breast cancer conveys estrogen resistance, and Tmem242 is a transmembrane protein with very little current annotation. We cannot detect any coordinated methylation changes (denoted by enriched annotations in our significant genes) in the frontal cortex and very few in the hippocampus. As mentioned above, we detect differential methylation of one CpG site within Ache, the gene coding for acetylcholinesterase. A number of the genes found to be differentially methylated can be linked to GWI. Examples include Lims1, which is known to be regulated by TNF and is involved in cell growth and survival ; Sesn1, Aplnr, Pxn and Actn1, which have been linked to ME/CFS [103, 104]; Col5a3, which was identified in a rat model of GWI in the hippocampus ; and Slc1a2, which has been linked to ALS, a disorder GW veterans are at increased risk of [105, 106]. This may indicate that networks of genes linked to GWI are only just beginning to be methylated, or that different groups of genes are methylated in different cell types. However, until these coordinated networks are elucidated, investigation of individual genes could lead to spurious associations.
We are able to show a range of changes in the transcriptome of this well-established mouse model, many of which reflect gene expression seen in veterans with GWI. Further, we see alteration in H3K27ac, showing potential chromatin configuration changes, which could lead to epigenetic effects with long-lasting implications. We also find differences in DNA methylation, although these are less easily interpretable than the transcriptome and H3K27ac changes.
Additional research is needed to assess whether effects of these epigenetic and transcriptional modifications on long-term health outcomes are cumulative and/or are potentiated by later exposures (e.g. infection). Notably, however, our findings reveal gene pathways known to be involved in long-term adverse health effects in GWI veterans. These results suggest that epigenetic and transcriptional regulation during the initial exposure period likely contribute to pathological outcomes in GWI. It would be important to examine these modifications in peripheral tissues from GWI veterans to ascertain whether biomarkers could be developed to predict future health outcomes.
Native chromatin immunoprecipitation sequencing
Differentially methylated cytosines
Differentially methylated regions
Disease Ontology Semantic and Enrichment
Gene Expression Omnibus
- GO BP:
Gene Ontology Biological Process
Gulf War illnesses
Kyoto Encyclopedia of Genes and Genomes
Myalgic encephalomyelitis/chronic fatigue syndrome
Oligodendrocyte precursor cells
Peripheral blood mononuclear cells
Ribonucleic acid sequencing
Reduced representation bisulfite sequencing
Tumor necrosis factor
Kerr KJ. Gulf War illness: an overview of events, most prevalent health outcomes, exposures, and clues as to pathogenesis. Rev Environ Health. 2015;30:273–86.
Committee on Gulf War and Health. Cory-Slechta D, Wedge R, editors. Gulf War and health: volume 10: update of health effects of serving in the Gulf War, 2016. Washington, D.C.: National Academies Press. p. 2016.
Dantzer R, O’Connor JC, Freund GG, Johnson RW, Kelley KW. From inflammation to sickness and depression: when the immune system subjugates the brain. Nat Rev Neurosci. 2008;9:46–56.
O’Callaghan JP, Kelly KA, Locker AR, Miller DB, Lasley SM. Corticosterone primes the neuroinflammatory response to DFP in mice: potential animal model of Gulf War illness. J Neurochem. 2015;133:708–21.
Gavaghan H. NIH panel rejects Persian Gulf syndrome. Nature. 1994;369:8.
Beale P. Gulf illness. BMJ. 1994;308:1574.
Robinson A. Veterans worry that unexplained medical problems a legacy of service during Gulf War. CMAJ. 1995;152:944–7.
White RF, Steele L, O’Callaghan JP, Sullivan K, Binns JH, Golomb BA, et al. Recent research on Gulf War illness and other health problems in veterans of the 1991 Gulf War: effects of toxicant exposures during deployment. Cortex. 2016;74:449–75.
Barth SK, Kang HK, Bullman TA, Wallin MT. Neurological mortality among U.S. veterans of the Persian Gulf War: 13-year follow-up. Am J Ind Med. 2009;52:663–70.
Zakirova Z, Crynen G, Hassan S, Abdullah L, Horne L, Mathura V, et al. A chronic longitudinal characterization of neurobehavioral and neuropathological cognitive impairment in a mouse model of Gulf War agent exposure. Front Integr Neurosci. 2015;9:71.
Abdullah L, Evans JE, Bishop A, Reed JM, Crynen G, Phillips J, et al. Lipidomic profiling of phosphocholine-containing brain lipids in mice with sensorimotor deficits and anxiety-like features after exposure to Gulf War agents. NeuroMolecular Med. 2012;14:349–61.
Pierce LM, Kurata WE, Matsumoto KW, Clark ME, Farmer DM. Long-term epigenetic alterations in a rat model of Gulf War illness. Neurotoxicology. 2016;55:20–32.
Research Advisory Committee (RAC) on Gulf War Veterans’ Illnesses. Gulf War illness and the health of Gulf War veterans: research update and recommendations, 2009-2013. Washington: U.S. Government Printing Office; 2014. https://www.va.gov/RAC-GWVI/RACReport2014Final.pdf
Steele L, Sastre A, Gerkovich MM, Cook MR. Complex factors in the etiology of Gulf War illness: wartime exposures and risk factors in veteran subgroups. Environ Health Perspect. 2012;120:112–8.
Locker AR, Michalovicz LT, Kelly KA, Miller JV, Miller DB, O’Callaghan JP. Corticosterone primes the neuroinflammatory response to Gulf War illness-relevant organophosphates independently of acetylcholinesterase inhibition. J Neurochem. 2017;142:444–55.
Koo B-B, Michalovicz LT, Calderazzo S, Kelly KA, Sullivan K, Killiany RJ, et al. Corticosterone potentiates DFP-induced neuroinflammation and affects high-order diffusion imaging in a rat model of Gulf War illness. Brain Behav Immun. 2018;67:42–6.
Golomb BA. Acetylcholinesterase inhibitors and Gulf War illnesses. Proc Natl Acad Sci U S A. 2008;105:4295–300.
Winkenwerder W. Environmental exposure report: pesticides final report. U.S. Dep. Defense, Off. Spec. Assist. to Undersecretary Def. (Personnel Readiness) Gulf War Illnesses Med. Readiness Mil. Deployments. Washington: US Department of Defense; 2003. https://gulflink.health.mil/pest_final/index.html
Sullivan K, Krengel M, Bradford W, Stone C, Thompson TA, Heeren T, et al. Neuropsychological functioning in military pesticide applicators from the Gulf War: effects on information processing speed, attention and visual memory. Neurotoxicol Teratol. 2017;65:1–13.
Research Advisory Committee (RAC) on Gulf War Veterans’ Illnesses. Gulf War illness and the health of Gulf War veterans: scientific findings and recommendations. Washington: U.S. Government Printing Office; 2008. https://www.va.gov/RAC-GWVI/docs/Committee_Documents/GWIandHealthofGWVeterans_RAC-GWVIReport_2008.pdf
Rao GV, Rao KS. Modulation in acetylcholinesterase of rat brain by pyrethroids in vivo and an in vitro kinetic study. J Neurochem. 1995;65:2259–66.
Corbel V, Stankiewicz M, Pennetier C, Fournier D, Stojan J, Girard E, et al. Evidence for inhibition of cholinesterases in insect and mammalian nervous systems by the insect repellent deet. BMC Biol. 2009;7:47.
Tuite JJ, Haley RW. Meteorological and intelligence evidence of long-distance transit of chemical weapons fallout from bombing early in the 1991 Persian Gulf War. Neuroepidemiology. 2013;40:160–77.
Haley RW, Tuite JJ. Epidemiologic evidence of health effects from long-distance transit of chemical weapons fallout from bombing early in the 1991 Persian Gulf War. Neuroepidemiology. 2013;40:178–89.
McGowan PO, Matthews SG. Prenatal stress, glucocorticoids, and developmental programming of the stress response. Endocrinology. 2018;159:69–82.
Klengel T, Binder EB. Epigenetics of stress-related psychiatric disorders and gene × environment interactions. Neuron. 2015;86:1343–57.
Faulk C, Dolinoy DC. Timing is everything: the when and how of environmentally induced changes in the epigenome of animals. Epigenetics. 2011;6:791–7.
Montrose L, Faulk C, Francis J, Dolinoy DC. Perinatal lead (Pb) exposure results in sex and tissue-dependent adult DNA methylation alterations in murine IAP transposons. Environ Mol Mutagen. 2017;58:540–50.
Colović MB, Krstić DZ, Lazarević-Pašti TD, Bondžić AM, Vasić VM. Acetylcholinesterase inhibitors: pharmacology and toxicology. Curr Neuropharmacol. 2013;11:315–35.
Zhang Q, Zhou XD, Denny T, Ottenweller JE, Lange G, LaManca JJ, et al. Changes in immune parameters seen in Gulf War veterans but not in civilians with chronic fatigue syndrome. Clin Diagn Lab Immunol. 1999;6:6–13.
Skowera A, Hotopf M, Sawicka E, Varela-Calvino R, Unwin C, Nikolaou V, et al. Cellular immune activation in Gulf War veterans. J Clin Immunol. 2004;24:66–73.
Vojdani A, Thrasher JD. Cellular and humoral immune abnormalities in Gulf War veterans. Environ Health Perspect. 2004;112:840–6.
Whistler T, Fletcher MA, Lonergan W, Zeng X-R, Lin J-M, Laperriere A, et al. Impaired immune function in Gulf War illness. BMC Med Genet. 2009;2:12.
Broderick G, Fletcher MA, Gallagher M, Barnes Z, Vernon SD, Klimas NG. Exploring the diagnostic potential of immune biomarker coexpression in Gulf War illness. Methods Mol Biol. 2012;934:145–64.
Smylie AL, Broderick G, Fernandes H, Razdan S, Barnes Z, Collado F, et al. A comparison of sex-specific immune signatures in Gulf War illness and chronic fatigue syndrome. BMC Immunol. 2013;14:29.
Khaiboullina SF, DeMeirleir KL, Rawat S, Berk GS, Gaynor-Berk RS, Mijatovic T, et al. Cytokine expression provides clues to the pathophysiology of Gulf War illness and myalgic encephalomyelitis. Cytokine. 2015;72:1–8.
Parkitny L, Middleton S, Baker K, Younger J. Evidence for abnormal cytokine expression in Gulf War illness: a preliminary analysis of daily immune monitoring data. BMC Immunol. 2015;16:57.
Johnson GJ, Slater BCS, Leis LA, Rector TS, Bach RR. Blood biomarkers of chronic inflammation in Gulf War illness. PLoS One. 2016;11:e0157855. Block ML, editor.
O’Callaghan JP. Quantification of glial fibrillary acidic protein: comparison of slot-immunobinding assays with a novel sandwich ELISA. Neurotoxicol Teratol. 1991;13:275–81.
Krueger F, The Babraham Institute. Trim Galore! [Internet]. [cited 2016 Aug 1]. Available from: http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/. Accessed 1 Aug 2016.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjournal. 2011;17:10.
Andrews S, The Babraham Institute. FastQC [Internet]. [cited 2017 Aug 1]. Available from: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 1 Aug 2016.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Mudge JM, Harrow J. Creating reference gene annotation for the mouse C57BL6/J genome assembly. Mamm Genome. 2015;26:366–78.
Harrow J, Denoeud F, Frankish A, Reymond A, Chen C-K, Chrast J, et al. GENCODE: producing a reference annotation for ENCODE. Genome Biol. 2006;7(Suppl 1):S4.1–9.
Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9:e1003118. Prlic A, editor
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Khang TF, Lau CY. Getting the most out of RNA-seq data analysis. PeerJ. 2015;3:e1360.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.
Witten DM. Classification and clustering of sequencing data using a Poisson model. Ann Appl Stat. 2011;5:2493–518.
Gong T, Szustakowski JD. DeconRNASeq: a statistical framework for deconvolution of heterogeneous tissue samples based on mRNA-Seq data. Bioinformatics. 2013;29:1083–5.
Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002;30:207–10.
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013;41:D991–5.
Zhang Y, Chen K, Sloan SA, Bennett ML, Scholze AR, O’Keeffe S, et al. An RNA-sequencing transcriptome and splicing database of glia, neurons, and vascular cells of the cerebral cortex. J Neurosci. 2014;34:11929–47.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Song Q, Decato B, Hong EE, Zhou M, Fang F, Qu J, et al. A reference methylome database and analysis pipeline to facilitate integrative and comparative epigenomics. PLoS One. 2013;8:e81148. El-Maarri O, editor
Dolzhenko E, Smith AD. Using beta-binomial regression for high-precision differential methylation analysis in multifactor whole-genome bisulfite sequencing experiments. BMC Bioinformatics. 2014;15:215.
Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A. 2010;107:21931–6.
Rye M, Sætrom P, Håndstad T, Drabløs F. Clustered ChIP-Seq-defined transcription factor binding sites and histone modifications map distinct classes of regulatory elements. BMC Biol. 2011;9:80.
Steinhauser S, Kurzawa N, Eils R, Herrmann C. A comprehensive comparison of tools for differential ChIP-seq analysis. Brief Bioinform. 2016;17:953–66.
Zhang Y, Lin Y-H, Johnson TD, Rozek LS, Sartor MA. PePr: a peak-calling prioritization pipeline to identify consistent or differential peaks from replicated ChIP-Seq data. Bioinformatics. 2014;30:2568–75.
Kikutake C, Yahara K. Identification of epigenetic biomarkers of lung adenocarcinoma through multi-omics data analysis. PLoS One. 2016;11:e0152918. Pradhan S, editor
Shen L, Shao N-Y, Liu X, Maze I, Feng J, Nestler EJ. diffReps: detecting differential chromatin modification sites from ChIP-seq data with biological replicates. PLoS One. 2013;8:e65598. Mantovani R, editor
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.
Gene Ontology Consortium. Gene Ontology Consortium: going forward. Nucleic Acids Res. 2015;43:D1049–56.
Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30.
Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40:D109–14.
Yu G, He Q-Y. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol BioSyst. 2016;12:477–9.
Yu G, Wang L-G, Yan G-R, He Q-Y. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. 2015;31:608–9.
Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017;33:2938–40.
Herculano-Houzel S, Lent R. Isotropic fractionator: a simple, rapid method for the quantification of total cell and neuron numbers in the brain. J Neurosci. 2005;25:2518–25121.
Herculano-Houzel S. The glia/neuron ratio: how it varies uniformly across brain structures and species and what that means for brain physiology and evolution. Glia. 2014;62:1377–91.
Darmanis S, Sloan SA, Zhang Y, Enge M, Caneda C, Shuer LM, et al. A survey of human brain transcriptome diversity at the single cell level. Proc Natl Acad Sci U S A. 2015;112:7285–90.
Lawson LJ, Perry VH, Dri P, Gordon S. Heterogeneity in the distribution and morphology of microglia in the normal adult mouse brain. Neuroscience. 1990;39:151–70.
Alonso G. Prolonged corticosterone treatment of adult rats inhibits the proliferation of oligodendrocyte progenitors present throughout white and gray matter regions of the brain. Glia. 2000;31:219–31.
Banasr M, Valentine GW, Li X-Y, Gourley SL, Taylor JR, Duman RS. Chronic unpredictable stress decreases cell proliferation in the cerebral cortex of the adult rat. Biol Psychiatry. 2007;62:496–504.
Kim S, Steelman AJ, Koito H, Li J. Astrocytes promote TNF-mediated toxicity to oligodendrocyte precursors. J Neurochem. 2011;116:53–66.
Henderson RF, Barr EB, Blackwell WB, Clark CR, Conn CA, Kalra R, et al. Response of rats to low levels of sarin. Toxicol Appl Pharmacol. 2002;184:67–76.
Henderson RF, Barr EB, Blackwell WB, Clark CR, Conn CA, Kalra R, et al. Response of F344 rats to inhalation of subclinical levels of sarin: exploring potential causes of Gulf War illness. Toxicol Ind Health. 2001;17:294–7.
Hubbard NA, Hutchison JL, Motes MA, Shokri-Kojori E, Bennett IJ, Brigante RM, et al. Central executive dysfunction and deferred prefrontal processing in veterans with Gulf War illness. Clin Psychol Sci a J Assoc Psychol Sci. 2014;2:319–27.
Liu P, Aslan S, Li X, Buhner DM, Spence JS, Briggs RW, et al. Perfusion deficit to cholinergic challenge in veterans with Gulf War illness. Neurotoxicology. 2011;32:242–6.
Haley RW, Charuvastra E, Shell WE, Buhner DM, Marshall WW, Biggs MM, et al. Cholinergic autonomic dysfunction in veterans with Gulf War illness: confirmation in a population-based sample. JAMA Neurol. 2013;70:191–200.
Concato J, Aslan M, Palmisano MM, Doebbeling CC, Peduzzi P, Ofek K, et al. Acetylcholinesterase activity in veterans of the first Gulf War. J Investig Med. 2007;55:360–7.
Edgar N, Sibille E. A putative functional role for oligodendrocytes in mood regulation. Transl Psychiatry. 2012;2:e109.
Heaton KJ, Palumbo CL, Proctor SP, Killiany RJ, Yurgelun-Todd DA, White RF. Quantitative magnetic resonance brain imaging in US army veterans of the 1991 Gulf War potentially exposed to sarin and cyclosarin. Neurotoxicology. 2007;28:761–9.
Chao LL, Zhang Y, Buckley S. Effects of low-level sarin and cyclosarin exposure on white matter integrity in Gulf War veterans. Neurotoxicology. 2015;48:239–48.
Rayhan RU, Stevens BW, Timbol CR, Adewuyi O, Walitt B, VanMeter JW, et al. Increased brain white matter axial diffusivity associated with fatigue, pain and hyperalgesia in Gulf War illness. PLoS One. 2013;8:e58493. Zang Y-F, editor.
Morris G, Berk M. The many roads to mitochondrial dysfunction in neuroimmune and neuropsychiatric disorders. BMC Med. 2015;13:68.
Golomb BA, Allison M, Koperski S, Koslik HJ, Devaraj S, Ritchie JB. Coenzyme Q10 benefits symptoms in Gulf War veterans: results of a randomized double-blind study. Neural Comput. 2014;26:2594–651.
Jin S, Kim JG, Park JW, Koch M, Horvath TL, Lee BJ. Hypothalamic TLR2 triggers sickness behavior via a microglia-neuronal axis. Sci Rep. 2016;6:29424.
Michalovicz L, Kelly K, O’Callaghan J, Locker A, Miller D. Corticosterone priming of the neuroinflammatory response to AChE inhibitors results in overexpression of Tlr2 and downstream targets, but not activation of the Nlrp3 inflammasome. Toxicol Suppl to Toxicol Sci. 2016;150:49–50.
Regueiro V, Moranta D, Campos MA, Margareto J, Garmendia J, Bengoechea JA. Klebsiella pneumoniae increases the levels of toll-like receptors 2 and 4 in human airway epithelial cells. Infect Immun. 2009;77:714–24.
Facci L, Barbierato M, Marinelli C, Argentini C, Skaper SD, Giusti P. Toll-like receptors 2, -3 and -4 prime microglia but not astrocytes across central nervous system regions for ATP-dependent interleukin-1β release. Sci Rep. 2014;4:6824.
Nakamura N, Lill JR, Phung Q, Jiang Z, Bakalarski C, de Mazière A, et al. Endosomes are specialized platforms for bacterial sensing and NOD2 signalling. Nature. 2014;509:240–4.
Broderick G, Ben-Hamo R, Vashishtha S, Efroni S, Nathanson L, Barnes Z, et al. Altered immune pathway activity under exercise challenge in Gulf War illness: an exploratory analysis. Brain Behav Immun. 2013;28:159–69.
Craddock TJA, Harvey JM, Nathanson L, Barnes ZM, Klimas NG, Fletcher MA, et al. Using gene expression signatures to identify novel treatment strategies in Gulf War illness. BMC Med Genet. 2015;8:36.
Abou-Donia MB, Conboy LA, Kokkotou E, Jacobson E, Elmasry EM, Elkafrawy P, et al. Screening for novel central nervous system biomarkers in veterans with Gulf War illness. Neurotoxicol Teratol. 2017;61:36–46.
Thiele JR, Zeller J, Bannasch H, Stark GB, Peter K, Eisenhardt SU. Targeting C-reactive protein in inflammatory disease by preventing conformational changes. Mediat Inflamm. 2015;2015:372432.
Jatiani A, Pannizzo P, Gualco E, Del-Valle L, Langford D. Neuronal PINCH is regulated by TNF-α and is required for neurite extension. J NeuroImmune Pharmacol. 2011;6:330–40.
Klimas NG, Broderick G, Fletcher MA. Biomarkers for chronic fatigue. Brain Behav Immun. 2012;26:1202–10.
Smith AK, Fang H, Whistler T, Unger ER, Rajeevan MS. Convergent genomic studies identify association of GRIK2 and NPAS2 with chronic fatigue syndrome. Neuropsychobiology. 2011;64:183–94.
Cloutier F, Marrero A, O’Connell C, Morin P. MicroRNAs as potential circulating biomarkers for amyotrophic lateral sclerosis. J Mol Neurosci. 2015;56:102–12.
Lee J, Hyeon SJ, Im H, Ryu H, Kim Y, Ryu H. Astrocytes and microglia as non-cell autonomous players in the pathogenesis of ALS. Exp Neurobiol. 2016;25:233–40.
Undergraduate assistants in the McGowan lab, Lisa Shao and Bryant Lim, assisted with the aspects of the lab work related to this project. We appreciate the excellent technical assistance provided by Brenda K. Billig, Christopher M. Felton and Ali Yilmaz.
This project was funded through the Congressionally Directed Medical Research Programs (CDMRP) from the US Department of Defence (DoD) (GW130053). This work was supported by the Assistant Secretary of Defense for Health Affairs, through the Gulf War Illness Research Program, GW130053. Opinions, interpretations, conclusions and recommendations are those of the authors and are not necessarily endorsed by the Department of Defense. The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of the National Institute for Occupational Safety and Health, Centers for Disease Control and Prevention.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. The datasets generated and/or analysed during the current study are not publicly available as they are being used for follow-up studies. They will be made publicly available after future manuscripts are published.
Ethics approval and consent to participate
All procedures were performed under protocols approved by the Institutional Animal Care and Use Committee of the Centers for Disease Control and Prevention, National Institute for Occupational Safety and Health and the US Army Medical Research and Materiel Command Animal Care and Use Review Office. The animal facility was certified by AAALAC International.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Poisson dissimilarity matrix for hippocampus RNA-seq. Heatmaps showing the Poisson dissimilarity matrix for our hippocampus RNA-seq samples. Darker blue indicates more similar samples. (CSV 52 kb)
Figure S2. Poisson dissimilarity matrix from frontal cortex RNA-seq. Heatmaps showing the Poisson dissimilarity matrix for our frontal cortex (B) RNA-seq samples. Darker blue indicates more similar samples. (CSV 2 kb)
Table S1. Genes identified as uniquely differentially expressed in the frontal cortex after CORT + DFP exposure, compared to all other groups. (TIFF 2391 kb)
Figure S3. Frontal cortex RNA-seq significantly enriched KEGG pathway annotations. KEGG pathways significantly enriched in genes that were differentially expressed in the frontal cortex of CORT + DFP exposed mice. (CSV 6 kb)
Table S2. KEGG pathways enriched in the genes identified as uniquely differentially expressed in the frontal cortex after CORT + DFP exposure, compared to all other groups. (CSV 4 kb)
Table S3. Gene Ontology (GO) Biological Process (BP) annotations enriched in the genes identified as uniquely differentially expressed in the frontal cortex after CORT + DFP exposure, compared to all other groups. (CSV 38 kb)
Table S4. Genes identified as uniquely differentially expressed in the hippocampus after CORT + DFP exposure, compared to all other groups. (CSV 1 kb)
Figure S4. Hippocampus RNA-seq significantly enriched KEGG pathway annotations. KEGG pathways significantly enriched in genes that were differentially expressed in hippocampus of CORT + DFP exposed mice. (CSV 3 kb)
Table S5. KEGG pathways enriched in the genes identified as uniquely differentially expressed in the hippocampus after CORT + DFP exposure, compared to all other groups. (CSV 103 kb)
Table S6. Gene Ontology (GO) Biological Process (BP) annotations enriched in the genes identified as uniquely differentially expressed in the hippocampus after CORT + DFP exposure, compared to all other groups. (CSV 55 kb)
Figure S5. Overlap of genes significantly differentially expressed with CORT + DFP exposure in the frontal cortex and hippocampus. (CSV 18 kb)
Table S10. Gene Ontology (GO) Biological Process (BP) annotations found to be enriched in genes identified as uniquely differentially expressed in both the frontal cortex (Additional file 3: Table S1) and the hippocampus (Additional file 7: Table S4) after CORT + DFP exposure, compared to all other groups. (CSV 12 kb)
Table S11. Genes identified as containing differentially methylated regions in the frontal cortex of combined CORT + DFP exposed animals. (TIFF 3010 kb)
Table S12. Genes identified as containing differentially methylated regions in the hippocampus of combined CORT + DFP exposed animals. (TIFF 3010 kb)
Table S13. Genes identified by PePr as having differential enrichment of H3K27ac in the frontal cortex of CORT + DFP exposed animals. (CSV 519 bytes)
Table S14. Genes identified by diffReps as having differential enrichment of H3K27ac in the frontal cortex of CORT + DFP exposed animals. (TIFF 2386 kb)
Table S15. Genes identified by both PePr and diffReps as having differential enrichment of H3K27ac in the frontal cortex of CORT + DFP exposed animals. (CSV 26 kb)
Figure S6. Frontal cortex H3K27ac ChIP-seq significantly enriched KEGG pathways. Top 50 KEGG pathways significantly enriched for differential enrichment of H3K27ac with CORT + DFP exposure. (TIFF 540 kb)
Table S16. Gene Ontology (GO) Biological Process (BP) annotations found to be enriched in genes identified as uniquely differentially enrichment H3K27ac in by both PePr and diffReps (Additional file 20: Table S15). (CSV 1 kb)
Table S17. KEGG pathways found to be enriched in genes identified as uniquely differentially enrichment H3K27ac in by both PePr and diffReps (Additional file 20: Table S15). (CSV 3 kb)
Figure S7. Comparison of significantly enriched gene ontology biological process annotations in each of our analyses. UpSetR diagram  of Gene Ontology (GO) biological process (BP) annotations enriched in genes significant in each of our differential analyses: frontal cortex RRBS (FC RRBS), frontal cortex H3K27ac ChIP-seq (FC ChIP), hippocampus RNA-seq (Hipp RNA), frontal cortex RNA-seq (FC RNA) and hippocampus RRBS (Hipp RRBS). (CSV 79 kb)
Figure S8. Comparison of significantly enriched KEGG pathway annotations in each of our analyses. UpSetR  diagram of KEGG pathways enriched in genes significant in each of our differential analyses: frontal cortex RRBS (FC RRBS), frontal cortex H3K27ac ChIP-seq (FC ChIP), hippocampus RNA-seq (Hipp RNA), frontal cortex RNA-seq (FC RNA) and hippocampus RRBS (Hipp RRBS). (TIFF 317 kb)
Table S18. Genes which are significantly differentially expressed uniquely in CORT + DFP and have a greater than 0.75-fold change in expression between saline and CORT + DFP in both frontal cortex and hippocampus. All six show an increase in expression. (CSV 3 kb)
About this article
- Gulf War illness
- Diisopropyl fluorophosphate
- Acetylcholinesterase inhibitors