Molecular Regulation of the Melatonin Biosynthesis Pathway in Unipolar and Bipolar Depression.

Melatonin is a neurohormone that maintains the circadian rhythms of the body. By regulating the secretion of other hormones and neurotransmitters, it acts as a pleiotropic modulator that affects, for example, reproductive, immune, cardiovascular, sleep, and wake systems and mood. Thus, synthetic melatonin has become an essential component in the treatment of depressive disorders. Although we know the pathway of melatonin action in the brain, we lack comprehensive cross-sectional studies on the periphery of depressed patients. This study aimed to comprehensively analyze the differences between healthy control subjects (n = 84) and unipolar and bipolar depression patients (n = 94), including an analysis of the melatonin pathway at the level of the genes and serum biomarkers. An innovative approach is a pilot study based on gene expression profiling carried out on clinical and cell culture models using agomelatine and melatonin. We confirmed the melatonin biosynthesis pathway's molecular regulation dysfunctions, with a specific pattern for unipolar and bipolar depression, at the AANAT gene, its polymorphisms (rs8150 and rs3760138), and examined the serum biomarkers (serotonin, AANAT, ASMT, and melatonin). The biological pathway analysis uncovered pathways and genes that were uniquely altered after agomelatine treatment in a clinical model and melatonin treatment in a cell culture model. In both models, we confirmed the immunomodulatory effect of melatonin agents in depression.

Introduction

Depression in its various forms, especially unipolar (UD) and bipolar (BD) depression episodes, is a common disorder in modern societies, one which severely impairs the quality of life and social functioning. The disease can gradually take a chronic form with a high relapse rate, leading to disability (World Health Organization, 2017). Depression can also contribute to suicide, the second most common cause of death among adolescents and young adults (Wasserman, 2016).

Based on the clinical picture alone, it is difficult to differentiate between depression episodes in terms of UD and BD. Therefore, there are no separate diagnostic criteria in the Diagnostic and Statistical Manual of Mental Disorders (DSM). For years, the search for laboratory markers that would facilitate the differentiation of types of depression has been carried out because incorrect diagnosis results in many years of inadequate, ineffective treatment associated with severe impairment of the patient's functioning and health. For example, the use of antidepressants in bipolar depression requires special care because of the risk of inducing mood changes or inducing mania (Baldessarini et al., 2019). Observations and clinical studies confirm the existence of differences between unipolar and bipolar depression on several levels: age of onset, relapse, family burden (Forty et al., 2008), or temperamental traits (Iasevoli et al., 2013; Fornaro et al., 2018). MRI neuroimaging also confirms that there are differences between the types of depression in terms of activation patterns in neural networks, including the amygdala, anterior cingulate gyrus (ACC), prefrontal cortex (PFC), and striatum during emotional, reward, or cognitive functions (Han et al., 2019). However, the cost and availability of the test limit its use in routine differential diagnoses.

Circadian rhythms form an integral part of the depressive symptom complex (Germain and Kupfer, 2008). Beyond disturbed sleep architecture (poor sleep quality, insomnia, or hypersomnia), a lower amplitude in the daily rhythms of locomotor activity, body temperature, and neurohormone secretion, such as norepinephrine, thyroid-stimulating hormone, and melatonin, is often observed in depressive patients (Mendoza, 2019). Based on these foundations, the chronobiological hypothesis of mood disorders was formulated (Malhi and Kuiper, 2013; Zaki et al., 2018), which became the basis for studying the use of exogenous melatonin in the treatment of depression (Boyce and Hopwood, 2013; Alston et al., 2019).

Melatonin (MEL), also known as the dark hormone, regulates the secretion of other hormones and maintains the body's circadian rhythm. Its pleiotropic effects include sleep modulation and antidepressant, anxiolytic, neuroprotective, anti-inflammatory, reproductive, antihypertensive, and antinociceptive traits, and disorders in each area may contribute to the severity of depression (Mahmood, 2019). In humans, the synthesis of MEL occurs primarily in the pineal gland, a small endocrine gland of specialized cells called pinealocytes. Small amounts of this hormone are also produced by the retina and lens of the eye, the lining of the gastrointestinal tract, and blood cells. The precursor of MEL synthesis is tryptophan (an exogenous amino acid), which is transformed into melatonin through an intermediate product, that is, serotonin. Two main rate-limiting enzymes, serotonin-N-acetyltransferase (AANAT) and acetyloserotonin-O-methyltransferase (ASMT), are involved in the process. The synthesis and release of melatonin are subject to fluctuations in the circadian cycle. The release of the hormone is regulated by postsynaptic receptors located in the hypothalamic suprachiasmatic nuclei (SCN), which receive stimuli from the retina and are considered the anatomical center of the biological clock associated with the diurnal light cycle. Melatonin achieves multiple biological effects by activating two membrane G-protein–coupled receptors, MT1 and MT2, and the retinoid orphan nuclear hormone receptor family (RZR/ROR). The primary mechanism of MEL synthesis and secretion in the brain is well understood and has been described in many previous studies (Singh and Jadhav, 2014; Dmitrzak-Weglarz and Reszka, 2017; Tordjman et al., 2017).

In both UD and BD depressive disorders, a decrease in nocturnal melatonin secretion and delayed peak secretion (phase shift) were observed (Khaleghipour et al., 2012; Melo et al., 2017). Treatment with exogenous MEL did not produce the expected antidepressant effects in humans (Carman et al., 1976). On the other hand, elevated levels of MEL in the blood plasma were noted after antidepressant treatment (Palazidou et al., 1992). In turn, a significant improvement in the state of depression and anxiety was observed after the administration of agomelatine (Salva and Hartley, 2012; Konstantakopoulos et al., 2020).

Agomelatine is a melatonergic agonist (MT1 and MT2 receptors) and an antagonist of 5-HT2C and 5-HT2B receptors. Agomelatine restores the circadian rhythm via the MT1 and MT2 receptors present in the suprachiasmatic nucleus. In turn, the antidepressant effect of agomelatine is related to its simultaneous influence on melatoninergic receptors and the 5-HT2C receptor. This action leads to an increase in dopaminergic and noradrenergic transmission in the prefrontal cortex, an increase in the 5-HT1A serotonergic receptor's activation without changes in its sensitivity, and normalization of disturbed circadian rhythms and sleep disorders. The antidepressant effect of agomelatine is also attributed to its neuroprotective properties and the reduction of glutamatergic activity resulting from its 5-HT2C receptor antagonism. It has been suggested that by acting on melatoninergic and serotonergic 5-HT2C receptors, agomelatine may also have an antianxiety effect and positively affect sexual behavior. The detailed mechanism of action, pharmacokinetics, and clinical characteristics of agomelatine have been characterized elsewhere (De Berardis et al., 2013; Fornaro et al., 2014; De Berardis et al., 2015).

Many melatonin derivatives are currently being used in the treatment of sleep disorders and depressive symptoms, many of which are available over the counter. However, we still do not fully understand the MEL pathway operation, although the available data have been updated in recent years (Hickie and Rogers, 2011).

