Transcriptomic and Network Analyses Reveal Immune Modulation by Endocannabinoids in Approach/Avoidance Traits

Approach and avoidance (A/A) tendencies are stable behavioral traits in responding to rewarding and fearful stimuli. They represent the superordinate division of emotion, and individual differences in such traits are associated with disease susceptibility. The neural circuitry underlying A/A traits is retained to be the cortico-limbic pathway including the amygdala, the central hub for the emotional processing. Furthermore, A/A-specific individual differences are associated with the activity of the endocannabinoid system (ECS) and especially of CB1 receptors whose density and functionality in amygdala differ according to A/A traits. ECS markedly interacts with the immune system (IS). However, how the interplay between ECS and IS is associated with A/A individual differences is still ill-defined. To fill this gap, here we analyzed the interaction between the gene expression of ECS and immune system (IS) in relation to individual differences. To unveil the deep architecture of ECS-IS interaction, we performed cell-specific transcriptomics analysis. Differential gene expression profiling, functional enrichment, and protein–protein interaction network analyses were performed in amygdala pyramidal neurons of mice showing different A/A behavioral tendencies. Several altered pro-inflammatory pathways were identified as associated with individual differences in A/A traits, indicating the chronic activation of the adaptive immune response sustained by the interplay between endocannabinoids and the IS. Furthermore, results showed that the interaction between the two systems modulates synaptic plasticity and neuronal metabolism in individual difference-specific manner. Deepening our knowledge about ECS/IS interaction may provide useful targets for treatment and prevention of psychopathology associated with A/A traits.

and especially CB1 receptors, with individual differences in A/A traits, and demonstrate regional and circuit-specific effects of EC signaling on emotional processes. Furthermore, it has been shown that individual differences, such as behavioral phenotypes or human personalities, are associated with disease susceptibility through immune regulation, even if the mechanisms do not converge on a clear pattern [26][27][28]. It has been reported that human individuals characterized by avoiding traits are prone to developing conditions that are broadly considered to reflect a T-helper 2 cell immune polarization [29][30][31], and to show elevated basal glucocorticoid production [32,33], associated with health and immune negative outcomes [34][35][36][37]. Accordingly, rodents characterized by avoiding traits, defined as having consistently slower-than-median approaching latencies in two novel conditions, have low-grade elevation in glucocorticoid production which is associated with shortened life span [9,10,38]. Furthermore, Michael et al. (2020) suggested that avoiding traits are associated with a glucocorticoid resistant state with poorly regulated innate inflammation and dampened cell-mediated immune responses, probably associated with exaggerated T-helper 2 responses [39].
To date, information on the relationship among A/A traits, ECS, and immune function is lacking. To fill this gap, in the present study the interaction between the gene expression of ECS and immune system (IS) was analyzed in relation to individual differences. To unveil the deep architecture of ECS-IS interaction, we used a cutting-edge tool: we performed cellspecific high-throughput stranded RNA-sequencing of amygdala pyramidal neuron from AV, BA, and AP mice and applied deep transcriptomic analyses consisting in differential gene expression profiling, functional enrichment, and protein-protein interaction network analyses. Results showed that interplay between the two systems is able to modulate synaptic plasticity and neuronal metabolism in an individual-difference-specific manner.

Altered Immune System Activity in the Pyramidal Neurons of Mice Characterized by A/A Traits
In Thy-1 YFP pyramidal neurons of amygdala, differential expression analysis identified 5480 significant DEGs in the comparison between BA and AV mice (2630 up-regulated; 2850 down-regulated), and 886 significant DEGs between BA and AP mice (531 upregulated; 355 down-regulated).
In the comparison between BA and AP mice, the AP mice showed a clear inflammatory pattern with the differentially expressed B-cell receptor (KEGG: mmu04662; p.adjusted < 0.01) and the chemokine signaling (KEGG:mmu04062; p.adjusted < 0.01) pathways. Therefore, we identified an alteration in the expression of immune-related genes in the amygdala pyramidal neurons of mice characterized by A/A traits, in comparison to BA mice.

