
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39300279
73042
10.1038/s41598-024-73042-2
Article
Brainstem transcriptomic changes in male Wistar rats after acute stress, comparing the use of duplex specific nuclease (DSN)
Lanshakov Dmitriy A. lanshakov@bionet.nsc.ru
dmitriylanshakov@yandex.com

12
Sukhareva Ekaterina V. 1
Bulygina Veta V. 3
Khozyainova Anna A. 4
Gerashchenko Tatiana S. 4
Denisov Evgeny V. 4
Kalinina Tatyana S. 23
1 grid.415877.8 0000 0001 2254 1834 Postgenomics Neurobiology Laboratory, Institute of Cytology and Genetics, Russian Academy of Science, Novosibirsk, Russian Federation
2 https://ror.org/04t2ss102 grid.4605.7 0000 0001 2189 6553 Natural Science Department, Novosibirsk State University, Novosibirsk, Russian Federation
3 grid.415877.8 0000 0001 2254 1834 Functional Neurogenomics Laboratory, Institute of Cytology and Genetics, Russian Academy of Science, Novosibirsk, Russian Federation
4 grid.415877.8 0000 0001 2254 1834 Laboratory of Cancer Progression Biology, Tomsk National Research Medical Center, Cancer Research Institute, Russian Academy of Sciences, Tomsk, Russian Federation
19 9 2024
19 9 2024
2024
14 2185625 12 2023
12 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
In this work, we have analyzed the transcriptomic changes in the brainstem of male Wistar rats 2 h after an acute stress exposure. We performed duplex-specific nuclease normalization of cDNA libraries and compared the results back-to-back for the first time. Based on our RNAseq data, we selected reference genes for RT-qPCR that are best suited for acute stress experiments. Most genes were upregulated. We detected a massive shift in neuropeptide Crh, Trh,Cga, Tshb, Uts2b, Tac4, Lep and neuropeptide receptor Hcrtr1, Sstr5, Bdkrb2, Crhr2 signaling, as well as glutamate Grin3b, Grm2 and GABA Gpr156, acetylcholine Chrm4,Chrne, adrenergic Adra2b receptors expression. A strong increase in the expression of intermediate filaments Krt83/Krt86/Krt80/Krt84/Krt87/Krt4/Krt76 and motor proteins Myo7a, Klc3 was detected. Remarkably, in the absence of astrocyte activation, we also observed signs of microglial activation at this time point. Both expression of anti-inflammatory cytokines Il13, Ccl24 and pro-inflammatory cytokine receptors Il9r, Il12rb1, Tnfrsf14, Tnfrsf13c, Tnfrsf25, Tnfrsf1b were increased. In the Wnt signaling pathway, we observed increased expression of ligands-receptors Wnt1, Wnt11, Ror2 and also negative regulators Notum, Sfrp5, Sost. RNAseq results after DSN treatment correlated at a high level with RNAseq results without DSN, but there was a proportion of genes that shifted their logFC values. They are mostly rare transcripts TPM 1–10 with higher 0.5–0.9 GC content.

Keywords

RNAseq
DSN
FST
Acute stress
Brainstem
Transcriptome
Subject terms

Molecular neuroscience
RNA sequencing
Neuroimmunology
Neural circuits
Gene expression profiling
http://dx.doi.org/10.13039/501100012190 Ministry of Science and Higher Education of the Russian Federation FWNR-2022-0002 issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Acute stress is an adaptive response by organisms to overcome adverse environmental conditions1,2. Over the course of evolution, the presence of predators has formed a stereotyped behavioral response, usually referred to as “fight or flight.” This behavioral response involves cooperative and fine-tuned neuronal activity in various brain structures. Notable among these is the brainstem—the structure responsible for autonomic control and reflexive responses to systemic and environmental stressors3. These functions are performed by several neuronal nuclei with different neurotransmitter systems in the brainstem2. Monoamine neuron nuclei are located in this brain structure. There is evidence that norepinephrine neurons in the locus coeruleus (LC) modulate cognitive functions such as memory and operant conditioning4. In acute stress conditions, the activity of LC norepinephrine neurons participates in the arousal response, along with the action of elevated levels of hormones5. Excessive increased activity of brainstem neurons causes oxidative stress, which may be a factor in the further development of various psychopathologies6.

Whole transcriptome analysis (RNA-seq) could be used to study gene expression changes in the brainstem after acute stress. The sequencing depth varies depending on the library preparation method. Reads from ribosomal RNA (rRNA) or highly expressed transcripts make up the majority of the total reads. Getting information about expression changes of rare transcripts requires deep sequencing, which can be expensive. Depletion of rRNA or purification of the polyA fraction is used to overcome this problem. Another approach that could be used is Duplex Specific Nuclease (DSN) treatment of the double-stranded cDNA library after denaturation and rehybridization7,8.

In this study, we performed brainstem RNA-seq analysis in Wistar rats 2 hours after acute stress, which in our case was a forced swim test. The rodent forced swim test has been used for decades and remains a classic model of inescapable stress. Reproducibility, robustness and widespread use by researchers are its strengths9,10. We chose the 2h time point to detect not only very rapid gene expression changes, but also secondary events that evolve over time, based on preliminary observations in the literature11. Major changes in gene expression were described along with functional annotations. DSN treatment of cDNA libraries was also performed, and we compared the RNAseq results of these two methods.

Results

Library characterization and consistency validation

As a model of acute stress, we used the classic forced swim test (FST) paradigm. The main experimental scheme was as follows (Fig. 1). Several stress models are used in rodent experiments. The forced swim test has been widely used since 1980 and is used not only to test antidepressants, but also to assess general stress response and level of adaptation12. We have previously shown that FST results in a rapid increase in plasma corticosterone levels, which remain elevated 2 h after the test when final test samples are collected for analysis13. The use of the forced swim test as a stressor is robust, reliable and reproducible. It allows us to compare the data obtained with our previous results.The approximate brainstem dissection scheme is shown on (Fig. 2) and described in the Methods section. In brief, the brainstem tissue block contains the caudal portions of the raphe nucleus and the nuclei of the norepinephrine neurons. We used template switching technology to generate double-stranded SMART-seq cDNA libraries (Fig. 1).Fig. 1 Schematic of the experiments. 2 h after FST, the brainstem region was isolated, cDNA libraries were prepared using SMART, template switching technology, and a fraction of the libraries were subjected to double-stranded nuclease normalization. All libraries were then sequenced on the NextSeq 500.

Fig. 2 Brainstem dissection schematic.

Duplex-specific nuclease from red king crab was used to normalize the cDNA libraries (trimming, i.e. DSN treatment). Here we used the term “normalization” as introduced by the inventors of the method8 and used by the kit manufacturer, but physically this process is more like rare transcript enrichment. After sequencing, you can clearly see the difference in the unnormalized library size and its reduction after DSN treatment (Fig. 3a). Interestingly, DSN untreated samples with FST had larger library sizes compared to the corresponding control, which could indicate a massive increase in transcription of a large number of genes (Fig. 3a). See the Methods section and Supplementary Materials for details on library preparations and concentrations. Principal component analysis revealed a sharp separation of the samples in the PCA dimensions (Fig. 3b)-stressed from unstressed, DSN treated from DSN untreated. Hierarchical clustering of the gene count data also showed a clear division of the samples into groups (Fig. 3c). Pearson’s correlation was lowest when comparing the sample without DSN treatment and the sample with treatment (R = 0.95, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{p} < 2.2{\textrm{e}}^{-16}$$\end{document}, n1_noFST vs tr1_noFST , Fig. 3d). Lin’s concordance coefficient of 0.05 confirms this more directly. In all n_vs_tr comparison pairs, Lin’s CC varied from 0.027 to 0.064. This is very low. On the graph (Fig. 3d), it is clearly shifted away from the trend line points. Comparing DSN untreated samples with and without stress also results in lower correlation coefficients (R = 0.96, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{p} <2.2{\textrm{e}}^{-16}$$\end{document}, n1_noFST vs n34_FST, Fig. 3d) compared to higher correlation coefficients of untreated samples from one experimental group (R = 0.98, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{p} < 2.2{\textrm{e}}^{-16}$$\end{document}, n1_noFST vs n34_FST, Fig. 3d). Lin’s CC were high in both cases, but for n1_vs_34 Lin’s CC = 0.99, for n1_vs_n3 Lin’s CC = 0.95.

Total RNA was used to prepare the libraries. Despite the fact that the oligo for the first strand synthesis contained oligo(dT) sequence, the majority of all reads were from rRNA in both DSN types of libraries (28S rRNA; mitochondrial rRNA AY172581.24, AY172581.9; 18S rRNA, Hsc-70-pseudogene and some others, Fig. 3e). Only 28S rRNA and AY172581.24 accounted for the majority of transcripts. Mitochondrial rRNA may be polyadenylated14. Non-abundant polyadenylated transcripts of 18S and 28S rRNA are also detected in human cell15. Although these are small rRNA species, their quantities per cell allow them to occupy the majority of the library. Among the protein-coding genes, cytochrome b, NADH dehydrogenase subunit 6, ATPase subunit 6, NDRG family member 4 were present in the 10 most expressed transcripts (Fig. 3e). Interestingly, the list of the 10 most abundant transcripts changed after DSN treatment, but still the majority consists of 28S rRNA, AY172581.24 with addition of Mt-cyb (Fig. 3e).

Functional annotation and transcriptomic alteration

