Unforgettable: optogenetic stimulation of prelimbic pyramidal neurons maintains fear memories by modulating amygdala pyramidal neuron transcriptome

Fear extinction requires coordinated neural activity within the amygdala and medial prefrontal cortex (mPFC). Any behavior has a transcriptomic signature that is modified by environmental experiences, and specific genes are involved in functional plasticity and synaptic wiring during fear extinction. Here, we investigated the effects of optogenetic manipulations of prelimbic (PrL) pyramidal neurons on amygdala gene expression to analyze the specific transcriptional pathways involved in adaptive and maladaptive fear extinction. To this aim, transgenic mice were (or not) fear-conditioned and during the extinction phase they received optogenetic (or sham) stimulations over PrL pyramidal neurons. At the end of behavioral testing, electrophysiological (neural cellular excitability and Excitatory Post-Synaptic Currents) and morphological (spinogenesis) correlates were evaluated in the PrL pyramidal neurons. Furthermore, transcriptomic cell-specific RNA-analyses (differential gene expression profiling and functional enrichment analyses) were performed in amygdala pyramidal neurons. Our results show that the optogenetic activation of PrL pyramidal neurons in fear-conditioned mice induces fear extinction deficits, reflected in an increase of cellular excitability, excitatory neurotransmission, and spinogenesis of PrL pyramidal neurons, and in strong modifications of the transcriptome of amygdala pyramidal neurons. Understanding the electrophysiological, morphological and transcriptomic architecture of fear extinction may facilitate the comprehension of fear-related disorders.


