Region- and time-dependent gene regulation in the amygdala and anterior cingulate cortex of a PTSD-like mouse model

Posttraumatic stress disorder is developed by exposure to a threatening and/or a horrifying event and characterized by the presence of anxiety, hyperarousal, avoidance, and sleep abnormality for a prolonged period of time. To elucidate the potential molecular mechanisms, we constructed a mouse model by electric foot shock followed by situational reminders and performed transcriptome analysis in brain tissues. The stressed mice acquired anxiety-like behavior after 2 weeks and exaggerated startle response after 4 weeks. Avoidance latency and freezing behavior were sustained up to 5 weeks post stress and abnormal static behavior was observed during the sleep period. RNA sequencing was performed in two of the emotional regulatory regions, anterior cingulate cortex and amygdala, at 2 and 5 weeks post stress. More than 1000 differentially expressed genes were identified at 2 weeks in both regions. The number of the regulated genes remained constant in amygdala at 5 weeks post stress, whereas those in anterior cingulate cortex were plummeted. Although synaptic remodeling and endocrine system were the most enriched signaling pathways in both anterior cingulate cortex and amygdala, the individual gene expression profile was regulated in a region- and time-dependent manner. In addition, several genes associated with PTSD involved in Hypothalamic-Pituitary-Adrenal axis were differentially regulated. These findings suggested that global gene expression profile was dynamically regulated in accordance with the disease development stage, and therefore targeting the distinct signaling molecules in different region and development stage might be critical for effective treatment to PTSD. Electronic supplementary material The online version of this article (10.1186/s13041-019-0449-0) contains supplementary material, which is available to authorized users.


Introduction
Post-traumatic stress disorder (PTSD) is a psychiatric disorder defined by profound disturbances in cognition, emotion, behavior and physiological function in response to a life threatening or a horrifying traumatic event. The symptoms of PTSD are characterized by three long-lasting symptomatic clusters (DSM-5): (1) re-experiencing (intrusive thoughts, recurrent nightmares, flashbacks and reactivity to reminders of the trauma); (2) avoidance and emotional numbness (avoiding stimuli associated with trauma and anhedonia); and (3) hyperarousal (hypervigilance, exaggerated startle response and sleep disturbance). PTSD patients may also exhibit impaired cognition, working memory, and reduced social activities. Another characteristic of PTSD is that those symptoms emerge long after traumatic events in some patients, suggesting that disease can develop after an extended asymptomatic period.
Anxiety disorders including PTSD encompass heterogeneous conditions that are associated with abnormalities in brain regions controlling fear response such as the medial prefrontal cortex, amygdala (AMY), hippocampus, and nucleus accumbens [1,2]. Changes in activity and functional connectivity of these brain regions lead to the development and maintenance of anxiety disorders [3][4][5]. AMY is a critical brain region involved in the detection of threat, fear learning and expression, and heightened emotional memory. The exaggerated activity in AMY following traumatic events is correlated with the severity of symptoms in PTSD patients. AMY is modulated by the cortical influence from medial prefrontal cortex. The altered activity in the anterior cingulate cortex (ACC), a subregion of medial prefrontal cortex, has been demonstrated in PTSD using functional MRI [6]. With extensive connection to limbic structures including AMY, ACC is thought to be the direct topdown regulator of PTSD susceptibility through modulating AMY activation to threatening stimuli [7,8]. In particular, dorsal ACC is thought to facilitate the expression of conditioned fear and increased neuronal activity in PTSD during fear conditioning [2,9]. The volume of ACC was consistently found to decrease in PTSD, however some studies have revealed hyperactivation of the ACC rather than hypoactivation [10]. Therefore, it is important to investigate the functional activity of ACC in relation to AMY hyperactivation in PTSD.
Recent genomic profiling and neuroimaging studies for human PTSD have uncovered several candidate genes for the pathogenesis, however detailed investigation on gene regulation and the underlying mechanisms are still limited [11][12][13]. Molecular and cellular encoding of traumatic events and behavioral responses can be reflected by the changes of synaptic plasticity. Specifically, signaling molecules associated with synaptic transmission and plasticity in connecting to AMY are implicated as primary substrates for fear learning and memory. Alteration in glutamatergic and GABAergic neurotransmission contributes to abnormal behavioral responses [14]. In addition, downstream signaling components, such as MAPK, cAMP and calcium mediated signaling molecules are involved in synaptic plasticity [15]. The alteration in several neuromodulatory systems such as endocannabinoid signaling, Wnt signaling and ErbB signaling can affect fear memory [16] and anxietylike behavior by synaptic remodeling in medial prefrontal cortex [17] and AMY [18,19].
Hypothalamic-pituitary-adrenal (HPA) axis plays a pivotal role in stress response upon exposure to emotional or physical stressors, and its feedback loop mediated by cortisol is altered in PTSD [20,21]. In fact, human genetic studies identified several PTSD-associated genes are involved in HPA axis and cortisol, for instance FKBP5 and corticotropin releasing hormone receptors [22]. Oxytocin and thyroid hormone signaling pathways are modulated by cortisol and the imbalance may cause emotional lability, impatience, anxiety and sleep disturbances [23,24].
As PTSD diagnostic criteria are not met at least 1 month after traumatic event according to DSM-5, sustained psychological abnormality is one of the prominent features. However, the molecular mechanisms underlying these psychological changes are not well-addressed. In this study we employed genome-wide unbiased gene analysis using RNA-seq to characterize the regulation of differentially expressed genes in a stress induced PTSD-like mouse model and found that several signaling pathways involved in synaptic remodeling and endocrine systems, as well as PTSD-associated genes were regulated in a region-and time-dependent manner.