For differentially expressed gene (DEG) analysis, genes with |log2FC| > 1 and Benjamini, Hochberg (“BH”) multiple comparison correction p.adj < 0.05 were selected. Acute stress generally causes an increase in gene expression, as seen in the volcano plots (q-value vs. log2 fold change) of untreated DSN samples (Fig. 4a, Supplementary Table 1). In the untreated samples (F.nT._vs_nF.nT) 1428 genes were upregulated (log2FC > 1) and only 45 genes were downregulated. In the samples with DSN treatment (F.T._vs_nF.T) observed a similar picture 1583 genes were upregulated and 179 were downregulated. An increased number of genes that passed the logFC selection criteria could be explained by the DSN treatment. It should be noted that 63% (1250) of genes were common between treated and untreated samples in the STRESS vs. CTRL comparison pair (Fig. 4b). Only 26% (512 genes) were unique for the DSN-treated STRESS vs. CTRL (F.T._vs_nF.T) comparison pair. Lower panel of volcano plots illustrates that DSN treatment depletes molecules from libraries and 60% (3031 genes) were common in this depletion. These images show many of the same genes with log2FC < 1 & p.adj “BH” < 0.05 in F.T_vs_F.nT and nF.T_vs_nF.nT. This means that these genes were present at higher levels in the untreated library and their levels decreased after DSN treatment. Pearson’s correlation of log2FC’s, i.e. data from treated and untreated samples, was about 80% (R = 0.78, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{p} < 2.2{\textrm{e}}^{-16}$$\end{document} Lin’s CC = 0.78 Fig. 4c). Comparing logFC in pairs with the same STRESS condition but with and without DSN treatment (nF.T_vs_nF.nT and F.T_vs_F.nT) gives a correlation of about 83% (R = 0.83, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{p} < 2,2{\textrm{e}}^{-16}$$\end{document}, Lin’s CC = 0.82 Fig. 4d). In general, log2FC data before and after DSN treatment correlated at a fairly high level (Fig. 4c), but there were a number of genes that changed their logFC after DSN treatment to completely opposite. We decide to try to clarify what is the cause of this shift and these distortion.

In comparison pairs with STRESS versus CTRL (F.nT_vs_nF.nT and F.T_vs_nF.T) we calculated log2FC rank before and after DSN treatment rank=||log2FC_DSN| − |log2FC_nDSN||. This parameter describes the shift of the obtained log2FC after DSN. Then we look at the dependencies of rank on transcript abundance in the library, GC content, and average cDNA length of the transcripts. At a glance, for the rare transcripts with TPM 1–15, log2FC is shifted more (Fig. 5a). GC content and mean cDNA length also influenced this process (Fig. 5b,c).It looked like the rare transcripts with a higher GC content of 0.5–0.9 were the most likely to be exposed to this effect. The cDNA length also affected this process (Fig. 5c) in the average range of 1000–5000 bp transcripts, but there were exceptions. Short transcripts 100–300 bp with high GC content were also mostly subject to log2FC shift. Further PCA analysis revealed that the variability is mostly described by GC content and mean transcript length, and to a lesser extent by TPM (Fig. 5d). Then we grouped genes according to log2FC rank, one rank unit is one group (Fig. 5e), and performed linear regression analysis, calculated estimated marginal means with “holm” adjusted multiple comparison correction (Fig. 5f–h). Only GC content and mean cDNA length were significantly different in the rank3 group compared to the rank1 group, p < 0.0001. At the same time, there is a shift, Lin’s CC = 0.34 (Fig. 5i), when we look at the transcript-TPM correlation without DSN and with DSN. Functional annotation and GO (Gene Ontology) enrichment analysis of samples without DSN treatment (comparison pair F.nT._vs_nF.nT) without log2FC cut-off condition (only p.adj < 0.05) showed the 5 most enriched GOs (Fig. 6, Table 1).Table 1 Top 10 enriched GO’s terms—biological process (BP).

GO.ID	Term	Annotated	Significant	Expected	p.adj	
GO:0043604	Amide biosynthetic process	898	571	433.88	1.8E−021	
GO:0044237	Cellular metabolic process	8489	4390	4101.53	2.8E−021	
GO:0006412	Translation	757	491	365.75	4.2E−21	
GO:0043043	Peptide biosynthetic process	779	503	376.38	5.4E−21	
GO:0006518	Peptide metabolic process	892	561	430.98	1.4E−19	
GO:0043603	Amide metabolic process	1161	707	560.95	2.5E−19	
GO:0006807	Nitrogen compound metabolic process	8580	4398	4145.5	9.1E−17	
GO:1901566	Organonitrogen compound biosynthetic process	1653	950	798.66	1.8E−15	
GO:0034641	Cellular nitrogen compound metabolic process	5456	2865	2636.11	5.8E−15	
GO:0044238	Primary metabolic process	9033	4596	4364.37	1.2E−14	

They were-GO:0043604 amide biosynthetic process, GO:0044237 cellular metabolic process, GO:0006412 translation, GO:0043043 peptide biosynthetic process, GO:0006518 peptide metabolic process. The scheme (Fig. 6) shows a hierarchy of the 5 most significant GO terms with upstream, more general terms that contain them. In general, these 5 terms reflect the intense cellular metabolic processes that take place in the brainstem after the arousal response. For subsequent GO analysis, genes were divided into UPregulated (logFC > 1, p.adj < 0.05) and DOWNregulated (logFC < − 1, p.adj < 0.05). GO enrichment in Biological Processes (BP) of DSN untreated samples comparison pairs showed that the majority of the upregulated genes were from GOs with overlapping sets of genes (Fig. 7, Supplementary Table 2)-GO:0007389 pattern specification process, GO:0003002 regionalization, GO:0008544 epidermis development, GO:0045165 cell fate commitment. It includes a diverse set of genes ranging from transcription factors (Dbx1 Tbx6 Lbx1 Pax2 Tcf7l1 Hoxa11 Hoxa1 Hoxa7) and Notch ligands (Dll3 Dll4) to integrin subunit alpha M (Itgam) and neurotrophin 4 (Ntf4). Mainly transcription factors and signaling molecules fulfill these GO terms (Figs. 7, 8a, Supplementary Table 2).

The genes enriched in GO:0030098 Lymphocyte Differentiation mostly contain genes strictly related to immune system (Ifnl1 Cd79a Ccr6 Ccr7 Clcf1 Cd19 Nfkbid Cd27 Lag3 Il9r Itk). Increased expression of these genes is evidence of microglial and possibly immune cell activation at some level 2 h after acute stress. So for Ifnl1 log2FC = 3.71 p.adj = 0.004492, for Ccr6 log2FC = 2.31 p.adj = 6.206451e−03, for Ccr7 log2FC = 2.48 p.adj = 0.0057, Clcf1 (ENSRNOG00000018752) log2FC = 2.31 p.adj = 0.048. Interestingly, the upregulated genes were enriched with the terms GO:0045109 intermediate filament organization, GO:0045104 intermediate filament cytoskeleton organization, GO:0045103 intermediate filament-based process (Figs. 7,  8). These terms consist of different keratin genes (Figs. 7,  8, Supplementary Table 2). This observation could be explained by massive intermediate filament and neuronal cytoskeletal reorganization after strong neurotransmitter release following acute stress16. The expression of motor proteins and vesicle cargo protein was also upregulated Myo7a log2FC = 3.06 p.adj = 0.0015, Klc3 log2FC = 2.50 p. adj = 0.0026, Tnni3 log2FC = 2.21 p.adj 0.031, Kif2c log2FC = 1.66 p.adj = 0.035, Tuba3a log2FC = 3.56 p.adj = 0.016 and some others. Enrichment of downregulated genes yields the following GO:0042445 hormone metabolic process (Igf2 Crabp2 Crhbp Cga Aldh1a2 Rbp1), GO:0046942 carboxylic acid transport (Slc22a2 Crabp2 Slc6a13 Trh Rbp1 Slc6a20), GO: 0016102 diterpenoid biosynthetic process (Crabp2 Aldh1a2 Rbp1), GO:0033189 response to vitamin A (Tshb Aldh1a2 Rbp1), GO:0051458 corticotropin secretion (Crh Crhbp) (Figs. 7,  8b, Supplementary Table 3). The GO enrichment of the comparison pairs (F.T._vs_nF.T) after DSN treatment gives similar results, but the Genes Ratio was higher in the same terms when compared without DSN treatment because more genes passed the logFC filter (Fig. 7, Supplementary Table 4–5).Fig. 3 Library characterization and data consistency validation. Initial analysis of the NGS data showed the size of the unnormalized libraries (a) and the sample distribution after normalization in pca dimensions (b, c) heatmap showing the distribution and hierarchical clustering of the samples (d) read count correlation in different samples n1, n3—untreated CTRL, n34—untreated FST, tr—respective DSN treated sample (e) TPM— transcripts per million, showing the 10 most abundant transcripts in the libraries.

Fig. 4 Differentially expressed genes (DEG) (a) volcano (q-value vs. log2 fold change) plots showing selected DEGs, cut-off filter |log2FC| > 1, p.adj < 0.05, in the names of the comparison pairs F = FST or noFST = nF, second letter indicates treatment T = DSN treated or nT = NON DSN treated (b) Venn diagrams showing the number of shared and separate genes from selected DEGs lists (c–e) Correlation of logFC data in different comparison pairs.

Fig. 5 (a–b) Dependence of logFC rank describing discrepancies in data by DSN from transcript representation, cDNA average length, and GC content (d) PCA analysis (e) grouping of genes according to log2FC rank. The box plots above show the number of genes in each group (f–h) Results of linear regression analysis with holm’s correction for multiple comparisons (i) correlation of gene’s TPM with and without DSN.

Fig. 6 Functional annotation of transcriptomic data using TopGO package in R. Five most significant GO terms (red rectangles) in Biological Processes from enrichment of all significant genes p.adj < 0.05 and hierarchy with more generalizing GO terms. Each box shows the GO term number, short description and p-value. The color of the box indicates statistical significance with yellow = 1 × \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-3}$$\end{document} orange = 1 × \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-5}$$\end{document} and red = 1 × \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-7}$$\end{document}, i.e. relative significance, ranging from dark red (most significant) to light yellow (least significant). Black arrows indicate is-a relationships and red arrows indicate part-of relationships.

Fig. 7 Enrichment analysis of DEGs based on clusterProfiler packages.Main Biological processes GO terms of |log2FC| > 1, p.adj “BH” < 0.05 genes with and without DSN treatment.

Fig. 8 Major genes and pathways altered after FST. (a) Gene concept network (cnetplot) of most enriched BP GO terms in UP regulated genes and overlapping set of genes consisting of these terms (b) gene concept network (cnetplot) of most enriched BP GO terms in DOWN regulated genes and overlapping set of genes consisting of these terms.

Alterations in neuroactive ligand–receptor interactions

To better understand the processes and pathways that were altered 2 h after acute stress in the brainstem, we performed KEGG17–19 pathway enrichment analysis (Supplementary Table 6). All KEGG pathway enrichment analysis supporting the entire data of this research with overlaid expression information could be founded at https://lda01.shinyapps.io/tr_app2/ as an interactive shiny application. Just wait for all the information on the page to load. Here we try to select KEGG pathways that cover a larger number of DEGs. One of the most interesting KEGG pathway terms was rno04080 Neuroactive ligand-receptor interaction, containing 33 DEGs (Fig. 9, Table 2).Fig. 9 Results of KEGG17–19 enrichment analysis for the pathway rno04080 “Neuroactive ligand-receptor interaction”. log2FC is displayed in color from green log2FC = − 5 to red log2FC = 5.

