Introduction

The use of selective serotonin reuptake inhibitor (SSRI) antidepressants during pregnancy has greatly increased [1,2,3,4]. Every year, hundreds of thousands of babies are born that have been exposed to SSRI medication in Europe and the US alone [5,6,7,8]. SSRIs cross the placental barrier [9] and target the fetal serotonin transporter (SERT) at a time during development when serotonin serves as a neurotrophic factor. SSRIs potentially affect brain circuit formation [10], evidenced by altered white and gray matter architecture [11] and connectivity [12,13,14] in babies. Although not all studies identify long-term associations between in utero SSRI exposure and cognitive and behavioral outcomes [15,16,17], many have reported higher anxiety [18], lower motor-, social- emotional- and adaptive behavior [19] and a higher risk of developing mental and behavioral disorders [20]. Some of these associations are sex specific [21,22,23]. A limitation of this field is the confounding factor of maternal depression [20, 24, 25], which has been associated with neurodevelopmental outcomes resembling those of SSRI exposure [16, 26,27,28,29], again with sex-specific effects observed [30,31,32,33]. Since prenatal exposure to SSRIs and maternal depression both affect the developing brain [34], in particular corticolimbic structures [12, 31, 32, 35, 36], it remains a challenge to assess the relative contribution of each.

Rodent experiments offer the ability to causally investigate the separate and combined effects of SSRI administration and aspects of maternal depressed mood. Neurodevelopmental patterns, including serotonin system-specific ones [37,38,39], are remarkably conserved between rodents and humans [40]. Meta-analysis of rodent behavioral outcomes after perinatal SSRI exposure revealed reduced activity and exploration behavior, a more passive stress-coping style, and less efficient sensory processing [41]. The neurobiological correlates of these behavioral changes likely include the serotonergic system [42], with roles for the prefrontal cortex (PFC) [43, 44], limbic structures [45], and the dorsal raphe nucleus [44, 45]. Rodent models recapitulate basic neurobiological changes seen in women with peripartum depression and their children [46]. In utero exposure to maternal stress alters behavioral outcomes in rodents, such as anxiety-like behavior and cognitive performance [47].

Rodent studies also enable investigation of molecular changes during fundamental neurodevelopmental periods that may trigger these long-term effects [48,49,50,51]. Whereas prenatal SSRI exposure and/or stress showed no effects in the amygdala, hippocampus or hypothalamus [52], early postnatal SSRI exposure altered transcriptomic status of the hippocampus [53]. Developmental screening of gene expression revealed that the hippocampus showed the largest paroxetine-induced differences in the first 2 postnatal weeks, and the amygdala around PND21 [54], suggesting that the transcriptomic response to SSRIs depends on brain region. Despite evidence for sex-specific effects of neurodevelopmental outcomes after both perinatal SSRI exposure [55,56,57] and prenatal stress [55, 58], transcriptome-wide analyses were only performed in males [52,53,54], and only one study considered prenatal stress [52].

To address these gaps in the literature, we aimed to investigate transcriptomic alterations in the corticolimbic circuitry of male and female juvenile rats exposed to maternal adversity and/or perinatal SSRIs. To this end, we used a genetic rat model of maternal vulnerability (MV). MV females exposed to early life stress (sMV) show a depressive-like phenotype in adulthood compared to control (cMV) females [59]. We treated sMV and cMV dams with the SSRI fluoxetine (FLX) or vehicle (Veh) throughout pregnancy and lactation. We collected brains from juvenile offspring for RNA sequencing of micropunched tissue from the medial PFC and the basolateral amygdala (BLA). To examine a potential role for epigenetic regulation, we quantified DNA methylation levels of differentially expressed genes in the same samples.

Materials and methods

Experimental animals

Heterozygous SERT knockout (SERT+/−, Slc6a4Hubr) female Wistar rats were used as MV females [60], because of their vulnerability to the effects of early life stress [59]. Animals were supplied with ad libitum chow (RMH-B, AB Diets; Woerden, the Netherlands) and water and kept on a 12:12 h light-dark cycle (lights off at 11:00 a.m., temperature 21 ± 2 °C, humidity 50 ± 5%). Cages were enriched with wooden gnawing sticks and nesting material (Enviro-dri™, Shepherd Specialty Papers, Richland, MI, USA), and cleaned weekly. Pregnant dams were housed individually in type III Makrolon cages. Pups were weaned at PND21 and group-housed in same-sex cages (type IV Makrolon) of 3–5 animals. All experimental were conducted in accordance with the governmental guidelines for care and use of laboratory animals (Centale Commissie Dierproeven) and were approved by the Institutional Animal Care and Use Committee of the University of Groningen.