Depressive disorders are strongly associated with genetic susceptibility (Serretti, 2017). Most of the associated studies have focused on the CLOCK genes that regulate the feedback loop of the master biological clock (Dmitrzak-Weglarz et al., 2015; Dmitrzak-Weglarz et al., 2016). Studies of single nucleotide polymorphisms (SNPs) in key genes of the melatonin biosynthesis and metabolism pathway may help clarify the pathological mechanisms of melatonin production associated with depressive disorders. However, the associated studies to date have focused on only two genes that encode the enzymes involved in the conversion of serotonin to melatonin: serotonin-N-acetyltransferase, AANAT (rs3760138, rs4238989, and rs8150) and acetyloserotonin-O-methyltransferase, ASMT (rs4446909) (Kripke et al., 2011). The associated studies on melatonin-related genes are insufficient and often not replicable in mental disorders; in depression, in some cases, there is a complete lack of them.

In recent years, expression microarrays have been used intensively as tools to search for biomarkers that change in disease or under the influence of applied treatment. Nonetheless, isolated reports refer to melatonin pathway profiling (Ha et al., 2006). Current modern microarray and bioinformatic methods may allow us to discover entirely new genes and pathways involved in MEL biosynthesis.

Although MEL was discovered more than 60 years ago, its relationship to depression is still not fully understood. Despite the well-established metabolism of MEL in the brain through animal models and postmortem studies, there is no comprehensive study of this pathway and its putative alterations in the peripheral blood in humans suffering depression.

Given the above, the purpose of this study is to investigate whether

- SNPs of selected candidate genes related to MEL biosynthesis pathway biomarkers (MBPBs) are associated with the risk of UD and BD depression

- different patterns of blood serum levels of key MBPBs feature in an episode of depression, depending on the type of depression, and after the treatment

- there are relationships between the studied polymorphisms and the levels of MBPBs in the blood serum

- screening with the use of expression microarrays allows determination of differentially expressed genes (DEGs) between agomelatine-treated vs. nontreated depressed patients (clinical model) and between MEL-treated vs. nontreated hippocampal neurons (cell culture model). We used primary cultures to analyze the specific cellular pathways after MEL treatment in a controlled environment. In this case, hippocampal neurons are used to study specific neuronal pathways that would otherwise prove difficult, if not impossible, to analyze in the intact brain (Seibenhener and Wooten, 2012).

This study is a pioneering study in the comprehensive analysis of peripheral blood MEL metabolism. In this project, we propose a comprehensive study of the MEL pathway at the levels of the genes, transcripts, and proteins through all MEL synthesis and metabolism stages.

Material and Methods

The study consisted of the following parts: genetic association analysis from peripheral blood, measurement of MEL biosynthesis pathway biomarkers (MBPBs) in serum, expression profiling analysis in PBMC after agomelatine treatment taken from depressed patients (clinical model), and after MEL treatment of mice primary hippocampal neuron culture (cell culture model).

Participants

All the participants were of Caucasian origin, and they were native Polish population. The study followed the rules of the Declaration of Helsinki and complied with Good Clinical Practice guidelines. The Bioethics Committee at the Poznan University of Medical Sciences (PUMS) (resolution no. 758/17) approved the study protocol. All participants gave written informed consent to participate in the study.

Recruitment was conducted in 2017–2019 at the Department of Adult Psychiatry, PUMS. Participants were recruited from among patients hospitalized for an episode of depression in the course of F.31 (bipolar disorder) or F.33 (recurrent depressive disorders) according to ICD-10 classification. A consensus lifetime diagnosis was made by two psychiatrists according to the ICD-10 and DSM-IV criteria, using SCID (Structured Clinical Interview for DSM Disorders) (First et al., 1996) and OPCRIT (the Operational Criteria Diagnostic Checklist) (McGuffin et al., 1991). The Hamilton Depression Rating Scale (HDRS17) was used to assess the severity of depression symptoms. A score >20 was required for inclusion in the study, indicating at least moderate severity of depression (pretreatment, Pre-T). A scoring <8 was required for achieving clinical remission or at least 50% reduction at HDRS, defined as treatment response (posttreatment, Post-T) (Hamilton, 1960). All patients after inpatient admission were treated with antidepressants and/or mood stabilizers according to doctors’ choice (Rybakowski, 2018). Patients enrolled in the study also had to meet the following inclusion/exclusion criteria (World Health Organization, 2017): age 18–65 years (Wasserman, 2016), at least the second episode in lifetime history (Baldessarini et al., 2019), and no chronic or acute somatic or neurological diseases. We offered participation in the study to a total of 110 patients. The diagnosis was not confirmed in 4 subjects, 3 subjects met the exclusion criteria, 2 subjects withdrew their consent to participate in the study, 4 subjects did not achieve minimal clinical improvement, and 3 subjects failed to obtain biological material for laboratory analyses.

The age-matched healthy control group was recruited from healthy volunteers recruited at the Nofer Institute of Occupational Medicine in Lodz. The control group was psychiatrically screened using a short diagnostic structured interview—the Polish version of the Mini International Neuropsychiatric Interview scale—to exclude those with any serious mental health problems (Lecrubier et al., 1997). Controls enrolled in the study also had to meet the following inclusion/exclusion criteria (World Health Organization, 2017): age 18–65 years (Wasserman, 2016), no family history of psychiatric disorders in first-degree relatives (Baldessarini et al., 2019), and no chronic or acute somatic or neurological diseases.

Eventually, the study was performed on 94 patients (UD, n = 41 and BD, n = 53) and 84 healthy controls (HC, n = 84) of both sexes.

Biological Material Collection

Blood samples from all subjects were restricted to those collected between 7.00 and 9.00 a.m. (after overnight fasting) to carry out all molecular analysis using the same procedures and under the same conditions to avoid daytime fluctuation and any pre-laboratory errors that could influence the final results.

Genotyping

The DNA was extracted from 5 ml of EDTA anticoagulated whole blood using the salting-out method (Miller et al., 1988). The SNP selection was based on the criteria described previously (Dmitrzak-Weglarz and Reszka, 2017). The 15 polymorphisms in the AANAT, ASMT, MTNR1B, MTNR1A, CYP1A2, and CYP1B1 genes have been identified with the use of KASP™ genotyping assays provided by LGC Genomic Laboratory (Hoddesdon, United Kingdom). In summary, KASP genotyping assays are based on competitive allele-specific PCR. Allelic discrimination is based on the competitive binding of two allele-specific primers labeled with FAM™ or HEX™ dye. The detailed described protocols can be found in the Biosearchtech Producers User Guide. , , Information about the analyzed SNPs and primers is shown in Table 1.

TABLE 1

Information about the studied polymorphisms.

Gene name ID gene	rs#	Alternative	Ancestral/reference allele	MAF b	Primer_AlleleX	Primer_AlleleY
AANAT–(serotonin-N-acetyltransferase) ID:15	rs8150	G > C, A a	G	C = 0.3201	AAT​ATT​ATT​TAA​TCC​ATG​AGA​CAT​CGT​CCA	ATT​TAA​TCC​ATG​AGA​CAT​CGT​CCG
	rs3760138	G > T, A a	G	G = 0.4781	GGG​GTG​CAG​GAT​GGG​GTG​TA	GGG​TGC​AGG​ATG​GGG​TGT​G
	rs4238989	C > G	C	C = 0.4423	CCC​TCG​GAG​TCT​TGG​ACA​AG	CCC​TCG​GAG​TCT​TGG​ACA​AC