Animals
Male, 8 to 10 weeks old C57BL/6J mice were purchased from Jackson Laboratories (Bar Harbor, ME). Animals were maintained under standard conditions of controlled temperature (23°C) and housed in a 12 h light/dark cycle (light on at 6 AM) with access to water and food ad libitum. They were group housed with 3 to 5 mice per cage and habituated to the vivarium for at least 3 days. Experimental procedures were carried out in accordance with NIH guidelines and approved by the Uniformed Services University Animal Care and Use Committee. One hundred twenty-four mice in total were used for performing various behavioral tests (Fig. 1) and molecular analyses.

Foot shock induction
Animal model for PTSD-like symptoms was constructed by foot shock induced fear conditioning and re-exposure to the shock chamber as previously described [25,26]. Briefly, trauma procedure consisted of a single exposure to an inescapable electric foot shock (Day 0) followed by re-exposure to the shock chamber every week until 4 weeks (Fig. 1). Mice received 15 s handling for 3 consecutive days before given foot shock. On the shock day, each animal was placed in the light compartment. After a 2 min adaptation period, the door to the dark shock apparatus made of plexiglass (17 cm × 17 cm × 25 cm) was open to let the animal entering the shock chamber. After another 5 min adaptation in the closed chamber, the animal received inescapable foot shock (1.5 mA × 2 s with 2 min interval for 8 times). The control group received the same procedures without stress. Animals were removed from the apparatus 1 min after the last shock.

Behavioral test
All the behavioral tests were performed during the light on phase from 8 AM to 3 PM. Mice were given one test a day to avoid the potential interference from the other tests ( Fig. 1). Animal cages were moved to a testing room at least 1 h prior to each test. After completion of the test session, the behavioral apparatus and chamber were cleansed with 70% ethanol, and then dried by hand-fan completely.

Avoidance latency/situational reminder and freezing tests
For measuring avoidance latency, animals were placed in a compartment adjacent to a foot shock chamber. After 30 s acclimation, the door was open and the time spent before entering foot shock chamber was measured within a period of 2 min. Test was terminated once mouse entered, and the mouse was immediately removed from the shock chamber to a new cage. For measuring freezing time, animals in 5 weeks post stress were placed inside chamber with the same context, and mouse behavior was monitored for 5 min with the digital video system. The data was analyzed by ANY-maze software (Stoelting Co, Wood Dale, IL) to determine freezing time and the number of episodes.

Open field test
Each mouse was placed in the center of a cubic chamber (40 cm × 40 cm × 35 cm) with a black Plexiglas wall under dim light in the testing room. The position and movement of animals were monitored for 10 min by video camera installed above the chamber. Data analyses were carried out by ANY maze software to determine the moving velocity, total travel distance, the time spent in the center zone, the entry number to the center zone. Center zone was assigned in a 20 cm × 20 cm square area in the center of arena.

Elevated plus maze test
The apparatus consisted of four arms (34 cm × 5 cm width), with two arms open and two closed by gray walls (15 cm height) arranged in the opposite side of the same type. The platform was located 50 cm above the floor of the testing room illuminated and four arms were connected in the center platform