Maternal adversity and fluoxetine treatment

SERT+/− female rats (MV) were exposed to early life stress (sMV) to induce anhedonia [59, 61] as previously described [62] (Fig. 1). In adulthood, sMV− and control (cMV) females were mated with wild-type males (gestational day (GD)0). Throughout pregnancy and lactation, dams received an oral gavage of 10 mg/kg fluoxetine (FLX, Fluoxetine 20 PCH, Pharmachemie BV, the Netherlands) or vehicle (Veh, Methylcellulose, Sigma-Aldrich Chemie BV, Zwijndrecht, the Netherlands) daily at 11:00 a.m. (Fig. 1). FLX (5 mg/mL) and MC (1%) solutions were prepared with autoclaved water. Oral treatment was given by gently picking up the animal without restraint, and using flexible PVC feeding tubes (40 cm length, Vygon, Valkenswaard, the Netherlands) in order to minimize discomfort. It should be noted that our lab has been confronted with thus far unexplained mortality in about one fourth of MV females as a result of fluoxetine treatment, a phenomenon we have described elsewhere. [63]. None of the females included in the current study showed any signs of toxicity as a consequence of FLX, except for one. The pups from this sMV-FLX female who died at PND17 were cross-fostered to a same-treated mother and included in the study.

Fig. 1: Overview of study design.
figure 1

Maternal vulnerability (MV) rats were either exposed to stress (sMV) or control handled (cMV) early in life. In adulthood, sMV and cMV females were bred. Throughout pregnancy and lactation, from gestational day (GD)1 until PND21, females received a daily oral injection of either 10 mg/kg fluoxetine (FLX) or methylcellulose (Veh). This resulted in eight offspring groups: cMV-Veh males, cMV-Veh females, sMV-Veh males, sMV-Veh females, cMV-FLX males, cMV-FLX females, sMV-FLX males, and sMV-FLX females. Offspring brains were taken at PND21 for molecular analysis (N = 8/group).

Full size image

Offspring brain collection and nucleic acid isolation

Wild-type pups from N = 31 nests were used (cMV-Veh N = 9, sMV-Veh N = 8, cMV-FLX N = 8, and sMV-FLX N = 6). Between PND14 and PND20, ear punches were taken for SERT genotyping as described previously [62]. On PND21, approximately 24 h after the last FLX or Veh injection of the mother, N = 8 males and N = 8 female offspring per group (Fig. 1, Supplementary File 1) were killed by rapid decapitation. Brains were isolated, snap-frozen in isopentane (Acros Organics) on dry ice, and then stored at −80 °C.

On a Leica CM3050 cryostat, 200 μM coronal brain sections were sliced. To locate brain areas of interest, a stereotaxic atlas of the PND21 rat brain was used [64]. The medial PFC was collected from Bregma 2.8 mm to Bregma 1.8 mm (Supplementary Fig. 1A), and the BLA from Bregma −1.6 mm to Bregma −3.0 mm (Supplementary Fig. 1B), using a 2.0 mm punch (Harris Uni-Core). Samples were stored in 2.0 mL safe-lock tubes (Sarstedt) at −80 °C.

The AllPrep DNA/RNA Mini Kit (Qiagen) was used to simultaneously isolate DNA and RNA according to the manufacturer’s protocol, using 350 μL Buffer RLT Plus, 3.5 μL β-mercaptoethanol (Sigma-Aldrich) and 1.75 μL Reagent DX (Qiagen). Samples were lysed using a TissueLyser II (Qiagen) and 5 mm stainless steel beeds, 2 times for 2:00 at 30 Hz. Nucleic acid isolation concentration and purity were checked using a NanoDrop 2000 (Thermo Scientific), and RNA concentration and integrity were quantified using an RNA 6000 Nano Kit 2100 on a Bioanalyzer Instrument (Agilent Technologies, CA, USA).