Altered ECS Activity in the Pyramidal Neurons of Mice Characterized by A/A Traits
In Thy-1 YFP pyramidal neurons of amygdala from AV mice, 68 DEGs belonging to the retrograde ECS signaling pathway (KEGG:mmu04723; p.adjusted < 0.001) were identified in comparison to BA mice ( Figure 2). Notably, FAAH and NAPEPLD genes were significantly down-regulated and the DAAG genes were over-expressed in AV mice when compared to BA mice ( Figure 3).
In Thy-1 YFP pyramidal neurons of amygdala from AP mice compared to BA mice, the obtained DEGs did not enrich any classical ECS pathway, as CB1 and CB2. However, we found that inflammatory mediator regulation of TRP channels was significantly altered (KEGG:mmu04750; p.adjusted < 0.05). In the comparison between BA and AP mice, the AP mice showed a clear inflammatory pattern with the differentially expressed B-cell receptor (KEGG: mmu04662; p. adjusted < 0.01) and the chemokine signaling (KEGG:mmu04062; p. adjusted < 0.01) pathways. Therefore, we identified an alteration in the expression of immune-related genes in the amygdala pyramidal neurons of mice characterized by A/A traits, in comparison to BA mice.

Altered ECS Activity in the Pyramidal Neurons of Mice Characterized by A/A Traits
In Thy-1 YFP pyramidal neurons of amygdala from AV mice, 68 DEGs belonging to the retrograde ECS signaling pathway (KEGG:mmu04723; p. adjusted < 0.001) were identified in comparison to BA mice ( Figure 2). Notably, FAAH and NAPEPLD genes were significantly down-regulated and the DAAG genes were over-expressed in AV mice when compared to BA mice ( Figure 3).
In Thy-1 YFP pyramidal neurons of amygdala from AP mice compared to BA mice, the obtained DEGs did not enrich any classical ECS pathway, as CB1 and CB2. However, we found that inflammatory mediator regulation of TRP channels was significantly altered (KEGG:mmu04750; p. adjusted < 0.05).   . Barplot representing differential expression of genes belonging to retrograde ECS signaling pathway along with genes from pro-inflammatory IL-2 and IL-12 production pathways seen in Section 2.1. These genes will be used in Section 2.3 to study the interplay between ECS and IS.

The Immune and ECS Interplay in the Pyramidal Neurons of Mice Characterized by the Avoiding Trait
Having identified significant alterations of both immune-and ECS-related genes in the comparison between AV and BA mice, we further investigated if these systems were interconnected. Thus, a PPI network was used to investigate the interplay between DEGs of ECS and IS, obtaining a significant PPI score (p < 0.0001).
Based on the network analysis, we identified four communities of functionally connected genes, demonstrating the interaction between DEGs of ECS and pro-inflammatory IL-2 and IL-12 production pathways in communities 1 and 2 ( Figure 4). . Barplot representing differential expression of genes belonging to retrograde ECS signaling pathway along with genes from pro-inflammatory IL-2 and IL-12 production pathways seen in Section 2.1. These genes will be used in Section 2.3 to study the interplay between ECS and IS.