Introduction
In variable and challenging environments with various contextual situations, the individuals face a number of approaching dangers. The knowledge of potential threats allows developing fear of threatening situations, choosing among the various behaviors the safest ones, and detecting future dangers. To develop adaptive fear responses, the brain has to discriminate different sensory cues and associate relevant stimuli with aversive events [1]. Thus, when a relevant stimulus (or a context) is associated with an aversive event, the fear associative learning and memory take place.
Learned fear has been widely studied using Classical Fear Conditioning (CFC), a very useful paradigm to analyze the neuronal and molecular bases of fear associative learning and memory [2,3]. In experimental models in which the CFC paradigm is implemented, the conditioned stimulus (CS), such as a specific cue or context, is associated with the unconditioned stimulus (US), such as foot-shock [3,4]. After the association has taken place, CS alone is able to induce the conditioned response (CR) of fear, such as freezing behavior. Interestingly, the CR to CS is gradually weakened by repeated exposure to unreinforced (presented without US) CS. In fact, the extinction of fear memory when the danger is gone is so crucial to implement other survival functions that the impairment of the coping mechanisms dampening fear memory results in maladaptive behaviors. Extinguishing a fear response does not simply involve the fading away of the previous learning. Rather, during extinction process the subject learns something newthe cue no longer predicts that the fearful event will occur [5]. Compromised fear extinction is the key clinical feature of several fear-related disorders, such as post-traumatic stress disorders (PTSD), generalized anxiety, major depression, and phobias [5]. The acquisition of CS-US associative memory and the acquisition and maintenance of extinction memory require coordinated neural activity within the amygdala, medial prefrontal cortex (mPFC), and hippocampus [6][7][8][9][10][11]. Specifically, the activation of pyramidal neurons of the baso-lateral amygdala (BLA) is necessary to associate sensory input with fear CR [12,13]. Interacting with the amygdala, the mPFC combines information from multiple inputs to exert top-down control allowing for appropriate responses. In this framework, it has been demonstrated that the balance between expression and extinction of fear CR is modulated by inputs from/to the amygdala to/from two sub-regions of the mPFC: the prelimbic cortex (PrL), which promotes fear responses, and the infralimbic cortex (IL), which promotes fear extinction [14,15].
To test the activation or deactivation of the BLA-mPFC network in transmitting learned fear CS-US association, in vivo optogenetic manipulations were performed by stimulating or inhibiting, with advanced spatial and temporal precision, specific neurons in the brain areas involved in fear processing [14,[16][17][18][19][20]. Fear extinction decreases the efficacy of excitatory synaptic transmission in the projections from mPFC to BLA, whereas inhibitory responses are not altered [21]. In parallel, projections from BLA to PrL exhibit cell-type-specific plasticity during states of high fear, whereas projections from BLA to IL are recruited during states of low fear [19]. Thus, in vivo optogenetics allows understanding the causal relationships between neuronal firing and behavior.
It has to be noted that any behavior has a specific genomic, transcriptomic, and epigenetic signature. Namely, transcriptomic changes are a crucial component of the neuronal modifications that underlie learning and memory [22][23][24], also after the fear exposure [25][26][27].
To date, specific genes involved in functional plasticity and synaptic wiring during fear memory are retained possible disease contributors and potential therapeutic targets for fear-related disorders [28][29][30][31][32]. The environmental experiences may modify the transcriptome [33][34][35][36], and these modifications, in turn, may explain the variability in resilience or predisposition to fearrelated disorders as well as the severity of their symptomatology [37]. Understanding the transcriptomic architecture of compromised fear inhibition may thus facilitate the comprehension of fear-related disorders and the identification of potential therapeutic targets.
In the present study, adult mice were submitted to CFC receiving (or not) the US, to test the synergy between the activity of the amygdala and mPFC in fear learning. In vivo optogenetic manipulations were delivered to maintain the US-related activation of the PrL pyramidal neurons and impair fear extinction. At the end of behavioral testing, electrophysiological, morphological, and transcriptomic correlates were evaluated. Namely, in PrL pyramidal neurons neural cellular excitability and excitatory neurotransmission as well as spinogenesis were evaluated. In parallel, amygdala pyramidal neurons were sorted to perform cell-specific RNA-analyses (differential gene expression profiling and functional enrichment analysis). These analyses permitted us to drive a genome-wide investigation of pyramidal neuron expression patterns without a priori selection of specific gene factors, thus preventing bias in gene expression brought by bulk tissue analysis and biological subject pooling. Notably, investigating the effects of mPFC optogenetic manipulations on amygdala pyramidal neurons gene expression revealed specific transcriptional pathways involved in impaired fear extinction.

Behavioral results: in vivo optogenetics of the PrL pyramidal neurons during contextual FC
The animals were (or not) fear-conditioned by using the CFC with repetitive sessions on day 1 (Conditioning phase), and 2, 3, 4, 7, and 14 (Extinction phase) ( Figure 1A). During the extinction phase, the mice received optogenetic (OPTO FEAR and OPTO NOT FEAR groups) or sham (SHAM FEAR and SHAM NOT FEAR groups) stimulations. A three-way ANOVA (stimulation x fear x day) on freezing behavior (measured during 0-3 min of CFC) ( Figure 1B showed consolidation processes of aversive memories on day 2, as revealed by their significantly increased freezing times between the Conditioning phase and day 2 (P=0.00004 for both OPTO FEAR and SHAM FEAR groups), and by the similar freezing times of the two groups (P=0.90). As expected, SHAM FEAR animals progressively extinguished fear memories over time and on day 14 their freezing behavior returned to a level similar (P=0.17) to that showed during the Conditioning phase. Interestingly, no extinction of fear memories was observed in OPTO FEAR group. On day 14, OPTO FEAR mice still displayed freezing times significantly longer than those of the Conditioning phase (P=0.00002). One-way ANOVA on freezing behavior (measured during 0-6 min of day 14) of all groups revealed a significant group effect (F3,36=11.05; P=0.00003) ( Figure 1C). As revealed by Newman-Keuls post-hoc comparisons, OPTO FEAR group showed the highest freezing times in comparison to the remaining groups (at least P=0.0005). (Conditioning phase) of CFC, each Thy1-COP4 mouse was allowed to explore the conditioning chamber for 3 min (Baseline). Afterward, half of the entire sample received three foot-shocks.
On days 2, 3, 4, 7, and 14 (Extinction phase), the fear-conditioned and not fear-conditioned mice were placed again in the conditioning chamber for 6 min. During the Extinction phase, no shock was delivered and the mice received three optogenetic or sham stimulations of PrL pyramidal neurons. B) Freezing behavior measured during 0-3 min of CFC. All animals showed similar responses in the Conditioning phase and only the fear-conditioned animals (OPTO FEAR and SHAM FEAR groups) showed increased freezing times on day 2. While SHAM FEAR group progressively extinguished fear memories over time, no extinction of fear memories was observed in OPTO FEAR group. C) Freezing times measured during 0-6 min of day 14. OPTO FEAR group showed the highest freezing times in comparison to the remaining groups (*** P=0.0005).