It is worth noting that the expression of neuropeptides Crh (corticotropin releasing hormone), Crhbp (corticotropin releasing hormone binding protein), Trh (thyrotropin releasing hormone), Tshb (thyroid stimulating hormone subunit beta), Cga (glycoprotein hormones, alpha polypeptide, on the Fig. 7 schema it was referred to FSH, LHB by the pathview package), Uts2b (urotensin 2B) decreased. With reduced expression of Crh and Crhbp, we observed increased expression of its receptor Crhr2, which could be evidenced by the compensatory regulatory loop between them.Table 2 log2FC for selected genes from rno04080.

Ensembl gene id	Rgd symbol	log2FC F.nT_vs_nF.nT	p value	p.adj	
ENSRNOG00000012703	Crh	− 1.28	0.0039	0.0124	
ENSRNOG00000017890	Crhbp	− 1.31	0.0002	0.0008	
ENSRNOG00000011824	Trh	− 1.41	0.001	0.006	
ENSRNOG00000016793	Tshb	− 2.4	0.01	0.027	
ENSRNOG00000009269	Cga	− 4.5	0.018	0.042	
ENSRNOG00000038512	Uts2b	− 1.26	0.0012	0.0047	
ENSRNOG00000011145	Crhr2	1.53	0.0081	0.022	
ENSRNOG00000006104	Tg	1.78	0.0144	0.035	
ENSRNOG00000004404	Tac4	1.63	0.00043	0.0019	
ENSRNOG00000045797	Lep	3.05	0.0077	0.021	
ENSRNOG00000009390	Edn2	3.84	0.0012	0.0048	
ENSRNOG00000068505	Insl3	2.94	0.00067	0.0028	
ENSRNOG00000029830	Adm2	2.76	0.0058	0.016	
ENSRNOG00000004488	Bdkrb1	3.95	0.00018	0.00096	
ENSRNOG00000047300	Bdkrb2	2.37	1.30E-05	0.0001	
ENSRNOG00000013838	Hcrtr1	1.05	0.00033	0.0016	
ENSRNOG00000018834	Sstr5	3.42	0.022	0.049	
ENSRNOG00000047457	Vipr1	4.09	0.014	0.035	
ENSRNOG00000004451	Mc3r	1.94	0.013	0.034	
ENSRNOG00000003880	Tph2	− 0.83	0.43	0.52	
ENSRNOG00000020410	Th	0.44	0.66	0.72	
ENSRNOG00000006641	Dbh	0.91	0.33	0.43	
ENSRNOG00000002797	Gpr156	2.46	0.00021	0.0010	
ENSRNOG00000013171	Grm2	1.68	6.74E−05	0.00041	
ENSRNOG00000012562	Grin3b	1.79	2.96E−06	2.98E−05	
ENSRNOG00000003777	Chrne	1.73	0.0053	0.015	
ENSRNOG00000017556	Chrm4	1.57	0.0009	0.0038	
ENSRNOG00000013887	Adra2b	3.33	5.65E–05	0.00035	
ENSRNOG00000049761	Htr6	2.01	0.016	0.039	
ENSRNOG00000061160	Htra4	2.50	0.019	0.044	
ENSRNOG00000002549	Htr5b	− 1.76	4.07E–06	3.9E–05	
ENSRNOG00000019270	P2ry6	1.26	0.0011	0.0043	
ENSRNOG00000017606	P2rx1	2.35	0.002	0.007	
ENSRNOG00000016756	Ptgir	3.29	0.019	0.045	
ENSRNOG00000004094	Ptger	1.15	4.32E–05	0.00028	
ENSRNOG00000018054	F2rl2	3.06	0.0062	0.0178	

Remarkably, increased expression of thyroglobulin Tg was observed. At the same time, the mRNA level of some other neuropeptide ligands was increased: Tac4 (tachykinin 4), Lep (leptin), Edn2 (endothelin 2), Insl3 (insulin-like 3, on the Fig. 7 scheme its expression was assigned to relaxin pathway), Adm2 (adrenomedullin 2, on the Fig. 7 scheme it was assigned to calcitonin pathway). 2 h after acute stress in the brainstem expression of neuropeptide receptors was increased. Thus, the mRNA level of bradykinin receptors (Bdkrb1, Bdkrb2), hypocretin receptor (Hcrtr1), somatostatin receptor (Sstr5), Vip receptor (vasoactive intestinal peptide receptor Vipr1), and melanocortin receptor (Mc3r) was increased. Less information was available on changes in the expression of receptors or enzymes that synthesize classical neurotransmitters. We did not detect expression changes of monoamines key synthesis enzymes (Dbh, Th, Tph2), possibly because the peak of its gene transcription after stress was earlier. At 2 h after stress we observed increased mRNA level of GABA, glutamate, acetylcholine and noradrenaline receptors—Gpr156, Grm2, Grin3b, Chrne, Chrm4, Adra2b. Expression shift of 5-HT receptors was controversial (because of this on Fig. 7 it expression shown in gray), with mRNA increase of Htr6 and Htra4 and mRNA decrease of Htr5b. It is also worth noting that pyrimidinergic receptor mRNA increases P2ry6, P2rx1. We also observed an increase in prostaglandin signaling, as reflected by an increase in Ptgir and Ptger1 mRNA. Pathview software also called F2rl2 coagulation factor II (thrombin) receptor-like 2 referred to proteinase-activated receptors, and its expression was also detected.

Changes in cytokine−cytokine receptor interactions

As mentioned after the Biological Processes GO enrichment analysis, we saw signs of microglial and immune system activation. Remarkably, at the same time there was no activation of acstrocytes at all Gfap, Vim, Aqp4 did not change mRNA levels (Table 3). To further elucidate this phenomenon, we performed KEGG pathway rno04060 Cytokine–cytokine receptor interaction (Fig. 10, Table 3) enrichment. It contained 30 DEGs.Fig. 10 Results of KEGG17–19 enrichment analysis of the pathway rno04060 “Cytokine–cytokine receptor interaction”. log2FC is color-coded from green log2FC = − 5 to red log2FC = 5.

If we look at the ligands, only one decreased expression Ccl19, another 10 increased their mRNA level—Ccl24, interleukin 13 (Il13). From the Il6/12-like family there were Clcf1, ENSRNOG00000018752 and Lif. From the class II helical cytokines IL10/28-like Ifnl3, interferon lambda 3 increased expression. From the prolactin family, it was colony stimulating factor 3 Csf3 and, as already mentioned, leptin. From the TGF-beta family they were growth differentiation factor 2 (Gdf2), bone morphogenetic protein 8b (Bmp8b), anti-Mullerian hormone (Amh). Of all the cytokine receptors, only Il13ra2 mRNA was reduced in our case. It could be a reaction to Il13’s increased expression. At the same time, the receptors of IL4-like family increased expression - colony stimulating factor 2 receptor subunit beta Csf2rb. The other chemokine and interleukin receptors were also upregulated Ccr6, Ccr7, Il9r, Il21r, Il22ra1, Il12rb1, Il27ra, Il17re. Of the TNF family receptors, only Tnfrsf11b decreased expression, Tnfrsf1b, Tnfrsf13c, Tnfrsf14, Tnfrsf25Cd27 increased its expression. Interestingly, from the prolactin family receptors Mpl - proto-oncogene thrombopoietin receptor mRNA level was increased.Table 3 log2FC for selected astrocytes activation marker genes and genes from rno04060.

Ensembl gene id	Rgd symbol	log2FC F.nT_vs_nF.nT	p value	p.adj	
ENSRNOG00000016043	Aqp4	− 0.26	0.35	0.44	
ENSRNOG00000002919	Gfap	0.18	2.10e−05	0.00015	
ENSRNOG00000018087	Vim	− 0.77	6.95E−06	6.08E−05	
ENSRNOG00000015668	Ccl19	− 1.18	0.0051	0.0154	
ENSRNOG00000031162	Ccl24	1.54	0.0131	0.0329	
ENSRNOG00000007652	Il13	3.6	0.0083	0.02259	
ENSRNOG00000018752	Clcf1	2.31	0.0211	0.048	
ENSRNOG00000007002	Lif	3.11	0.0016	0.0061	
ENSRNOG00000071068	Ifnl3	3.16	0.0026	0.0087	
ENSRNOG00000008525	Csf3	4.28	0.0102	0.0269	
ENSRNOG00000057751	Gdf2	6.49	0.0165	0.0395	
ENSRNOG00000033209	Bmp8b	3.07	0.0022	0.0078	
ENSRNOG00000019377	Amh	2.86	0.0003	0.0017	
ENSRNOG00000032973	Il13ra2	− 1.5	0.0032	0.0105	
ENSRNOG00000000187	Csf2rb	2.08	0.00077	0.0032	
ENSRNOG00000012964	Ccr6	2.31	0.0017	0.0062	
ENSRNOG00000010665	Ccr7	2.48	0.0015	0.0057	
ENSRNOG00000020630	Il9r	3.26	0.0122	0.031	
ENSRNOG00000015773	Il21r	2.12	0.0209	0.047	
ENSRNOG00000066331	Il22ra1	5.26	0.0018	0.0065	
ENSRNOG00000019216	Il12rb1	3.27	0.0005	0.0025	
ENSRNOG00000005747	Il27ra	3.73	0.00019	0.000099	
ENSRNOG00000009204	Il17re	2.63	0.0063	0.018	
ENSRNOG00000008336	Tnfrsf11b	− 1.07	3.69E−06	3.59E−05	
ENSRNOG00000016575	Tnfrsf1b	1.06	0.0039	0.0124	
ENSRNOG00000007664	Tnfrsf13c	3.72	0.00098	0.0039	
ENSRNOG00000013820	Tnfrsf14	3.35	0.0098	0.0259	
ENSRNOG00000021814	Tnfrsf25	2.14	0.011	0.030	
ENSRNOG00000027466	Cd27	2.59	0.016	0.039	
ENSRNOG00000028377	Mpl	2.23	0.0027	0.0091	

Wnt signaling pathway