RNA sequencing

RNA isolated from the PFC and BLA from N = 5 male and N = 5 female offspring per treatment group (N = 80 in total, Supplementary File 1) was selected for sequencing based on RNA quantity and quality. RNA integrity scores ranged from 8.2 to 9.6. Only one animal per litter was used. RNAseq, including library preparation, quality control and mapping of reads was performed by Novogene, Hong Kong.

Library preparation and sequencing

RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). Sequencing libraries were generated from 3 μg RNA per sample using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using HiSeq PE Cluster Kit cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform and 125 bp/150 bp paired-end reads were generated.

Quality control

Raw FASTQ reads were first processed through in-house Perl scripts. Clean reads were obtained by removing reads containing adapter, reads containing poly-N, and low quality reads. Final clean bases ranged from 5.7G to 7.8G per sample. To ascertain the quality of the reads, the Q20, Q30 and GC content were calculated. Q30 ranged from 88 to 92%, indicating a high quality.

Reads mapping

To map the reads to their respective genes, the Ensembl reference genome release 79 for Rattus Norvegicus was used [65]. Index of the reference genome was built using Bowtie v2.2.3 and paired-end reads were aligned to the reference genome using TopHat v2.0.12. The number of reads that could be mapped to the reference genome ranged from 80 to 87%. HTSeq v0.6.1 was used to count the number of reads mapped to each gene (Supplementary File 2). The number of reads that could be uniquely mapped to the reference genome ranged from 78 to 84%, and the number of non-splice reads ranged from 51 to 61%.

RNAseq data exploration

The Novogene read count file was used or further analyses using DESeq2 v1.22.2, Bioconductor release 3.8 [66] in R v3.5. First, the data was preprocessed with a variance stabilizing transformation and explored using heat maps of the highest expressing genes and the sample-to-sample distances, and through principal component analysis (PCA). The PCA plot showed that PFC and BLA samples cluster separately. One BLA sample did not cluster with any other samples, suggesting its gene expression pattern is entirely different (Supplementary Fig. 2A). It was excluded, leaving the cMV-Veh male BLA group with N = 4.

Gene set enrichment analysis

To examine changes in sets of related genes, Gene Set Enrichment Analysis (GSEA) v4.0.1 was used through the Broad Institute desktop application [67]. We ran GSEA on the control group (cMV-Veh) versus 1. the maternal adversity group (sMV-Veh) and 2. the FLX-exposed group (cMV-FLX) and 3. the maternal adversity- and FLX-exposed group (sMV-FLX) per combination of sex and brain area. The gene sets database was generated by generating a list of Gene Ontology (GO) terms for all genes in our data using the Biomart database system of Ensembl genes 95, dataset Rnor_6.0 [65]. The GSEA results yield a list of gene sets that differentiate between the two phenotypes, and also which genes within these sets contribute the most to the observed difference. The Normalized Enrichment Score (NES) is the degree to which the gene set is overrepresented in a phenotype. Significance cutoff was set at False Discovery Rate (FDR) <0.25.

Differential expression analysis

DESeq2 was used to identify the genes that were differentially expressed depending on maternal adversity exposure, FLX treatment, or a combination of both. Although male and female gene expression did not differ based on PCA (Supplementary Fig. 2B), we decided to analyze them separately based on literature and previous findings from our lab [63]. A sMV + FLX + sMV*FLX model was applied to every combination of sex and brain area (PFC males, BLA males, PFC females, BLA females). Significance cutoff was set at FDR < 0.1.

qPCR validation of RNAseq results

cDNA synthesis and qPCR

To validate key RNAseq results, we performed quantitative reverse transcriptase PCR (qRT-PCR) validation of four genetic targets in independent biological samples. RNA isolated from the PFC and BLA from N = 3 male and N = 3 female offspring per treatment group (N = 48 in total, Supplementary File 1) were used as input for complementary DNA (cDNA) synthesis using oligo(dT)18 primers and RevertAid H Minus Reverse Transcriptase (Thermo Scientific) according to the manufacturer’s protocol using a Veriti 96 well thermal cycler (Applied Biosystems).