ASMT–acetyloserotonin-O-methyltransferase ID:438	rs4446909	G > A	G	A = 0.286	GAA​CTG​ATC​CCA​AAC​CAA​TAT​GCT​C	AGA​ACT​GAT​CCC​AAA​CCA​ATA​TGC​TT
	rs5989681	G > C, A a	G	C = 0.302	AGG​GGA​TCA​CAA​AGC​GTG​TGG​T	GGG​ATC​ACA​AAG​CGT​GTG​GC
MTNR1B–melatonin receptor 1B ID:4544	rs4753426	C > T, A a	T	T = 0.4980	AAT​ACA​ACA​TAT​TTG​TGA​TTA​ATC​CAT​GCC	CAA​TAC​AAC​ATA​TTT​GTG​ATT​AAT​CCA​TGC​T
	rs10830963	C > G,T a	C	G = 0.2883	GGC​AGT​TAC​TGG​TTC​TGG​ATA​GG	GGC​AGT​TAC​TGG​TTC​TGG​ATA​GC
	rs156244 c	C > A,G a ,T a	C	C = 0.0010	ATT​AAC​GTC​ATC​ATT​AAA​ATA​ATA​AAA​CCC​AAA	AAC​GTC​ATC​ATT​AAA​ATA​ATA​AAA​CCC​AAC
	rs6483221	C > T,A a ,G a	C	C = 0.1968	AAA​AGA​TCA​GAT​ATG​TTA​ACA​TTC​TTG​ATA​AG	CAA​AAG​ATC​AGA​TAT​GTT​AAC​ATT​CTT​GAT​AAA
MTNR1A- melatonin receptor 1A ID:4543	rs2119882	C > T		T = 0.3837	CCC​CAA​TCC​CAT​TTC​GCA​TTT​GG	CCC​CAA​TCC​CAT​TTC​GCA​TTT​GA
	rs12506228	C > A	C	A = 0.2913	GAA​ACT​AAG​ATG​GGT​AGA​AAG​TCA​GAT	AAA​CTA​AGA​TGG​GTA​GAA​AGT​CAG​AG
CYP1A2–cytochrome P450 1A2 ID:1544	rs2470890	T > C	T	C = 0.35424	CAG​AAT​GGT​GGT​GTC​TTC​TTC​AA	CAG​AAT​GGT​GGT​GTC​TTC​TTC​AG
	rs762551	C > A, G a	C	C = 0.35424	CCA​TCT​ACC​ATG​CGT​CCT​GT	CCA​TCT​ACC​ATG​CGT​CCT​GG
	rs2472304	A > G	G	G = 0.4006	GTA​GAC​TGA​ACA​AAC​AAC​CTG​GGT​T	AGA​CTG​AAC​AAA​CAA​CCT​GGG​TC
CYP1B1–cytochrome P450 1B1 ID:1545	rs1056836	G > C, (exon) Leu432Val	G	C = 0.4434	GTG​GTC​TGT​GAA​TCA​TGA​CCC​AC	GTG​GTC​TGT​GAA​TCA​TGA​CCC​AG

Allele not identified in studied population.

1000GENOMES Europe population based on GRCh38. p12 version.

Variant not polymorphic in studied population.

Serum Samples

5 ml of venous blood was withdrawn into anticoagulant-free tubes. After 30 min of incubation, serum was separated by centrifugation, aliquoted, and after two-stage freezing, at −20 and then −80°C, stored until further analyses. For measurement of the MBPB enzyme-linked immunosorbent assays were performed using the commercially available tests. The detailed description is presented in Supplementary Excel File_1–ELISA Tests.

All ELISA tests were performed according to the manufacturer's instructions, without any modifications. The status of all used tests was for research use only (RUO). All samples and standards were run in duplicates, and the mean value of the two assays was used for statistical evaluation. Optical density was read with a spectrophotometric plate reader (Asys UVM 340 Microplate Reader from Biochrom Ltd., Cambridge, United Kingdom) for a wavelength of 450 nm ± 10 nm. A four-parameter logistic ELISA curve fitting was used to assay concentration in the tested samples.

RNA Sample Processing

Clinical Sample

To perform pilot profile gene expression on microarrays between agomelatine-treated and non–agomelatine-treated patients, we selected only women who did not smoke, did not take hormonal drugs, and did not work in shifts. The detailed description is presented in Supplementary Excel File_1–Microarray Clinical Samples.

4 ml of venous blood was withdrawn into an EDTA anticoagulated tube and immediately drawn through a LeukoLOCK™ filter, washed with 1xPBS (Life Technologies, United States), and saturated with RNAlater® (Life Technologies, United States). Then, the inlet ports of the filter were sealed, and the filter was stored at −80°C until processing (Gautam et al., 2019).

RNA isolation and purification were performed according to alternative protocol for extraction of RNA from PBMC captured on LeukoLOCK™ filters.

Cell Culture

Primary cultures of dissociated hippocampal neurons (HNs) were prepared from FVB mice at postnatal day 0 (P0), according to the standard protocol (Seibenhener and Wooten, 2012). The neurons were plated at a density of 1.6 × 105 cells per well in 12-well plates. Cell cultures were allowed to settle at 37°C in a 5% CO2 atmosphere. Cultured neurons from two independent cultures were stimulated for 24h with 10 µM MEL dissolved in ethanol at DIV19. The nontoxic concentration of MEL was defined based on previously published data (Park et al., 2007). Total RNA was extracted from cell culture samples using the RNeasy® Mini Kit (QIAGEN). Cell line testing is not an experimental procedure. The animals used to obtain the primary lines were treated under Directive 2010/63/EU (Appendix no. VI) by qualified staff.

Microarray-Based Gene Expression Profiling

RNA concentration and purity were evaluated using a NanoDrop 2000c UV-Vis spectrophotometer (Thermo Scientific, United States). The RNA integrity number (RIN) was measured for each sample using Agilent 2100 Expert software. Samples with a RIN value of 6 and higher were further processed (Schroeder et al., 2006). The SurePrint G3 Human Gene Expression v3 8 × 60K Microarray (Agilent, United States) was utilized in this study. The arrays were processed according to the One-Color Microarray Based Gene Expression Analysis protocol v. 6.9.1. The slides were scanned using a SureScan Dx Microarray Scanner (Agilent, United States). The images obtained after scanning were analyzed using Agilent Feature Extraction software v. 12.0.3.1. The analysis included quality control metrics (filtering of outlier spots, background subtraction from features, and dye normalization) and report. An image and a detailed description of the quality control metrics are in Supplementary Excel File_1-QC Metrics.

The principal component analysis (PCA) was conducted on microarray data (submitted to GEO NCBI resources) using filtered flags [detected, not detected] (Figure 1).

FIGURE 1