Light/dark box test
The light-dark box was made of white and black opaque Plexiglas (40 × 20 × 30 cm light chamber, 23 × 20 × 30 cm dark chamber). The chambers were connected by a 6 × 6 cm door in the middle of the wall separating the two chambers. Animals were placed in the middle of the dark chamber with the door closed. Video monitoring was started when the door was open and continued for 10 min to measure transition number into the lit area, retention time in lit, and travel distance inside lit using ANY-maze software.

Acoustic startle test
Startle response was measured using ventilated startle chambers (SR-LAB system, San Diego, CA). Sound levels within each chamber were measured using a sound-level meter to ensure consistent presentation. Animals were placed in a Plexiglas cylinder (∅4.7 cm × 10 cm long) inside a chamber (29 cm × 29 cm × 29 cm). The cylinder resting on a platform was connected to a piezoelectric accelerometer to monitor movement. Each test session was started with a 3 min acclimation period to background white noise of 70 dB, followed by 10 times acoustic stimulus (120 dB × 20 ms) with a 1 min interval. Animals were returned to home cage immediately after test session. Data analyses were carried out by SD software (SR-LAB system) and summarized with Microsoft Excel. The latency to first peak, maximum amplitude, and average startle amplitude from 1 ms to 100 ms were analyzed.

Home-cage immobility test
Video, magnet, photobeam, or infrared sensor system has been used to measure home cage activity [27]. This system has been applied for monitoring sleep-wake behavior  Fig. 1 Timeline of PTSD-like model construction and behavior tests. C57BL/6J mice at 8-10 weeks of age were subjected to acclimation followed by electric foot shock. The control mice underwent the same acclimation procedure without foot shock. Avoidance test and situational reminder (AV/SR) were weekly performed, and the freezing test was done at 5 weeks post stress (PS). Open field test (OF) and elevated plus maze test (EPM) were performed at 2 weeks and 4 weeks PS. Light/dark box test (LDB) was performed around 4 weeks PS. Acoustic startle response test (ASR) was performed around 2 weeks and 4 weeks PS, and locomotive activity test were carried out from 4 weeks to 5 weeks PS instead of polysomnography [28,29]. Group housed animals were separated individually into a cage (20 cm × 38 cm) and nourished with water gel and food ad libitum for 3 days. One day prior to the test, the cages were placed on a test platform before 6 PM to adapt the environment. The room light was turned off at 6 PM and turned on at 6 AM the following day, in accordance with the previous light-dark cycle. Locomotive activity was monitored immediately after 6 AM with the video camera system and continued for 11.5 h. Data was analyzed numerically by ANY maze software to determine the immobile period. Locomotive activity during the entire test period was dissected with a 2 min interval to assess immobile period which showed no locomotion within each 2 min-section. Total travel distance and total immobile time were measured by ANY maze. Short immobile episode was determined as the number of immobile periods of 1 or 2 sequential immobile sections.

RNA extraction
Mice were sacrificed after deep anesthetization from 10 AM to 2 PM. Brains were immediately fresh frozen on dry ice and stored at − 80°C. Coronal section of the brain at 300 μm was cut at − 20°C. Three adjacent brain sections of bilateral amygdala (AMY) (Bregma − 0.94 mm to − 1.84 mm) and those of bilateral anterior cingulate cortex (ACC) (Bregma + 1.34 mm to + 0.44 mm) [30] were collected by punching with 1 mm-diameter brain punch based on the method described by Palkovits [31]. Total RNAs from AMY and ACC were isolated by Trizol (Invitrogen, Carlsbad, CA) and purified with the RNeasy Mini Kit (Qiagen, Germantown, MD). The RNA concentration was measured by Qubit 3.0 Fluorometer (Invitrogen).

RNA-seq
Library preparation for RNA-seq was performed by True-Seq stranded mRNA sample preparation kit (Illumina, San Diego, CA) starting from 600 ng of total RNA. Sequence data were generated by Illumina 3000 and Next-Seq platforms. We converted raw BCL files to FASTQ files and performed multiplexing using Illumina bcl2fastq2 software (version 2.17). Sequencing quality control was done by Fastqc program. We aligned the resulting FASTQ files to mm 10 mouse reference genome by STAR aligner. Quantification to read counts, FPKM and TPM was done by RSEM program. With raw read counts matrix as input, we identified genes that were differentially expressed by DESeq2 R package, using unpaired two-class significance analysis and a false discovery rate threshold of 0.05. Gene expression visualization was plotted by heatmap R package. KEGG Pathway and GO analysis were performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 [32] or "DOSE" R Bioconductor package [33].