qRT-PCR reactions were performed using one custom designed TaqMan® gene expression assay—Actb—and five pre-mixed assays—Gapdh, Cldn11, Cnp, Mag, Mbp—(Supplementary File 3, Thermo Scientific). Each reaction volume was 20 µL, including 10 µL iTaq Universal Probes Supermix (Bio-Rad) and 20 ng cDNA. Reactions were run in triplicate on 96 fast half skirt qPCR plates (Sarstedt). An Applied Biosystems 7500 Fast Real-Time PCR System (Life Technologies, the Netherlands) was used with the following conditions: 2 min 95 °C; 40 cycles of 15 s 95 °C, 1 min 60 °C.

qPCR analysis

PCR efficiencies and Cq values were calculated using LinRegPCR v. 2021.1 [68]. If the standard deviation (SD) of the mean Cq value of a triplicate was ≥0.3, the most deviating measurement was excluded. If the two remaining measurements still had an SD ≥ 0.3, we repeated the qPCR for that triplicate and used the new measurements. Because of the low yield of the starting material, we were unable to perform repeat qPCR reactions for some samples (Supplementary File 8). To calculate target gene expression, we used the mean normalized expression (MNE), based on the ratio between the Cq values of the target and the reference genes (Actb and Gapdh), and the PCR efficiency [69]. For BLA samples, RNA yield was especially limited and therefore we used only one reference gene (Actb) for these samples. Results are visualized as log(MNE).

R v4.0.2 was used for statistical analysis. Log(MNE) results were visually inspected using a Q–Q plot and deemed to be normally distributed. Two-way ANOVAs were used to analyze each target gene per combination of sex and brain region, with maternal adversity sMV and FLX as factors. Statistical significance was set at p < 0.05. Figures were produced using GraphPad prism v5.

DNA methylation analysis

Pyrosequencing was used to examine epigenetic regulation of genes of interest in the same samples that were used for RNA sequencing. Genomic DNA (240 ng) was bisulfite-treated using the EZ DNA methylation gold kit (Zymo Research, Leiden, The Netherlands) according to manufacturer’s protocol. This treatment converts cytosine but not 5-methylcytosine residues to uracil. Bisulfite-specific primers for the promoter regions of Cldn11, Cnp, Mag and Mbp were designed using PyroMark Assay Design Software 2.0 (Qiagen) (Supplementary Fig. 3; Supplementary File 4). HotStarTaq master mix (Qiagen, Hilden, Germany) was used for amplification of 1 μL bisulfite-treated DNA using the following steps: DNA polymerase activation (95 °C, 15 min), 3-step denaturation cycle (94 °C, 30 s), annealing (variable temperatures, 30 s), and extension (72 °C, 30 s) repeated for variable number of cycles (Supplementary Fig. 4). The final extension was performed at 72 °C for 7 min. Three samples (2 cMV-Veh male, 1 cMV-FLX female) failed to amplify during the PCR reaction, and were replaced by samples from the same treatment group but not previously used for RNAseq. Percentage methylation of selected CpG positions was determined using a PyroMark Q48 Autoprep Instrument (Qiagen) and PyroMark Q48 software (Qiagen). CpG positions 2 and 7 for Cnp and CpG position 3 for Mbp were excluded due to high peak height deviation. Since the DNA methylation of CpGs of the same gene correlated highly, we used the average DNA methylation per gene for further analyses.

DNA methylation percentages for every combination of sex and brain area were visually inspected using a Q–Q plot in R v3.5 and deemed to be normally distributed. GraphPad prism v8 was used for linear regression analyses and two-way ANOVAs, with statistical significance set at p < 0.05. Error bars in graphs represent SEM.

Results

Brain region- and sex-specific effects of perinatal fluoxetine and maternal adversity on myelin-related gene sets

First, we used GSEA to investigate the effects of maternal adversity and FLX exposure on the expression of gene sets. To this end, we compared the cMV-Veh group to the cMV-FLX group, the sMV-Veh group, and the sMV-FLX group for every combination of sex and brain area, yielding 12 comparisons (Supplementary File 5). The top 30 gene sets with the largest enrichment score for the cMV-Veh group versus the cMV-FLX group (effect of FLX) and versus the sMV-Veh group (effect of sMV) were plotted (Fig. 2A–D).