Surprisingly, we observed significant changes in the Wnt signaling pathway (Fig. 11, Table 4). The mRNA of the canonical Wnt1 ligand was highly elevated, and the expression of the planar cell polarity ligand Wnt11 was also increased. At the same time, the expression of almost all Wnt signaling negative regulators Notum, Rspo1, Sfrp5, Apcdd1l was increased except for Serpinf1. The mRNA level of the noncanonical Wnt receptor Ror2 was also increased. Despite the increase in expression of both ligands and negative regulators, the expression of genes downstream of the cascade were also increased Tle7, Fosl1, Plcb2, Nfatc4.Fig. 11 Results of KEGG17–19 pathway rno04310 “Wnt signaling pathway”. log2FC shown in color from green log2FC = − 5 to red log2FC = 5.

Table 4 log2FC for selected genes from rno04310.

Ensembl gene id	Rgd symbol	log2FC F.nT_vs_nF.nT	p value	p.adj	
ENSRNOG00000061818	Wnt1	2.74	0.0027	0.0091	
ENSRNOG00000015982	Wnt11	1.82	0.0026	0.0088	
ENSRNOG00000036680	Notum	1.16	0.0028	0.0094	
ENSRNOG00000009656	Rspo1	1.54	0.019	0.044	
ENSRNOG00000014940	Sfrp5	1.86	0.0012	0.0047	
ENSRNOG00000028440	Apcdd1l	2.44	0.0099	0.026	
ENSRNOG00000003172	Serpinf1	− 1.03	0.0043	0.013	
ENSRNOG00000053232	Ror2	1.86	0.011	0.028	
ENSRNOG00000038741	Tle7	3.35	0.012	0.031	
ENSRNOG00000020552	Fosl1	1.11	0.0001	0.0006	
ENSRNOG00000058337	Plcb2	1.98	0.0064	0.018	
ENSRNOG00000020482	Nfatc4	1.64	0.0052	0.015	

Calcium signaling, cAMP response and MAPK pathways

As we expected, the arousal response results in a massive shift in a large number of signaling pathways. Enhancement of Calcium Signaling Shown in Voltage Gated Calcium Channels mRNA Elevation Cacna1s, Cacna1h, Cacna2d4 (Fig. 12, Table 5). It is worthy of note that the expression of another channel Kcnj4 has been increased. Expression of inner cell ER calcium release channel Ryr1 and inositol 1,4,5-trisphosphate receptor, type 3 Itpr3 was also increased. Increased expression of Mst1 and Erbb2 upstream MAPK kinases was observed. The mRNA level of intracellular MAPKs was increased - Mapk13, Map3k6, Map3k21, Map4k1 (Fig. 13, Table 5). An increased expression of the phosphatases Ptpn7 and Dusp4 was also observed. It could be explained by the triggering of compensatory mechanisms after massive neurotransmitter release and excitation. If we look at the cAMP signaling pathway, then the mRNA level of cAMP-dependent phosphodiesterase 4C Pde4c, troponin I3 Tnni3, which activates cAMP-dependent PKA, Creb3l3, Atp2a3 was also increased.Fig. 12 Results of KEGG17–19 enrichment analysis for the pathway rno04010 “MAPK signaling pathway”. log2FC is displayed in colors from green log2FC = − 5 to red log2FC = 5.

Fig. 13 Results of KEGG17–19 enrichment analysis of the pathway rno04024 “cAMP signaling pathway”. log2FC is displayed in colors from green log2FC = − 5 to red log2FC = 5.

Table 5 log2FC for selected genes from cAMP, Calcium and MAPK signaling pathways.

Ensembl gene id	Rgd symbol	log2FC F.nT_vs_nF.nT	p value	p.adj	
ENSRNOG00000046231	Cacna1s	3.09	0.0026	0.0088	
ENSRNOG00000033893	Cacna1h	1.72	2.06E−07	3.50E−06	
ENSRNOG00000008031	Cacna2d4	2.51	0.0038	0.0121	
ENSRNOG00000013869	Kcnj4	3.31	0.0075	0.020	
ENSRNOG00000020557	Ryr1	1.76	1.40E−06	1.61E−05	
ENSRNOG00000052795	Itpr3	1.77	1.71E−08	5.32E−07	
ENSRNOG00000019680	Mst1	4.25	1.48E−05	0.00011	
ENSRNOG00000006450	Erbb2	1.92	0.0009	0.0036	
ENSRNOG00000000515	Mapk13	1.83	0.0031	0.010	
ENSRNOG00000008936	Map3k6	3.27	1.44E−06	1.64E−05	
ENSRNOG00000019931	Map3k21	3.12	0.0034	0.011	
ENSRNOG00000020505	Map4k1	1.30	0.0050	0.015	
ENSRNOG00000005807	Ptpn7	1	0.012	0.03	
ENSRNOG00000011921	Dusp4	1.17	0.012	0.03	
ENSRNOG00000043249	Pde4c	3.17	0.0049	0.014	
ENSRNOG00000018250	Tnni3	2.21	0.012	0.031	
ENSRNOG00000032202	Creb3l3	3.71	0.0001	0.00061	
ENSRNOG00000017912	Atp2a3	1.29	0.0008	0.0034	

Selection of RT-qPCR reference genes

Prior to RT-qPCR validation, we decided to use our RNAseq data to select the most appropriate housekeeping reference genes for PCR expression estimation after acute stress. For this purpose, we used a method based on disjoint subgraphs20. First, we limit sets of selected genes with |log2FC| < 0.04 that were common in STRESS vs. CTRL comparison pairs (F.nT._vs_nF.nT and F.T._vs_nF.T) without and with DSN treatment. Then we add to this list commonly used reference genes such as GAPDH (G6pdx), Actb (Actx), Hprt (Hgprtase), Sdha and ribosomal proteins Rpl8, Rpl30, Rps13a, Rps16, Rps17. Due to discrepancies in different species databases for orthologous genes, there may be different symbols, which is why different gene symbols are shown on the images. We then performed analysis on the counts per million data. Analysis using the equivalence test procedure revealed that most of the nodes are not connected to edges (Fig. 14a). Genes that formed a maximum clique without DSN treatment were Frmd3, Sox18, Kctd10, Zfp617, Zfp498. Similar maximal clique were observed after DSN treatment, but branching on the dendrogram starts outside the p-value cutoff (Fig. 14b). Since we already preselected genes with |log2FC| < 0.04, we also performed an analysis with the ANOVA test (anva1.fpc), which detects which set of genes show different behavior between conditions (Fig. 15a,b). On the complement of the graph without DSN, you can see nodes of genes that behave more similarly—Trfr, TFIID, Sox18, Rps17, B3Galt4, Hgprtase, Rpl13a, Rps16, Frmd3, Kctd10. In libraries with DSN treatment, we observed a similar situation. Scatterplots (Fig. 15c) of CPM data generally confirm conclusions from graphs and dendrograms. For further pcr validation we selected Rpl13a, Rps16, Rps17, Rpl8, Rpl30, Hgprtase, B3galt4, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-actin.Fig. 14 Selection of reference genes for RT-qPCR using a method based on disjoint subgraphs. Complementary graphs and trees showing the analysis for selected genes in the SARP.compo package with equivalence function, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta$$\end{document} = 0.6 and p = 0.2 (a) graph complements (b) dendrogram representation of the graphs.

Real-time RT-qPCR analysis

The following set of reference genes was selected for RT-qPCR verification Rpl13a, Rps16, Rps17, Rpl8, Rpl30, Hgprtase, B3galt4, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-actin (Supplementary Table 8). PCR validation was performed on RNA samples from all experimental animals (Supplementary Tables 7, 8). Reverse transcription was performed in the usual manner using OligodT and MMLV reverse transcriptase. As expected, the Ct values of the reference genes correlate differently according to the selection method. The correlation matrix is shown in (Fig. 16a). As target genes, we selected genes with low expression levels or genes expressed in specific cell populations that did not meet the selection criteria for log2FC, p.adj - Esyt, Tcf7l1, Dcn (Fig. 16b). RNAseq data confirm that TPM values without DSN are at low levels (Fig. 16c).For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \Delta$$\end{document}Ct calculations, normalization was performed to the root mean square of all selected reference gene Ct’s (Supplementary Table 9). Welch two-sample t-test with “Bonferroni” adjustment for multiple comparisons revealed a 1.38-fold increase in Esyt1 mRNA (t = − 2.5556, df = 7.8702, p = 0. 03433) and approximately 1.61 and 1.38 fold decrease in expression of Dcn (t = 2.296, df = 9.4661, p = 0.04595) and Tcf7l1 (t = 3.5648, df = 7.4323, p = 0.008294). These observations confirm the log2FC data from our RNAseq experiments, except for Tcf7l1 (Fig. 16d–e). The choice of reference genes obtained with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \Delta$$\end{document}Ct was also verified with the SARP.compo graph-based approach and equivalence function (Fig. 16f–g, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta$$\end{document} = 0.6 and p = 0.2)21. On the graph complement (Fig. 16f) the maximum clique of nodes-set of genes that vary equivalently Act, B3galt4, Rps17, Rpl13a could be seen and confirmed on the dendrogamm (Fig. 16g).Fig. 15 Selection of reference genes for RT-qPCR using a method based on disjoint subgraphs. (a) analysis with the ANOVA test (anva1.fpc function)—the resulting graphs and their complements, on the graphs the nodes represent the RNA amounts of the covariating gene (b) dendrogram representation of the graphs (c) scatterplots of CPM data for selected genes.

Fig. 16 Results of RT-qPCR (a) Reference genes Ct’s Pearson correlation matrix (b) Expression of selected target genes in the brain stem, Image credit: Allen Institute for Brain Science (c) transcript abundance of target genes, TPM - transcripts per million (d) Fold change expression after RT-qPCR, Welch two-sample t-test *p < 0.05, **p < 0.01 with Bonferroni adjustment for multiple comparisons (e) log2FC data of selected target genes from RNAseq data (f–g) plots showing analysis of PCR experimental data in SARP.compo package, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta$$\end{document} = 0.6 and p = 0.2.

Discussion