Statistics
The GraphPad Prism 7 (GraphPad Software Inc., San Diego, CA) was used for statistical analysis. The weekly avoidance latency test was analyzed by two-way ANOVA followed by Bonferroni's multiple comparison test. The statistical comparison between the two groups was determined by Welch's unpaired two tailed t-test. All data are expressed as mean ± SEM. Statistical significance was set at p < 0.05. *, **, and *** depicting p < 0.05, p < 0.01, and p < 0.001, respectively.

Sustained avoidance and fear memory in foot shock stressed mice
For assessment of avoidance behavior, latency till moving into the same foot shock chamber from the light compartment was measured every week. The latency of stressed mice was longer than that of control mice starting from the first week and sustained for at least 4 weeks (F 4,104 = 2.782, p = 0.0305) (Fig. 2a). Notably, almost all stressed mice did not move out from light compartment to the stressed chamber except for the first trial (Fig. 2a). When placed inside the same shock chamber at 5 weeks after stress, the stressed mice possessed a longer freezing time compared to the control mice (p = 0.0038) (Fig. 2b), although the frequency of freezing episode was not significantly altered.

Anxiety-like behavior emerged at 2 weeks post stress
As anxiety is a typical behavioral feature in PTSD, we examined if stressed mice acquired and sustained anxiety using several behavioral tests. Open field test was performed in a black walled squire space for 10 min. At 2 weeks PS, the stressed mice showed decreased retention time in the center zone and entry number to the center zone (p = 0.0056) (Fig. 3a, b). Total travel distance was significantly shorter in the stressed mice than control mice (p = 0.0006) (Fig. 3c). At 4 weeks PS, the stressed mice showed significant reduction in retention time in the center zone (Fig. 3d, p = 0.0126), entry number to the center zone (Fig. 3e, p = 0.0304), and total travel distance (Fig. 3f, p = 0.0047). Shorter travel distance was unlikely caused by motor deficit of the stressed mice, since maximum running speed was not impaired (Fig. 3g). In the elevated plus maze test, the stressed mice spent less time in open arms than control mice at 2 weeks PS, despite the significant difference was not reached (Fig. 4a). The entry number to open arms was significantly lower in the stressed mice (p = 0.0095) (Fig. 4b). When examined at 4 week PS, there were no significant differences in both time spent and entry number in open arms between the stressed mice and control mice ( Fig. 4c and d). To further examine the anxiogenic behavior at these time points, we also performed a light/dark box test. It was shown that the number of transitions to light area and travel distance in light area during 10 min test period were significantly reduced in the stressed group than the control group ( Fig. 5a and c). Although time spent in light area was not significantly different (Fig. 5b), the stressed mice preferred dark protected area and less explored in an unprotected area compare to control mice at 4 weeks PS. Overall, these tests demonstrated that the stressed animals started to show anxiety-like behavior at 2 weeks PS and continued till 4 weeks PS.

Hyperarousal symptoms in stressed mice
Acoustic startle response against 120 dB stimulus repeated for ten times were monitored by a motion sensor system, and the maximum peak amplitude, latency to the first peak, and average amplitude were measured. As shown in Fig. 6a, both control and stressed mice showed a similar pattern to startle response. The average amplitude of the maximum peak out of 10 trials was not different between control and stressed mice at 2 weeks PS. The average amplitude from 1 ms to 100 ms was also similar despite a slight increase observed in the stressed mice (Fig. 6c). However, at 4 weeks PS the peak amplitude was significantly higher in the stressed mice (p = 0.0048) (Fig. 6d), and the average amplitude was also greater in stressed mice than control mice (p = 0.0440) (Fig. 6e). There was no difference in startle response latency between the stressed and control mice at either test period. These results suggested that startle response was exaggerated in stressed animals at relatively later time point.