Fig. 2: Brain region- and sex-specific effects of perinatal fluoxetine and maternal adversity on gene expression.
figure 2figure 2

A The effect of fluoxetine exposure (cMV-FLX vs. cMV-Veh) in males. B The effect of maternal adversity (sMV-Veh vs. cMV-Veh) in males. C The effect of fluoxetine exposure (cMV-FLX vs. cMV-Veh) in females. D The effect of maternal adversity (sMV-Veh vs. cMV-Veh) in females. Myelin-related gene sets are highlighted in yellow. Asterisk (*) indicates FDR <0.25. Related to Supplementary File 5.

Full size image

FLX-exposed males

In males, several gene sets related to oligodendrocytes and myelin formation were positively associated with the cMV-FLX phenotype in the PFC, but negatively associated with the cMV-FLX phenotype in the BLA (Fig. 2A). More specifically, FLX exposure was associated with an upregulation of the gene sets Oligodendrocyte differentiation (NES = 2.38, FDR < 0.001), Oligodendrocyte development (NES = 2.25, FDR < 0.05), Myelin sheath (NES = 2.23, FDR < 0.01) and Myelination (NES = 2.12, FDR < 0.05) in the PFC, but with a downregulation of the gene sets Myelin sheath (NES = −2.66, FDR < 0.001), Oligodendrocyte differentiation (NES = −2.19, FDR < 0.05) and Myelination (NES = −2.18, FDR < 0.05) in the BLA (Figs. 2A, 3A).

Fig. 3: Gene expression related to “Myelin sheath” is associated with perinatal fluoxetine and maternal adversity in a brain region- and sex-specific manner.
figure 3

The black vertical lines correspond to individual genes in this gene set, and their horizontal position corresponds to how strongly they associate with either of the two phenotypes that are being compared. A The effect of fluoxetine exposure (cMV-FLX vs. cMV-Veh) in males. B The effect of maternal adversity (sMV-Veh vs. cMV-Veh) in males. C The effect of fluoxetine exposure (cMV-FLX vs. cMV-Veh) in females. D The effect of maternal adversity (sMV-Veh vs. cMV-Veh) in females. NES normalized enrichment score, FDR false discovery rate. N/A indicates that the “Myelin sheath” gene set was not in the top 100 gene sets for this comparison. Related to Supplementary File 7.

Full size image

sMV-exposed males

Maternal adversity had a similar effect on gene set enrichment in males; myelin-related gene sets were among the top gene sets positively associated with sMV in the PFC, while they were negatively associated with sMV in the BLA (Fig. 2B). In the PFC, maternal adversity was associated—although not with an FDR lower than 0.25—with an upregulation of the gene sets Myelination (NES = 1.83, FDR = 0.47), Oligodendrocyte development (NES = 1.83, FDR = 0.39) and Myelin sheath (NES = 1.61, FDR = 0.60). In the BLA, however, maternal adversity was associated with a downregulation of the gene sets Myelin sheath (NES = −2.57, FDR < 0.001), Oligodendrocyte differentiation (NES = −2.13, FDR < 0.01), Myelination in peripheral nervous system (NES = −2.11, FDR < 0.05) and Myelination (NES = −1.94, FDR < 0.05) (Figs. 2B,  3B). Maternal adversity was also negatively associated with several gene sets related to neuronal communication in the male BLA, such as Learning (NES = −2.39, FDR < 0.001) and Synapse (NES = −2.22, FDR < 0.01) (Fig. 2B).

FLX-exposed females

In females, FLX exposure was negatively associated with an array of gene sets related to general cell maintenance and proliferation in the PFC, such as Ribosome (NES = −2.40, FDR < 0.001) and Translation (NES = −2.29, FDR < 0.001) (Fig. 2C). In the BLA, there was only one significantly (FDR < 0.25) downregulated pathway in the cMV-FLX females, which was Cellular oxidant detoxification (NES = −1.97, FDR < 0.25) (Fig. 2C). Although not significantly, the gene set Myelin sheath in the female BLA showed the same negative association with FLX as in the male BLA (NES = −1.55, FDR = 0.47) (Fig. 3C).

