Transcriptomic Studies of Antidepressant Action in Rodent Models of Depression: A First Meta-Analysis

Antidepressants (ADs) are, for now, the best everyday treatment we have for moderate to severe major depressive episodes (MDEs). ADs are among the most prescribed drugs in the Western Hemisphere; however, the trial-and-error prescription strategy and side-effects leave a lot to be desired. More than 60% of patients suffering from major depression fail to respond to the first AD they are prescribed. For those who respond, full response is only observed after several weeks of treatment. In addition, there are no biomarkers that could help with therapeutic decisions; meanwhile, this is already true in cancer and other fields of medicine. For years, many investigators have been working to decipher the underlying mechanisms of AD response. Here, we provide the first systematic review of animal models. We thoroughly searched all the studies involving rodents, profiling transcriptomic alterations consecutive to AD treatment in naïve animals or in animals subjected to stress-induced models of depression. We have been confronted by an important heterogeneity regarding the drugs and the experimental settings. Thus, we perform a meta-analysis of the AD signature of fluoxetine (FLX) in the hippocampus, the most studied target. Among genes and pathways consistently modulated across species, we identify both old players of AD action and novel transcriptional biomarker candidates that warrant further investigation. We discuss the most prominent transcripts (immediate early genes and activity-dependent synaptic plasticity pathways). We also stress the need for systematic studies of AD action in animal models that span across sex, peripheral and central tissues, and pharmacological classes.