Increased frequency of short immobile episodes during the sleep period
Sleep disturbance is a common feature and impacts the pathophysiology of PTSD and other psychiatric diseases [34]. We characterized the sleep-wake behavior of stressed mice with long-term video monitoring in their home cages. There was no significant difference between control and stressed mice in total immobile time during the entire inactive phase (Fig. 7a). However, the number of short immobile episodes was significantly larger in the stressed mice than the control mice (p = 0.0353) (Fig. 7b), suggesting that the sleep/wake cycle was disturbed in stressed mice. The total travel distance showed no significant difference between the control and stressed mice, indicating the locomotive activity was not impaired in stressed animals (Fig. 7c).

Differentially expressed genes in ACC and AMY
We performed genome wide RNA expression analysis in the anterior cingulate cortex (ACC), which is involved in the modulation of emotional and fear encoding and consolidation. Of the 3233 differentially expressed genes (DEGs) in ACC at 2 weeks post stress (PS), 1335 were upregulated, and 1898 were downregulated (Additional file 1: Figure S1A). However, at 5 weeks PS the numbers plunged to 43 (Additional file 1: Figure S1B) and only 16 of them were overlapped between these time points (Additional file 1: Figure S1C). In AMY 2 weeks PS, there were 1063 DEGs, in which 873 genes were predominantly down-regulated (Additional file 2: Figure  S2A). Different from ACC, there were a substantial number of DEGs in 5 weeks PS (Additional file 2: Figure  S2B). Of more than 1000 DEGs, approximately 9% (94 DEGs) were overlapped between 2 weeks and 5 weeks PS as shown in the Venn diagram (Additional file 2: Figure S2C).

Regulation of PTSD associated genes in the mouse model
Several human genomic studies have been conducted recently to identify candidate genes for PTSD patients [35][36][37][38]. We found that ten genes associated with PTSD were regulated in either AMY or ACC ( Table 1). Several of these genes were involved in HPA axis and glucocorticoid regulation, which include Crhr1 and 2, Fkbp5, Rora, Sgk1, and Stat5b. Other genes including Adcy8, Cacna1c, Klhl1, and Shank1 were involved in synapse plasticity and regulation. Gene expression regulation in our animal model was correlated, at least to some extent, with that identified in different brain areas or blood samples from patients with PTSD.

Differentially expressed genes in AMY
Those DEGs in AMY were applied for Gene Ontology (GO) analysis using DAVID to determine enriched functional annotation in terms of molecular functions (Additional file 3: Table S1) or biological processes (Additional file 4: Table S2). First, to characterize the regulation of global mRNA expression, we searched for the GO terms related to transcriptional regulation and found comparable numbers of DEGs in the Transcription and Negative regulation of transcription between 2 weeks and 5 weeks PS, whereas the Positive regulation of transcription was enriched only in 5 weeks PS (Additional file 5: Figure S3). Of note, the Covalent chromatin modification was enriched only in 2 weeks PS. Consistently, there were four more enrichments related to the chromatin modifications: unmethylated CpG binding, histone deacetylase binding, histone H3-K4 methylation and histone H4-K5 acetylation. There were several gene regulations in terms of cellular alteration in the central nervous system (CNS) (Additional file 6: Figure S4). Nervous system development was enriched both in 2 weeks and 5 weeks PS, while neuron projection development was only seen in 2 weeks PS. Neurogenesis and Neural precursor cell proliferation were enriched in 5 weeks PS. GOs related to neuron differentiation were enriched both in 2 weeks PS and in 5 weeks PS. For non-neuronal cell differentiation, Oligodendrocyte differentiation and Epithelial cell differentiation were enriched in 2 and 5 weeks PS, respectively.