The Immune and ECS Interplay in the Pyramidal Neurons of Mice Characterized by the Avoiding Trait
Having identified significant alterations of both immune-and ECS-related genes in the comparison between AV and BA mice, we further investigated if these systems were interconnected. Thus, a PPI network was used to investigate the interplay between DEGs of ECS and IS, obtaining a significant PPI score (p < 0.0001).
Based on the network analysis, we identified four communities of functionally connected genes, demonstrating the interaction between DEGs of ECS and pro-inflammatory IL-2 and IL-12 production pathways in communities 1 and 2 ( Figure 4). In particular, community 1 enriched the glutamatergic synapse pathway and included all DEGs of IL-2 and IL-12 production systems (p < 0.00001) and of ECS with MAPK1, MAPK3, and MAPK10, representing the connection hubs between the IS and ECS components (Community 3: anandamide metabolism; p < 0.0001) in the glutamatergic and GABAergic transmission (Community 4: GABA receptor complex; p < 0.0001) ( Figure  3). Furthermore, community 2 enriched the mitochondrial respiratory chain complex I (KEGG:mmu00190; p < 0.0001). Specifically, community 2 showed a dense overexpression pattern of cytochrome c oxidase DEGs that regulate the inflammation in association with several genes coding for NADH subunits ( Figure 5). In the comparison between BA and AV phenotypes, the neurotransmission and synaptic plasticity were the most altered functions. In fact, another network analysis on GO terms in the Biological Process component of Gene Ontology showed three functional modules linking together neurotransmission regulation (GO:0099177, regulation of transsynaptic signaling; p < 0.0001), electrophysiological activity (GO:0043269, regulation of ion  In particular, community 1 enriched the glutamatergic synapse pathway and included all DEGs of IL-2 and IL-12 production systems (p < 0.00001) and of ECS with MAPK1, MAPK3, and MAPK10, representing the connection hubs between the IS and ECS components (Community 3: anandamide metabolism; p < 0.0001) in the glutamatergic and GABAergic transmission (Community 4: GABA receptor complex; p < 0.0001) ( Figure 3). Furthermore, community 2 enriched the mitochondrial respiratory chain complex I (KEGG:mmu00190; p < 0.0001). Specifically, community 2 showed a dense overexpression pattern of cytochrome c oxidase DEGs that regulate the inflammation in association with several genes coding for NADH subunits ( Figure 5).
In the comparison between BA and AV phenotypes, the neurotransmission and synaptic plasticity were the most altered functions. In fact, another network analysis on GO terms in the Biological Process component of Gene Ontology showed three functional modules linking together neurotransmission regulation (GO:0099177, regulation of trans-synaptic signaling; p < 0.0001), electrophysiological activity (GO:0043269, regulation of ion transport; p < 0.0001), and neuronal plasticity (GO:0048858, cell projection morphogenesis; p < 0.0001). These three modules converge on the synapse organization (GO:0050808; p < 0.0001) and regulation of membrane potential (GO:0042391; p < 0.0001) pathways, resulting thus in a modulation of behavior (GO:0007610; p < 0.0001) ( Figure 6).

Discussion
Individual differences in A/A traits originate from differences in predisposition toward rewarding or punitive stimuli, respectively. While most studies addressed the differences in the conflict paradigms by manipulating the amygdala-frontal network [4,5,7,[47][48][49][50][51][52], in the present study we faced the A/A traits investigating the molecular differences. Interestingly, the spontaneous individual differences that characterize the three sub-populations of AV, BA, and AP mice are associated with differences in CB1 density in the amygdala [3][4][5]. As described by the extensive literature on the ECS modulation of

Discussion
Individual differences in A/A traits originate from differences in predisposition toward rewarding or punitive stimuli, respectively. While most studies addressed the differ-ences in the conflict paradigms by manipulating the amygdala-frontal network [4,5,7,[47][48][49][50][51][52], in the present study we faced the A/A traits investigating the molecular differences. Interestingly, the spontaneous individual differences that characterize the three sub-populations of AV, BA, and AP mice are associated with differences in CB1 density in the amygdala [3][4][5]. As described by the extensive literature on the ECS modulation of IS, the EC receptors are expressed also on immune cells [53,54]. Furthermore, molecules previously thought to be specific for the immune-cells, such as pro-and anti-inflammatory cytokines, can be also produced and released from neurons.
In the present research, to investigate the interaction between ECS and IS, we performed both RNA-Seq on the pyramidal neurons of the amygdala and bioinformatics analyses in the AV, BA, AP phenotypes.
We identified several altered pro-inflammatory pathways in the pyramidal neurons of the amygdala of AV when compared to the BA mice. In particular, we found in AV mice the alteration of the pathways leading to the production of TNFα, IL-2, and IL-12. The production of these pro-inflammatory cytokines is tightly regulated by the superfamily of transient receptor potential (TRP) and by the B-cell receptor signaling pathways controlling cell differentiation and cytotoxicity. Further, TRP action is facilitated by pro-inflammatory chemokines, responsible for the recruitment of immune cells [55][56][57]. Such findings indicate the chronic activation of adaptive immune response [58].
Although we did not identify any direct difference in the activation of the ECS in AP mice, in comparison to BA mice, the same pattern of inflammation was found in AP and AV mice when compared to BA mice. In fact, the B-cell receptor and chemokine signaling pathways as well as modulation of the inflammatory response by TRP channels were altered in both AP and AV mice. This evidence suggests that the opposite approach and avoidance behaviors may be related to a unique molecular background. In other words, we are advancing that the up-or down-regulation of the genes involved in the interplay between ECS and IS may be a building block of the switch between approach and avoidance.
Since EC interact with and are expressed by almost all immune cells [46], ECS contributes to maintaining the immune homeostasis, functioning as a gate-keeper of IS [40][41][42][43][44] through the activation of the CB1 and CB2 receptors as well as of TRPV1 and TRPV2 receptors. The TRP receptors, mainly TRPV1, are EC targets implied in the functionality of CB1 and in EC degradation [57,59,60]. We found that AV mice displayed an altered retrograde ECS signaling, with NAPEPLD and FAAH coding genes significantly downregulated. It has been proposed that the inhibition of FAAH increases CB1 activity [18]. However, FAAH inhibition also results in increased activation of other receptors, such as TRPV1 and peroxisome proliferator-activated receptor-α (PPARα) [18,61]. These non-classical receptors often have roles opposite to classical EC receptors, hindering the anti-inflammatory role of CB1 in immune modulation [53]. Thus, the downregulation of FAAH expression might sustain inflammatory processes in AV mice. Subsequently, we applied network analysis to further analyze the pattern of co-expression of IS and ECS genes in the regulation of avoidance behavior. Interestingly, in AV mice, ECS and IS genes were mostly co-expressed in the altered glutamatergic synapse pathway of the amygdala, in which the activity-dependent flow of glutamate and EC controls the secretion of pro-inflammatory cytokines through the activity of CB1 and TRP receptors in general [62], and likely of TRPV1 in particular.
This study was performed with limited sample size, as the study of individual differences needs highly curated and homogeneous, and thus inevitably restricted groups. When coming to high-throughput sequencing for gene expression data, a small sample size results in HDLSS (high dimension, low sample size) data problems. This hindered possibilities with the network analysis, forcing us to study gene interactions by using STRINGdb. Although our study was limited to bioinformatics analyses of gene expression from RNA-Seq data, and no functional experiments were conducted, we were able to study gene networks and to map genes to known biological pathways, which contributed to the understanding of ECS/IS interaction.
To further investigate the role of such an interaction in conflict behaviors, the chemical blockade of enzymes linked to ECS metabolism (such as FAAH or NAPEPLD) could be performed. Moreover, anti-inflammatory drugs could be administered to investigate whether AV or AP phenotypes shift toward a balanced approach to conflicting stimuli. Given that the A/A traits represent psychological risk factors [63], shedding light on ECS/IS interaction in such behaviors may provide useful targets for their treatment and prevention.

Subjects and Experimental Procedures
B6.Cg-Tg(Thy1-YFP)HJrs/J (Thy1-YFP; Jackson Laboratories, Bar Harbor, Maine, USA) transgenic mice were used in the present study. Thy1-YFP mice express yellow fluorescent protein (YFP) transgene driven by the thymus cell antigen 1 (Thy1) promoter that targets the pyramidal neurons of the neocortex (layer 5), amygdala, and hippocampus, and to a lesser extent the cerebellar mossy fibers, neurons in the thalamus, midbrain and brainstem, and olfactory bulb mitral cells. Transgene-expressing neurons are morphologically and physiologically comparable to non-mutant neurons.
The animals were group-housed (4 mice/cage) with standard food (Mucedola, Milan, Italy) and water ad libitum, and kept under a 12-h light/dark cycle (with the light on at 07:00 h), controlled temperature (22-23 • C), and constant humidity (60 ± 5%). The behavioral testing took place during the light phase. All efforts were made to minimize animal suffering and to reduce their number, in accordance with the European Directive (Directive 2010/63/EU). The animals assigned to the same experimental group were never siblings.
The first step of the present research was the categorization of mice according to their tendency to A/A responses. To this aim, a sample of transgenic mice (Thy1-YFP; N = 50) was submitted to the A/A Conflict Task. By means of this task, we selected transgenic mice which spontaneously responded with balanced behavior (BA; N = 3), and transgenic mice which spontaneously responded with avoiding (AV; N = 3) or approaching (AP; N = 3) behaviors in the presence of the same conflicting stimuli [3,4]. One month after the end of behavioral task, the AV, BA, and AP mice were sacrificed to sort the amygdala Thy1-YFP pyramidal neurons in order to perform the transcriptomic cell-specific RNA-analyses, tool for unveiling the deep architecture of ECS-IS interaction. To verify whether and which transcriptomic elements (especially linked to endocannabinoid system and immune system) were associated with the specific and stable predisposition to A/A behavior, the transcriptomic analyses were evaluated in the three phenotypes at a time-point distant from any behavioral testing, to rule out any acute effect of the behavioral performance on gene expression levels.

A/A Conflict Task
As previously described [3,5], the apparatus consisted of a Plexiglas Y-maze with a starting gray arm from which two arms (8 cm wide, 30 cm long, and 15 cm high) stemmed, arranged at an angle of 90 • to each other. A T-guillotine door was placed at the end of the starting arm to prevent backward movements of the animal. An arm entry was defined as four legs entering one of the arms. The two choice arms differed in both color and brightness. One of the two arms had a black and opaque floor and walls, and no light inside, while the other had a white floor and walls, and was lit by a 16-W neon lamp. Notably, the colored "furniture" and the neon lamp were exchangeable between arms to alternate the spatial positions of the white and black arms. The apparatus was placed in a room that was slightly lit by a red light (40 W). It was always thoroughly cleaned with 70% ethanol and dried after each trial to remove scent cues. At the end of each choice arm, there was a blue chemically inert tube cap (3 cm in diameter, 1 cm deep) used as a food tray that prevented mice from seeing the reward at a distance but allowed easy access to the reward. Because appetites for palatable foods have to be learned, a week before testing the animals were exposed in their home cages for three days to the novel palatable food (Fonzies, KP Snack Foods, Munchen, Germany). At the end of this phase and during successive testing, to increase the motivation to search for the reward, the animals were slightly food deprived by limiting food access to 12 h/day. About 24 h after the habituation to the apparatus, the slightly food-deprived mice were submitted to two 10-trial sessions of the testing. In Session 1 (S1), the animal was placed in the starting stem and could choose to enter one of the two arms, both containing the same standard food reward. After eating, the animal was allowed to stand in its cage for a 1 min-inter-trial interval. In Session 2 (S2), which started 24 h after S1, the white arm was rewarded with the highly palatable food (Fonzies), while the black arm remained rewarded with the standard food. Notably, the A/A Conflict Task requires to choose between two conflicting drives: reaching a palatable reward placed in an aversive (white and lighted) environment, or reaching a standard food placed in a reassuring (black and dark) environment. The A/A Conflict Index represents the difference in the number of white choices between S1 and S2. Given that this index was normally distributed (mean ∆ = 1, SD = ±1.7), it allowed identifying the specific phenotypes of the mice [4]. In particular, in the presence of conflicting inputs, we identified 3 BA mice that reacted with balanced responses between A/A traits and their values in the A/A Conflict Index corresponded to the mean of the distribution. Furthermore, we identified 3 AV mice that exhibited avoiding responses to the conflicting stimuli and that had values of the A/A Conflict Index corresponding to minus two standard deviations of the mean. We identified also 3 AP mice that exhibited approaching responses to the conflicting stimuli and that had values of the A/A Conflict Index corresponding to minus two standard deviations of the mean.

Dissociation of Amygdala Tissue for Fluorescence-Activated Cell Sorting (FACS)
One month after the end of A/A behavioral task, the brains of AV, BA, and AP mice were cut to take bilateral amygdala 1-mm punches. Manual and enzymatic dissociations were performed using the Neural Tissue Dissociation Kit (P) (Miltenyi Biotec, Bergisch Gladbach, Germany) with some modifications. Each solution was kept on ice to minimize RNA degradation. Pipette tips were pre-coated in a 0.2 µM filtered 1× PBS-0.5% BSA solution (DPBS without Mg 2+ and Ca 2+ , Gibco by Life Technologies, Grand Island, NY, USA; BSA Fraction V (pH 7.0), PanReac AppliChem GmbH, Darmstadt, Germany). Briefly, the amygdala punches were placed on a 35-mm diameter Petri dish, cut into small pieces using a scalpel, and 1 mL of cold Hanks' Balanced Salt Solution without Mg 2+ and Ca 2+ (HBBS w/o) (Sigma-Aldrich, St. Louis, MO, USA) was added. The tissue was transferred into a 1.5 mL protein LoBind tube. Additional 1 mL HBBS w/o was used to rinse the dish and added to the 1.5 mL tube. Tissue was centrifuged at 300× g for 2 min at room temperature, and the supernatant was carefully aspirated. Then, 975 µL of pre-heated enzyme mix 1 (enzyme P 25 µL, buffer × 950 µL) was added to the tissue, and the 1.5 mL tube was incubated for 15 min at 37 • C under slow, continuous rotation using the MACSmix Tube Rotator (Miltenyi Biotec, Bergisch Gladbach, Germany). Then, 15 µL enzyme mix 2 (enzyme A 5 µL, buffer Y 10 µL) was added to the sample. The sample was gently inverted to mix and mechanically dissociated using the wide-tipped fire-polished Pasteur pipette by pipetting up and down 10 times slowly, followed by a further incubation in the rotator for 10 min at 37 • C under slow rotation. The second round of mechanical dissociation was performed using serially fire-polished filtered-glass Pasteur pipettes with gradual diameter diminution, and pipetting slowly up and down 10 times with each pipette, or until no tissue pieces were observable. The sample was again incubated at 37 • C for 10 min using rotator under slow rotation. The sample was strained through a MACS Smart Strainer (70 µm) (Miltenyi Biotec, Bergisch Gladbach, Germany), placed on a 15 mL tube, pre-coated with 0.2 µM filtered 1× PBS-0.5% BSA, using 8 mL of HBBS with Mg 2+ and Ca 2+ . Cells were pelleted by centrifugation at 300× g for 10 min at room temperature. In order to increase the material, the supernatant was transferred into a new 15 mL tube, and centrifuged again at 300× g for 10 min at room temperature. The second supernatant was discarded, and the pellets obtained from these two centrifugations were pooled into a 1× PBS-0.5% BSA pre-coated SNAP-cap tube containing 1 mL of PBS. Finally, 20U Superase-Inhibitor (Ambion, Invitrogen, ThermoFisher Scientific, Walthem, MA, USA) was added and samples were stored on ice up to sorting.

Cell Sorting and Isolation of Purified Pyramidal Neurons
For the instrument set-up, the samples collected from the amygdala of wild-type YFP-negative mice were used to gate YFP-positive neurons based on forward scatter (FSC) and side scatter (SSC) light scattering and to set YFP negativity. Afterwards, amygdala samples were collected from the Thy1-YFP AV, BA, and AP mice and stained with 1 µL of propidium iodide (PI) in order to identify dead cells. Pyramidal neurons were then sorted by using the MoFlo Astrios EQ (Beckman Coulter, Brea, CA, USA) and the pyramidal neurons characterized by YFP were collected on the basis of their physical parameters, singlets, PI negative (live cells), and YFP intensity. For initial characterization, samples were collected in PBS and samples were examined under a fluorescent microscope to verify correct sorting. Thereafter, cells were sorted directly into ice-cold lysis buffer (Reliaprep RNA Cell Miniprep System, Promega, Fitchburg, WI, USA), mixed by vortexing, kept on ice, and then stored at −80 • C until RNA extraction.

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 was tested by Agilent 2100 Bioanalyzer RNA assay (Agilent Technologies, Santa Clara, CA, USA) or Caliper (PerkinElmer, Waltham, MA, USA).
Library preparation and sequencing were performed at IGATechnology (Udine, Italy). At least three independent biological replicates were used for each group. Each replicate corresponds to the amygdala of a single Thy1-YFP mouse.
Libraries were generated from each sample individually, starting from 0.06-4.41 ng of total RNA, using the Ovation SoLo RNA-Seq kit for Ultra-low input (NuGEN, Tecan Genomics, Redwood City, CA, USA), following the manufacturer's instructions (library type: fr-second strand). Final libraries were checked with both Qubit 2.0 Fluorometer (Invitrogen by Life technologies, Carlsbad, CA, USA) and Agilent Bioanalyzer DNA assay (Agilent Technologies, Santa Clara, CA, USA) or Caliper (PerkinElmer, Waltham, MA, USA). Libraries were then prepared for sequencing and sequenced on paired-end 2 × 75 bp mode on NextSeq500 (Illumina, San Diego, CA, USA) producing 33.9 MR on average (min 29.3 MR, max 40.7 MR). For the processing of raw data (format conversion and de-multiplexing), Bcl2Fastq 2.20 version of the Illumina pipeline was used.

Differential Expression Analysis
After TMM normalization and low counts filtering, the resulting genes (AMY = 9872) underwent the downstream analysis. Batch effect correction was applied with ARSyN and a principal component analysis (PCA) was performed to assess sample clustering based on their expression profiles. After PCA, 1 AV mouse was identified as outlier and removed from downstream analyses. Differentially expressed genes (DEGs) were identified using nonparametric analysis for biological replicates included in the NOISeq library. Significant differentially expressed genes were identified for a q > 0.95, equivalent to a p < 0.05 after FDR correction [64]. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Over-Representation Analyses (ORAs) were performed using Clusterprofiler to identify enriched pathways. The obtained GO terms were further clustered using GOSemSim and visualized using pathview [65] and the enrich map method to visualize and interpret results [66,67]. All biostatistical analyses were performed in R v.4.1 [68]. DEGs from retrograde endocannabinoid signaling and immune system were enriched on STRINGdb [69] to assess shared functions and known co-expression patterns, along with their protein-protein interaction (PPI) enrichment scores.