Quality control for analyzed A-cell culture and B-clinical model (normalized intensity value; PCA filtered on flags [detected, not detected]; hierarchical clustering).

Statistical and Bioinformatics Analysis

Results were expressed as numbers and percentages, mean ± standard deviation, and median, as appropriate. Statistical analyses were conducted in Statistica v13.3 (StatSoft, Poland), VassarStats (http://vassarstats.net/), and G*Power 3.1 (University Dusseldorf) software. Normality of distribution was assessed with the Shapiro–Wilk test, and equality of variances was checked using Levene's test. Due to the lack of normality for most of the variables, more restrictive nonparametric tests were used. Comparison of two unpaired groups was performed using the U Mann–Whitney test, while comparison of paired groups was performed using the Wilcoxon pair test. Multiple linear regression analysis was used to assess the influence of age and sex on other analyzed clinical and molecular variables. Spearman's rank correlation was used to detect the relationship between variables. Genotype and allele frequencies were compared using the χ2-test and Fisher exact test, respectively. All analyzed polymorphisms were at Hardy–Weinberg equilibrium. All tests were two-tailed, and α was set at 0.05 (as we compared two groups simultaneously). To determine a priori minimal sample size, we used G*Power 3 software (Faul et al., 2007; Faul et al., 2009). Given α = 0.05, power (1−β) = 0.8, a medium effect size, and the total desired sample size was as follows:

− n = 108 for genetic analyses and the χ2-test (w = 0.3),

− n = 134 for the Mann–Whitney test (d = 0.5),

− n = 35 for the Wilcoxon pair test (dz = 0.5),

− n = 68 for multiple linear regression analysis (f 2 = 0.15).

Therefore, the study included adequate sample size, and the statistical power was appropriate to detect significant differences in the studied groups.

Differentially Expressed Genes

The data obtained after extraction were further analyzed using GeneSpring 14.9.1 (Agilent, United States). The analysis aimed to determine differentially expressed genes (DEGs) in the PBMC of the patients treated with agomelatine and primary HN culture treated with MEL.

To identify up- and downregulated gene lists, we used restricted cutoff with fold change >2 and a significant p-value of <0.05 according to previous indications (Dalman et al., 2012). Significant differences in gene expression were determined via moderated Student's t-test. The differences were considered statistically significant at p < 0.05. According to the Benjamini–Hochberg procedure, the obtained p-values were corrected using the false discovery rate (FDR) as multiple testing correction methods (Cheng and Pounds, 2007).

Pathway Analysis

Lists of differentially expressed genes, with the restricted cutoff described above, were uploaded to the DAVID 6.8—Available Online Database for Annotation, Visualization, and Integrated Discovery Classification System (DAVID; https://david.ncifcrf.gov/home.jsp) to analyze their functional affiliations and participation in curated pathways (Huang et al., 2009). We used the Expression Analysis Systematic Explorer (EASE) tool to automate biological theme determination for a downloaded list of genes. We used the analysis inclusion threshold as the EASE Score—a modified Fisher exact p < 0.05 (Kleber et al., 2015). Ontologies and pathways were assigned independently to upregulated and downregulated gene lists. Annotations were limited to HGNC gene symbols and the Homo sapiens genome as a background (Wain et al., 2002). The Kyoto Encyclopedia of Genes and Genomes (KEGG) was used as a reference database (Kanehisa et al., 2019). To present significant protein–protein association networks, we used STRING v11 software (Szklarczyk et al., 2019).

Results

Sample Description

The clinical characteristics of patients and control subjects are presented in Table 2.

TABLE 2

Clinical description of the study population.

	BD (n = 53)			UD (n = 41)			HC (n = 84)			M–W test
	Mean	SD	Median	Mean	SD	Median	Mean	SD	Median	BD vs. UD	BD vs. HC	UD vs. HC
Age	41.91	13.87	43.00	43.49	15.99	47.00	43.11	10.91	40.00	0.616	0.845	0.716
Years of education	13.22	2.87	12.00	13.72	2.83	12.00	15.86	2.06	17.00	0.488	1.047E−06	1.712E−04
BDI	34.47	9.44	35.50	30.23	10.83	29.00	4.04	2.95	4.00	0.037	1.936E−29	6.280E−27
BMI	24.46	4.88	23.10	24.89	4.96	24.16	23.80	3.20	23.62	0.679	0.764	0.428
Onset of illness	27.57	9.80	26.00	35.16	14.49	33.00				0.013
Number of hospitalizations	7.35	15.76	3.00	1.89	2.16	1.00				2.470E−04
Duration of current hospitalization (days)	50.64	24.92	50.00	43.32	16.02	43.00				0.283
HDRS Pre-T	25.93	6.60	24.00	25.50	4.76	25.00				0.783
HDRS Post-T	3.17	2.50	3.00	4.41	4.19	3.00				0.316
W-P test (HDRS Pre-T vs. Post-T)			7.623E−09			5.392E−07

BD, bipolar disorder; UD, unipolar disorder; HC, healthy control; n, number; SD, standard deviation; M–W test, Mann–Whitney U test; W-P test, Wilcoxon pair test.

Significant p-values are indicated in bold (significance considered p < 0.05).

In our group, both the BD and UD patients had significantly shorter education (p < 0.001) and higher BDI scores (p < 0.001) than the HC subjects. Bipolar patients demonstrated a higher score of BDI (p = 0.037), an earlier age of onset (p = 0.013), and more hospitalizations (p < 0.001) than UD patients. All patients achieved a significant and required posttreatment improvement in depression symptoms on the HDRS scale (p < 0.001 for BD and UD, respectively). The duration of current hospitalizations did not differ significantly between bipolar and unipolar depression (p = 0.283).

During recruitment, male patients more often refused to take part in the study, resulting in the sex ratio disproportion: BD: n = 7 (13.2%), UD: n = 4 (9.8%), HC: n = 13 (15.5%). To control sex and age as confounding factors, we performed multiple linear regression analysis. This analysis showed no significant effects of age (p = 0.821) or sex (p = 0.472) on clinical and molecular variables (F (2,178) = 0.293; p = 0.746; R 2 = 0.003). Thus, in further analyses, groups were not subdivided according to age and sex of participants.

Out of 15 analyzed polymorphisms related to the MEL pathway, only rs156244 in the MTNR1B gene turned out to be non-polymorphic in the analyzed population. No deviation from the Hardy–Weinberg equilibrium law was observed for the rest of the analyzed polymorphisms. Case–control comparisons between patients with bipolar depression and control subjects showed an association of two AANAT polymorphisms (rs8150; p = 0.028 and rs3760138; p = 0.043). For the other analyzed polymorphisms, we did not find any significant differences in genotype and allele frequencies between the compared groups of patients and controls (Supplementary Excel File_1-Genotypes & Alleles).

Serum Expression of MEL Biosynthesis Pathway Biomarkers (MBPBs)

The serum concentrations of MBPBs are presented in Table 3 and Figure 2.

TABLE 3

Comparisons of melatonin biosynthesis pathway biomarkers in the serum level.

		BD				UD				HC				Comparisons
		SERT ng/ml	AANAT ng/ml	ASMT ng/ml	MEL pg/ml	SERT ng/ml	AANAT ng/ml	ASMT ng/ml	MEL pg/ml	SERT ng/ml	AANAT ng/ml	ASMT ng/ml	MEL pg/ml	M–W test p-value	SERT ng/ml	AANAT ng/ml	ASMT ng/ml	MEL pg/ml
Pre-T	Mean	26.03	2.07	5.97	1267.85	31.28	2.29	5.19	1508.37	36.68	1.60	4.01	2329.27	BD vs. HC	0.0101 (↓)	0.0334 (↑)	0.0142 (↑)	<0.0001 (↓)
	SD	9.74	2.07	4.11	1018.38	16.30	2.45	6.31	1266.95	17.22	2.56	2.44	1200.01	UP vs HC	0.4275	0.1479	0.9265	<0.0001 (↓)
	Median	25.02	1.50	7.09	965.70	26.81	1.40	3.69	1282.52	17.82	1.17	3.30	1938.24	BD vs. UD	0.2259	0.8534	0.1421	0.2364
Post-T	Mean	27.39	2.23	7.70	1873.95	31.01	2.70	6.82	2265.82					BD vs. HC	0.0965	0.0643	0.0955	0.0001 (↓)
	SD	8.43	1.99	10.69	1928.45	12.79	2.65	6.13	1684.21					UP vs HC	0.8961	0.0079 (↑)	0.0194 (↑)	0.0159 (↓)
	Median	25.30	1.42	4.57	1519.71	27.68	1.93	6.08	1674.58					BD vs. UD	0.1759	0.4982	0.5802	0.1160
W-P test p-value		0.0173 (↑)	0.1992	0.5313	0.0018 (↑)	0.2137	0.2845	0.9250	0.0006 (↑)

UD, unipolar disorder; BD, bipolar disorder; Mann–Whitney U test; Wilcoxon pair test; (↑)—the arrow shows the direction of changes in the concentration of the tested metabolite in relation to the reference group.

Significant p-values are indicated in bold (significance considered p < 0.05).

FIGURE 2

Melatonin biosynthesis pathway biomarkers’ serum level pre- and posttreatment (box and whisker plots showing median values, interquartile ranges (box), non-outlying values (whiskers), and outlier points (dot)).

We observed significant differences in the concentration of all the measured MBPBs between pretreatment BD and HC. In BD patients, serotonin and MEL concentration were significantly decreased (SERT p = 0.0101; MEL p < 0.0001), while both measured enzymes were increased (AANAT p = 0.0334; ASMT p = 0.0142). After treatment (comparison between pre- and posttreatment), we observed an increase in all biomarkers’ concentration, although it was significant only for serotonin (SERT p = 0.0173) and MEL (MEL p = 0.0018). A recomparison of posttreatment BD vs. HC has revealed the existence of a still significantly lower level of MEL concentration (MEL p = 0.0001).

Comparison between pretreatment UD and HC showed that similarly to BD, the concentrations of serotonin and MEL were decreased, and the enzymes increased, although the significant difference was only in the concentration of MEL (MEL p < 0.0001). After treatment (pre- vs. posttreatment), we observed a significant increase in MEL concentration (MEL p = 0.0006). Recomparison of posttreatment UD vs. HC revealed that the differences in the concentrations of both enzymes worsened in UP patients (AANAT p = 0.0079; ASMT p = 0.0194). Also, a significantly lower level of MEL concentration (MEL p = 0.0159) was still observed.

Serum Concentration of MBPBs Dependent on Genotype Carriers

In the patients, with regard to the type of depression and the disease's state, we found several significant dependencies: rs2119882 and ASMT (p = 0.0012), rs3760138 and MEL (p = 0.0330), rs762551 and AANAT (p = 0.0022), rs3760138 and AANAT (p = 0.0094), rs12506228 and AANAT (p = 0.0269), and rs1056836 and ASMT (p = 0.0267) (Supplementary Excel File_1-Spearman's Rank Correlation).

For significantly correlated pairs, we provided comparisons of serum concentration of MBPBs between genotypes. Only the rs3760138 polymorphism of the AANAT gene, whose GG genotype, is associated with the highest concentrations of MEL in the pretreatment state (p = 0.0281), while the TT genotype is associated with the highest concentrations of AANAT in the posttreatment state (p = 0.0085), turned out to differentiate the type of depression and the disease state (Figure 3).

FIGURE 3

Differences of MEL and AANAT serum level between AANAT rs3760138 genotypes.

Differentially Expressed Genes (DEGs)

Clinical Model

We found numerous differentially expressed genes (DEGs) in the clinical model. For example, in agomelatine (+), we discovered 2217 upregulated and 2431 downregulated genes with minimal FC > 2. After gene ontology term enrichment, only 542 upregulated and 733 downregulated genes were functionally known by the DAVID/KEGG database. We used the lists of the upregulated and downregulated genes in agomelatine (+) and (-) as a set of genetic markers. The comparison of these collections allowed for identifying unique and shared gene expression profiles between the compared groups of patients. For comparisons, Venn diagrams were constructed (Figure 4). The list of unique and shared genes was processed again in the DAVID 6.8 database to obtain renewed data.

FIGURE 4

Venn diagrams representing unique and shared differential gene expression profiles of agomelatine (+) & (−) patients, numbers of genes upregulated (UP) and downregulated (DOWN), blue ring—number of genes uniquely up- or downregulated in agomelatine (+), and red ring—number of genes uniquely up- or downregulated in agomelatine (−). Overlapping parts of the diagram—the number of shared genes uniquely up- or downregulated between groups of patients. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) *functionally known by the DAVID/KEGG database.