sMV-exposed females

Maternal adversity was associated with changes in many pathways in the female PFC. For example, gene sets related to endocrine signaling such as Thyroid gland development (NES = −2.08, FDR < 0.05) and Response to glucocorticoid (NES = −2.03, FDR < 0.05) were downregulated in the PFC in sMV females (Fig. 2D). In the BLA, gene sets related to neuronal cell proliferation and synaptic plasticity were upregulated in sMV females, such as Postsynaptic density (NES = 2.17, FDR < 0.05) and Regulation of postsynaptic membrane potential (NES = 2.06, FDR < 0.05) (Fig. 2D). In addition, like in the male sMV BLA, the gene set Myelin sheath was downregulated in the female sMV BLA (NES = −1.90, FDR < 0.05) (Figs. 2D,  3D).

sMV-FLX-exposed offspring

To examine the effect of the combination of sMV and FLX on enrichment of gene sets, we also performed GSEA on the cMV-Veh vs. the sMV-FLX groups for each combination of sex and brain area (Supplementary File 5). Gene sets associated with exposure to both sMV and FLX are highly similar to those exposed to either treatment. In males, Myelin Sheath is the top gene set upregulated in sMV-FLX-exposed animals (NES = 2.20, FDR < 0.05), while in the male BLA it is downregulated (NES = −1.79, FDR < 0.25).

Differential expression analysis reveals that maternal adversity and perinatal fluoxetine exposure interact to affect myelin-related gene expression in the BLA

Interactions on the single-gene level we examined constructing a FLX * maternal adversity interaction model using DESeq2 for every combination of sex and brain area, yielding 4 models (Supplementary File 6). Overall, FLX was associated with the most significantly differentially expressed genes in the male BLA (17 genes at FDR < 0.1, Supplementary Figs. 4B, 5B), and maternal adversity was associated with the most significantly differentially expressed genes in the female PFC (29 genes at FDR < 0.1, Supplementary Figs. 4C, 5C). In the BLA of both sexes, the top genes that showed a FLX * maternal adversity interaction effect were related to myelin, although not significantly (FDR < 0.1) in females (Supplementary File 6). Some of these genes were also the main contributors to the effects found on the gene set Myelin sheath (Supplementary File 7, Fig. 4). Specifically, claudin-11 (Cldn11), 2′,3′-cyclic-nucleotide 3′-phosphodiesterase (Cnp), myelin−associated glycoprotein (Mag) and myelin basic protein (Mbp) showed a FLX * maternal adversity interaction effect in the BLA (Fig. 4B, D). These genes were downregulated by FLX in cMV animals, but upregulated by FLX in sMV animals in the BLA. In the PFC, no significant effects of FLX and maternal adversity on the expression of Cldn11, Cnp, Mag and Mbp were found (Supplementary File 6, Fig. 4A, C).

Fig. 4: Maternal adversity and perinatal fluoxetine exposure interact to affect myelin-related gene expression in the BLA.
figure 4figure 4

Heat maps of the expression of all individual genes in the gene set “Myelin sheath”, and bar plots of the gene counts of 4 selected genes: Cldn11, Cnp, Mag, and Mbp for: A The male PFC. B The male BLA. C The female PFC. D The female BLA. Related to Supplementary Files 2, 6, 7.

Full size image

In order to validate these single-gene results from the RNAseq analysis, we used N = 3 independent biological replicates per combination of sex and brain area (N = 48 in total) to measure gene expression of Cldn11, Cnp, Mag, and Mbp using qRT-PCR (Supplementary File 8). Although small sample sizes limited statistical power, an overall pattern of lower myelin-related gene expression in the BLA in FLX-treated groups was observed in both males (Supplementary Fig. 6B) and females (Supplementary Fig. 6D). This effect was statistically significant for Mag in the male BLA (FLX main effect p < 0.01, Supplementary Fig. 6B) and for Mbp in the female BLA (FLX main effect p < 0.05, Supplementary Fig. 6D).

DNA methylation of Mag and Mbp correlates to gene expression