The arousal reaction in response to acute stress primarily involves activation of brainstem and midbrain structures, as it is explicitly described in previous research works2,5,22. Interestingly, similar transcriptomic changes were observed in the hippocampus after acute stress, with mostly upregulation of gene transcription that was not detectable after 4h11. Acute stress also increased the phosphoproteome in the hippocampus. In our work, we took a time point between only the closest, rapid effects-minutes-and time point 4 h, where the effects are absent11. We did not detect significant changes in the key enzymes of monoamine synthesis (Dbh,Th, Tph), probably because the peak of its gene transcription was earlier, and they are more rapidly affected by FST. We observed increased expression of acetylcholine Chrm4, Adra2b, glutamate Grin3b, Gpr156 and GABA Gpr156. We have seen massive changes in the expression of neuropeptides and their receptors after acute stress. Roughly speaking, increase in receptor mRNA levels and decrease in peptide mRNA levels. For example, maintaining wakefulness is critical for the fight-or-flight response. Orexin signaling is responsible for maintaining wakefulness23,24, and we observed Hcrtr1 (orexin receptor) increased mRNA levels. A shift in this signaling balance could lead to sleep disturbances in PTSD or chronic stress25. Another observation that could be explained from an evolutionary point of view is the increase in tachykinin 4-substance P/enkephalin expression. Tac4 processing could produce hemokinin-1, a peptide with antinociceptive effects26. This could be an evolutionary adaptation to the likelihood of encountering predators under stressful conditions. In addition to nociception, Tac4 has been implicated in a variety of biological actions, including vasodilation, smooth muscle contraction, neurogenic inflammation, and immune system activation27.

Increased expression of bradykinin receptors serves to ensure local blood and oxygen supply to organs28. Classical concept of stress adaptation based on two major regulatory pathways: the hypothalamic-pituitary-adrenocortical (HPA) axis and the sympathetic adrenomedullary axis. The role of corticotropin-releasing hormone in the regulation of the HPA axis is well described. In addition, there are Crh-expressing neurons in various brain structures-cortex, amygdala, and especially midbrain and brainstem2. A clear understanding of the role of these neurons in stress response and adaptation to stress is in progress and requires further investigation. In our case, we observed an intensification of the Crhr2 signaling only by increasing its expression. Remarkably, in addition to Crh, significantly reduced expression of Cga, Trh, Tshb and Utsb2b reduced expression. These types of brainstem neurons are even less understood. As was mentioned, mainly decrease in peptides i.e. ligands and increase in receptors expression could be explained that wave of elevated transcription of ligands mRNA was early. Significant downregulation of retinoic acid and Igf2 pathways was observed, confirming previous results that retinoic acid signaling affects rats emotional behavior29.

Massive neurotransmitter release is observed in the arousal response. We observed widespread increased expression of intermediate filaments, keratins, cytoskeleton, and motor proteins at 2 h after FST. Increased expression of transcription factors that control the regulation of various cell programs was also observed at this time. Interestingly, we have seen and detected some signs of microglial or perhaps immune cell activation 2 hours after acute stress. This is primarily due to microglial activation, as we only detected increased Il13 mRNA levels, without IL4 or IL6 upregulation. During LPS-induced microglial activation, upregulation of proinflammatory cytokines Il4,Il6, Il1b,Tnfsf13b is observed30,31. Interleukin 13 is an anti-inflammatory cytokine that regulates microglia/macrophage polarization toward an anti-inflammatory phenotype and is expressed exclusively in activated microglia32. Another anti-inflammatory chemokine that was upregulated was Ccl24. This chemokine may be produced by anti-inflammatory microglia and can induce T cell differentiation of Th2 and Treg cells33. Microglia is a resident macrophage of the brain and a very unstable cell population. It has been shown that mice lacking Bdnf precisely in microglial cells have impaired memory and synaptic plasticity34. Microglial cells may also be involved in synaptic pruning35. In the hypothalamus and hippocampus, stress-induced microglial activation occurs through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}1-AR and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}2-AR adrenergic receptors36. We observed microglial activation without astrocyte activation. Notably, together with evidence of anti-inflammatory microglial activation, we observed increased expression of receptors for proinflammatory cytokines Il9r, Il12rb1 and TNF-receptors Tnfrsf14, Tnfrsf13c, Tnfrsf25,Tnfrsf1b and pleiotropic Il21r receptor. Such duality in activation of a pro- and anti-inflammatory pathway could be explained by compensatory activation of anti-inflammatory pathways in response to a predominantly inflammatory previous impulse. Interestingly, increased expression of IFN-lambda Ifnl3, one of the key cytokines in the innate antiviral defense, was observed after stress. This could not be explained at this time

Increased expression of the Tgf-beta family ligands Gdf2 (Bmp9) and Bmp8b could be explained by the need for high energy balance during arousal. Increased leptin levels could also be explained by this reason. Bmp8b has a thermogenic effect mediated by the inhibition of AMP-activated protein kinase (AMPK) in the ventromedial nucleus of the hypothalamus (VMH) and the subsequent increase in orexin OX signaling via the OX receptor 1 (OX1R) in the lateral hypothalamic area (LHA)37. It is possible that similar signaling occurs in the brainstem. Emotional stress is recognized as a major risk factor for cardiovascular disease and hypertension, which may be associated with anxiety and neurological disorders38. The link between anxiety, emotional stress, and hypertension may result from an altered neuronal firing rate in central pathways that control sympathetic activity38. In the results of our functional annotation analysis, we did not find genes associated with cell death, necrosis or apoptosis that could directly link acute stress to neurological damage. However, we did find promotion of neuroinflammation and a massive shift in neuropeptide signaling. Some recent clinical and preclinical evidence suggests that neuroinflammation is a key factor involved in the pathogenesis of major depressive disorder39, and in the pathogenesis of Parkinson’s disease40 and Alzheimer’s disease41. On the other hand, neuropeptides are powerful regulators of multiple physiological functions and behaviors. It appears that highly repetitive stress in particular could lead to neurological disorders resulting from the imbalance of neuropeptide signaling systems that also control blood pressure, and that increasing neuroinflammation accompanies or actively participates in this process. For example, melatonin is involved in sleep-wake cycle regulation, but it also has anti-inflammatory properties and is implicated in the pathophysiology of major depressive disorder (MDD)42. Agonists for Rxfr3 - receptor for another neuropeptide relaxin now considered a promising target for MDD pharmacology43.

Another signaling pathway activated after stress is Wnt. Wnt signaling regulates various aspects of brain development such as neural stem cell proliferation, neuronal differentiation, axon guidance, dendritogenesis, and synaptogenesis44. Recent studies have shown that the activity of the Wnt pathway can be critically regulated in response to neuronal activity45. Wnt signaling has been implicated in postsynaptic NMDA and GABA receptor clustering and synaptic plasticity. Our results regarding other cAMP, calcium signaling, and MAPK pathways are consistent with expectations and previous works.

Our RNAseq results were verified and confirmed by RT-qPCR. As target genes, we select genes with low expression levels or expressed in restricted cell populations. Decorin (Dcn) and Tcf7l1 are genes related to Wnt/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-catenin signaling. Tcf7l1 - transcription factor, is downstream on the pathway, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-catenin as a coactivator, it binds to HMG high mobility group domain proteins and affects transcription46. Decorin, an ECM molecule regulated by Wnt7a, could be characterized as a neurogenic factor and was downregulated after FST. Both Wnt pathway genes were downregulated according to our RT-qPCR results. Esyt1 is a candidate sensor protein for asynchronous neurotransmitter release47, and was upregulated. For the RT-qPCR, we performed a screen for suitable reference genes. According to the approach of Curist et. al. 201921, it is Act,B3galt4,Rps17,Rpl13a.

In general, results after DSN treatment correlated at a high level with RNAseq results without DSN, but there was a proportion of genes that shifted their logFC values. It is mostly transcripts with higher GC content and shorter average length. Representativeness of transcripts in library did not significantly affect log2FC shift. This method is more suitable for the normalization of screening libraries, i.e. two hybrid or expression libraries. If you just need to increase the representation of rare transcripts in it. In conclusion, in this work we analyzed transcriptomic changes and main pathways affected by acute stress in the male rat brainstem. Historically, the characterization of arousal as a response with a set of mutually exclusive and opposite states5 finds its reflection in the activation of opposite molecular cascades and pathways.

Methods

Animals and experimental design

All animal procedures were approved by the Institute of Cytology and Genetics Institutional Animal Care and Use Committee (Protocol 151 dated 28.04.2023) and were performed in accordance with the European Communities Council Directive 63/2010/EU and ARRIVE guidelines. Adult P90 male Wistar rats were obtained from the conventional vivarium of the Institute of Cytology and Genetics and housed in groups of 4 per cage at 22–24 °C, 12/12 light/dark cycle with free access to food and water.

Force swim test and sample collection

Twenty-four hours prior to any manipulation, the animals were housed in a single cage with free access to food and water. The day before, test animals (n = 6 per group) underwent a 15-min pretest in a 30 × 60 cm glass swimming cylinder filled with water t = 25 °C. The following day, a 5-min test session was performed using the same setup. After stress, the animals were placed in their single cage for 2 h. Control animals were left in their individual cages without pre-testing for the entire period prior to sacrifice. Animals were euthanized by rapid decapitation. Animals were quickly decapitated and regions dissected on ice as described48 and snap frozen in liquid nitrogen. Midbrain samples that were not included in the current study included the block of tissue from the rostral border of the superior colliculus to the rostral border of the pons and contained the DRN to approximately − 8.7 mm bregma (Fig. 2). The brainstem that was caudal to the midbrain region included the MRN and the remaining raphe complex located in the pons and medulla oblongata (Fig. 2). Brainstem samples, excluding the midbrain region and cerebellum, were then snap frozen in liquid nitrogen

Sequencing libraries preparation and DSN treatment