Cell Culture Model

A fold change of more than 2 was observed for 762 upregulated and 17729 downregulated genes in MEL (+) vs. MEL (-) treated cell cultures. The 762 up- and 5343 downregulated genes were functionally known by the DAVID/KEGG database.

Gene Ontology Term Enrichment

In Table 5, we present a summary of gene set enrichment and pathway analysis.

TABLE 5

Gene set enrichment.

	Up				Down
	Gene ontologies			KEGG_P	Gene ontologies			KEGG_P
	BP	CC	MF		BP	CC	MF
Agomelatine (+) vs. HC	67	10	43	6	667	65	45	37
Agomelatine (−) vs. HC	39	14	16	1	305	41	70	20
Unique for agomelatine	76	6	30	3	639	52	30	45
Melatonin (+) vs. (−)	9	4	2	—	213	28	63	10

BP, biological process; CC, cellular component; MF, molecular function; KEGG_P, KEGG pathway; analysis inclusion threshold, EASE score, a modified Fisher exact p < 0.05.

Genes upregulated in agomelatine (+) significantly enriched 120 Gene Ontologies: 67 from the Biological Process (BP) category, 10 from the Cellular Component (CC) category, and 43 from the Molecular Function (MF) category. Upregulated genes significantly enriched 6 KEGG pathways.