Gene regulation in AMY
A number of KEGG pathways were involved in synaptic remodeling (Fig. 8) from the enriched pathway analyses (Additional file 7: Table S3). Pathways for neurotransmitters including acetylcholine, dopamine, GABA, glutamate, and serotonin were enriched in 5 weeks PS. We found gene regulation of glutamatergic kainite receptor subunits, Grik3, 4 and 5, metabotropic receptor subunits, Grm4 and   In GABAergic synapse, type A GABAergic receptor subunits (Gabrg1, Gabrg3, and Gabrq) and GABA transporter (Slc38a1) were regulated after 5 weeks. Moreover, acetylcholine esterase, Ache, and muscarinic cholinergic receptor, Chrm3, dopamine receptor1, Drd1, and serotonergic receptors, Htr1b, 1d, 2c, and 7 were regulated. In terms of synapse mediated signaling, LTP and LTD pathways and several downstream signaling pathways crucial for synapse plasticity such as MAPK, Ras, and PI3K/ Akt were enriched. Ras signaling pathways including Mras, Kras, as well as Ras GTPase regulatory factors, and Rap1 signaling pathways including Rap1 GTPase regulatory factors were overlapped in both 2 weeks and 5 weeks PS. And the key molecules in calcium mediated signaling, calmodulin (Calm2), CAM kinase (Camk2a), and CREB binding protein (Crebbp) which are involved in LTP and LTD were also altered. The key components of Wnt signaling, four receptor genes and three ligands including β-catenin, were regulated. Notably, the retrograde endocannabinoid signaling was highly regulated at 5 weeks PS. Moreover, pathways related to synaptic formation and activity, including Adherens junction, Axon guidance, Focal adhesion, Endocytosis, and Regulation of actin cytoskeleton pathways were enriched.
To further characterize the gene regulation of synapse remodeling, the pathways overlapped in 2 weeks and 5 weeks PSs were characterized by individual DEGs. Most of DEGs were regulated either after 2 weeks or 5 weeks PS, but there were a few DEGs overlapped in both time periods, particularly involved in the Endocytosis pathway ( Table 2). These results suggest that individual DEG was differentially regulated at 2 weeks and 5 weeks PS. Consistent with the differential regulation of synapse remodeling genes in AMY, DEGs related to transcription showed only a few overlapped genes between 2 and 5 weeks PS, suggesting a diverse transcriptional regulation.

Gene regulation in ACC
There were substantial number of GOs and pathways enriched in 2 weeks PS. As seen in AMY, we found several GOs and KEGG pathways related to hormone and neuroendocrine system and synapse remodeling. Figure 10 shows enriched KEGG pathways including Estrogen signaling, GnRH signaling, Insulin signaling, Oxytocin signaling, Prolactin signaling, and Thyroid hormone signaling. Substantial number of KEGG pathways were also involved in synapse remodeling. MAPK signaling pathway was the only pathway seen in 5 weeks PS among them. Thus, synapse remodeling and hormone and neuroendocrine system were the mainly regulated pathways in 2 weeks PS, and they were mostly terminated in 5 weeks PS.

Gene regulation related to PTSD-like behavior in AMY and ACC
GOs associated with PTSD-like behavior such as fear response, startle response, social behavior were predominantly seen after 5 weeks in AMY and 2 weeks in ACC (Fig. 11)

Discussion
Behavioral characterization of the foot shock induced stress model Electric foot shock has been widely used as fear conditioning, because it can provide consistent stress to animals and is easy to assess fear memory. Severe foot shock has also been used as a stressor for PTSD-like animal model [39]. In previous reports, animals receiving electric foot shock was subsequently re-exposed to the same contextual environment [40]. This procedure which is called situational reminder seems critical to acquire long-term PTSD-like abnormal behaviors including long-lasting hyperarousal [26], sustained freezing behavior [25] as well as higher evoked field potential in basolateral AMY [41]. In this study, animals were given the weekly avoidance test and we observed that, except for the first week, all the stressed animals avoid to enter the previously foot-shock chamber, suggesting these stressed animals were threatened under this situation. Together with freezing test at 5 weeks PS, we found that the stressed animals acquired sustained fear memory for more than 1 month. Since anxiety-like behavior is a prominent characteristic of PTSD, we carried out three different behavioral tests (Figs. 3, 4, and 5) and found that the stressed mice acquired anxiety-like behavior, which emerged after 2 weeks and sustained until 4 weeks PS. In contrast, exaggerated acoustic startling response of the stressed mice was only acquired after 4 weeks PS (Fig. 6). Thus, two abnormal behavioral traits in our model were differentially evolved. Moreover, short  immobile episodes were more frequent in stressed mice than control animals, suggesting a sleep disturbance in the stressed animals and a reminiscence of the abnormal sleep-wake cycle in PTSD patients. Differential susceptibility of animals to stress has been recognized as an important trait of the animal model because of the existence of human resilience or vulnerability to trauma [42,43]. In fact, we observed a wider sample variance in stressed group than in control group in the acoustic startle test at 4 weeks PS (Fig. 6c and e), but not at 2 weeks PS ( Fig. 6b and d), suggesting that differential susceptibility was created when stress induced hyperarousal was prominent. Cohen's group has evaluated individual rodent susceptibility based on the acoustic startle response [44]. Currently it is unknown whether the wide variation was derived from individual difference occurred during the development of PTSD-like abnormality, since there is lack of consistency for the stressed animals to behave as vulnerable or resilient in different behavioral assays.