Introduction
Major depression, which affects nearly 5% of the world's adult population, has become the leading cause of disability, and thus represents a heavy burden for our societies [1][2][3]. However, the clinical management of major depression remains "artisanal", involving many "trial-and-error" attempts based on the intuition of the medical practitioner [4].
Moderate to severe depressive episodes are treated with antidepressant (AD) drugs including selective serotonin reuptake inhibitors (SSRIs), serotonin and norepinephrine reuptake inhibitors (SNRIs), and older, less selective compounds. More than 30 different ADs are available today at an ambulatory setting, and there is no clear difference in overall efficacy between the different AD treatments [5]. Despite their known activity toward targets involved in the pathophysiology of major depression [6] (all of these ADs but the newly approved ketamine target the monoaminergic system, acting on a more or less selective combination of monoamine transporters and/or receptors), they do not provide an effective response for many patients. More than 60% of patients suffering from a major depressive episode (MDE) fail to achieve remission following the AD they are prescribed. Even more, one-third of patients will not reach remission after up to four AD trials [7]. For those who respond, full response is only observed after several weeks of treatment. This delay negatively affects patient compliance, increases overall pharmacological burden, and potentially dangerous side effects.
The discovery of novel targets for better and safer Ads, as well as the identification of biomarkers related to therapeutic efficacy, is urgently needed [8]. For the moment, despite efforts in this direction, there is no reliable biomarker that can predict the therapeutic response in a patient suffering from an MDE, and there is no absolute predictor to guide the choice of the therapeutic approach [9]. For this, we believe it is necessary to know exactly what the effects of AD drugs are, in particular the changes in the gene expression program that they induce over time [10], in the parts of our body directly involved in the processing of emotions and mood.
Indeed, RNA, as an immediate product of gene expression and epigenetic programs, is a perfect reflection of the functional status of an individual. For this reason, many investigators have compared peripheral RNA expression profiles between subjects suffering from major depression and control subjects or between patients at different times during their AD treatment [11,12]. Unfortunately, the complexity of the conditions of human subjects in terms of personal history (stress at different times of life, age, sex, obesity, cardiovascular problems, infectious history, tobacco and alcohol consumption, hereditary factors, and dietary habits) and the different undergoing pharmacological treatments at the time of study require very large cohorts to achieve sufficient statistical power and allow the identification of universal signatures. Furthermore, it would be necessary to study the signatures induced by AD drugs in a context of healthy subjects to compare them to those of subjects suffering from an MDE if we want to be able to effectively select new therapeutic targets. However, from this point of view, the data available in humans are very limited.
To better understand and manage affective disorders, major efforts have been made in recent years to characterize the transcriptional signatures of major depression and AD action that distinguish patients from healthy subjects or their equivalent in animal models, particularly in mice and rats [13]. Animal models provide a relatively homogeneous study population in which the events triggering the onset of depressive-type symptoms can be precisely controlled, as well as the timing of the administration of a particular drug at a specific dose and for a defined duration [14]. In addition, animal models allow access to all tissues that may be relevant for both pathophysiology and biomarker development, which is not always the case in humans for evident ethical and technical reasons. Furthermore, it should be noted that in humans, for the time being, most of the transcriptional data have been obtained from post-mortem brain tissue and that it is complicated in these circumstances to disentangle the signature of treatments from that of the long disease process and the terminal effect of death itself. This is something that can be much better controlled in an animal model where, also, new technologies may rapidly advance our understanding of psychiatric disorders at the cellular, molecular, and brain circuit levels [15].
Much progress has been made on the mechanisms leading to the phenomena of vulnerability or, on the contrary, resilience to stress using animal models and the modalities of the response to AD treatment, in particular with regard to transcriptional signatures [16][17][18][19][20]. As in humans, animal models reveal absence of homogeneous response, which offers the possibility of selecting animals that respond well to AD treatment to study the biological pathways involved in the recovery of normal social, emotional, hedonic, affective, mnemonic, exploratory, and feeding behaviors [21]. Going a step beyond, recent conver-gence studies have been undertaken from human and animal transcriptomes and should advance our knowledge of the pathophysiology of major depression [16,[22][23][24][25]. Concerning the identification of signatures of AD treatments, several teams have tried to gather common features of genome-wide transcriptional variations induced by the action of ADs in patients suffering from affective disorders [11,12,26], but to our knowledge, studies carried out in rodent models are very scarce and fragmented.
Despite the privileged use of SSRIs in the pharmacological armamentarium to treat major depression, and the widespread use of this pharmacological class in omics studies in humans and animal models [27], we found only one review regarding studies conducted in animals to evaluate the transcriptional effects of these drugs [28]. However, this work is necessary to establish the targets common to humans and animals that would deserve special attention to improve therapeutic efficacy, which remains too limited today.
Since, to our knowledge, no study has established common patterns of gene expression program variation in rodent models of AD drug exposure, we searched all available transcriptomic data in rats and mice with AD drug exposure and sought to define recurrent elements of a biological signature with already known biological processes involved in mood disorder pathophysiology, as well as other processes that should merit further investigations.

Preponderance of Studies Investigating the Effect of Fluoxetine
We first listed all studies reporting effects of a confirmed or putative antidepressant intervention on the variation of transcriptional expression in the mouse or the rat in (i) naïve animals and (ii) in animals subjected to a paradigm modeling negative affects and depression-like behaviors. Focusing on studies using sufficiently mature technologies and considering the whole genome, we identified 108 published articles and 6 publicly available but unpublished transcriptomic datasets ( Figure 1, Table S1). We noticed, however, that in some cases the same animal cohorts have been used in multiple publications. We inventoried 107 independent animal cohorts.
Nature of AD interventions: When we look at the class of the AD treatments tested ( Figure 2D), we see that almost half of the cohorts used SSRIs, in particular fluoxetine (FLX), which was used in more than a third of the studies. The second most represented class is tricyclics, especially imipramine, followed by the NMDA receptor (NMDAR) antagonists, essentially ketamine. Only after that do we find the SNRIs, followed by mood stabilizers (lithium), atypical antipsychotics, electrical stimulation of the brain (essentially by electroconvulsive therapy, ECT), monoamine oxidase inhibitors (MAOIs), and histone deacetylase (HDAC) inhibitors. Finally, while there is a wide variety of AD interventions (both drugs and non-drugs), of the 65 ADs we have listed (Table S2), more than half (N = 40) have only been tested for pan-genomic activity once.
Animal models: Numerous paradigms have been developed in rodents to model symptoms related to human mood disorders and this diversity is reflected in the transcriptomic studies that we surveyed ( Figure 2E). Thus, while no fewer than 13 different paradigms have been used (Table S3), the most popular is unpredictable chronic moderate stress (UCMS) used in 40% of studies, while electric shocks, restraint stress, social defeat, and maternal separation each concerned 8-10% of the studies.
Region of interest: From an anatomical point of view ( Figure 2F), the brain is the most studied organ, overwhelmingly with 105 out of 107 cohorts, while for the periphery, just three studies profiled blood [13,21,39,40], and only one study examined the adrenal gland [41], liver and kidney [42,43], or mammary glands [31]. Concerning the brain, there is a strong heterogeneity in the size of the cerebral areas or the type of cells profiled ranging from the whole brain to the paraventricular nucleus of the hypothalamus [44] to a particular type of cortical neurons [37,45,46], but, in fact, the most studied area is the hippocampus (as a whole or more precisely the dentate gyri) that we found in almost half of the studies, far ahead of the prefrontal cortex, the frontal cortex, the amygdala and the nucleus accumbens. Finally, although animal models offer the possibility of easily recovering several different anatomical areas and thus determine the possible regionalization of expression profiles, most studies only sampled one anatomical area and only four studies distinguished more than three anatomical areas from the same animals [30,[47][48][49].
three studies profiled blood [13,21,39,40], and only one study examined the adrenal gland [41], liver and kidney [42,43], or mammary glands [31]. Concerning the brain, there is a strong heterogeneity in the size of the cerebral areas or the type of cells profiled ranging from the whole brain to the paraventricular nucleus of the hypothalamus [44] to a particular type of cortical neurons [37,45,46], but, in fact, the most studied area is the hippocampus (as a whole or more precisely the dentate gyri) that we found in almost half of the studies, far ahead of the prefrontal cortex, the frontal cortex, the amygdala and the nucleus accumbens. Finally, although animal models offer the possibility of easily recovering several different anatomical areas and thus determine the possible regionalization of expression profiles, most studies only sampled one anatomical area and only four studies distinguished more than three anatomical areas from the same animals [30,[47][48][49].
Given the great variability in terms of AD intervention, paradigms used, species, sex, brain regions, and assay methods, it is difficult to make meaningful comparison across treatments. This is even more difficult considering the fact that a little more than half of the expression data are not publicly available ( Figure 1). Therefore, we decided to focus on the effect of FLX on the hippocampus, which nevertheless represents 12 studies (Table 1).   Given the great variability in terms of AD intervention, paradigms used, species, sex, brain regions, and assay methods, it is difficult to make meaningful comparison across treatments. This is even more difficult considering the fact that a little more than half of the expression data are not publicly available ( Figure 1). Therefore, we decided to focus on the effect of FLX on the hippocampus, which nevertheless represents 12 studies (Table 1).

Signature of FLX Response in Stressed Rodents
After retrieving the complete data for each of the 12 selected studies, we established for each dataset the lists of differentially expressed genes between the control condition and the FLX treatment condition. First, we grouped studies in mice and rats where FLX treatment restored behavior in stressed animals, which represents four studies in mice and three in rats (Table 1). We felt that it would not be appropriate to merge all the raw data (after different calibration procedures) into one large dataset because of the great heterogeneity in the specifics of each study. In addition to the differences between species and strains, sample sizes, transcriptomic platforms, and behavioral protocols also differed between studies. Instead, we kept from each individual study the list of the top 300 genes with nominal p values below the arbitrary threshold of 0.05. Then, we combined the individual lists into a single graph to allow identification of genes with repeated dysregulated expression in multiple experiments ( Figure S1). In this way, we extracted 1973 genes that were modulated by FLX in at least one experiment in stressed animals. Of these genes, 98 were dysregulated in at least four different experiments (Table S4). We also defined for each gene the global trend of the direction of variation of the transcriptional expression with a positive consensus score for overexpression and negative for underexpression. Of note, although nine genes (Arc, Ddah1, Egr1, Hmgcs1, Kcng2, Klhl5, Nr4a1, Oxtr, and Zfhx2) show significant expression variation in at least five of the experiments, only four show a clear directional pattern (Ddah1 and Oxtr are predominantly upregulated, whereas Klhl5 and Zfhx2 are predominantly downregulated). Interestingly, at best, only 15 genes show a significant direction of variation that is consistent across most experiments (Table 2).
In the literature there is no consensus way to establish a signature based on a collection of different datasets. As an alternative approach, we considered defining the FLX action portrait based on the methodology recently proposed by Stephen Gammie to establish a signature of major depression in the human brain [23]. Applying this methodology, we ranked 39,629 genes and obtained scores ranging from +4.56 for the most overexpressed consensually in the seven studies (immediate early response gene, Ier5) to −4.46 for the most underexpressed (Sh3d19 , Table S5). We listed in Table 3 the signature of effect of FLX with genes showing absolute values of portrait scores equal or above 4. When we compared the integration ( Table 2, N = 22) with the portrait score (Table 3, N = 12) signature, we identified two genes (Hmgcs1 and Prkar1b) that are common to both signatures. The advantage of having a sign indicating the overall direction of expression variation, whether with the consensus or the portrait score, allows us to assess whether there is any distortion in the directionality of FLX-induced transcriptional expression variation. We can see that this is not the case because with both scoring methods a little less than 15% of the genes show a directionality of their expression variation in one direction or another and thus for nearly 70% of the genes no trend is apparent. To extract the most salient variations with the portrait score, we separated the genes with scores whose absolute value is equal or greater than 3 and thus obtained 82 overexpressed genes and 80 underexpressed genes (Table S5). As a more general and consensual signature of FLX in a stress paradigm, we combined both scoring procedures with absolute values equal or above 2, resulting in 412 upregulated genes and 411 downregulated genes in an almost perfect balance (Table S5).
To identify the biological processes targeted by the signatures described above, we conducted ontological analysis with the DAVID algorithm. Table S6 shows ontological results for the list of the 98 most commonly deregulated genes from Table S4 (integration method). Two overrepresented processes emerged: the activity of transcription factors and the synapse. Then, we performed ontological analysis on the 823 genes with consensus and portrait scores above 2 after separating up- (Table 4) and downregulated genes ( Table 5). Overexpressed genes were strongly associated with the synapse, postsynaptic density, and cell junction, as well as with glutamatergic and oxytocin signaling pathways. We also found recurrent modules classically associated with AD function. On the contrary, underexpressed genes were largely enriched for the translation machinery, the ribosome.

Signature of FLX Response in Naive Rodents
Regarding the signature of the effect of FLX on naïve animals, we surveyed seven different datasets, five in mice and two in rats; 1763 genes showed a variation in expression in at least one experiment ( Figure S2) and 117 in at least four experiments (Table S7). Sixteen genes were significantly affected by FLX in five out of seven experiments, and in fact, 36 genes showed significant variation in the same direction in most experiments ( Table 6).

Shared Signature of FLX in Stressed and Naïve Rodents
To establish the core of FLX-induced transcriptional activity in the hippocampus, regardless of the paradigm employed, we first examined the convergence between the stress and non-stress signatures previously described and obtained by either integration or scoring methods. Figure 3 shows that only one gene is common to all comparisons, Zfhx2. Next, we combined in the same dataset, the 14 comparisons made previously (7 in a stress context and 7 in a naïve context, Figure S3). With the integration method, among 3360 genes modulated in at least one experiment following FLX treatment, 78 were altered in at least 7 experiments, 15 in at least 8 experiments, while Pdlim5 and Zfhx2 were altered in 9 and 10 comparisons, respectively (Table S11). Second, by applying the portrait method to the 14 comparisons, we could rank 40,113 genes with top scorer genes, Sel1l3 (+7.88) and Nfia (−8.13) (Table S12). When we examine the best scoring genes, with absolute values of both integration and portrait scores equal or above 6, we obtain nine upregulated genes: Ddr1, Ier5, Igfbp6, Nptx2, Prkar1b, Ptpn5, Sel1l3, Tyro3, and Zfp703, and seven downregulated genes: Akt3, Fat4, Nfia, Pcdh19, Rab27a, Scn3b, and Zfhx2.
Without surprise, the ontological analysis of the processes enriched from the list of 78 genes most affected by FLX action (Table S11 integration method) reveals mainly a synaptic activity (Table S13). When performing ontological analysis on the 113 upregulated genes with consensus and portrait scores above or equal to 4, the axon and the MAPK signaling pathways were significantly enriched (Table 9). For the 151 downregulated genes, the MAPK signaling pathway was only process significantly affected (Table 10). A secondary analysis of the genes involved in the MAPK signaling pathway identified distinct components for up and downregulated genes. Upregulated genes were associated with neurotrophic pathways while downregulated genes were restricted to protein phosphorylation.
The ontological analysis of the processes enriched from the list of 78 genes deregulated after FLX action confirms with both ontological tools that neuronal structure, development, and signaling are the main processes concerned (Tables 8 and 9). Genes belonging to the core signatures (Tables 2 and 3) are in bold. Genes belonging to the core signatures (Tables 2 and 3) are in bold.      By evaluating lists of equal size obtained by integration method (Table S11) versus consensus and portrait score (absolute values equal or above 5, Table S12), we defined 15 FLX-modulated genes: Baalc, Igfbp6, Itga4, Nptx2, Prkar1b, Rasgrf1, S100a6, Sel1l3, Slc4a4, Sorcs1, Tmem47, Trpm3, Zfhx2, and Zfp316. Among these genes, five are also listed in Figure 3 (Baalc, Igfbp6, Prkar1b, Sel1l3, and Zfhx2).
Without surprise, the ontological analysis of the processes enriched from the list of 78 genes most affected by FLX action (Table S11 integration method) reveals mainly a synaptic activity (Table S13). When performing ontological analysis on the 113 upregulated genes with consensus and portrait scores above or equal to 4, the axon and the MAPK signaling pathways were significantly enriched (Table 9). For the 151 downregulated genes, the MAPK signaling pathway was only process significantly affected (Table 10). A secondary analysis of the genes involved in the MAPK signaling pathway identified distinct components for up and downregulated genes. Upregulated genes were associated with neurotrophic pathways while downregulated genes were restricted to protein phosphorylation.   Genes belonging to the core signature (Tables S11 and S12: genes with absolute scores ≥6, Figure 3) are in bold. Genes belonging to the core signature (TableS S11 and S12: genes with absolute scores ≥6, Figure 3) are in bold.

Discussion
In this study we have voluntarily limited our field of research to post-2006 publications on transcriptional signatures of the effect of AD treatments in animal models and have identified 114 studies. Nevertheless, we were confronted with the difficulty of accessing complete transcriptional data for more than 60 studies for which data were not available in a public manner. This unfortunately reflects, at the minimum, a lack of reflexes on the part of investigators, or perhaps also a lack of sufficient incentive on the part of the authors and publishers, to share data [59]. Indeed, omics studies come at a substantial cost, and it is conceivable that some authors would prefer to fully capitalize on future analyses before sharing the data publicly.

Current Picture: A Fragmented Landscape
A detailed review of the available studies shows a fragmented landscape regarding the design of the studies and an important imbalance in many experimental factors regarding sex, region, and/or ADs studied (Figure 1).
Although women are twice as affected by major depression, preclinical study models of this disease and its treatment are still too often focused on male animals. Thus, among all the studies studying the transcriptome following the action of AD molecules, less than 10% concern females (Figure 1). In addition, recent studies suggest that the transcriptional signatures of depression are different in men and women [60]. It is therefore essential that more preclinical studies seriously address the changes induced in females by the major therapeutic agents that are or will be offered to patients with mood disorders.
A similar imbalance exists for brain regions targeted in transcriptomic studies of AD effects (Figure 1). This is an important caveat because, while stress (acute and prolonged) causes molecular and functional changes in several brain areas such as the medial prefrontal cortex, amygdala, nucleus accumbens, and hippocampus, among the most studied, the stress-induced molecular signatures differ across regions. For instance, activation of the CREB-BDNF axis in the hippocampus is antidepressant, but the same activation in the nucleus accumbens is pro-depressant [61]. Region-dependent effects of stress on neuronal activation and BDNF signaling could also be sex-dependent [62,63]. Comparing AD transcriptomic signatures across regions in male and female mice could provide potential insights into how brain pharmacology can be modulated in a region-dependent and sexspecific manner; however, we could not undertake this analysis because of insufficient power.
Another big discrepancy that we observed across studies concerns pharmacological classes tested and compounds inside the same class. Although SSRIs remain the most prescribed ADs, it is important to revisit older, less selective classes of ADs. Indeed, many open questions remain in terms of both the delay of onset [64] and of the mechanism of action [65]. Again, there was not enough power (numbers) in the available studies to consider for our analysis. Therefore, our review of the literature underlines the need for systematic studies designed to include male and female mice, different regions, and more than one class of ADs.

Methodological Choices
It is not a surprise that the majority of the available studies concerns the most prescribed AD in the world, Prozac, approved by the Food and Drug Administration (FDA) in 1987, whose active ingredient is FLX [66]. Moreover, the preferred study area for the transcriptional effect induced by ADs is the hippocampus, a region consistently implicated in major depression and AD action [67]. Thus, we decided to focus on the transcriptomic signature produced by FLX in the hippocampus. In the different studies analyzed for some is the whole hippocampus that has been collected; in other cases, it is the dentate gyrus and/or the cornu ammonis 1 (CA1) and sometimes it is specified whether it is the ventral or dorsal part ( Table 1). The hippocampus is functionally divided into a dorsal region that is primarily engaged in cognitive functions and a ventral region that regulates emotion [68,69]. Although it has been observed for one of the experiments we included that similar gene expression responses to FLX were found in dorsal and ventral dentate gyrus [51], more recent investigations show that the ventral hippocampus was particularly sensitive to the effects of stress. Therefore, it was proposed to consider the dorsal and ventral hippocampus separately when conducting high-throughput molecular analyses [24,70]. We are aware of this caveat, but we had no choice if we wanted to keep enough power in our analysis. Stephen Gammie has recently drawn a picture of depression based on publicly available transcriptomic data from human brains [23]. His comparison of the human signature with the signatures of nearly 200 data sets, mainly from mice or rats treated with ADs and other drug classes, showed that the effect of FLX assessed at the hippocampus level provided the closest inverted signature to that of human depression [23].
Another choice we made, for the sake of homogeneity, was to select, at least for one experiment, samples from mice responding to FLX (GSE84185). Indeed, if it seems that most investigators retained for transcriptomic analysis only samples from animals showing an improvement in behavior following FLX administration, in at least one case the overall response was profiled [21]. It may be noted that at least two studies had highlighted interindividual variability in behavioral responses to FLX [21,51]. It was even proposed that it would be interesting to evaluate ambiguous responders, as they could be useful in dissociating anti-depression-like effects from anti-anxiety-like effects [51]. Later, such behavioral heterogeneity, especially concerning anxiety, following chronic FLX treatment was confirmed and may reflect a specific gene expression profile [46].
A limitation of our work is that we almost only processed protein-coding mRNAs. However, emerging literature, especially in terms of biomarker identification for mood disorders, concerns microRNAs (miRNAs) [71,72]. Thus, among the studies that we identified as potentially relevant for our analysis, nine concerned miRNAs, two of which were interested in the effect of FLX in the hippocampus [38,73]. It will therefore be important in the future to compare miRNA and mRNA signatures to study whether they are related [74,75], which will help to better understand the mechanism of action of ADs and undoubtedly improve their effectiveness [76]. In addition to miRNAs, a large family of non-coding RNAs could not be considered in our analysis. These are the myriad of long non-coding RNAs (lncRNAs) [77]. However, not only are they very present for a while on microarrays, and a fortiori are widely detected by high-throughput sequencing, but it is often in this class of RNA that we find the strongest expression variations [78,79]. Efforts to standardize annotations and a better knowledge of their exact role should in the future help us to better integrate them into a more global signature of the effect of ADs [80,81]. Circular RNAs are another category that is nowadays emerging in the field of psychiatry transcriptomics [82]. However, the available data in animal models of depression and antidepressant action are still limited, albeit promising [83].

Serotonin System and BDNF
The primary target of FLX is the serotonin transporter 5-HTT, encoded by the SLC6A4 gene. Although there may be transient effects on variations in Slc6a4 mRNA expression levels in rodent models [84,85], there is little evidence to suggest a transcriptional effect of FLX on either SLC6A4 [86] or the limiting enzyme for serotonin synthesis in brain, tryptophan hydroxylase 2 (TPH2) [28], and we did not observe robust modifications for these two genes. On the other hand, in the large family of serotonin receptors, several studies have implicated transcriptional variations in the gene encoding the 1B receptor (HTR1B). We confirm in our meta-analysis this effect for Htr1b but also for Htr5b (Table 6), for which there is no functional equivalent in humans [87]. One target related to serotonin is the overexpressed Vip messenger RNA (Tables 7, S8 and S12) encoding vasoactive intestinal peptide, as it has been known for a long time that the neurotransmitter effect of serotonin in the hippocampus can be modulated by Vip through regulation of cyclic adenosine 3':5'-monophosphate (cAMP) levels and serotonin receptors [88][89][90]. It was also found that FLX, while improving depressed behavior in a rat model of chronic stress-induced depression, increased Vip expression [91]. Other expected targets for FLX include brainderived neurotrophic factor (BDNF), which is widely implicated in depression and the mechanisms of its treatment and would be expected to be increased [28], which is confirmed both in a stress paradigm (Tables 4 and S5), or in naive animals (Table S8) and thus also in the overall meta-analysis (Tables 9 and S12).

Immediate Early Genes
In the present meta-analysis, among the genes significantly modulated by FLX appear the immediate early genes Nr4a1, Nr4a2, Arc, Egr1, Fosb, Fosl2, and Junb (Tables 2, S4 and S6). These are genes that have long been identified as particularly important in the program of gene expression modification during the induction of depression in different animal models [92][93][94]. EGR1, in particular, is the major downstream partner of the ERK/Elk-1 cascade that we have recently proposed as a novel target for AD development [67]. Interestingly, such immediate early genes are regulated not only in the hippocampus but also in other brain areas such as frontal and prefrontal cortex, lateral amygdala, and nucleus accumbens [95][96][97][98][99], and they are the indirect target of other AD drugs than FLX such as agomelatine, duloxetine, imipramine, ketamine, paroxetine, or vortioxetine [100][101][102][103][104][105] but also of other AD-related interventions such as sleep deprivation and deep brain stimulation [106][107][108]. In addition to this set of genes, we also noted the immediate early response 5 gene (Ier5), which appears upregulated in our meta-analysis (Tables 3, S5, S8 and S12), was previously found dysregulated in peripheral blood mononuclear cells of unmedicated mood-disorder patients compared to healthy controls [109].

Signal Transduction Pathways
It is not a surprise that our analysis highlights universal signal transduction pathways. Following receptor activation, second messengers fine-tune neuronal activity, synaptic remodeling, long-term potentiation, and neurogenesis. They thus regulate synaptic plasticity and cellular resilience, processes by which the brain perceives, adapts, and responds to a variety of internal and external stimuli (including stress and other depressogenic factors [110]. Over the last decades, evidence from different research groups showed that kinase-phosphatase pathways are important mediators of AD action, and have used translational models, namely genetically modified animals to phenocopy AD-sensitive behaviors and to highlight AD-sensitive brain circuits [111][112][113][114].

a MAP kinases
In line with the neurotrophic hypothesis of depression, genes coding for members of the neurotrophic tyrosine receptor kinase family such as Ntrk2 and Ntrk3 were normalized by FLX action (Tables 4, 6, 9 and S12) in a reverse direction compared to observation in rodent stress models of depression [115]. These genes encode kinases that, upon neurotrophin binding, phosphorylate members of the MAPK pathway. Thus, it makes sense to find the MAPK signaling pathway as one of the pathways most consistently affected by FLX (Tables 9 and 10). Among the multiple connections between BDNF, MAP kinases, and AD action, we found an upregulation of the Ptpn5 gene expression (Tables S5, S8, S12 and 9 and Figure 3). It encodes a striatal-enriched tyrosine protein phosphatase (STEP) that inactivates key neuronal signaling proteins such as MAP kinases, tyrosine kinases NMDA, and AMPA receptors and whose inactivation decreases the expression of BDNF [116]. Therefore, STEP inhibitors have been proposed as a novel target for AD drugs of the new generation [117].

b WNT/catenin
Another potential player in this activation of the MAPK pathway is insulin-like growth factor binding protein 6 encoded by Igfbp6 mRNA, which is overexpressed in the different paradigms in response to FLX (Tables 6, 7 and 9 and Figure 3) including studies not included in the current meta-analysis [118]. There is even another member of this gene family with Igfbp4 (Figure 3), which can be linked to the activation of another biological pathway important for the resolution of depressive symptoms in connection with hippocampal neurogenesis, the Wnt/β-catenin pathway [119,120]. That same pathway is also activated by the neuronal pentraxin 2 [121], encoded by the Nptx2 gene that is upregulated by FLX (Tables 3, 9, S5, S8 and S12) as already demonstrated in a previous study focused on a shared mechanism of AD effect of chronic FLX and exercise [118].
The response to ADs involves several overlapping mechanisms and, in addition to the signaling mentioned above, the PI3K/AKT pathway can also be mentioned [122]. Thus, a dysregulation of Akt3 expression is observed (Tables 7 and 10), knowing that the invalidation of this gene in mice revealed an endophenotype reminiscent of psychiatric manifestations such as schizophrenia, anxiety, and depression [123]. The gene Tyro3, which we found upregulated after FLX action (Tables 4, 7, 9 and S12), encodes the most widely expressed receptor protein tyrosine kinase and is linked to the PI3K/AKT pathway [124]. d Cyclic AMP It is noteworthy that among the genes common to the signatures of FLX activity in a stress or naive context is the Pde4b gene (Tables 6 and S6, Figure 3). It encodes a cyclic phosphodiesterase, regulating concentrations of cyclic nucleotides, and thereby plays a role in signal transduction. Not only have several studies specifically shown variations in the expression of transcripts and protein encoded by Pde4b in preclinical models and humans [125][126][127][128][129], but this gene is in fact the target of a drug with AD properties, rolipram. A link between neurotrophin synthesis and the production of second messengers like cAMP is protein kinase A. It turns out that one of the genes most reliably affected by FLX codes for a subunit of this kinase, Prkar1b (Tables 2-4  Perhaps a pivot of the effect of antidepressants on major depression through several previously mentioned signaling pathways, cAMP/PKA, phosphodiesterase, and neurotrophins, materializes through the activation of nuclear receptors, the glucocorticoid receptor (GR) and the mineralocorticoid (MR) [130][131][132][133][134]. Two genes share the task of regulating, among others, inflammatory responses, proliferation and differentiation processes, NR3C1 and NR3C2. In our data, if it seems that the tendency is rather to underexpress these two transcription factors, Nr3c2 is more consensually deregulated in the different signatures that we obtain (Tables S5, S8 and S12), without being one of the key markers either. On the other hand, considering the partners of the GR, we notice that AHI1, which regulates the nuclear translocation of the GR, and which has already been associated with stressinduced depressive behavior in mice [135], has its mRNA significantly underexpressed in the different signatures we have generated (Tables S5, S8 and S12). f Synaptic plasticity In a way, in fine, the beneficial action of SSRIs can be linked to a capacity to achieve synaptic plasticity and from this point of view it seems important to us to mention the decrease in expression of two mRNAs coding for transcription factors involved in brain development [136], Nfia and Nfib (Tables 2, 7, S5, S6, S8 and S11 and Figure 3). For the latter, a study in rats had shown that its expression variation was during the molecular signature at the level of the frontal cortex characterizing the AD action of quetiapine in a chronic stress model in rats [137]. Recently, a study on lymphoblastoid cell lines from depression patients treated with citalopram reported a significant association of NFIB expression with improvement in depression scale [138]. The NMDA subtype of glutamate receptors (NMDAR) plays a key role in synaptic plasticity in the context of depression [139], and it has been shown that a Ca 2+ /calmodulin-dependent Ras-guanine-nucleotide-releasing factor (RasGRF1) served as an NMDAR-dependent regulator of the ERK kinase pathway. In our meta-analysis, we noted that Rasgrf1 was one of the most consistent downregulated genes by FLX activity (Tables 6-8 and 10). Consistently, it has been shown that Rasgrf1 knockdown reversed the effect of UCMS on mice [140], and it has been shown in humans that RASGRF1 may be a potential specific biomarker of treatment response for bipolar disorder [141]. Concerning the NMDAR-positive modulation underlying AD action, one identified initial trigger of such an effect is Drd1-pyramidal cell signaling [142,143]. We uncovered a Drd1 gene upregulation induced by FLX in our meta-analysis (Tables 4, 7 and 9) in agreement with these observations. Another gene common to the different signatures of FLX activity is Pdlim5 (Table 6, Figure 3), which encodes a protein that tethers kinases to the Z-discs in striated muscles but also restrains the postsynaptic growth of excitatory synapses. Several investigations suggested that the human gene plays a key role in the pathophysiology of mood disorders but is also likely involved in AD response [144][145][146].
Overall, whatever the paradigm explored (stress vs. naïve), it appears that FLX has a marked effect on synaptic activity in the hippocampus (Tables 4, 8, S9 and S13). Thus, many genes characterizing the transcriptional signature of FLX code for proteins playing an important role at the synaptic level. This is the case for the Sorcs1 gene [147], encoding the sortilin-related VPS10 domain containing receptor 1, upregulated by FLX (Tables 6 and S11) as previously observed [148]. Another example is the Itsn1 gene (Tables 6, 7, 8, S9 and S12), encoding the intersectin 1 that regulates synaptic vesicle recycling [149]. Interestingly, we can even point out a link, perhaps not the most expected, between synapse and transcriptional activators. For example, the Pcdh19 gene (Tables 6, 7, S7 and S12), which codes for a protocadherin, bridges neuronal activity with gene expression by regulating immediateearly genes expression to favor maintenance of neuronal homeostasis [150]. In addition, Sox11 (Tables 7, S4, S6, S7 and S11) encodes a transcriptional factor playing a major role in brain development and activated in an activity-dependent fashion specific to the dentate gyrus of the hippocampus [151], which had already been shown to be dysregulated by FLX [118], as well as by paroxetine, another SSRI [100].

Emerging Pathways
A novelty in this meta-analysis is the zinc finger homeobox 2 gene, Zfhx2 (Tables 2 and 6, Figure 3), for which there is generally very little literature. This putative transcriptional regulator is highly expressed in the developing brain and, very interestingly, its global knockout has been reported to show several behavioral abnormalities, namely hyperactivity, enhanced depression-like behaviors, and an aberrantly altered anxiety-like phenotype [152], as well as hyposensitivity to noxious mechanical stimuli [153]. A human variant of ZFHX2 has been associated with a pain-insensitive phenotype [153]. Another gene related to pain and nociception is Trpm3 (Tables 7 and S12, Figure 3), encoding a tetrameric cation channel. It had been shown that FLX but also imipramine and the antipsychotic drug chlorpromazine inhibit the channel activity [154]. Among the genes whose relevance in terms of belonging to a biological pathway is not at first obvious to relate to the pathophysiology of depression or the action of ADs is Slc4a4 (Tables 7, S8 and S12). It encodes a sodium bicarbonate cotransporter (NBC) involved in the regulation of bicarbonate secretion and absorption and intracellular pH, but, interestingly, it has been proposed as one of the best biomarkers for predicting suicidal ideation [155,156]. Conversely, our analysis pinpoints to the oxytocin pathway altered in animals under stress exposure ( Table 3) that has received renewed attention in regard to psychiatric disorders [157].
Last but not least, components of the ribosome that we saw as particularly affected in stressed animals responding to FLX (Table 5) are also affected in the study of Hori et al., who explored the blood transcriptome of three homogeneous groups, namely "resilient", "vulnerable", and "resistant" [158]. Remarkably, among genes dysregulated between the different groups are RPL17, RPL21, RPL34, RPS15A, and RPS27, which were all dysregulated in rodent models (Table 5). In another study, Guilloux et al. were able to identify a group of six genes with predictive value for response to citalopram treatment for major depression [159]. Half of the genes were ribosomal units including again RPL17, as well as RPL24, that we found dysregulated (Table 5). . The filter "English language" was also activated. In addition, references from eligible studies were examined to include additional studies that would not be retrieved by the PubMed search. To compile a complete list of eligible studies (Table S1), we excluded the studies that did not use an AD or profile the action of an AD genome-wide. Because we focused on the in vivo action of AD in rodents as a model of major depression in humans, we also excluded studies performed in systems other than rodent. Namely, studies on other species, cells in culture, and on bacterial infection or brain injury models were excluded. We searched all eligible studies for either an ID number linked to the public download of full transcriptomic data or a table containing the complete data as supplementary information. We also searched NCBI resources, the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/, accessed on 19 April 2022), the Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra, accessed on 19 April 2022), and the European repository ArrayExpress (https://www.ebi.ac.uk/arrayexpress/, accessed on 19 April 2022) for microarray and high-throughput sequencing data related to published and unpublished studies by typing "antidepressant AND (mouse OR rat)". Finally, when full transcriptomic data were not publicly available, we emailed the corresponding authors to request raw data.

Individual Reprocessing of Published Experiments-Microarrays
Normalized microarrays from seven studies (GSE84185, GSE43261, GSE56028, GSE118669, GSE54307, GSE6476, and GSE42940) were retrieved from GEO using GEOquery [160]. One last microarray was retrieved directly from the authors [50] and was normalized using RMA from the oligo R package [161]. The microarrays' latest annotations were retrieved from Ensembl Biomart for Affymetrix GeneChip Mouse Genome 430 2.0, Affymetrix GeneChip Mouse 430A 2.0, Affymetrix Rat Gene 1.0 ST arrays, and Agilent SurePrint G3 Mouse Gene Expression 8 × 60 K microarrays. When unavailable, annotations from the original authors were completed using megablast queries (GSE42940). Probes with mapping to multiple genes were discarded. Differential gene expression analyses were produced with Limma [162].

Integration Analyses
To compare differential expression analyses between different experiments, we kept only genes with p values below 0.05 and in the top 300 of any compared designs. For RNA-Seq datasets, the log2 fold change values have been calculated using the "normal" shrinkage type and not the apeglm from Deseq2 because the latter was considered too stringent for this kind of integration analysis. Plots were produced using the R ggplot2 package [168]. To be able to specify the general trend of the direction of expression variation of a gene, we were inspired by the procedure defined in the MetaVolcano tool, the vote counting approach [169]. This consists of counting the number of comparisons with p values below 0.05 showing overexpression and subtracting the number of comparisons with p values below 0.05 corresponding to underexpression. The resulting integer, whether positive, zero, or negative, is called the consensus score.

Portraits of FLX Action
As an alternative procedure to extract the transcriptional signature of FLX treatment, we applied the method proposed by Stephen Gammie to generate a portrait of depression [23]. The p value for each gene and for each individual microarray or RNA-Seq dataset used in integration analyses was −log10 transformed and then multiplied by the sign of the direction of change. Then, for each dataset, both upregulated and downregulated genes were ranked from low to high. With the 1000-increment cutoff, for each gene we counted how many times it received a ranking of 1000 or less among the upregulated genes. The same approach was applied to count 2000 or less, but that time the count number was multiplied by 0.1. We used this method up to 8000 by multiplying the sum by 10 −(k−1) , where k is the number of the thousands of ranks considered. The same calculation was applied to downregulated genes. Finally, we calculated the sum of the above values for a given gene. Genes were sorted by the absolute value so that the final list contains all genes ordered by magnitude of the difference from controls with also information on direction of expression change.

Ontological Analysis
Gene lists were uploaded on DAVID (database for annotation, visualization and integrated discovery) Bioinformatics Resources (2021 Update) [170] for identifying statistically relevant biological processes. The default settings were kept unchanged and only functional annotations related to the gene ontology (GO) biological process, molecular process, and cellular component, as well as the KEGG pathway, with corrected p-value (Bonferroni) < 0.05, were retrieved. We limited the search to Mus musculus species annotation. To simplify the lists of biological processes obtained, we only kept, for identical gene lists, the process associated with the lowest p-value.

Conclusions
We provided the first meta-analysis of a transcriptional signature of FLX in animal models and highlighted the impact and pitfalls of gene expression variation studies. Future studies aiming at extending such signature to other widely prescribed ADs and considering other tissues and brain regions in both males and females are warranted to provide better insights into optimal targets for AD response in human.
Despite all the limitations mentioned before, our study shows how murine models are a valuable source for understanding the mechanisms involved in the response to antidepressants. Indeed, with the benefit of hindsight on the main processes regulated by FLX, we realize that there are solid points of convergence with what is described in the context of the follow-up of patients suffering from major depression. It would obviously be important to know whether the changes induced by FLX, in the context of an animal model recapitulating the response to AD treatment, are similar to the biological program that is set up in a depressed patient who responds favorably to AD treatment. The major concern is that all human studies based on post-mortem tissues cannot distinguish the death-induced signature from those induced by the treatment and those corresponding to the depressive illness. In contrast, studies that allow for clinical interviews at the same time the biological samples are taken, such as studies based on blood samples, can define the AD treatmentrelated signature in relation to response or non-response to treatment. In fact, studies on human cohorts have already underlined the importance of some biological processes significantly enriched from ontological analyses of genome-wide transcriptional variations in FLX-treated rodent models of major depression. That includes synaptic plasticity [171], even when the variations in the corticolimbic regions of the brain are of opposite directions between men and women [60,172], and in particular at the blood level with glutamatergic signaling, the intervention of neurotrophin [173,174], and MAPK signaling [74].
Overall, our results highlight activity-dependent gene transcription as the major AD target. This is a major finding with an important potential for future biomarker clinical research. At the center of cellular plasticity and resilience, in the CNS and the periphery, such cell signaling cascades are likely early-stage mediators of AD response body-wide. The findings of this meta-analysis allow us to therefore hypothesize that ADs do induce early changes in cellular biology and processes at the interface of CNS and periphery that could constitute clinically relevant biomarkers and suggest that more focus is needed in the earliest point of biomarker variations. Funding: This research was funded by ERA-NET NEURON ANTaRES (https://www.neuron-eranet. eu/projects/ANTaRES/) with reference ANR-19-NEUR-0005 to E.T.T., R.B., G.T., E.C.I., V.G. and P.O.-T., and by Gebra (Imaging GEne and BRAin networks in stress susceptibility: novel pathways and biomarkers in mice and men) with reference ANR-19-CE37-0017 to E.T.T., R.B., G.T., E.C.I. and V.G.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The source code for all the analyses is available at https://github.com/ guillaumecharbonnier/mw-ibrahim2022 and a compiled version of the Bookdown report is available at https://guillaumecharbonnier.github.io/mw-ibrahim2022.
Conflicts of Interest: E.S. is founder and acting chief officer of Damona Pharmaceuticals, a drug development company with small molecules in the pipeline for treatment of cognitive deficits across brain disorders and aging. E.T.T. is co-founder of Melkin Pharmaceuticals, a biotech company. Neither company had any role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. All other authors report no biomedical financial interests or potential conflicts of interest.