Genes downregulated in agomelatine (+) enriched 777 Gene Ontologies: 667 from the Biological Process, 65 from the Cellular Component, and 45 from the Molecular Function categories. Downregulated genes significantly enriched 37 KEGG pathways.

The factsheet for Gene Ontology (GO) term enrichment and pathway analysis provided in all individual analyzed types of models is presented in Table 5.

The complete list of significantly altered GO terms and pathways was included in Supplementary Excel File_2 with the appropriate Excel spreadsheet.

Pathway Analysis and Gene Interaction Network

In Supplementary Excel File_1-KEGG Pathways & Genes, we present a significant pathway and genes in both analyzed clinical and cell culture models. The pathway and genes that pass FDR correction are presented in Table 6. Interaction analyses of genes creating unique pathways for agomelatine treatment showed 61 nodes, 256 edges, 8.39 average node degree, and 0.497 average local clustering coefficient (PPI enrichment p-value< 1.0e−16) (Figure 5). Approximately 60 proteins that had interacted together were found.

TABLE 6

KEGG pathways pass after FDR correction.

	KEGG pathway	Genes in pathway	p-value	FDR
	Agomelatine (+) vs. HC	Downregulated genes
	hsa04666: Fc gamma R-mediated phagocytosis	AKT2, LIMK2, DOCK2, WASF2, HCK, CRK, MAPK3, RAC2, PLPP2, GAB2, ARPC1B, LIMK2, WAS, RAF1, PIK3CD	0.0001	0.0161
	hsa04062: chemokine signaling pathway	AKT2, DOCK2, ADCY1, GNB2, IKBKG, SHC1, HCK, CRK, MAPK3, CCL2, PRKACA, CCR4, RAC2, CXCR1, NFKBIB, CCR1, WAS, RAF1, PIK3CD, ARRB1, STAT3	0.0002	0.0161
	hsa04380: osteoclast differentiation	AKT2, LILRB3, TYK2, IKBKG, IL1B, LILRA1, MAPK3, NFKB2, LILRA4, FOS, GAB2, LILRB4, LILRA3, PIK3CD, TGFB1, JUNB, LILRA2	0.0002	0.0161
	hsa04722: neurotrophin signaling pathway	AKT2, NTRK1, MAP3K5, SHC1, CRK, MAPK3, NTRK2, MAPKAPK2, BAX, RAPGEF1, NFKBIB, RAF1, PIK3CD, NGFR, SH2B1	0.0008	0.0478
	Agomelatine (-) vs. HC	Downregulated genes
	hsa04080: neuroactive ligand–receptor interaction	AVPR2, GABRQ, GRIN2C, GLP2R, TACR1, GPR50, MC5R, PRSS3, VIPR2, LPAR2, GRIK5, NPY5R, CHRNB4, CRHR2, RXFP2, ADORA1, SCTR, GRIA1, MC3R, RXFP3, TRHR, APLNR, PTGER1, GABRE, GABRR3, GHRHR, VIPR2, GRIN2C, NMUR2, GLRA1, GRIN3A	6.8960E−08	1.5378E−05
	hsa04972: pancreatic secretion	ADCY6, PLA2G5, CELA2A, PLA2G12B, GNAS, ATP1B2, PLA2G2A, ATP2B3, ATP2A1, SCTR, ATP2B2, PRSS3	0.0002	0.0221
	Unique pathways for agomelatine treatment	Downregulated genes
	hsa04062: chemokine signaling pathway	AKT2, DOCK2, ADCY1, GNB2, IKBKG, SHC1, HCK, CRK, MAPK3, CCL2, PRKACA, CCR4, RAC2, CXCR1, NFKBIB, CCR1, WAS, RAF1, PIK3CD, ARRB1, STAT3	3.3486E−06	0.0007
	hsa04666: Fc gamma R-mediated phagocytosis	AKT2, LIMK2, DOCK2, HCK, CRK, MAPK3, RAC2, PLPP2, GAB2, ARPC1B, LIMK2, WAS, RAF1, PIK3CD	1.9636E−05	0.0020
	hsa04380: osteoclast differentiation	AKT2, LILRB3, IKBKG, IL1B, LILRA1, MAPK3, NFKB2, LILRA4, GAB2, LILRB4, LILRA3, PIK3CD, TGFB1, JUNB	0.0004	0.0277
	hsa05220: chronic myeloid leukemia	AKT2, CDKN1A, GAB2, RAF1, PIK3CD, IKBKG, SHC1, CRK, MAPK3, TGFB1	0.0006	0.0309
	hsa05020: prion diseases	ELK1, BAX, NOTCH1, IL1B, MAPK3, STIP1, PRKACA	0.0008	0.0327
	hsa05145: toxoplasmosis	AKT2, HSPA6, BIRC8, IKBKG, IL10RB, LAMA2, MAPK3, MYD88, NFKBIB, TLR4, STAT3, TGFB1	0.0011	0.0361
	hsa05166: HTLV-I infection	AKT2, ITGB2, ELK1, ADCY1, ZFP36, IKBKG, ATF3, LTBR, NFKB2, PRKACA, CDKN1A, BAX, ANAPC2, PIK3CD, DVL3, PDGFB, TGFB1, DVL1, VAC14	0.0020	0.0488
	hsa04722: neurotrophin signaling pathway	AKT2, BAX, NFKBIB, RAPGEF1, MAP3K5, RAF1, PIK3CD, SHC1, CRK, MAPK3, SH2B1, MAPKAPK2	0.0022	0.0488
	hsa05142: Chagas disease (American trypanosomiasis)	AKT2, MYD88, ADCY1, TLR4, CFLAR, PIK3CD, IKBKG, IL1B, MAPK3, TGFB1, CCL2	0.0024	0.0488
	hsa05140: leishmaniasis	ELK1, ITGB2, PTPN6, MYD88, NFKBIB, TLR4, IL1B, MAPK3, TGFB1	0.0024	0.0488
	hsa04668: TNF signaling pathway	AKT2, CREB5, MAP3K5, CFLAR, PIK3CD, IKBKG, TNFRSF1B, IL1B, MAPK3, CCL2, JUNB	0.0029	0.0547
	hsa05161: hepatitis B	AKT2, CREB5, ELK1, IKBKG, MAPK3, CDKN1A, BAX, MYD88, TLR4, RAF1, PIK3CD, STAT3, TGFB1	0.0033	0.0553

FIGURE 5

Interaction network of unique genes for agomelatine treatment.

Unique Genes for MBPBs in Analyzed Models

We compared the lists of genes in altered KEGG pathways (Supplementary excel file_1 - Most Frequent Genes in Pathways). We found that five downregulated genes (STAT3, CDKN1A, PDGFB, SHC1, and ELK1) were most frequently repeated in the KEGG pathways unique for agomelatine treatment.

Discussion

Genes and Polymorphisms