Total cellular RNA was isolated using a one-step acidic phenol extraction as described previously49. Only RNA samples with OD 260/280 ratios close to 2.0 were selected for further processing. RNA from two different animals was pooled for a final amount of 1 μg of each sample according to Supplementary Tables 7, 9 and selected for RNA-Seq and subsequent DSN treatment. Each sample was a different biological replicate. Samples for DSN treatment were derived from the same single-stranded cDNA, just another replicate of the same biological experimental point amplified the same number of cycles as for DSN untreated samples. cDNA synthesis and amplification to generate double-stranded libraries were performed using the Mint cDNA Synthesis Kit (Evrogen, SK001). Number of cycles required-10, for amplification was determined using sybr-green real-time PCR (Supplementary Material). DSN treatment and cDNA normalization were performed using the Trimmer-2 cDNA normalization kit (Evrogen, NK003). In brief, 1.2 μg of the double-stranded library was used as the starting material in a reaction volume of 16 μl. Denaturation was performed at 98 °C for 2 min. The samples were then hybridized at 68 °C for 5h. DSN treatment was performed with 0.25U of enzyme for 25 min at 68 °C. The double-stranded library was then synthesized with 11 amplification cycles. In the Supplementary Materials, concentrations and banding profiles of double-stranded libraries and libraries after DSN treatment are presented. The concentration of DSN libraries was much lower even after 11 cycles of amplification. To obtain sufficient DNA after DSN treatment for sequencing libraries, 2 reactions were set up, pooled, and precipitated. For NGS, libraries were ultrasonically sheared to 200 bp fragments. Adapters and indices were ligated using Kappa HyperPrep Kit (Roche). The concentration of cDNA libraries was measured using the dsDNA High Sensitivity Kit on a Qubit 4.0 fluorometer (Thermo Fisher Scientific). The quality of cDNA libraries was assessed using High Sensitivity D1000 ScreenTape on a 4150 TapeStation (Agilent). Libraries were sequenced on a NextSeq 500 instrument (Illumina) using single end 75 bp reads with 20M reads per sample after DSN treatment and 80M reads per sample without DSN treatment.

Bioinformatics analysis

Quality of reads was assessed using FASTQC v0.11.8 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Adapter sequences were trimmed with Trimmomatic50 v0.36 http://www.usadellab.org/cms/?page=trimmomatic. Raw reads were aligned to the Rnor7.2 reference genome using Hisat251 v2.2.1 http://daehwankimlab.github.io/hisat2. The number of mapped reads was counted using the featureCounts tool from the Rsubread package52 2.16.1 https://bioconductor.org/packages/3.18/bioc/html/Rsubread.html. Library counts, normalization, and GLM differential gene expression were calculated using Edger53 4.0.16 https://bioinf.wehi.edu.au/edgeR/. Functional annotation and visualization was done using clusterprofiler54 4.10.1 https://github.com/GuangchuangYu/clusterProfiler/issues and pathview55 1.42.0 https://github.com/datapplab/pathview packages. Reference genes were selected based on counts per million data using the SARP.compo package20,21 0.1.8 https://CRAN.R-project.org/package=SARP.compo and p = 0.05 cutoff for graphs plotting. For the analysis of Ct’s pcr data p = 0.23 cutoff were selected based on simulation proposed by Curis et.al. The Rscript used in the study can be found in the Supplement.

Real-time RT-qPCR analysis

A single-step acidic phenol extraction was used to isolate total cellular RNA as described previously1,49,56. Reverse transcription was performed using MMLV Reverse Transcriptase (Sibenzyme). 1 μg of total RNA was reverse transcribed using 100U MMLV Reverse Transcriptase (Sibenzyme), 1 mM dNTP, 2 mM DTT, 1 μM OligodT primer (Evrogen, SB001) and standard thermocycler temperature conditions for MMLV. All real-time PCR reactions were performed on the ABI ViiA7 System (Thermo). Amplification was performed using the real-time 5X qPCRmix-HS SYBR+LowROX (Evrogen, PK156S) with primers and Taqman probes (actin) from Table 6. The cycling conditions were as follows: 92 °C −3 min, 94 °C−10 s, Tm °C (from Table 6)−15 s, 72 °C−10 s for 40 cycles. Product sizes were checked by agarose gel electrophoresis (Fig. 16, original photo is in the Supplementary Material). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \Delta$$\end{document}Ct was calculated with the root mean square of all selected reference genes Ct’s in Excel, then in R tested for normality with the Shapiro-Wilk test. Group effect was tested with Welch Two Sample t-test with “bonferroni” multiple-testing correction. Significant considered values with p < 0.05. Ct’s and fold change data could be found in Supplementary Tables 7, 8 respectively.Fig. 17 Product sizes were checked with agarose gel electrophoresis.

Table 6 Oligo sequence information.

Name	Seqence 5′–3′	Product size, bp	Tm ,°C	
Rps17F	AGACTGTGAAGAAGGCGGC	81	64	
Rps17R	CGCGCTTGTTGGTGTGAAA	81	64	
Rps16F	AGAGCTGTGCCTGCTTCTG	180	64	
Rps16R	CTACATGCCTGGGGTTGGG	180	64	
Rpl8F	CTGTTTGGGGGCTCTCTGG	163	62	
Rpl8R	TAGGATGGTCACCTGCGGA	163	62	
Rpl13aF	AGGTGGTGGTTGTACGCTG	109	62	
Rpl13aR	CGAGACGGGTTGGTGTTCA	109	62	
B3galt4F	TGTGGGAGTAAGCGCAAGG	240	62	
B3galt4R	ATGAACCGGCACCTCAAGG	240	62	
Rpl30F	TGCATGCATGGTCCCCATT	214	64	
Rpl30R	GTTCCTGGGACCCAAGAGC	214	64	
HgprtaseF	ATGTCGACCCTCAGTCCCA	209	62	
HgprtaseR	CCCTTCAGCACACAGAGGG	209	62	
Esyt1F	CTGACTTCTCGCCTTGCCTC	299	61.7	
Esyt1R	GCTGTATTTGGCATGAGCTACG	299	61.7	
Tcf7l1F	CCCCTTGTGCTGTGGTAGG	75	61.7	
Tcf7l1R	GGAGGGCAGAGAGTCAGGA	75	61.7	
DcnF	TGGAGCCTTGCAGGGAATG	148	63.8	
DcnR	TTTCAGGCTGGCTGCATCA	148	63.8	
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-actin	Rn00667869_m1(Thermo)	–	60	

Supplementary Information

Below is the link to the electronic supplementary material. Supplementary Legends.

Supplementary Information.

Supplementary Figures.

Supplementary Table 1.

Supplementary Table 2.

Supplementary Table 3.

Supplementary Table 4.

Supplementary Table 5.

Supplementary Table 6.

Supplementary Table 7.

Supplementary Table 8.

Supplementary Table 9.

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-73042-2.

Acknowledgements

Studies was supported by joint research project FWNR-2022-0002. RNA-sequencing was carried out in The Core Facility «Medical genomics»(Tomsk NRMC) and the Tomsk Regional Common Use Center. Bioinformatics analysis was carried out at the Information and Computing Center of Novosibirsk State University.

Author contributions

D.A.L. conceived the project, designed, and performed the experiments, and evaluated the data; E.V.S., T.S.K.—RNA isolation, cDNA synthesis for RT-qPCR, qPCR, Mint cDNA libraries preparation, E.V.S., T.S.K., D.A.L., V.V.B.—animal testing, samples collection, D.A.L.—Mint libraries preparations, DSN treatment, D.A.L.—Bioinformatics analysis, A.A.K., T.S.G., E.V.D.—KAPA cDNA library preparation, pooling, and sequencing.

Data availibility

Raw reads and sorted bam files were deposed to CNGB https://db.cngb.org/ under the project number CNP0004878. Rscript and other findings of the paper are included as supplementary tables. PCR results are also included as supplementary tables. KEGG pathway graphs with overlaid expression information could be founded at https://lda01.shinyapps.io/tr_app2/ as interactive shiny application. Just wait until all information will be loaded on the page

Declarations

Competing interests