We next examined DNA methylation levels around the transcription start sites of Cldn11, Cnp, Mag and Mbp (Supplementary Figs. 3, 9) in the same samples. Linear regression analyses showed there was a significant negative correlation between DNA methylation and Mag gene expression in both the PFC (R2 = 0.1689, p < 0.05, Fig. 5A) and the BLA (R2 = 0.2299, p < 0.01, Fig. 5B). There were no significant interaction effects of FLX * maternal adversity on Mag DNA methylation, but a strong trend in the female PFC (p = 0.059) (Fig. 5A, B). For Mbp, there was a negative correlation between DNA methylation and gene expression in the PFC (R2 = 0.2525, p < 0.01, Fig. 5C) but not in the BLA (R2 = 0.0000, p = 0.9597, Fig. 5D). There was a significant main effect of FLX on PFC Mbp methylation in females (Fig. 5C, D). There were no correlations between DNA methylation and gene expression of Cldn11 and Cnp (Supplementary Fig. 7), so these were not further examined for FLX * maternal adversity interactions.

Fig. 5: DNA methylation of Mag and Mbp correlates to gene expression.
figure 5

Scatter plots of DNA methylation percentages and gene counts, and bar plots of DNA methylation percentages in males and females for: A Mag in the PFC. B Mag in the BLA. C Mbp in the PFC. D Mbp in the BLA. Related to Supplementary Files 2, 8.

Full size image

Discussion

We found consistent alterations in gene sets related to myelination in juvenile rats exposed to perinatal fluoxetine and maternal adversity. In fluoxetine-exposed males, myelin-related gene expression was upregulated in the PFC, and downregulated in the BLA. In fluoxetine-exposed females, myelin-related genes were not affected in the PFC, but the BLA showed a similar response to the male BLA. Myelination of axons by oligodendrocytes plays a critical role in brain development and functioning; it influences neuronal circuit formation, allows for fast nerve conduction, and provides metabolic support to neurons [70]. Myelination starts in mid-to-late gestation in humans [71] and shortly after birth in rodents [72]. At PND21, the limbic system is among the last regions to become fully myelinated [72]. Our results suggest that developmental exposure to SSRIs and maternal adversity interfere with myelination at a critical time in development.

This aligns with earlier research showing effects on myelin-related genes after early SSRI exposure [73]. Moreover, constitutive SERT knockout decreased Cldn11, Mag and Mbp expression at PND8, increased it at PND14, and decreased it again at PND21 versus wildtypes in the male medial PFC [74]. This is contrary to our findings of enhanced PFC myelin-related gene expression after fluoxetine exposure at PND21, highlighting that total lack of the SERT is not equal to pharmacological SERT inhibition. Another study linked perinatal citalopram exposure to hypo- and hypermyelination in the adult corpus callosum, and higher levels of abnormal axons, with more severe deficits in males than in females [45]. Work in SERT knockout versus wild-type rats also indicated lower connectivity in the corpus callosum [75]. Overall, the effect of developmental SSRI exposure or constitutive SERT knockout on myelination appears to be highly dependent on brain region, age and sex.

SSRI exposure may be directly linked to alterations in myelination, as in vitro evidence showed that high serotonin levels damage immature oligodendrocytes [76]. Alternatively, high serotonin levels might indirectly affect myelination by affecting axons [76]. That is, myelination is linked to neuronal differentiation [77] and activity [71], while serotonin availability during neurodevelopment modulates the maturation of thalamocortical axons [78] and key brain circuits [43, 44]. It has been suggested that increased perinatal serotonin levels would lead to an accelerated rate of neurodevelopment [79]. In line with this, developmental citalopram treatment is associated with earlier onset of synaptogenesis [80] and in vitro evidence suggests that human embryo development is accelerated after treatment with fluoxetine [81]. The current findings may also be reflective of accelerated brain maturation.

Further support of a fairly unspecific effect of perinatal fluoxetine exposure on the rate of brain maturation is provided by our finding that the effects of maternal adversity are remarkably similar. This was true especially in the BLA, where myelin-related genes were also downregulated in the maternal adversity offspring. In line with this, the developing amygdala in both humans and animals is vulnerable to the effects of prenatal stress on functional and structural connectivity [82]. Although we did not find significant effects of maternal adversity on myelin-related gene expression in the medial PFC, other studies have identified an acute increase in myelination-related genes including Mag, Mbp, Cldn11, and Cnp after early postnatal stress at PND15 [83], but a decrease in these genes and associated proteins in adulthood [84]. This again suggests an effect on brain maturation, evidenced by a precocious differentiation of oligodendrocytes in early life [83].