In the present study, we analyzed genes involved in the MEL metabolism pathway and selected polymorphisms associated with symptoms of mental disorders (Dmitrzak-Weglarz and Reszka, 2017). The studies of polymorphisms in the MEL pathway genes are few, ambiguous, and outdated (Kripke et al., 2011). Only in rs4446909 ASMT is the association confirmed with UD, BD, reduced ASMT mRNA expression, and lower enzymatic activity of ASMT (Etain et al., 2012; Geoffroy et al., 2014). Our study found an association between two polymorphisms (rs8150 and rs3760138) in the AANAT gene and bipolar depression. We also observed that variants of the rs3760138 polymorphism correlate with the concentration of MEL measured in the blood serum of patients experiencing aggravation of depressive symptoms and the concentration of the enzyme itself in remission. The results of our research are consistent with those of the analysis presented by Soria et al. (2010). We discovered that the TT variant genotype is associated with both the risk of depression and a higher concentration of the encoded enzyme. Interestingly, in a case–control study of Norwegian nurses working night shifts, the GG carriers of rs3760138 in the AANAT gene showed an increased risk of breast cancer (Zienolddiny et al., 2013). The present results provide evidence that AANAT gene variants play a role in the pathophysiology of depression and support the hypothesis of the MEL signaling pathway in depression. The lack of significant correlations in the remaining polymorphisms may result from a weak association that is elusive in a limited clinical trial and does not exclude the influence of these genes in the metabolism of MEL and the depressive disorders related to it. The overall evidence from our study and other studies encourages the continuation of more detailed analyses of this pathway.

Serum Biomarkers

Measurements of the four critical biomarker concentrations in the MEL pathway in the blood serum confirmed the disturbances of this pathway in depressive disorders and their partial normalization after treatment. In both UD and BD, patients revealed different patterns of concentration of individual biomarkers before and after treatment and compared to the control group.

Our study found that the concentration of global serotonin in the blood serum was decreased nominally in UD and statistically significantly in BD. These results are consistent with the serotonin concentration measurements in the plasma of patients with depression (Saldanha et al., 2009). It is indicated that due to the rapid mitochondrial metabolism of serotonin by the MAO enzyme to 5-hydroxyindole acetic acid (5HIAA), the latter forms a better biomarker of SERT biotransformation (Jayamohananan et al., 2019). A significantly higher concentration of 5HIAA in the plasma of patients with depression suggests a rapid breakdown of serotonin in the brain (T P et al., 2008). The serotonin hypothesis of depression indicates that impaired serotonin function may cause clinical depression, but it is unnecessary or sufficient (Cowen and Browning, 2015).

We observed that the concentrations of AANAT and ASMT enzymes were nominally elevated in UD and statistically significantly increased in BD. These results identify and confirm changes in the expression level of the enzymes depending on the type of depression, although earlier in depression (Talarowska et al., 2014) and the autism spectrum, the decreased activity of these enzymes was demonstrated (Pagan et al., 2017). Of particular interest was the increase in the concentration of both the AANAT and ASMT enzymes observed in UD after treatment. Antidepressants, in particular SSRIs, inhibit serotonin reuptake by increasing its concentration in the inter-synaptic space, which leads to increased stimulation of the postsynaptic serotonin receptors.

The expression of enzymes that use serotonin to synthesize MEL can potentially be stimulated by pharmacologically disturbed serotonin levels. On the other hand, faster serotonin turnover may prevent the capturing of subtle changes in peripheral serotonin concentration (Stapelberg et al., 2018).

In our study, the MEL concentration in the blood serum was significantly lower in patients with depression regardless of the type of depression, which is consistent with previous reports (Khaleghipour et al., 2012; Oglodek et al., 2016). This observation also partially overlaps with the results of Bumb et al. (2016). Additionally, they also confirmed altered MEL levels in depressed patients. Like us, they selected a group of patients with UD and BD and also measured the MEL concentration in the cerebrospinal fluid. In the serum, they observed a significant decrease in MEL for UD but not for BD, while we observed a decrease in both. The discrepancies in the results may be due to our inclusion of a homogeneous BD group only in the depressive state, while Bumb et al. included patients in manic and depressive states. On the other hand, in the CSF, they observed a significantly reduced MEL level in BD, unlike in the serum, but not in UD. Our study, in turn, is in contradiction with the results of the Chinese population. The serum MEL levels in the first-episode and recurrent groups were significantly higher than those in the control group. The authors stated that the results may have been caused by the measurement technique used, the treatment applied, and the stratification effect (Tao et al., 2020). We have observed that MEL secretion disturbance is more severe in BD and more challenging to normalize with treatment. Persistent disturbance of MEL may be associated with the maintenance of sleep disorders (mainly reduced need for sleep), and these, in turn, contribute to more frequent relapses and greater intensity of depression symptoms, which we observed in the studied BD group. The results of our research provide molecular evidence for the validity of the methods used for the regulation and normalization of sleep and wake rhythm in BD therapy (Miklowitz et al., 2012; Alston et al., 2019).

To summarize, the results suggest the existence of diverse, specific patterns of serum concentration of MBPBs for types of depression, which do not exclude the existence of tissue-specific patterns. Longitudinal studies are needed to investigate the phase–disease relationship in MEL changes and the causes and consequences of site-specific changes.

Clinical Model

We performed a screening analysis of gene expression profiling in PBMC obtained from patients hospitalized due to a depressive episode to detect altered genes in patients treated with agomelatine. Agomelatine is a synthetic analog of MEL and acts as an agonist on the MT1 and MT2 receptors. Agomelatine also acts as an antagonist on the 5-HT2C receptors and has a low affinity for most other receptors, including adrenoceptors, dopamine, GABA, muscarinic, histamine, benzodiazepine, and sigma receptors, as well as ion channels (Norman and Olver, 2019). We selected a list of genes specific to agomelatine action that formed twelve significant downregulated pathways. Interaction analyses showed that approximately 60 proteins of these pathways had interacted together.

The common denominator of the identified pathways was participation in inflammatory processes (i.e., chemokines, neurotrophins, and cytokines). In particular, neurotrophic factors regulate the survival, development, and function of nervous tissue. Inhibition of inflammation is essential for neuroprotection and antiapoptotic activities (Skaper, 2018). We know that MEL acts as a pro-inflammatory modulator when the body is attacked by pathogens (Hardeland et al., 2015). However, when the inflammation increases to the extent that there is a risk of cell damage, it becomes an anti-inflammatory modulator (Carrillo-Vico et al., 2013). In this way, it protects the body against the development of chronic inflammatory diseases (Radogna et al., 2010). The results confirmed that agomelatine, as an analog of MEL, may act as a "buffer" of the immune system and improve depressive symptoms. Thus, attempts have been made to use exogenous MEL to treat COVID-19 caused by the SARS-CoV-2 virus, especially in the elderly with a reduced level of MEL (Kleszczynski et al., 2020; Shneider et al., 2020).