The authors declare no conflict of interest.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Lanshakov DA Single neonatal dexamethasone administration has long-lasting outcome on depressive-like behaviour, Bdnf, Nt-3, p75ngfr and sorting receptors (SorCS1-3) stress reactive expression Sci. Rep. 2021 11 8092 10.1038/s41598-021-87652-7 33854153
Lanshakov, D. A. et al. Single neonatal dexamethasone administration has long-lasting outcome on depressive-like behaviour, Bdnf, Nt-3, p75ngfr and sorting receptors (SorCS1-3) stress reactive expression. Sci. Rep. 11, 8092. 10.1038/s41598-021-87652-7 (2021).33854153
2. Chaves T Stress Adaptation and the Brainstem with Focus on Corticotropin-Releasing Hormone Int. J. Mol. Sci. 2021 22 9090 10.3390/ijms22169090 34445795
Chaves, T. et al. Stress Adaptation and the Brainstem with Focus on Corticotropin-Releasing Hormone. Int. J. Mol. Sci. 22, 9090. 10.3390/ijms22169090 (2021).34445795
3. Sattin D Leonardi M Picozzi M The autonomic nervous system and the brainstem: A fundamental role or the background actors for consciousness generation? Hypothesis, evidence, and future directions for rehabilitation and theoretical approaches Brain Behav. 2020 10 e01474 10.1002/brb3.1474 31782916
Sattin, D., Leonardi, M. & Picozzi, M. The autonomic nervous system and the brainstem: A fundamental role or the background actors for consciousness generation? Hypothesis, evidence, and future directions for rehabilitation and theoretical approaches. Brain Behav. 10, e01474. 10.1002/brb3.1474 (2020).31782916
4. Giustino TF Maren S Noradrenergic modulation of fear conditioning and extinction Front. Behav. Neurosci. 2018 12 43 10.3389/fnbeh.2018.00043 29593511
Giustino, T. F. & Maren, S. Noradrenergic modulation of fear conditioning and extinction. Front. Behav. Neurosci. 12, 43. 10.3389/fnbeh.2018.00043 (2018).29593511
5. Ross JA Van Bockstaele EJ The locus coeruleus- norepinephrine system in stress and arousal: Unraveling historical, current, and future perspectives Front. Psych. 2021 11 601519 10.3389/fpsyt.2020.601519
Ross, J. A. & Van Bockstaele, E. J. The locus coeruleus- norepinephrine system in stress and arousal: Unraveling historical, current, and future perspectives. Front. Psych. 11, 601519. 10.3389/fpsyt.2020.601519 (2021).
6. Chaoui N Long lasting effect of acute restraint stress on behavior and brain anti-oxidative status AIMS Neurosci. 2022 9 57 75 10.3934/Neuroscience.2022005 35434276
Chaoui, N. et al. Long lasting effect of acute restraint stress on behavior and brain anti-oxidative status. AIMS Neurosci. 9, 57–75. 10.3934/Neuroscience.2022005 (2022).35434276
7. Yi H Duplex-specific nuclease efficiently removes rRNA for prokaryotic RNA-seq Nucleic Acids Res. 2011 39 e140 10.1093/nar/gkr617 21880599
Yi, H. et al. Duplex-specific nuclease efficiently removes rRNA for prokaryotic RNA-seq. Nucleic Acids Res. 39, e140. 10.1093/nar/gkr617 (2011).21880599
8. Zhulidov PA Simple cDNA normalization using kamchatka crab duplex-specific nuclease Nucleic Acids Res. 2004 32 e37 10.1093/nar/gnh031 14973331
Zhulidov, P. A. et al. Simple cDNA normalization using kamchatka crab duplex-specific nuclease. Nucleic Acids Res. 32, e37. 10.1093/nar/gnh031 (2004).14973331
9. Commons KG Cholanians AB Babb JA Ehlinger DG The rodent forced swim test measures stress-coping strategy. Not depression-like behavior. ACS Chem. Neurosci. 2017 8 955 960 10.1021/acschemneuro.7b00042 28287253
Commons, K. G., Cholanians, A. B., Babb, J. A. & Ehlinger, D. G. The rodent forced swim test measures stress-coping strategy. Not depression-like behavior.. ACS Chem. Neurosci. 8, 955–960. 10.1021/acschemneuro.7b00042 (2017).28287253
10. Musazzi L Acute inescapable stress rapidly increases synaptic energy metabolism in prefrontal cortex and alters working memory performance Cereb. Cortex 2019 29 4948 4957 10.1093/cercor/bhz034 30877789
Musazzi, L. et al. Acute inescapable stress rapidly increases synaptic energy metabolism in prefrontal cortex and alters working memory performance. Cereb. Cortex 29, 4948–4957. 10.1093/cercor/bhz034 (2019).30877789
11. Von Ziegler LM Multiomic profiling of the acute stress response in the mouse hippocampus Nat. Commun. 2022 13 1824 10.1038/s41467-022-29367-5 35383160
Von Ziegler, L. M. et al. Multiomic profiling of the acute stress response in the mouse hippocampus. Nat. Commun. 13, 1824. 10.1038/s41467-022-29367-5 (2022).35383160
12. Molendijk ML De Kloet ER Forced swim stressor: Trends in usage and mechanistic consideration Eur. J. Neurosci. 2022 55 2813 2831 10.1111/ejn.15139 33548153
Molendijk, M. L. & De Kloet, E. R. Forced swim stressor: Trends in usage and mechanistic consideration. Eur. J. Neurosci. 55, 2813–2831. 10.1111/ejn.15139 (2022).33548153
13. Shishkina GT Kalinina TS Berezova IV Bulygina VV Dygalo NN Resistance to the development of stress-induced behavioral despair in the forced swim test associated with elevated hippocampal Bcl-xl expression Behav. Brain Res. 2010 213 218 224 10.1016/j.bbr.2010.05.003 20457187
Shishkina, G. T., Kalinina, T. S., Berezova, I. V., Bulygina, V. V. & Dygalo, N. N. Resistance to the development of stress-induced behavioral despair in the forced swim test associated with elevated hippocampal Bcl-xl expression. Behav. Brain Res. 213, 218–224. 10.1016/j.bbr.2010.05.003 (2010).20457187
14. Baserga SJ Polyadenylation of a human mitochondrial ribosomal RNA transcript detected by molecular cloning Gene 1985 35 305 312 10.1016/0378-1119(85)90009-5 4043734
Baserga, S. J. et al. Polyadenylation of a human mitochondrial ribosomal RNA transcript detected by molecular cloning. Gene 35, 305–312. 10.1016/0378-1119(85)90009-5 (1985).4043734
15. Slomovic S Laufer D Geiger D Schuster G Polyadenylation of ribosomal RNA in human cells Nucleic Acids Res. 2006 34 2966 2975 10.1093/nar/gkl357 16738135
Slomovic, S., Laufer, D., Geiger, D. & Schuster, G. Polyadenylation of ribosomal RNA in human cells. Nucleic Acids Res. 34, 2966–2975. 10.1093/nar/gkl357 (2006).16738135
16. Margiotta A Bucci C Role of Intermediate Filaments in Vesicular Traffic Cells 2016 5 20 10.3390/cells5020020 27120621
Margiotta, A. & Bucci, C. Role of Intermediate Filaments in Vesicular Traffic. Cells 5, 20. 10.3390/cells5020020 (2016).27120621
17. Kanehisa M KEGG: Kyoto encyclopedia of genes and genomes Nucleic Acids Res. 2000 28 27 30 10.1093/nar/28.1.27 10592173
Kanehisa, M. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. 10.1093/nar/28.1.27 (2000).10592173
18. Kanehisa M Toward understanding the origin and evolution of cellular organisms Protein Sci. 2019 28 1947 1951 10.1002/pro.3715 31441146
Kanehisa, M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 28, 1947–1951. 10.1002/pro.3715 (2019).31441146
19. Kanehisa M Furumichi M Sato Y Kawashima M Ishiguro-Watanabe M KEGG for taxonomy-based analysis of pathways and genomes Nucleic Acids Res. 2023 51 D587 D592 10.1093/nar/gkac963 36300620
Kanehisa, M., Furumichi, M., Sato, Y., Kawashima, M. & Ishiguro-Watanabe, M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 51, D587–D592. 10.1093/nar/gkac963 (2023).36300620
20. Curis E Determination of sets of covariating gene expression using graph analysis on pairwise expression ratios Bioinformatics 2019 35 258 265 10.1093/bioinformatics/bty629 30010788
Curis, E. et al. Determination of sets of covariating gene expression using graph analysis on pairwise expression ratios. Bioinformatics 35, 258–265. 10.1093/bioinformatics/bty629 (2019).30010788
21. Curis E Selecting reference genes in RT-qPCR based on equivalence tests: A network based approach Sci. Rep. 2019 9 16231 10.1038/s41598-019-52217-2 31700128
Curis, E. et al. Selecting reference genes in RT-qPCR based on equivalence tests: A network based approach. Sci. Rep. 9, 16231. 10.1038/s41598-019-52217-2 (2019).31700128
22. Libersat F Pflueger H-J Monoamines and the orchestration of behavior Bioscience 2004 54 17 10.1641/0006-3568(2004)054[0017:MATOOB]2.0.CO;2
Libersat, F. & Pflueger, H.-J. Monoamines and the orchestration of behavior. Bioscience 54, 17. 10.1641/0006-3568(2004)054[0017:MATOOB]2.0.CO;2 (2004).
23. Alexandre C Andermann ML Scammell TE Control of arousal by the orexin neurons Curr. Opin. Neurobiol. 2013 23 752 759 10.1016/j.conb.2013.04.008 23683477
Alexandre, C., Andermann, M. L. & Scammell, T. E. Control of arousal by the orexin neurons. Curr. Opin. Neurobiol. 23, 752–759. 10.1016/j.conb.2013.04.008 (2013).23683477
24. Inutsuka A Yamanaka A The physiological role of orexin/hypocretin neurons in the regulation of sleep/wakefulness and neuroendocrine functions Front. Endocrinol. 2013 10.3389/fendo.2013.00018
Inutsuka, A. & Yamanaka, A. The physiological role of orexin/hypocretin neurons in the regulation of sleep/wakefulness and neuroendocrine functions. Front. Endocrinol.[SPACE]10.3389/fendo.2013.00018 (2013).
25. Kaplan GB Lakis GA Zhoba H Sleep-wake and arousal dysfunctions in post-traumatic stress disorder: Role of orexin systems Brain Res. Bull. 2022 186 106 122 10.1016/j.brainresbull.2022.05.006 35618150
Kaplan, G. B., Lakis, G. A. & Zhoba, H. Sleep-wake and arousal dysfunctions in post-traumatic stress disorder: Role of orexin systems. Brain Res. Bull. 186, 106–122. 10.1016/j.brainresbull.2022.05.006 (2022).35618150
26. Fu C-Y Tang X-L Yang Q Chen Q Wang R Effects of rat/mouse hemokinin-1, a mammalian tachykinin peptide, on the antinociceptive activity of pethidine administered at the peripheral and supraspinal level Behav. Brain Res. 2007 184 39 46 10.1016/j.bbr.2007.06.019 17675256
Fu, C.-Y., Tang, X.-L., Yang, Q., Chen, Q. & Wang, R. Effects of rat/mouse hemokinin-1, a mammalian tachykinin peptide, on the antinociceptive activity of pethidine administered at the peripheral and supraspinal level. Behav. Brain Res. 184, 39–46. 10.1016/j.bbr.2007.06.019 (2007).17675256
27. Page NM Characterization of the endokinins: Human tachykinins with cardiovascular activity Proc. Natl. Acad. Sci. 2003 100 6245 6250 10.1073/pnas.0931458100 12716968
Page, N. M. et al. Characterization of the endokinins: Human tachykinins with cardiovascular activity. Proc. Natl. Acad. Sci. 100, 6245–6250. 10.1073/pnas.0931458100 (2003).12716968
28. Golias C Charalabopoulos A Stagikas D Charalabopoulos K Batistatou A The kinin system–bradykinin: Biological effects and clinical implications. Multiple role of the kinin system–bradykinin Hippokratia 2007 11 124 128 19582206
Golias, C., Charalabopoulos, A., Stagikas, D., Charalabopoulos, K. & Batistatou, A. The kinin system–bradykinin: Biological effects and clinical implications. Multiple role of the kinin system–bradykinin. Hippokratia 11, 124–128 (2007).19582206
29. Zhang Y Manipulation of retinoic acid signaling in the nucleus accumbens shell alters rat emotional behavior Behav. Brain Res. 2019 376 112177 10.1016/j.bbr.2019.112177 31449909
Zhang, Y. et al. Manipulation of retinoic acid signaling in the nucleus accumbens shell alters rat emotional behavior. Behav. Brain Res. 376, 112177. 10.1016/j.bbr.2019.112177 (2019).31449909
30. Babenko VN Shishkina GT Lanshakov DA Sukhareva EV Dygalo NN LPS administration impacts glial immune programs by alternative splicing Biomolecules 2022 12 277 10.3390/biom12020277 35204777
Babenko, V. N., Shishkina, G. T., Lanshakov, D. A., Sukhareva, E. V. & Dygalo, N. N. LPS administration impacts glial immune programs by alternative splicing. Biomolecules 12, 277. 10.3390/biom12020277 (2022).35204777
31. Shishkina GT Genes involved by dexamethasone in prevention of long-term memory impairment caused by lipopolysaccharide-induced neuroinflammation Biomedicines 2023 11 2595 10.3390/biomedicines11102595 37892969
Shishkina, G. T. et al. Genes involved by dexamethasone in prevention of long-term memory impairment caused by lipopolysaccharide-induced neuroinflammation. Biomedicines 11, 2595. 10.3390/biomedicines11102595 (2023).37892969
32. Chen D Interleukin 13 promotes long-term recovery after ischemic stroke by inhibiting the activation of STAT3 J. Neuroinflammation 2022 19 112 10.1186/s12974-022-02471-5 35578342
Chen, D. et al. Interleukin 13 promotes long-term recovery after ischemic stroke by inhibiting the activation of STAT3. J. Neuroinflammation 19, 112. 10.1186/s12974-022-02471-5 (2022).35578342
33. Xu Y The reciprocal interactions between microglia and T cells in Parkinson’s disease: A double-edged sword J. Neuroinflammation 2023 20 33 10.1186/s12974-023-02723-y 36774485
Xu, Y. et al. The reciprocal interactions between microglia and T cells in Parkinson’s disease: A double-edged sword. J. Neuroinflammation 20, 33. 10.1186/s12974-023-02723-y (2023).36774485
34. Parkhurst CN Microglia promote learning-dependent synapse formation through brain-derived neurotrophic factor Cell 2013 155 1596 1609 10.1016/j.cell.2013.11.030 24360280
Parkhurst, C. N. et al. Microglia promote learning-dependent synapse formation through brain-derived neurotrophic factor. Cell 155, 1596–1609. 10.1016/j.cell.2013.11.030 (2013).24360280
35. Paolicelli RC Synaptic pruning by microglia is necessary for normal brain development Science 2011 333 1456 1458 10.1126/science.1202529 21778362
Paolicelli, R. C. et al. Synaptic pruning by microglia is necessary for normal brain development. Science 333, 1456–1458. 10.1126/science.1202529 (2011).21778362
36. Sugama S Stress-induced microglial activation occurs through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}-adrenergic receptor: Noradrenaline as a key neurotransmitter in microglial activation J. Neuroinflamm. 2019 16 266 10.1186/s12974-019-1632-z
Sugama, S. et al. Stress-induced microglial activation occurs through -adrenergic receptor: Noradrenaline as a key neurotransmitter in microglial activation. J. Neuroinflamm. 16, 266. 10.1186/s12974-019-1632-z (2019).
37. Martins L A Functional Link between AMPK and Orexin Mediates the Effect of BMP8B on Energy Balance Cell Rep. 2016 16 2231 2242 10.1016/j.celrep.2016.07.045 27524625
Martins, L. et al. A Functional Link between AMPK and Orexin Mediates the Effect of BMP8B on Energy Balance. Cell Rep. 16, 2231–2242. 10.1016/j.celrep.2016.07.045 (2016).27524625
38. Fontes MAP Neurogenic Background for Emotional Stress-Associated Hypertension Curr. Hypertens. Rep. 2023 25 107 116 10.1007/s11906-023-01235-7 37058193
Fontes, M. A. P. et al. Neurogenic Background for Emotional Stress-Associated Hypertension. Curr. Hypertens. Rep. 25, 107–116. 10.1007/s11906-023-01235-7 (2023).37058193
39. Troubat R Neuroinflammation and depression: A review Eur. J. Neurosci. 2021 53 151 171 10.1111/ejn.14720 32150310
Troubat, R. et al. Neuroinflammation and depression: A review. Eur. J. Neurosci. 53, 151–171. 10.1111/ejn.14720 (2021).32150310
40. Han T Xu Y Sun L Hashimoto M Wei J Microglial response to aging and neuroinflammation in the development of neurodegenerative diseases Neural Regen. Res. 2024 19 1241 1248 10.4103/1673-5374.385845 37905870
Han, T., Xu, Y., Sun, L., Hashimoto, M. & Wei, J. Microglial response to aging and neuroinflammation in the development of neurodegenerative diseases. Neural Regen. Res. 19, 1241–1248. 10.4103/1673-5374.385845 (2024).37905870
41. Asamu MO Oladipo OO Abayomi OA Adebayo AA Alzheimer’s disease: The role of T lymphocytes in neuroinflammation and neurodegeneration Brain Res. 2023 1821 148589 10.1016/j.brainres.2023.148589 37734576
Asamu, M. O., Oladipo, O. O., Abayomi, O. A. & Adebayo, A. A. Alzheimer’s disease: The role of T lymphocytes in neuroinflammation and neurodegeneration. Brain Res. 1821, 148589. 10.1016/j.brainres.2023.148589 (2023).37734576
42. Won E Na K-S Kim Y-K Associations between Melatonin, Neuroinflammation, and Brain Alterations in Depression Int. J. Mol. Sci. 2021 23 305 10.3390/ijms23010305 35008730
Won, E., Na, K.-S. & Kim, Y.-K. Associations between Melatonin, Neuroinflammation, and Brain Alterations in Depression. Int. J. Mol. Sci. 23, 305. 10.3390/ijms23010305 (2021).35008730
43. Smith CM Relaxin-3/RXFP3 networks: An emerging target for the treatment of depression and other neuropsychiatric diseases? Front. Pharmacol. 2014 10.3389/fphar.2014.00046 24711793
Smith, C. M. et al. Relaxin-3/RXFP3 networks: An emerging target for the treatment of depression and other neuropsychiatric diseases?. Front. Pharmacol.[SPACE]10.3389/fphar.2014.00046 (2014).24711793
44. Teo S Salinas PC Wnt-Frizzled Signaling Regulates Activity-Mediated Synapse Formation Front. Mol. Neurosci. 2021 14 683035 10.3389/fnmol.2021.683035 34194299
Teo, S. & Salinas, P. C. Wnt-Frizzled Signaling Regulates Activity-Mediated Synapse Formation. Front. Mol. Neurosci. 14, 683035. 10.3389/fnmol.2021.683035 (2021).34194299
45. Tang S-J Synaptic Activity-Regulated Wnt Signaling in Synaptic Plasticity, Glial Function and Chronic Pain CNS & Neurol. Disorders - Drug Targets 2014 13 737 744 10.2174/1871527312666131223114457
Tang, S.-J. Synaptic Activity-Regulated Wnt Signaling in Synaptic Plasticity, Glial Function and Chronic Pain. CNS & Neurol. Disorders - Drug Targets 13, 737–744. 10.2174/1871527312666131223114457 (2014).
46. Bem J Wnt/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B$$\end{document}-catenin signaling in brain development and mental disorders: Keeping TCF7L2 in mind FEBS Lett. 2019 593 1654 1674 10.1002/1873-3468.13502 31218672
Bem, J. et al. Wnt/-catenin signaling in brain development and mental disorders: Keeping TCF7L2 in mind. FEBS Lett. 593, 1654–1674. 10.1002/1873-3468.13502 (2019).31218672
47. Balle, F. & Kaiserslautern, T. U. (eds) Tagungsband / Young Researcher Symposium (YRS) 2013 (Fraunhofer Verlag, Stuttgart, 2013).
48. Shishkina G Kalinina T Dygalo N Up-regulation of tryptophan hydroxylase-2 mRNA in the rat brain by chronic fluoxetine treatment correlates with its antidepressant effect Neuroscience 2007 150 404 412 10.1016/j.neuroscience.2007.09.017 17950541
Shishkina, G., Kalinina, T. & Dygalo, N. Up-regulation of tryptophan hydroxylase-2 mRNA in the rat brain by chronic fluoxetine treatment correlates with its antidepressant effect. Neuroscience 150, 404–412. 10.1016/j.neuroscience.2007.09.017 (2007).17950541
49. Lanshakov DA Sukhareva EV Kalinina TS Dygalo NN Dexamethasone-induced acute excitotoxic cell death in the developing brain Neurobiol. Dis. 2016 91 1 9 10.1016/j.nbd.2016.02.009 26873551
Lanshakov, D. A., Sukhareva, E. V., Kalinina, T. S. & Dygalo, N. N. Dexamethasone-induced acute excitotoxic cell death in the developing brain. Neurobiol. Dis. 91, 1–9. 10.1016/j.nbd.2016.02.009 (2016).26873551
50. Bolger AM Lohse M Usadel B Trimmomatic: A flexible trimmer for Illumina sequence data Bioinformatics 2014 30 2114 2120 10.1093/bioinformatics/btu170 24695404
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. 10.1093/bioinformatics/btu170 (2014).24695404
51. Kim D Langmead B Salzberg SL HISAT: A fast spliced aligner with low memory requirements Nat. Methods 2015 12 357 360 10.1038/nmeth.3317 25751142
Kim, D., Langmead, B. & Salzberg, S. L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. 10.1038/nmeth.3317 (2015).25751142
52. Liao Y Smyth GK Shi W The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads Nucleic Acids Res. 2019 47 e47 e47 10.1093/nar/gkz114 30783653
Liao, Y., Smyth, G. K. & Shi, W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 47, e47–e47. 10.1093/nar/gkz114 (2019).30783653
53. McCarthy DJ Chen Y Smyth GK Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation Nucleic Acids Res. 2012 40 4288 4297 10.1093/nar/gks042 22287627
McCarthy, D. J., Chen, Y. & Smyth, G. K. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40, 4288–4297. 10.1093/nar/gks042 (2012).22287627
54. Wu T clusterProfiler 4.0: A universal enrichment tool for interpreting omics data Innov. 2021 2 100141 10.1016/j.xinn.2021.100141
Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov. 2, 100141. 10.1016/j.xinn.2021.100141 (2021).
55. Luo W Brouwer C Pathview: An R/Bioconductor package for pathway-based data integration and visualization Bioinformatics 2013 29 1830 1831 10.1093/bioinformatics/btt285 23740750
Luo, W. & Brouwer, C. Pathview: An R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics 29, 1830–1831. 10.1093/bioinformatics/btt285 (2013).23740750
56. Shaburova EV Lanshakov DA Effective Transduction of Brain Neurons with Lentiviral Vectors Purified via Ion-Exchange Chromatography Appl. Biochem. Microbiol. 2021 57 890 898 10.1134/S0003683821080044
Shaburova, E. V. & Lanshakov, D. A. Effective Transduction of Brain Neurons with Lentiviral Vectors Purified via Ion-Exchange Chromatography. Appl. Biochem. Microbiol. 57, 890–898. 10.1134/S0003683821080044 (2021).