Electrophysiological results: Cellular excitability of PrL pyramidal neurons
To test the ability of the optogenetics to modulate the activity of PrL cortex in fear-conditioned or not fear-conditioned animals, an extensive characterization of the cellular excitability of pyramidal layer 5 neurons after sham or optogenetic stimulation was performed.
In OPTO FEAR group, the optogenetic stimulation induced a robust increase in the evoked firing activities triggered by growing pulses of depolarizing current (0 to 400 pA Conversely, in both groups of not fear-conditioned animals (SHAM NOT FEAR and OPTO NOT FEAR groups) the optogenetic stimulation did not produce any significant difference in cellular excitability. Neither the evoked firing nor the rheobase nor EPSC frequency showed any significant difference between SHAM NOT FEAR and OPTO NOT FEAR groups ( Figure 2D).
Taken together these data demonstrate that optogenetic stimulation was able to modify the intrinsic cellular excitability of PrL pyramidal neurons supporting the assumption that this area is involved in modulation of the fear responses during extinction phase.      The top twenty significant KEGG terms resulted associated to signaling processes related to several neurodegenerative diseases and metabolic pathways ( Figure 5D) (namely: mmu05022 Pathways of neurodegeneration -multiple diseases, mmu04144 Endocytosis, mmu05010 Alzheimer Disease, mmu05012 Parkinson Disease, mmu04140 Autophagy -animal, mmu05014 Amyotrophic Lateral Sclerosis, mmu04137 Mitophagy -animal, mmu05132 Salmonella infection, mmu04714 Thermogenesis, mmu04360 Axon guidance, mmu04919 Thyroid hormone signaling pathway, mmu04141 Protein processing in endoplasmic reticulum, mmu05231 Choline metabolism in cancer, mmu04722 Neurotrophin signaling pathway, mmu03010 Ribosome, mmu04072 Phospholipase D signaling pathway, mmu01212 Fatty acid metabolism, mmu01200 Carbon metabolism, mmu05020 Prion disease, and mmu00190 Oxidative phosphorylation).

Comparison between OPTO NOT FEAR vs. SHAM NOT FEAR groups
Differential Expression analysis showed 1661 significant DEGs (308 up; 1353 down, with reference level set on the SHAM NOT FEAR condition) (Figure 5E-F; Table 2). The ORA performed on the obtained DEGs resulted only in 36 significantly enriched KEGG terms and 0 significantly enriched GO terms. Here, the genes universe did not allow enriching many GO terms because of sampling bias correction [38]. Again, the top twenty over-represented KEGG terms were related to neurodegenerative diseases and metabolic pathways ( Figure 5G) (namely: mmu04144 Endocytosis, mmu00190 Oxidative phosphorylation, mmu04140 Autophagy -animal, mmu05014 Amyotrophic Lateral Sclerosis, mmu04714 Thermogenesis, mmu05010 Alzheimer Disease, mmu05012 Parkinson Disease, mmu05022 Pathways of neurodegeneration -multiple diseases, mmu05016 Huntington Disease, mmu01212 Fatty acid metabolism, mmu05020 Prion Disease, mmu04141 Protein processing in endoplasmic reticulum, mmu03040 Spliceosome, mmu04530 Tight junction, mmu04137 Mitophagyanimal, mmu01040 Biosynthesis of unsaturated fatty acids, mmu04932 Non-alcoholic fatty liver disease, mmu04142 Lysosome, mmu01200 Carbon metabolism, and mmu04071 Sphingolipid signaling pathway).  Component, and Molecular Function) were sorted by q-value before plotting them together.