Gene regulation associated with PTSD
Genome-wide association studies of human PTSD have identified the disease specific polymorphism of several associated genes in their gene loci. Some of the associated genes are involved in HPA axis and glucocorticoids to regulate stress response and emotional control [13]. In our model, several genes were regulated in either AMY or ACC, including Crhr1 and 2 [45][46][47], Rora [48], Sgk1 [49], and Stat5b [50] (Table 1). Among them, FKBP5 has been intensively investigated since it is a co-chaperone for glucocorticoid (GC) receptor to suppress its transcriptional activity [51]. In addition, we found genes involved in synaptic plasticity: Adcy8, whose deficiency led to stress-induced anxiety [52] implying a role in long-term potentiation and synaptic plasticity for fear related memory [38], Cacna1c (Cav1.2 subunit of L-type calcium channel), which mediates long-term synaptic plasticity through calcium signaling, Klhl1, the actin binding protein modulating calcium current [53] and neurite outgrowth [54], and Shank1, the anchoring protein playing a role in glutamatergic synaptic plasticity [55]. These results suggest that those PTSD-associated genes can play an important role on the disease development.

Time-dependent gene regulation in AMY
Numerical gene regulation was comparable between 2 weeks and 5 weeks PS since there were more than 1000 DEGs in each condition (Additional file 2: Figure S2). However, DEGs in the overlapped pathways were predominantly in either 2 weeks or 5 weeks PS, but only a few DEGs overlapped, consistent with the transcription related DEGs ( Table 2). Alteration of drastic transcriptional regulatory genes resulted in dynamic changes of global gene expression including synapse remodeling or hormone/neuroendocrine system. In addition, histone mediated epigenetic regulation may contribute to time-dependent regulation (Additional file 5: Figure S3), as recently postulated that epigenetic alteration is deeply involved in the pathogenesis as well as in the disease susceptibility [43,56,57]. Of note, the glutamatergic, GABAergic and endocannabinoid signaling were all highly enriched in 5 weeks PS (Fig. 8), consistent with the notion that the excitatory and inhibitory neurotransmission is tightly regulated by the endocannabinoid signaling under stressful conditions [58]. In addition, cholinergic and serotonergic input can modulate neurotransmission at 5 weeks PS. Ras, Rap1, and MAPK pathways which play an important role in synaptic transmission, including fear memory encoding and expression [59,60] and calcium mediated signaling, including LTP were also predominantly regulated (Fig. 8). Moreover, hormone and neuroendocrine system were also predominantly enriched in 5 weeks PS (Fig. 9). However, we did not find significantly regulated expression of the core HPA axis genes including corticotropin releasing hormone, its receptor, or GC receptor, which was different from a previous study [61]. It is possible that the differences might be derived from the different animal species or experimental procedures. PTSD-like model constructed by foot shock followed by situational reminder may not clearly provoke HPA axis activation in male animal as reported previously [62]. Instead, we found that oxytocin signaling, which is known to exert anxiolytic and antidepressant effect [63,64], and also as an important negative modulator of HPA axis and GC signaling [65], was highly enriched in 5 weeks PS (Fig.  9). In fact, oxytocin was upregulated in 2 weeks PS, then plummeted in 5 weeks PS. Oxytocin receptor gene was also upregulated in 2 weeks PS, while returned to basal level in 5 weeks PS. Thus, there is a possibility that extra-hypothalamic GC signaling is modulated by regulatory neuroendocrine system such as oxytocin, although further investigation is necessary to clarify the regulation of HPA axis. Moreover, several GOs that are related to the PTSD-like behavioral traits such as startle response, response to drug, and sensory perception to pain were enriched at 5 weeks PS (Fig. 11), suggesting that those DEGs play a key role in the molecular mechanism to disease development. Taken together, gene regulation involved in synapse remodeling and hormone/neuroendocrine system in AMY is likely culminated in 5 weeks PS, which might be correlated to the development of PTSDlike abnormal behaviors.

Termination of gene regulation after 5 weeks in ACC
ACC is reported to alter the neuronal activity in PTSD patients during fear conditioning test as well as in the resting condition [66]. ACC has extensive afferent connections to emotion regulatory limbic regions including