Our lab previously identified interactions between maternal adversity and fluoxetine on a behavioral level [63]. In the current study, maternal adversity and fluoxetine interacted to affect Mag, Mbp, Cldn11, and Cnp gene expression, especially in the BLA. The group receiving both treatments resembled the control group, suggesting that fluoxetine may normalize some effects of in utero exposure to maternal adversity. This is in line with earlier work showing that in utero citalopram exposure prevented the increased fetal frontal cortex serotonin levels induced by chronic prenatal stress [85]. Similarly, perinatal fluoxetine exposure largely reversed the effects of pre-gestational maternal stress on increased serotonin levels in the PFC [55], and stress-coping behavior and hippocampal neurogenesis in adolescence [86].

Preclinical evidence suggests that males are more likely to develop long-term behavioral effects after perinatal SSRI exposure than females [41]. The current results also point to a stronger phenotype in males. In females, the dampening effect of fluoxetine on myelin-related genes in the BLA is present but non-significant, whereas myelin-related gene sets in the PFC are unaffected. Connectivity studies in humans suggest that white matter development may occur earlier in girls than in boys [71], potentially modulating the ability of SSRIs to alter this process. In contrast, human studies have shown that girls usually are more affected by exposure to maternal depressive symptoms during pregnancy than boys on neuroimaging outcomes [31,32,33]. Accordingly, our results show that the transcriptomic state of the male PFC at PND21 is rather insensitive to the effects of maternal adversity, while a substantial number of genes and gene sets are affected in the female PFC.

Epigenetic regulation of gene expression, in particular by DNA methylation at CpG dinucleotides, is thought to link early life experiences to later-life health and behavior [87]. SSRI exposure did not affect whole-genome DNA methylation in human cord blood [88], and effects of maternal depression on offspring DNA methylation are difficult to replicate [89]. We investigated the brain directly, and identified negative correlations between gene expression and DNA methylation for Mag in the PFC and BLA, and for Mbp in the PFC. No such correlations were found for Cldn11 or Cnp, potentially because of low overall DNA methylation. DNA methylation is suggested to contribute to the differentiation of oligodendrocyte precursor cells and the process of myelination [90]. In agreement with the current results, a negative correlation between Mag expression and DNA methylation in oligodendroglial-lineage cells has been reported [90]. Overall, we provide the first indication that fluoxetine-induced changes in myelination may be mediated by an epigenetic mechanism.

The most important strengths of our study are the investigation of both the PFC and BLA, and both sexes, revealing clear brain region- and sex-dependent effects. In addition, we directly correlated gene expression and DNA methylation. Limitations are the inclusion of only one time-point, precluding inferences about the long-term effects of perinatal SSRI exposure and maternal adversity. In addition, we did not validate our results at the protein level. Moreover, we used SERT+/− rats as a vulnerable model for early life stress. These animals have reduced SERT expression which may alter the responses to fluoxetine. We previously showed that fluoxetine and the metabolite norfluoxetine are within normal ranges during pregnancy and lactation in SERT+/− dams, suggesting that they metabolize FLX similarly to wild-type animals [63], but further research is necessary to confirm this. Finally, we did not assess maternal care behavior. FLX may influence maternal care patterns and thereby influence myelination patterns. Future research should include cross-fostering experiments to disentangle the effects of FLX treatment from maternal care influences.

Human studies often observe similar neurodevelopmental outcomes in children exposed to SSRIs and children exposed to unmedicated depression, evoking questions about the source of the effects. Here, we observe striking similarities between the effects of SSRI exposure and maternal adversity on rat gene expression, suggesting that they both contribute. The next step is to elucidate whether these similarities are due to increased brain serotonin levels. Furthermore, the time course of changes in myelination and white matter development after perinatal insults needs to be delineated. Both the causes of altered myelination, such as differences in neuronal activity [83], and its consequences, such as altered brain circuit formation [91], deserve to be further explored.