Comparison between OPTO FEAR vs. OPTO NOT FEAR groups
Differential Expression analysis showed 3506 significant DEGs (2104 up; 1402 down, with reference level set on the OPTO NOT FEAR condition) ( Figure 6A-B, Table 3  Protein processing in endoplasmic reticulum, mmu05017 Spinocerebellar ataxia, mmu05211 Renal cell carcinoma, mmu04928 Parathyroid hormone synthesis, secretion and action, mmu01212 Fatty acid metabolism, and mmu04714 Thermogenesis).

Comparison between SHAM FEAR vs. SHAM NOT FEAR groups
In line with PCA highlighting poor segregation between SHAM FEAR vs. SHAM NOT FEAR groups, differential Expression analysis showed only 107 significant DEGs (81 up; 26 down, with reference level set on SHAM NOT FEAR condition) ( Figure 6D-E, Table 4). KEGG analysis showed modulation of the pathways associated to cAMP signaling (mmu04024), Prostate cancer (mmu05215), and Thyroid hormone signaling (mmu04919).  Component, and Molecular Function) were sorted by q-value before plotting them together.

DEGs involved in learning/memory and fear response in the comparison between OPTO FEAR vs. SHAM FEAR groups
Despite the relative pathways were not identified as significantly enriched, we wondered whether genes involved in learning/memory and fear response were significantly modulated in the case of impaired fear extinction caused by the optogenetic stimulation on PrL pyramidal neurons. Finally, we also noted the strong down regulation of Thy1 (log2FC=-2.77) gene. Numbers into brackets indicate the number of DEGs related to each indicated BP over the gene universe. Light blue dots indicate genes described in literature as associated to fear extinction.
When DEGs are part of the top twenty DEGs list and/or part of the considered BP the gene name is indicated.