Melatonin and its ability to modulate inflammation in multiple sclerosis have been the subject of intense research in recent years. The high coexistence of depression is also observed in this disease with one common denominator—immunological changes (Feinstein et al., 2014; Solaro et al., 2018). Multiple sclerosis (MS) is an inflammatory demyelinating disease. An increased risk of this disease has been observed in residents in high-latitude countries where sunlight is limited and populations are deficient in vitamin D and high levels of melatonin (Ghareghani et al., 2018). It has been suggested that oxidative stress could be responsible for this disease. Thus, melatonin, as a natural antioxidant in the CNS, may constitute helpful therapy. Ghareghani et al. (2018), using a mouse model of MS, showed that melatonin promoted remyelination, decreased inflammation, and stimulated the activity of antioxidant enzymes (Ghareghani et al., 2017; Ghareghani et al., 2018). The research team also showed that MEL administration in the acute phase of MS might exacerbate encephalomyelitis in line with the aforementioned dualistic MEL function as a pro- or anti-inflammatory modulator depending on inflammation intensity (Ghareghani et al., 2017). In MS treatment, a significant impact is caused by time and dose of MEL and bilateral interaction between used drugs. For example, MEL improved the baclofen effect on spasticity (Ghareghani et al., 2018). In turn, corticosteroids decreased MEL levels and caused sleep disturbance (Dokoohaki et al., 2017).

The analysis of the most frequently repeated genes in the identified pathways, which changed explicitly under the influence of agomelatine, revealed the following genes:

STAT3 (signal transducer and activator of transcription 3) acts as a transcription activator in response to cytokines and growth factors

CDKN1A (cyclin-dependent kinase inhibitor 1A) functions as a regulator of cell cycle progression at G1

PDGFB (platelet-derived growth factor subunit B) activates PDGF receptor tyrosine kinases, which play a role in a wide range of developmental processes

SHC1 (SHC adaptor protein 1) is an adapter protein in the signal transduction pathways

ELK1 (ETS transcription factor ELK1) is a transcription factor that binds to purine-rich DNA sequences and induces target gene transcription upon JNK-signaling pathway stimulation.

The genes listed above belong to the transcription factor group that influences gene expression regulation and, thus, the etiology of many mental disorders or symptoms. Therefore, their further detailed study is necessary (Semenza, 2005; Fuxman Bass et al., 2015; Papavassiliou and Papavassiliou, 2016).

Cell Culture Model

To check the changes taking place in the brain under the influence of MEL, we conducted a pilot experiment using cell cultures based on mouse hippocampal neurons. We identified ten downregulated pathways (hsa03018: RNA degradation, hsa05200: pathways in cancer, hsa03460: Fanconi anemia pathway, hsa01100: metabolic pathways, hsa04713: circadian entrainment, hsa01130: biosynthesis of antibiotics, hsa05134: legionellosis, hsa04810: regulation of actin cytoskeleton, hsa04726: serotonergic synapse, and hsa04610: complement and coagulation cascades). These results demonstrate the influence of MEL in the CNS on the generally understood neuronal metabolism, including the serotoninergic system and circadian rhythms. Additionally, they indicate the participation of disturbances of the MEL level in the etiology of depression.

Wiechmann (2002), for the first time, presented the results of DEG analyses in two types of MEL-treated rat eye nerve cells using the microarray technique. Although identified, the genes with altered expression in both types of examined cells did not coincide, indicating not only the tissue but also the cellular specificity of MEL action (Wiechmann, 2002).

In turn, Ha et al. (2006), for the first time, presented the results of the microarray analysis of MEL-treated PBMC cells obtained from healthy volunteers. Forty-six upregulated and twenty-three downregulated genes were identified. Most genes formed a cellular physiological response pathway. In this pathway, the most remarkable upregulation change was shown by the CFTR gene (cystic fibrosis transmembrane conductance regulator) and the downregulation by gene UBE2M (ubiquitin-conjugating enzyme E2M) (Ha et al., 2006). In 2007, the same team performed a microarray analysis with a modified protocol (Park et al., 2007). MEL reduced the expression of genes associated with inflammation. In particular, chemokine genes were downregulated, which we also observed in patients taking agomelatine. These results confirmed that MEL exerts anti-inflammatory properties by regulating the chemokines, which act as crucial mediators in various inflammatory processes (Ban et al., 2011).

In the cell culture model under study, we observed a relatively large number of downregulated genes. The distribution of melatonin receptors in the hippocampal cells is lower than that of the cerebellum or occipital cortex (Mazzucchelli et al., 1996; Pandiperumal et al., 2008; Klosen et al., 2019). Despite the selection of the experimental dose of melatonin based on the available literature data, we cannot exclude the possibility that the unbound melatonin exerts a subtoxic effect on NH cells, expressed by an increased number of downexpressed genes.

To summarize, we obtained a joint conclusion confirming that MEL is involved in modulating the immune system. The effect on the expression of the number and type of genes depends on the model tested (human/mouse), type of cells tested (PBMC/macrophages), protocol (incubation time/concentrations of mediators used), and type of microarray (global/dedicated). Hence, it is not easy to make a detailed comparison, as each of the discussed studies showed different sets of genes either up- or down-regulated by MEL.

Limitations

Our screening microarray results include female patients only, so the study should be repeated in an independent and homogeneous group of men. The expression data obtained from PBMC reflect only about 85% of the changes in the brain, disregarding the tissue-specific ones. Measurement of the general SERT and MEL levels in serum may not be the optimal method. Although the NH cell line is currently the best research model for cell pathways in neuroscience (Seibenhener and Wooten, 2012), it is not necessarily an equally good model for observing melatonin effects. Compared to other brain structures, the hippocampus has a lower expression of melatonin receptors (Mazzucchelli et al., 1996; Klosen et al., 2019). Finally, the total sample size was small, although comparable to other available studies, and achieved 80% power.

Conclusion

The microarray analysis results should be treated as preliminary because they are not without limitations, mainly of a small test sample. However, they are a valuable complement to clinical observations at the molecular level consistent in the modulation of the immune system by melatonin. These results may help determine further research directions, such as searching for more adequate cell models for the melatonin pathway or functional research focused on transcription factors that are challenging to interpret due to multivariate determinants. It may also be essential to develop a uniform pipeline of bioinformatic analyses in the future (critical cutoff thresholds, applied corrections, and algorithms to predict gene ontology annotations), which currently may give partially different or complementary results in different centers.

In summary, we confirm the dysfunctions in the molecular regulation of the MEL biosynthesis pathway with a specific pattern for unipolar and bipolar depression at the tier of genes, their polymorphisms, and examined serum biomarkers. Using a systematic search of DEGs with restricted thresholds, we confirm that microarray profiling is a useful biological marker search tool. The biological pathway analysis uncovered pathways and genes uniquely altered after agomelatine treatment in a clinical sample and MEL treatment in cell culture. In both models, we confirmed the immunomodulatory effect of MEL agents in depression.

TABLE 4

Differentially expressed genes (DEGs) in studied models.

Model	Total DEGs, p < 0.05, FC > 2	Up	Down	DAVID/KEGG a	Up	Down
Clinical model
Agomelatine (+) vs. HC	4648	2217	2431	1275	542	733
Agomelatine (−) vs. HC	3999	567	3432	803	152	651
Common and different genes (Venn diagram)	3224	1688	1536	947	402	545
Cell culture model
Melatonin (+) vs. melatonin (−)	18,489	762	17727	13,105	762	5343

FC, fold change; PBMC, peripheral blood mononuclear cells; HN, hippocampal neurons.

Functionally known by the DAVID/KEGG database.