Discussion
Starting from the assumption that fear learning and extinction are adaptive processes caused Optogenetic experiments have confirmed that modified synaptic transmission in the amygdala-PrL network during fear learning interferes with long-term fear consolidation [14]. By using ex vivo electrophysiology, combined with optogenetic techniques, it has been demonstrated that fear extinction decreases the efficacy of glutamatergic transmission from mPFC to the amygdala, whereas inhibitory responses are not altered [21]. In parallel, it has been reported that the amygdala neurons projecting to the PrL cortex are active during states of high-fear, whereas those projecting to IL cortex are recruited and exhibit cell-type-specific plasticity during fear extinction [19].
Cell-type-specific transcriptome analysis represents cutting-edge tool to reveal targets useful for understanding and treating fear-related disorders. Differential gene expression and coexpression network analyses identified diverse networks activated or inhibited by fear learning vs. extinction, and upstream regulator analysis and viral vector-manipulations demonstrated that fear extinction is associated with reduced cAMP/Ca 2+ responsive element binding (CREB) expression [39].
Our study identifies cell-type amygdala pyramidal neuron-specific gene networks (by means of the RNA-sequencing) disclosing pathways mediating adaptive or maladaptive fear extinction and opening innovative possibilities to understand deeper the underlying mechanisms of fear process and its impairment. Here we discuss some of the genes and pathways that result more relevant.
Optogenetic stimulation of the PrL pyramidal neurons resulted in strong modifications of the amygdala pyramidal neuron transcriptome ( Figure 5A-B, 3), which are associated with synaptic transmission. Namely, shunting inhibition via GABAA receptors reduces activation of N-methyl-D-aspartate receptors, and impairs long-term potentiation [40], as well as voltage-gated potassium channels control cellular excitability by regulating a variety of neuronal properties, such as inter-spike membrane potential, action potential waveform, and firing frequency [41]. Within the brain, Kcnh3 is expressed in the cerebral cortex, amygdala, hippocampus, and striatal regions, with specific expression in pyramidal neurons [42].
Other DEGs that we identified as associated with impaired extinction were linked to inflammation processes, such as Csmd1 [43] and Bcl2 [44]. Notably, overexpression of Bcl2 blocks the apoptotic death of the pro-B-lymphocyte cells and neurons [45,46]. Furthermore, in the same comparison OPTO FEAR vs. SHAM FEAR, the top twenty GO terms associated with impaired extinction highlighted the differential involvement of pathways associated with neuronal plasticity and glutamatergic synaptic signaling, resulting in modulation of overall preand post-synaptic organization. In parallel, the top twenty KEGG terms were associated with processes related to several neurodegenerative diseases, metabolic pathways, production of lipids and proteins, and thyroid hormone and neurotrophin signaling. Remarkably, fearassociated enrichments have been related with dendritic and post-synaptic processes, while extinction-associated enrichments have been related with cellular metabolism and proliferation [39].
By looking at genes related to Learning/memory or Fear response processes, several DEGs were scored in the comparison between OPTO FEAR vs. SHAM FEAR groups (Figure 7). Among these, The photo-stimulation of PrL pyramidal neurons in the presence or absence of fear differently modulated amygdala pyramidal neuron transcriptome ( Figure 6A-B, Table 3). Among the top twenty DEGs, some relevant genes resulted down regulated in OPTO FEAR mice when compared with OPTO NOT FEAR mice. Among these genes, there was Ngef (Neuronal guanine nucleotide exchange factor, also known as Ephexin1), an ephrin (Eph) receptor-interacting exchange protein that promotes EphA4 binding and leads to cell morphology changes [53] and reduces spine density [54]. Accordingly, the OPTO FEAR mice (in which Ngef is down regulated) exhibited increased number and density of dendritic spines in the apical arborizations.
Furthermore, the binding of Eph with its receptors constitutes a molecular link between Eph receptors and actin cytoskeleton and modulates pre-synaptic calcium channel activity [55]. To Accordingly, long-term fear memory is impaired in mice lacking EphB2 forward signaling [56].
Among top twenty DEGs, another gene associated with actin-binding protein was Enc1 (ectodermal-neural cortex 1) that resulted down regulated in OPTO FEAR mice. Kim et al. (1998) showed that expression of Enc1 induces neuronal process formation, whereas antisense treatment inhibits neurite development [57].  Table 2). Among DEGs (down regulated in OPTO NOT FEAR mice when compared with SHAM NOT FEAR mice), some genes were implicated in cell proliferation, synaptic activation, and long term potentiation (such as Nab1 and Sos1), others in inflammation processes (such as Gpx4 and Nkap), with roles in transcriptional repression and RNA splicing and processing (Nkap) [61]. Furthermore, in the same comparison the gene universe did not allow enriching many GO terms because of sampling bias correction [38]. Again, the top twenty over-represented KEGG terms were related to neurodegenerative diseases, metabolic pathways, and biosynthesis of protein and unsaturated fatty acids.
Finally, we observed that fear conditioning per se, without optogenetic stimulation, did not greatly impact gene expression in amygdala pyramidal neurons ( Figure 6D-E, separation between the transcriptomes of fear-and not fear-conditioned animals [39]. The authors discussed their findings as result of stress induced translational changes due to the handling of animals and CS exposure. This hypothesis tested in a separate cohort of mice was confirmed, demonstrating that stress-related genes were similarly regulated in both groups.
Thus, it is likely that the signature of associative fear learning was obscured by generalized stress-related changes. However, the limited segregation between SHAM FEAR and SHAM NOT FEAR gene profiles might also be caused even by the time point of the transcriptomic analyses.
In fact, gene profiling of the SHAM FEAR group was performed once the animals had extinguished the fear memories, and fear response was over, although the fear engram was likely still stored in the amygdala-hippocampus-mPFC circuit [62].
Actually, among DEGs present in the comparison between SHAM FEAR vs. SHAM NOT FEAR groups, some genes (such as Slx4ip and Kif6) were up regulated in SHAM FEAR group and are implicated in DNA/RNA regulation, maintenance and repair [63,64]. Conversely, the Ino80 gene resulted down regulated in the same comparison. Notably, it has been reported that the deletion of Ino80 results in defective cellular proliferation and premature entry into cellular senescence, due to activation of the DNA damage response [65]. Furthermore, other genes associated to metabolism (such as the down regulated Smcr8 and the up regulated Clstn3), neurotransmission (such as the up regulated Hcn2), and inflammation (such as the up regulated Sumf1 and Plaa) were found differentially expressed in SHAM FEAR group when compared with SHAM NOT FEAR group. Interestingly, also Adcyap1r1 (adenylate cyclase activating polypeptide 1 receptor 1) resulted up regulated in SHAM FEAR group. The pituitary adenylate cyclaseactivating polypeptide is a hormone that stimulates the secretion of growth hormone, adrenocorticotrophic hormone, catecholamines, and insulin, by its interaction with specific receptors. Interestingly, the methylation of Adcyap1r1 in peripheral blood has been associated with PTSD, and Adcyap1r1 mRNA is induced by fear learning [66]. Since DNA methylation of regulatory elements usually acts to repress gene transcription, the findings by Ressler et at.
(2011) indicating that methylation of Adcyap1r1 is associated with impaired fear extinction (as typically occurring in PTSD) together with the present ones indicating that the up regulation of this gene is associated with efficient fear extinction converge in highlighting that the system pituitary adenylate cyclase-activating polypeptide with its receptors might be an important mediator of abnormal fear responses following trauma [66].
In the same comparison (SHAM FEAR vs. SHAM NOT FEAR), the top twenty GO terms were related to morphogenesis and synaptic plasticity, with the important specificity of the excitatory synapse. KEGG analysis showed modulation of the pathways associated to cAMP signaling, cancer, and thyroid hormone signaling.
Overall specifically the inhibitory GABAergic signaling, as well as by differential involvement of pathways associated with neuronal plasticity and glutamatergic signaling.
Our findings demonstrate that impaired extinction is featured by specific changes of transcriptome that validate previous findings [39], and provide targets for future translational research into cell type-specific control of fear extinction, emphasizing the key role of pyramidal neurons belonging to amygdala-mPFC fear matrix. This type of comprehensive cell-type specific analysis produces an important array of targets potentially useful for diagnosis, treatment, and prevention of fear-related disorders.

Subjects
Male adult (2. The animals assigned to the same experimental group were never siblings.

Experimental procedure
Thy1-COP4 mice were unilaterally implanted with a guide cannula on PrL sub-region of right   Franklin and Paxinos (1997) and Mouse Brain Atlases (The Mouse Brain Library; www.nervenet.org). Ferrule-terminated implanted optical fibers were secured to the skull using dental acrylic.

Stereotaxic surgery and fiber optic implantations
Mice were allowed to recover from surgery for 1 week before behavioral testing. During the recovery period, they were habituated to handling and connection of the optic fiber with the optogenetic ferrule. Locations of implanted optical fibers were validated using histology in all experimental mice.

CFC and in vivo optogenetic stimulations of the PrL pyramidal neurons
As previously described [67,68], the CFC was carried out in a soundproof conditioning chamber (50 cm long, 24.5 cm wide, 26.5 cm high) (Ugo Basile, Varese, Italy) made of gray Perspex with a metal grid floor. A video camera placed above the conditioning chamber allowed observing animal behavior. Before the behavioral testing, the mice were handled to connect the optic fiber with the optogenetic ferrule.
As depicted in Figure 1A Freezing was recorded by an experimenter blind to the group the animal belonged to and freezing times during the first 3 min for Conditioning (Baseline) and Extinction phases as well as during the entire 6 min for the last day of Extinction phase were compared among groups.

Slice preparation and electrophysiological recordings of PrL pyramidal neurons
On day 14 of CFC, half of the entire sample of mice (n=5/group) was anesthetized with an overdose of halothane (Sigma-Aldrich, Missouri, USA) and decapitated. Once removed, the brain was attached with cyanoacrylate glue to a tray and then sectioned into 275 µM-thick coronal slices by means of a vibratome (Leica VT1200s, Wetzlar, Germany). During slicing, brain was maintained in ice-cold artificial cerebrospinal fluid solution (ACSF) containing (in mM): NaCl

RNA-seq library preparation
After thawing on ice in presence of additional proteinase K, RNA was isolated according to manufacturer's instructions including on-column DNase treatment. RNA samples were quantified and the quality tested by Agilent 2100 Bioanalyzer RNA assay (Agilent Technologies, California, USA) or Caliper (PerkinElmer, Massachusetts, USA) ( Table 5).
Library preparation and sequencing were performed at IGATechnology (Udine, Italy). At least 3 independent biological replicates were used for each group (  did not pass quality control and were excluded from further analyses.

Analysis of RNA sequencing data
GRCm38.p6 genome was used to map the reads, and transcript abundances were estimated using Salmon v1.2 [70]. To obtain gene-level count matrices the quantification data were imported using tximport [71]. All further analyses based on these count matrices were performed with the free software R v4.0.2, Bioconductor v3.11 [72], and the package NOISeq v2.31.0 [73]. Differences in RNA composition between samples were corrected by the Trimmed mean of M-values (TMM) normalization [74], and filtered for low counts based on a count per million reads (CPM) criteria. Subsequently, ARSyNseq was used to remove the technical batch effect and NOISeqBIO was used to assess differential gene expression (q>0.95, equivalent to FDR-corrected P<0.05) [73,75]. GO and KEGG pathway analyses were performed by using ClusterProfiler v3.16.1 [76]. The enrichment map method was used to identify functional modules of mutually overlapping gene sets [76,77].

Statistical Analysis
As regard the behavioral results, a three-factor ANOVA (stimulation x fear x day) on freezing behavior (measured during 0-3 min of contextual FC) and one-factor ANOVA on freezing behavior (measured during 0-6 min of day 14) were used. Newman-Keuls post-hoc comparisons were applied when permitted.
As regard the morphological results, one-factor ANOVAs on spine number and density were used. Newman-Keuls post-hoc comparisons were applied when permitted.
As regard the electrophysiological results, Pearson r Correlation test (for evoked firing activities and EPSC values) and Mann-Whitney U Test (for rheobase and EPSC values) were used.
As for the behavioral, morphological, and electrophysiological results values of P=0.05 were considered statistically significant.
As regard the transcriptomic results, after TMM normalization and low counts filtering, the resulting genes underwent the downstream analysis. Batch effect correction was applied with ARSyN and a PCA was performed to assess sample clustering based on their expression profiles.
Differential expression analysis was performed on Group x Condition design and DEGs were identified using NOISeqBIO, a non-parametric analysis for biological replicates. Significant differentially expressed genes were identified for a q>0.95, equivalent to an FDR-corrected P<0.05. Subsequently, GO and KEGG over-representation analyses were performed and significant pathways were represented by means of enrichment map method to visualize and interpret results.

Conclusions
Given the critical role the amygdala-PrL pyramidal neurons played in fear processing, the characterization of the structural, neurophysiological and molecular changes of this neuronal