A protein-coding gene expression atlas from the brain of pregnant and non-pregnant goats

Background: The brain is an extraordinarily complex organ with multiple anatomical structures involved in highly specialized functions related with behavior and physiological homeostasis. Our goal was to build an atlas of protein-coding gene expression in the goat brain by sequencing the transcriptomes of 12 brain regions in seven female Murciano-Granadina goats, from which three of them were 1-month pregnant.

Results: Between 14,889 (cerebellar hemisphere) and 15,592 (pineal gland) protein-coding genes were expressed in goat brain regions, and most of them displayed ubiquitous or broad patterns of expression across tissues. Principal component analysis and hierarchical clustering based on the patterns of mRNA expression revealed that samples from certain brain regions tend to group according to their position in the anterior-posterior axis of the neural tube, i.e., hindbrain (pons and medulla oblongata), midbrain (rostral colliculus) and forebrain (frontal neocortex, olfactory bulb, hypothalamus, and hippocampus). Exceptions to this observation were cerebellum and glandular tissues (pineal gland and hypophysis), which showed highly divergent mRNA expression profiles. Differential expression analysis between pregnant and non-pregnant goats revealed moderate changes of mRNA expression in the frontal neocortex, hippocampus, adenohypophysis and pons, and very dramatic changes in the olfactory bulb. Many genes showing differential expression in this organ are related to olfactory function and behavior in humans.

Conclusion: With the exception of cerebellum and glandular tissues, there is a relationship between the cellular origin of sampled regions along the anterior-posterior axis of the neural tube and their mRNA expression patterns in the goat adult brain. Gestation induces substantial changes in the mRNA expression of the olfactory bulb, a finding consistent with the key role of this anatomical structure on the development of maternal behavior.

1 Introduction

The mammalian brain is an extraordinarily complex organ integrating multiple highly specialized structures involved in the regulation of memory, behavior, learning, sensory function, motor skills and body homeostasis, amongst others. During embryo development, three vesicles emerge from the neural tube (Chhetri and Das, 2021), which roughly correspond to the forebrain (prosencephalon), midbrain (mesencephalon) and hindbrain (rhombencephalon). Further subdivisions take place in the forebrain (telencephalon and diencephalon) and hindbrain (metencephalon and myelencephalon), resulting in the formation of five encephalic vesicles (Chhetri and Das, 2021). In the bovine embryo, these five encephalic vesicles are visible at 24 days after conception, but they do not become fully developed until 110 days of gestation (Ferreira et al., 2018). Multiple anatomical structures with specialized functions (Chhetri and Das, 2021) are subsequently derived from the telencephalon (cerebral hemisphere, basal ganglia, hippocampus, olfactory bulb, lateral ventricles, etc.), diencephalon (thalamus, hypothalamus, pineal body, neurohypophysis, infundibulum, third ventricle, etc.), mesencephalon (rostral colliculus, tegmentum, crus cerebri, cerebral aqueduct, etc.), metencephalon (pons, cerebellum, upper part of the fourth ventricle, etc.) and myelencephalon (medulla oblongata, spinal cord lower part of the fourth ventricle, etc.).

Transcriptomic analyses may hold the key to significantly advance our knowledge about the biological functions of brain regions. In mice, patterns of gene expression have been used to establish a molecular atlas of the adult brain (Ortiz et al., 2020), and as much as 737 brain structures have been identified with microscopy techniques complemented with other approaches, providing a comprehensive view about the high functional complexity of this organ (Erö et al., 2018). An important feature of the brain is its high plasticity, so any atlas of gene expression is necessarily dynamic, not only in space but also in time. For instance, pregnancy in mice is associated with extensive changes in gene expression in the neocortex, cerebellum, hippocampus and hypothalamus, and there is evidence that several of such modifications might be long lasting (Ray et al., 2016; Arbeitman, 2019).

Very few atlases of brain gene expression have been generated in domestic animals. Recently, Sjöstedt et al. (2020) investigated the gene expression profiles of 10 major mammalian brain regions in humans, mice, and pigs. They found that global transcriptomic profiles are, in general, well conserved in these three species, with cerebrum and brainstem regions clustering apart in hierarchical trees, and the cerebellum showing a highly divergent profile of gene expression (Sjöstedt et al., 2020). The mRNA expression of several brain regions has also been reported in sheep (Clark et al., 2017) and cattle (Harhay et al., 2010), although not comprehensively or systematically. In goats, mRNA expression of two neural tissues, frontal lobe cortex and cerebellum, has also been characterized by RNA-Seq (Muriuki et al., 2019). The main goal of the current work was to establish an atlas of protein-coding gene expression of the caprine brain by sequencing the transcriptomes of 12 encephalic regions in seven Murciano-Granadina goats. Given that three of these goats were pregnant at the time of slaughter, we have also investigated whether gestation affects the transcriptomic profiles of the 12 brain regions under study.

2 Materials and methods

2.1 Sample collection

Seven adult non-lactating Murciano-Granadina goats raised in the experimental farm of the Faculty of Veterinary Sciences at the Universitat Autònoma de Barcelona (UAB) were slaughtered because of reasons unrelated with this project. Sampled goats were kept under the same management and environmental conditions and had similar ages (6.28 ± 1.38 years). Besides, three of these goats were 1-month pregnant at the time of slaughtering. To minimize pain, goats were administered pentobarbital (150 mg/kg) in the jugular vein. A lethal intravenous dose of sodium pentobarbital causes a mammal to lose consciousness within seconds and results in clinical death within just minutes, but we cannot completely rule out the possibility that the euthanasia procedure used in the current work might have altered the brain mRNA expression of certain genes. Since animals were killed due to the routine culling process implemented in this experimental farm, no permission from the Ethics Committee on Animal and Human Experimentation at UAB was required.

After slaughtering, goats were transported to the Necropsy room of the Faculty of Veterinary Sciences at UAB. The cranial vault of each goat was opened with a bone saw and twelve brain anatomical structures were carefully dissected and biopsied by an expert anatomist. Samples were drawn from the adenohypophysis, cerebellar hemisphere, cerebellar trunk or peduncle, frontal neocortex, hippocampus, hypothalamus, medulla oblongata, neurohypophysis, pineal gland, pons, olfactory bulb and rostral colliculus. Biopsies were immediately submerged in RNAlater (Thermofisher Scientific, Barcelona, Spain) to be stored at -80°C until processing.

To purify total RNA, tissue samples were mixed with 1 ml QIAzol (QIAGEN Inc., Barcelona, Spain) and homogenized using the Lysing Matrix D reagent (MP Biomedicals, Santa Ana, CA) in a Precellys 24 tissue homogenizer (Bertin Instruments, Rockville, MD). The extraction of total RNA was performed using the RNeasy lipid tissue mini kit (QIAGEN Inc., Barcelona, Spain) following the protocol described by the manufacturer. The concentration and purity of extracted RNA molecules were analyzed using the Nanodrop ND-1000 spectrophotometer (Thermofisher Scientific, Barcelona, Spain) and RNA integrity was assessed with a Bioanalyzer-2100 equipment (Agilent Technologies, Santa Clara, CA) using the RNA 6000 Nano Kit 4.2 (Agilent Technologies, Santa Clara, CA).

2.2 Sequencing of total RNA

Paired-end sequencing (2 x 50 bp) of total RNA was carried out at the Centre Nacional de Anàlisi Genòmica (CNAG). Sequencing methods have been reported by Guan et al. (2020). Briefly, the RNA-Seq library was prepared with the KAPA Stranded mRNA-Seq Illumina Platforms Kit (Roche, Sant Cugat, Spain) by using 500 ng total RNA as template. Oligo-dT magnetic beads were used to enrich the poly-A fraction and subsequently, RNA was fragmented. Strand cDNA synthesis was performed in the presence of dUTP to enforce strand-specificity. The blunt-ended double stranded cDNA was 3′-adenylated and ligated to Illumina adaptors with unique dual indexes and unique molecular identifiers (Integrated DNA Technologies, Coralville, IA). Enrichment of the ligation product was ensured by performing 15 cycles of polymerase chain reaction amplification. An Agilent 2100 Bioanalyzer equipment was employed to verify the quality of the final library by using the DNA 7500 assay (Agilent Technologies, Inc., Santa Clara, CA). Library sequencing was carried out with a HiSeq 4000 equipment (Illumina, San Diego, CA) in accordance with the protocol for dual indexing advised by the manufacturer. Image analysis, base calling and quality scoring of the sequencing run were checked with the Real-Time Analysis (RTA 2.7.7) tool (Illumina, San Diego, CA) and FASTQ sequence files were subsequently generated.

2.3 Quality control, alignment and quantification

The quality of the sequences was assessed with the FastQC software v.0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and adapters were trimmed using the TrimGalore v.0.6.6 tool (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). In addition, reads with more than five ambiguous bases (Ns) were removed. Filtered sequences were aligned to the goat ARS1 reference genome (Bickhart et al., 2017) with the HISAT2 v.2.2.1 aligner using paired end option (-1 and -2) and default parameters (Kim et al., 2015). Gene assembly and expression quantification were performed with StringTie v.2.1.0 by extracting, with the -e option, those genes annotated in the reference GTF file as reported in a previous protocol (Pertea et al., 2016). The estimated count matrix was obtained from coverage values of each feature (genes in our case) using the dedicated script “prepDE.py” from the StringTie pipeline and based on the following formula (https://github.com/gpertea/stringtie/blob/master/prepDE.py):

$Readspergene=\frac{coverage\times genelength}{readlength}$

2.4 Construction of a brain atlas based on the expression of protein-coding genes

Genes expressed in each tissue were cataloged according to their Ensembl biotype. To remove systematic technical effects (e.g., library size), count data from protein-coding genes were normalized to CPM with the edgeR package (Robinson et al., 2009). Finally, data were log₂ transformed adding a pseudo-count of 1. Each tissue was normalized separately. Normalization factors were estimated by jointly considering pregnant and non-pregnant goats because library size features were quite similar in both groups. To perform downstream analyses, only genes with a CPM over 0.5 in at least two samples were considered.

In order to detect errors during the sampling procedure, which might alter our results (accidental cross-contaminations when retrieving samples, errors when annotating the anatomical origin of samples etc), the quality of the regional sampling of the brain was evaluated in two ways. First, Euclidean distances were estimated in a pairwise manner based on the corrected level of gene expression (log₂(CPM + 1)) considering all protein-coding genes (17,054 genes) to identify outlier samples. The distribution of distances between samples from the same tissue was analyzed using the boxplot.stats function from the grDevices package in R (https://r-universe.dev/manuals/grDevices.html), and samples deviating more than 1.5 times from the interquartile range were considered as outliers and, in consequence, eliminated. Additionally, we built a heatmap based on the gene expression level of 61 genes selected by their high tissue-specificity in a set of 12 human tissues (https://www.proteinatlas.org/, version 22) equivalent to those investigated in goats herewith. This analysis was performed to make sure that the expression patterns of the goat samples match what is expected based on human transcriptomic data.

A PCA was constructed by using the prcomp function from the stats R package (R Core Team, 2015). Additionally, a hierarchical clustering dendrogram was built based on Euclidean distances and applying the Ward D clustering method with the hclust function from the R stats package (Murtagh and Legendre, 2014). For each brain region, the top 1,000 genes with the highest expression levels were selected to perform pathway enrichment analyses with the Enrichr R package (Chen et al., 2013). The KEGG release 99.0 human database (Kanehisa et al., 2016) was used as reference to annotate pathways enriched in the 12 sets of genes (one for each brain region). In order to identify such enriched pathways, a combined score (c) was estimated from the output of a traditional Fisher exact test and a Z-score was estimated from expected ranks and variances of the set of genes. In accordance with Chen et al. (2013) the following formula was applied:

$c=\log \left(p\right)\times z$

where c is the combined score, p corresponds to p-values estimated with the Fisher exact test, and z is the Z-score from the expected rank.

The tissue-specificity expression of protein-coding genes was assessed by calculating the tissue specificity index tau (τ) defined by Yanai et al. (2005):

$\mathrm{\tau }=\frac{{\sum }_{i=1}^{\mathrm{n}}\left(1-x_i\right)}{\mathrm{n}-1}$

where n is the number of brain regions, and x_i corresponds to the expression profile component normalized by the maximal component value (Yanai et al., 2005). This calculation was performed using the tspex v.0.6.1 program (Camargo et al., 2020). A τ-value of 0 would correspond to a housekeeping gene with ubiquitous tissue expression. In contrast, genes with a restricted regional expression show τ-values above 0.85, while genes with τ = 1 can be assumed to be expressed in a single brain region.

2.5 Differential expression and pathway enrichment analyses

Before doing the differential expression analysis, the level of variability across replicates from the same condition (pregnant or non-pregnant) was assessed by calculating the coefficient of variation (CV) of the expressed genes and genes differentially expressed between conditions:

$CV=\frac{\sigma }{\mu }$

where μ is the average expression of each gene and σ is the corresponding standard deviation.

The software DEseq2 (Love et al., 2014) was used to perform a differential expression analysis comparing the transcriptomic profiles of pregnant vs. non-pregnant goats. Genes with a number of estimated counts below 10 were removed. Correction for multiple testing was applied using the false discovery rate (FDR) method (Benjamini and Hochberg, 1995). We considered that a gene is differentially expressed when two conditions are met: 1) absolute logarithm of the fold change |log₂FC| > 0.58, and 2) Statistical significance threshold: q-value < 0.05. Differentially expressed genes were used to perform pathway enrichment analyses using the Enrichr R package (Chen et al., 2013). The non-pregnant goats were considered as the baseline reference, meaning that any upregulation or downregulation of gene expression is referred to changes in the pregnant goats with regard to the non-pregnant. The annotation of differentially expressed genes (DEGs) was performed as explained in Section 2.4, while the enrichment pathway analyses were carried out independently in the sets of upregulated and downregulated genes. The level of enrichment was estimated by calculating a c score (Chen et al., 2013), as previously defined.

3 Results

3.1 Sequencing of total RNA from 12 goat brain regions

The average RNA integrity number (RIN) of the 84 RNA extractions was 7.53 ± 0.62, ranging from 6.2 to 9 (Table 1). These values are in agreement with the conservative threshold of RNA degradation (RIN = 6.4–7.9) defined by Gallego Romero et al. (2014) for the in-depth analysis of RNA transcripts. As shown in Table 1, the average sequencing depth ranged from 37.78 to 41.71 million paired-end reads per sample across all 12 brain regions, and the average alignment rate ranged from 84.44% (adenohypophysis) to 94.35% (pons). The read counts per gene per sample, including pregnant and non-pregnant goats are available in Supplementary Table S1.

TABLE 1

TABLE 1. Quality of RNA samples and number of genes expressed in 12 brain goat regions. Average RNA integrity number (RIN), alignment rate and number of expressed genes (CPM > 0.5 in at least two samples) per brain region in seven Murciano-Granadina female goats.

3.2 Clustering of brain regions according to their mRNA expression profile

The total number of protein-coding loci expressed (counts per million, CPM > 0.5 in at least two samples) in each tissue ranged from 14,889 genes for the cerebellar hemisphere to 15,592 genes for the pineal gland. Considering all brain regions, 17,054 protein-coding genes passed the quality control (QC) filters and were expressed in at least two samples. One pineal gland sample deviated more than 1.5 times from the interquartile range, so it was identified as an outlier and excluded from any further analysis. As a complementary measure to detect sampling errors, we built a heatmap based on the expression levels of 61 marker genes which in humans show high tissue specificity for each one of the brain regions considered in this research. As shown in Supplementary Figure S1, the expression patterns of the samples based on the selected 61 genes is concordant across replicates from the same brain region, thus evidencing that dissection and sampling were correctly performed. In general, all brain regions display high expression levels of the marker genes in close correspondence to expectations based on human data, with the only exception of pons and medulla oblongata. Samples from the medulla oblongata display high mRNA levels of genes PHOX2A and SLC6A2, which are considered markers for the pons based on human data. In congruence with our results, the expression of these genes is reported to be high in the medulla oblongata of mice (Sunkin et al., 2013). Furthermore, the gene EXD1 selected as a marker for pons was highly expressed in the pineal gland of goats. This gene is also expressed at high levels in the pituitary gland of pigs, and at moderate levels in multiple brain structures of pigs and mice (Sjöstedt et al., 2020). We also evaluated the variability of the mRNA expression patterns across replicates by calculating the corresponding coefficient of variation per group. In general, the coefficient of variation was lower within the pregnant and not pregnant groups than when all samples were pooled together in a single group (Supplementary Figure S2). This result was even more evident when analyzing the sets of differentially expressed genes in brain regions moderately or highly affected by pregnancy (Supplementary Figure S2B).

In close agreement with data reported above, the PCAs, with (Figure 1A) and without (Figure 1B) cerebellum and glandular samples, and the dendrogram (Figure 2) showed that most samples clustered according to their region of origin. We also observed that four regions, namely adenohypophysis, neurohypophysis, pineal gland, and cerebellum (hemisphere and trunk), displayed highly divergent gene expression patterns when compared to the remaining encephalic regions (Figure 1A; Figure 2). The dendrogram indicated that pons and medulla oblongata, which are hindbrain structures, tended to group together. Moreover, frontal neocortex, olfactory bulb, and hippocampus, which derive from the forebrain telencephalon vesicle, also tended to group together. In strong contrast, cerebellum trunk and hemisphere clustered far apart from any other brain region (including pons despite its anatomical proximity), and adenohypophysis (derived from the oral ectoderm) and neurohypophysis (derived from the diencephalon), which have completely different embryonic origins, also grouped jointly and clearly separated from the remaining brain structures.

FIGURE 1

FIGURE 1. Principal component analysis (PCA) of samples from 12 goat brain regions. (A) Principal component analysis (PCA) of samples from 12 brain regions (adenohypophysis, cerebellar hemisphere, cerebellar trunk, frontal neocortex, hippocampus, hypothalamus, medulla oblongata, neurohypophysis, pineal gland, pons, olfactory bulb and rostral colliculus) retrieved from seven Murciano-Granadina female goats. (B) Same PCA without samples from adenohypophysis, cerebellar hemisphere, cerebellar trunk/peduncle, neurohypophysis, and pineal gland. Samples obtained from pregnant goats have been marked with a black dot.

FIGURE 2

FIGURE 2. Hierarchical clustering of samples from 12 brain regions retrieved from seven Murciano-Granadina goats based on their mRNA expression profiles. Dendrogram displaying the clustering patterns of samples from 12 brain regions. In most cases, samples cluster according to their region of origin and, moreover, regional clusters tend to group according to the forebrain, midbrain or hindbrain embryonic vesicle they originate from. Exceptions to this general observation are hypothalamus, that groups with rostral colliculus, and two main groups of outliers: cerebellum and glandular tissues (pineal gland and hypophysis). In the right part of the figure, we indicate the embryonic vesicle from which each region is derived. In an initial developmental stage of mammals, three forebrain (prosencephalon), midbrain (mesencephalon) and hindbrain (rhombencephalon) vesicles are formed. Subsequently, prosencephalon is subdivided into two further vesicles (telencephalon and diencephalon) and so does rhombencephalon (metencephalon and myelencephalon). Adenohypophysis does not originate from any of these five vesicles, but from the oral ectoderm. Samples obtained from pregnant goats have been marked with a black triangle on the side.

3.3 Regional specificity of mRNA expression

The distribution of the tissue-specificity index τ-values was highly skewed to the left (Figure 3A). Indeed, 5,655 (33.16%) protein-coding genes showed τ-values below 0.15, indicative of a highly ubiquitous pattern of expression across all brain regions (Supplementary Table S2A). In contrast, 10,036 (58.84%) protein-coding genes displayed intermediate τ-values (τ = 0.15–0.85) and 1,363 protein-coding genes (7.99%) had a highly region-specific profile (τ > 0.85) of mRNA expression (Figures 3A, B). A heatmap based on the level of expression of tissue-specific genes (τ ≥ 0.85) identified for each tissue is shown in Supplementary Figure S3.

FIGURE 3

FIGURE 3. Tissue specificity of protein-coding genes expressed in 12 goat brain regions. (A) Histogram representing the tissue specificity of 17,054 protein-coding genes expressed in 12 goat brain regions (in at least one sample) and surpassing QC filters. The number of genes is indicated in the y-axis, while Tau specificity scores (τ, see Materials and Methods for details) are shown in the x-axis. The dashed line indicates = 0.85 threshold that defines tissue-specificity. (B) Number of genes with tissue-specific expression (τ > 0.85) in each of the 12 brain regions under study.

As shown in Supplementary Table S2B, the gene pathway enrichment analysis for loci with an ubiquitous pattern of expression (τ = 0–0.15) revealed that they are mainly involved in general biological processes such as RNA transport (combined score, c = 110.230, q-value = 1.16 E-13), autophagy (c = 105.214, q-value = 5.66 E-12), oxidative phosphorylation (c = 110.941, q-value = 3.34 E-12), spliceosome (c = 124.950, q-value = 1.16 E-13), proteasome (c = 185.773, q-value = 5.52 E-10) and ubiquitin mediated proteolysis (c = 127.938, q-value = 1.56 E-13), as well as in the development of neurological diseases such as amyotrophic lateral sclerosis (c = 213.777, q-value = 5.57 E-26), Parkinson disease (c = 195.845, q-value = 1.86 E-21), and Huntington disease (c = 142.384, q-value = 3.91 E-19), to mention a few.

As indicated in Supplementary Table S3 and Supplementary Figures S4–S15, we have also investigated whether genes with highly region-specific gene expression patterns (τ = 0.85–1, listed in Supplementary Table S3A) are enriched in specific functional pathways in all analyzed regions (Supplementary Table S3B–S3M). The neuroactive ligand-receptor interaction pathway was highly enriched in genes with τ = 0.85–1 in several of the analyzed regions, i.e., frontal neocortex, hippocampus, hypothalamus, and olfactory bulb. Other pathways associated with neuronal transduction were also enriched in genes with a tissue-restricted expression, e.g., serotonergic synapse (cerebellar hemisphere and cerebellar trunk with q-values < 0.05 and c of 74.868, and 65.095, respectively) and GABAergic synapse (cerebellar hemisphere and trunk, with q-values < 0.05 and c of 69.725 and 60.823, respectively). Besides, the cholinergic synapse pathway showed enrichment in the rostral colliculus (c = 43.169, q-value = 0.05). Several brain regions showed an enrichment pattern that agreed well with their biological functions (Supplementary Tables S3B–S3M). For instance, in the pineal gland the phototransduction pathway showed a c of 1101.90 (q-value = 1.06 E-13).

3.4 KEGG pathway enrichment analysis of highly expressed genes

The results of the Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis for the 1,000 mRNA genes with the highest expression in each region (Supplementary Tables S4A–S4K) revealed several pathways broadly shared across brain regions and also reaching strong statistical significance (q-value < 0.05), e.g., synaptic vesicle cycle, endocrine and other factor-regulated calcium reabsorption, endocytosis, phagosome, long-term potentiation, dopaminergic synapse, glutamatergic synapse, circadian entrainment, gap junction, gastric acid secretion and adrenergic signaling in cardiomyocytes. We also detected many pathways related to neurological conditions that attained high statistical significance (q-value < 0.05) and were represented in a broad array of encephalic regions, e.g., Alzheimer disease, Huntington disease, Parkinson disease, prion disease, amyotrophic lateral sclerosis, and, to a much lesser extent, spinocerebellar ataxia. Pathways related with Vibrio cholerae and Salmonella infection were also commonly found. Protein processing in endoplasmic reticulum (q-value < 0.05), RNA transport (q-value < 0.05) and ribosome (q-value < 0.05) were amongst the most significant pathways in the adenohypophysis (Supplementary Table S4A) and neurohypophysis (Supplementary Table S4H), but not in the other regions. As said before, in the pineal gland (Supplementary Table S4J) phototransduction was also a highly significant pathway (q-value = 0.006) while in the remaining brain regions it did not reach statistical significance.

3.5 Analysis of differential gene expression in pregnant and non-pregnant goats

Samples from the pineal gland were not included in the differential expression analysis, since the removal of an outlier sample decreased the number of replicates from the pregnant group to two, which is insufficient to make a sound statistical analysis. When comparing pregnant vs. non-pregnant goats, brain regions showed strong differences in terms of the number of differentially expressed genes. Six regions displayed little changes in their expression levels in response to 1 month-pregnancy: in the cerebellar hemisphere, cerebellar trunk, hypothalamus, medulla oblongata, neurohypophysis and rostral colliculus, only 2, 1, 1, 12, 4 and 1 DEGs were identified, respectively (Figure 4B, Supplementary Table S5). In strong contrast, we observed remarkable changes in the expression profiles of five brain regions (Figure 4A, Supplementary Table S5): adenohypophysis (201 DEGs, 13 downregulated and 188 upregulated in pregnant goats, Supplementary Table S5A), frontal neocortex (82 DEGs, 37 downregulated and 45 upregulated in pregnant goats, Supplementary Table S5D), hippocampus (70 DEGs, 15 downregulated and 55 upregulated in pregnant goats, Supplementary Table S5E), pons (190 DEGs, 57 downregulated and 133 upregulated in pregnant goats, Supplementary Table S5J) and, most remarkably, olfactory bulb (1207 DEGs, 381 downregulated and 826 upregulated in pregnant goats, Supplementary Table S5I). Pathways significantly associated with DEGs in each tissue are shown in Supplementary Tables S6A–S6E and Supplementary Figures S16–S20. In the olfactory bulb (Supplementary Table S6D), the majority of biological functions enriched in the set of upregulated DEGs in pregnant goats were associated with axon guidance (c = 89.246; q-value = 7.29 E-07), dopaminergic synapse (c = 79.526; q-value = 6.07 E-06), neuroactive ligand-receptor interaction (c = 53.217; q-value = 5.27 E-06), cholinergic synapse (c = 55.712; q-value = 1.32 E-04) and phenylalanine, tyrosine and tryptophan biosynthesis (c = 53.087; q-value = 8.51 E-02). In contrast, downregulated DEGs in the olfactory bulb of pregnant goats were enriched in pathways related to nervous development and control of cellular activities, including migration, differentiation, and proliferation (Supplementary Table S6D). A heatmap based on the expression levels of genes differentially expressed in five tissues moderately or highly affected by pregnancy (i.e., adenohypophysis, frontal neocortex, hippocampus, olfactory bulb, and pons) is shown in Supplementary Figure S21. It can be seen that the main determinant of gene expression patterns is the region-of-origin of samples rather than pregnancy.

FIGURE 4

FIGURE 4. Volcano plots of differentially expressed genes in eleven brain regions when comparing pregnant and non-pregnant goats. (A) Brain regions little affected by pregnancy, i.e., cerebellar hemisphere, cerebellar trunk/peduncle, hypothalamus, medulla oblongata, neurohypophysis and rostral colliculus. (B) Brain regions moderately or strongly affected by pregnancy, i.e., adenohypophysis, frontal neocortex, hippocampus, olfactory bulb, and pons. Genes with a fold change below -1.5 and a q-value below 0.05 are depicted in blue, while genes with a fold change over 1.5 and a q-value < 0.05 are depicted in red. Grey dots represent genes that do not display differential expression.

By making a literature search for each one of the genes showing differential expression in the olfactory bulb of pregnant vs. non-pregnant goats, we observed that many of them are related to human behavioral traits (Supplementary Table S7), including, but not limited to: maternal behavior (e.g., downregulated: DDC; upregulated: DBH, DRD1, NTS and HTR2A), affective behavior, sociability and exploration (e.g., downregulated: FGFR4, GRIP2, NTRK2 and RYR3; upregulated: HTR2C, FBX O 45, GPR3, KCNQ2 and PLXNA2), anxiety and depression (e.g., downregulated: NOS2, PAN2 and H3-3B; upregulated: BSCL2, HRH1, SIK2, PDYN, GLRB and NRN1), autism (e.g., downregulated: MOCOS, CACNA1D, MBD6 and AUTS2; upregulated: RAB39B, BTBD11, KCNQ3, CDH9, CADPS2, DOCK4 and EXT1), aggression (e.g., downregulated: NOS1; upregulated: HRH3, PRNP, HNMT and GRIA3), cognition, memory and learning (e.g., downregulated: EPHA10, STAT5 and MMP28; upregulated: EPHA6, SYP, SORBS2, ARHGEF4, HCN1, CAMK2N2, MMP17, BTBD9, CLSTN3, STAU2RIMKLA, PAK6, SLC22A4 and NEURL1), response to stress (e.g., downregulated: IFIT1; upregulated: PDYN, EPOP, DPYSL2 and HCN2), feeding behavior (e.g., upregulated: NELL2, CXCL14, GPR45, GPR162, NPY and ACBD7) and diverse neuropsychiatric disorders (e.g., downregulated: SLITRK6, IL1RAPL1, BAHCC1, WDR62 and SLC6A1; upregulated: SCN1A, GABBR2, RTN4R, ACOT7, KIF5A. HECW2, NRG1, SNCA, ATP8A2, WFS1, HTR5A, STX1A and JPH3).

As shown in Supplementary Table S8, our literature search also revealed several DEGs associated with the migration and maturation of olfactory bulb interneurons (e.g., PROKR2 and SALL3), olfactory bulb angiogenesis (e.g., ANOS1), olfactory bulb morphogenesis and establishment of a functional olfactory neural circuitry (e.g., ARX, TSHZ1, OLFM2 and FGFR1), axonal growth and guidance (e.g., RAP1GAP2, SHH, GLI3 and DSCAM) and sensitivity of olfactory sensory neurons to food cues (e.g., NPY). A heatmap displaying the gene expression patterns of 207 genes related to behavior (Supplementary Table S7) and 53 genes involved in olfaction (Supplementary Table S8) that show differential expression in the olfactory bulb of pregnant vs. non-pregnant goats is shown in Supplementary Figure S22.

4 Discussion

4.1 The clustering of samples from certain brain regions follows an anterior-posterior pattern consistent with embryonic origin

The majority of the 17,054 protein-coding genes expressed in the goat brain had τ-values below 0.25 (55.57% of the genes), implying that they have a ubiquitous pattern of expression. Indeed, only 469 genes (2.75%) were expressed in just one of the 12 sampled brain regions (τ = 1). Results obtained in the Genotype-Tissue Expression (GTEx) project were consistent with this finding (Melé et al., 2015), since only 200 genes showed tissue-specific expression (95% of these genes were exclusively expressed in the testis).

Brain region of origin, rather than pregnancy, was the major factor explaining the clustering of samples. By using microarrays, Ferraz et al. (2008) investigated the profiles of expression of 16 porcine tissues and observed that the factor “tissue of origin” accounted for ∼11 times more variability than sex or breed. In humans, the GTEx Consortium retrieved 1641 post-mortem samples covering 54 body sites from 175 individuals and reported that tissue type was the primary factor explaining differences in gene expression (Ardlie et al., 2015), a result that is fully consistent with ours. Similar findings have been obtained when building atlases of gene expression in cattle (Harhay et al., 2010), pigs (Freeman et al., 2012) and sheep (Clark et al., 2017). We have observed that a few samples do not cluster with the remaining samples from the same region of origin (see for instance the two hippocampus samples in the upper left quadrant of the PCA shown in Figure 1B). This lack of consistent clustering has been observed in other atlases of gene expression (Harhay et al., 2010; Ardlie et al., 2015; Clark et al., 2017). Rather than mix-ups or contaminations, unexpected patterns of clustering could be due to the expression of thousands of housekeeping genes, that are transcribed at similar rates in different tissues, and, additionally, to the existence of “communication pathways” between tissues which might contribute to generate a high correlation of expression profiles between them, as evidenced in previous studies (Wang et al., 2019).

We have observed that, in goats, frontal neocortex, olfactory bulb and hippocampus samples tend to group together, and the same trend is observed for rostral colliculus, medulla oblongata and pons samples, while hypothalamus samples would be placed in the middle of these two clusters. In strong contrast, cerebellum, pineal gland and adeno/neurohypophysis show highly differentiated patterns of gene expression (Figure 1A; Figure 2). It could be argued that, with few exceptions, brain samples cluster together because of anatomical proximity. However, it should be noticed that the olfactory bulb is located closer to hypothalamus than to the hippocampus, or that the distance between the rostral colliculus and the hippocampus is shorter than that between the olfactory bulb and the hippocampus. The common denominator between frontal neocortex, olfactory bulb and hippocampus is that these three organs derive from the telencephalon, while pons and medulla oblongata come from the rhombencephalon. In the comprehensive analysis of porcine tissue expression carried out by Ferraz et al. (2008), tissues clustered according to the germ layer (ectoderm, mesoderm, or endoderm) they derive from. Measurement of transcriptome profiles in 24 murine neural tissues also highlighted the existence of a relationship between the cellular position along the anterior-posterior axis of the neural tube and gene expression in different regions from the adult brain (Zapala et al., 2005), as also published by Ortiz et al. (2020). These published results and those obtained by us agree with the hypothesis that embryogenesis might leave a durable footprint in the profile of mRNA expression of mammalian brain regions.

4.2 Evolutionary conservation of the genome-wide patterns of expression across brain regions

Sjöstedt et al. (2020) described the patterns of genome-wide expression of protein-coding genes in a number of brain regions from humans, mice and pigs. They found that, in general, transcriptomic patterns are evolutionarily conserved in these three species, with the three forebrain structures (cerebral cortex, hippocampus, and amygdala) grouping together, and the midbrain, thalamus, and pons and medulla samples forming another cluster close to the hypothalamus (Sjöstedt et al., 2020). They also found that the cerebellum behaves as an outlier, a finding fully coherent with ours.

Although transcriptomic profiles of brain regions seem to be well conserved across species, relevant differences also exist. For instance, Sjöstedt et al. (2020) indicated that in humans the olfactory bulb clusters with other forebrain structures (cerebral cortex, hippocampus, amygdala, etc.), while in pigs and mice this brain structure almost behaves as an outlier. They reasoned that this might be due to the fact that these two latter species have more evolved olfactory systems than humans (Sjöstedt et al., 2020). In goats, the olfactory bulb is also much more developed than in humans (Kavoi and Jameela, 2011), a feature consistent with the fundamental role of the olfactory system in the maternal and social behaviors of ungulates (Keller and Lévy, 2012). However, our data indicate that in goats the pattern of mRNA expression of the olfactory bulb is closely aligned with that of other forebrain regions (Figure 1B; Figure 2). This result agrees much more with the PCA reported by Sjöstedt et al. (2020) for human brain structures than for those corresponding to pigs or mice. A higher resemblance of goats to humans, rather than to pigs, is unexpected because goats and pigs are ungulates (Benton et al., 2015). Sjöstedt et al. (2020) stated that the functional importance of olfaction might be a key factor explaining the differential patterns of olfactory bulb mRNA expression observed in pigs/mice vs. humans. Such hypothesis is not fully consistent with our results and should be interpreted cautiously until further data are available in other mammalian species.

4.3 Cerebellum, pineal gland and hypophysis have highly differentiated patterns of gene expression

Cerebellum displayed a highly differentiated pattern of mRNA expression when compared to other brain structures (Figure 1A; Figure 2), a finding that, as mentioned before, is consistent with previous reports (Ortiz et al., 2020; Sjöstedt et al., 2020). Moreover, cerebellum trunk and hemisphere mRNA profiles were quite similar. Cerebellum is strongly specialized in the learning and coordination of motor activities as well as in the triggering of reflex responses, and it has been suggested that it is scarcely connected with cognitive areas of the brain cortex (Glickstein and Doron, 2008). One of the most distinctive features of the cerebellum is its extraordinarily high cell density: while this organ represents 10% of brain volume, it encompasses 42% of all brain cells and 59.8% of all excitatory neurons, mainly due to the conspicuous abundance of tightly packed granular cells (Erö et al., 2018). Another relevant feature of the cerebellum is that the same circuit, composed by mossy fibers that excite granule cells that, in turn, excite Purkinje cells, forms a fundamental unit that is replicated thousands of times (Kozareva et al., 2020). These biological particularities might have contributed to transform the cerebellum into an organ with a very specific transcriptomic profile, not only in goats but also in other mammals.

We also observed highly divergent patterns of gene expression in three glandular tissues: the pineal gland and the adeno/neurohypophysis (Figure 1A), a finding consistent with results reported by Harhay et al. (2010). The strong functional specialization of these three anatomical structures in hormonal secretion might explain these findings. In response to light, the pineal gland, which develops from an evagination of neuroepithelium in the dorsal midline of the diencephalon, synthesizes and releases melatonin, which is a key regulator of the circadian sleep-wake cycle and seasonal rhythms (Patel et al., 2020). Remarkably, in our study the pineal gland was the second organ showing the highest number of genes with region-specific expression (τ > 0.85, Figure 3B and Supplementary Table S3A). Neurohypophysis and adenohypophysis are also specialized in the secretion of molecules with key physiological roles, and data collected in sheep and cattle indicate that their profiles of mRNA expression are highly differentiated from those of other brain structures (Harhay et al., 2010; Clark et al., 2017). The main hormones produced by the adenohypophysis are prolactin, adrenocorticotropic hormone, luteinizing hormone, follicle-stimulating hormone, growth hormone and thyroid-stimulating hormone, which regulate a very diverse set of biological processes including growth, metabolism, lactation, stress and reproduction (Le Tissier et al., 2012). One particular feature of the adenohypophysis is that it does not develop from any of the five neural vesicles but from the oral ectoderm (Le Tissier et al., 2012). In contrast, the neurohypophysis has a diencephalic origin and stores vasopressin and oxytocin (both are synthesized in the hypothalamus), two hormones regulating diuresis and a broad array of reproduction and behavioral processes, respectively (Santos Fontanez and De Jesus, 2021). Despite having completely different embryological origins, histological structure and biological functions, our data indicate that neurohypophysis and adenohypophysis share similar profiles of mRNA expression (Figure 1A; Figure 2). Our interpretation is that high functional specialization in hormonal secretion might erase, at least partially, the transcriptomic footprint associated with embryogenesis, a hypothesis supported by the distinctive patterns of gene expression observed in both pituitary structures when compared with the remainder encephalic regions. Further studies are needed to assess whether such interpretation is correct.

4.4 Heterogeneous effects of pregnancy on the expression profiles of eleven goat brain regions

Amongst the seven individuals sampled in our study, there were three pregnant goats providing the opportunity to investigate the effect of pregnancy on the brain transcriptome. Since the number of replicates for each pregnant (N = 3) and non-pregnant (N = 4) category is low, sensitivity to detect DEGs is expected to be quite limited (Schurch et al., 2016). To circumvent this difficulty, we have used the total number of DEGs as the main criterion to identify which brain regions are mostly affected by 1-month pregnancy. In principle, low number of replicates should affect the number of detected DEGs to a similar extent in all 11 brain regions (pineal gland was excluded from the differential expression analysis), so the number of DEGs seems to be an appropriate indicator of which regions are more affected by pregnancy.

In comparisons involving cerebellar hemisphere, cerebellar trunk, hypothalamus, medulla oblongata, neurohypophysis and rostral colliculus, the total number of DEGs when comparing pregnant vs. non-pregnant goats was very low or inexistent (Figure 4B). In contrast, between 70 and 201 DEGs were detected in the frontal neocortex, hippocampus, adenohypophysis and pons (Figure 4A). By far, the organ which displayed the largest number of DEGs was the olfactory bulb (826 upregulated and 381 downregulated genes in pregnant goats). From these data, we conclude that 1-month pregnancy does not have the same effect on all goat brain regions. Ray et al. (2016) investigated changes in the mRNA expression of four brain structures (hypothalamus, neocortex, hippocampus and cerebellum) in virgin, pregnant, and postpartum mice, and they found that in the virgin vs. pregnant comparison the number of DEGs was much higher in the hippocampus than in the cerebellum, while hypothalamus and neocortex showed intermediate values. These results support the notion that pregnancy does not affect all brain regions to the same extent.

Brain gene expression during pregnancy is dynamic and it might change depending on the time point under consideration. By using magnetic resonance imaging, it has been shown that the brain of primiparous women experiences substantial morphological changes during gestation which mostly affect the right middle temporal gyrus, inferior frontal gyrus and posterior cingulate cortex (Hoekzema et al., 2017). A longitudinal morphometric study in mice also provided evidence of transient hypertrophy associated with gestation and/or lactation in the medial preoptic area, bed nucleus of the stria terminalis, amygdala, caudate nucleus, and hippocampus (Barrière et al., 2021). These findings indicate that certain brain structures undergo significant alterations in their profiles of expression and morphology because of gestation, and that biological changes are highly dynamic. Thus, the lists of DEGs detected for 11 goat brain regions at 1-month of gestation might not be representative of the whole gestation, but just of this specific time point.

In the goat frontal neocortex and the hippocampus, 82 DEGs (45 upregulated and 37 downregulated in pregnant goats) and 70 DEGs (55 upregulated and 15 downregulated in pregnant goats) were detected, respectively (Supplementary Tables S5D and S5E). The hippocampus and the frontal neocortex play important roles in modulating memory and learning as well as social behavior (Zemla and Basu, 2017; Reinert et al., 2021), and there is evidence that hippocampal neurogenesis is affected by pregnancy (Rolls et al., 2008). In this way, genes included in pathways involved in the immune response, including NFATC3 and CCL2, and apoptosis (CTSH) were upregulated in these tissues in pregnant goats (Supplementary Table S6C). Interestingly, these genes have been associated with neurodegenerative processes in the hippocampus and memory loss (Abdul et al., 2009; Wang et al., 2010). Besides, hippocampal lesion in the female rat results in deficient nest construction and reduced pup survival, suggesting that this organ has a key influence on the development of maternal behavior (Murphy and Wildeman, 1986).

In the adenohypophysis 201 DEGs (188 upregulated and 13 downregulated in pregnant goats) were identified (Supplementary Table S5A). Adenohypophysis provides hormonal signals that are fundamental for the maintenance of gestation. As listed in Supplementary Table S6A, genes included in hormonal pathways like oxytocin signaling pathway (i.e., RYR2, ROCK2, GNAQ, NFATC3 and ADCY7) and cortisol synthesis and secretion (i.e., NCEH1, GNAQ, CREB3L2 and ADCY7) were upregulated in pregnant goats. The anterior pituitary is strongly enlarged during gestation in response to hormones produced by the placenta and ovaries, and the secretion of prolactin increases gradually to prepare the female for lactation (Carlson, 2011).

Differential mRNA expression in the pons, with 190 DEGs (133 upregulated and 57 downregulated in pregnant goats) is harder to interpret (Supplementary Table S5J). This brain region is mostly involved in the regulation of breathing and sleep, and in relaying information to or from the cerebellum to other brain regions (Venkatraman et al., 2017). However, there is evidence that the pons has a basic role in the generation and experience of emotions through the integration of arousal, autonomic function, motor control, and somatosensory signals (Venkatraman et al., 2017). Moreover, it has been reported that preoptic area projections to lower brainstem regions affect maternal behavior in postpartum rats (Numan and Numan, 1991).

4.5 Several genes are differentially expressed in two or more brain regions in pregnant vs. non-pregnant goats

We have detected 94 genes that are differentially expressed in more than one tissue when comparing pregnant and non-pregnant goats (90 genes are differentially expressed in two tissues and four genes are differentially expressed in three tissues). The most significant pathways detected in this list of 94 genes are Wnt signalling pathway, adherens junctions signalling pathways, and FoxO signalling pathway (Supplementary Table S9). These pathways are related to cellular growth and proliferation and they are involved in the formation and modulation of neuronal circuits, neurogenesis, neuronal migration and axon guidance (Inestrosa and Varela-Nallar, 2015; Renault et al., 2009; Stocker and Chenn, 2015).

Amongst this list of 94 loci, one gene that is particularly interesting is that encoding the dopa decarboxylase (DDC) mRNA, which was consistently downregulated in the pons and olfactory bulb (log₂FC -2.08 and -1.3, respectively). This gene is responsible for the synthesis of dopamine, essential for a healthy pregnancy (Gratz et al., 2018). Dopamine is a key modulator of maternal behavior in rodents, partly because it inhibits the pituitary release of prolactin which modulates lactation, corpora lutea function (in certain species), and maternal behavior (Nevard et al., 2022).

Notably, the direction of differential gene expression was not always consistent across tissues. For instance, the vesicular associated membrane protein 1 (VAMP1) gene was significantly upregulated in the olfactory bulb of pregnant goats (log₂FC = 1.15), while in the frontal neocortex was downregulated in pregnant goats (log₂FC = -0.67). The VAMP1 gene encodes a protein that plays an important role in synapsis and synapsis plasticity and forms part of the soluble N-ethylmaleimide-sensitive factor attachment protein receptor complex (Dalmaz et al., 2021). These proteins are involved in neurocognitive processes and associated with the Wnt pathways that, among other things, are important for neurogenesis and cell proliferation (Dalmaz et al., 2021).

4.6 Dramatic changes in the mRNA expression of the olfactory bulb in response to early gestation

The olfactory bulb was the brain region that displayed, by far, the highest number of DEGs (Supplementary Table S5I), implying that this organ is strongly affected by 1 month-pregnancy. The olfactory bulb processes smell information transmitted from olfactory sensory neurons, which express odorant receptors in the olfactory epithelium, and relays it to the olfactory cortex (Nagayama et al., 2014). Importantly, the sense of olfaction is strongly connected to a broad array of behavioral responses related with aversion to food, avoidance of predators, sex arousal, partner preference, and aggression, amongst others (Sullivan et al., 2015). For instance, bilateral bulbectomy of male rats is associated with indifference to the sexual status of females and suppression of mating (Edwards et al., 1990), and it also decreases aggressive behavior between males (Rowe and Edwards, 1971).

Of course, the strong changes in the mRNA expression profile of the olfactory bulb observed in pregnant goats, when compared to their non-pregnant counterparts, are not triggered by odorant signals delivered by the offspring. More likely, this remarkable change of gene expression might be promoted by other brain regions and/or endocrine glands delivering chemical signals to the olfactory bulb. As expected, several DEGs are related with the sense of olfaction (Supplementary Table S8). For instance, the inactivation of the zinc finger homeodomain factor teashirt zinc finger family member 1 (TSHZ1, |log₂FC| 1.36, q-value < 0.05) or of the prokineticin receptor 2 (PROKR2, |log₂FC| 1.28, q-value < 0.05) genes results in olfactory bulb hypoplasia (Prosser et al., 2007; Ragancokova et al., 2014); and the loss of the Spalt like transcription factor 3 (SALL3, |log₂FC| 1.02, q-value < 0.05) inhibits the formation of the glomerular layer in the olfactory bulb (Harrison et al., 2008). The aristaless related homeobox (ARX, |log₂FC| 0.97, q-value < 0.05) gene is fundamental for the proper development of mouse olfactory system (Yoshihara et al., 2005), and loss of sonic hedgehog signaling molecule (SHH, |log₂FC| 1.02, q-value < 0.05) or olfactomedin 2 (OLFM2, |log₂FC| 1.09, q-value < 0.05) impairs olfactory function (Balmer and Lamantia, 2004; Sultana et al., 2014). Moreover, the formation of the olfactory bulb requires fibroblast growth factor (FGF) signaling mediated by the FGF receptor one FGFR1 (|log₂FC| 0.59, q-value < 0.05) (Hébert et al., 2003), and anosmin 1 (ANOS1, |log₂FC| 1.68, q-value < 0.05), which is essential for blood vessel formation in the olfactory bulb (Matsushima et al., 2020). So, modifications in the expression of these and other genes might result in changes of smell perception and/or of behavioral responses associated with it as a transition to maternity.

More enigmatic are changes in the mRNA expression of a large number of genes known to influence behavioral traits in humans (Supplementary Table S7). For instance, we have observed an upregulation, in the olfactory bulb of pregnant goats, of the dopamine β-hydroxylase gene (DBH, |log₂FC| 0.93, q-value < 0.05), which encodes an oxidoreductase that catalyzes the conversion of dopamine to norepinephrine. In mice, the inactivation of the DBH gene causes a profound impairment of maternal behavior characterized by the abandonment of pups which die a few days after birth (Thomas and Palmiter, 1997). The dopamine D1 receptor (DRD1, |log₂FC| 1.91, q-value < 0.05) and the 5-hydroxytryptamine receptor 2A (HTR2A, |log₂FC| 0.91, q-value < 0.05) show an increased expression in the olfactory bulb of pregnant goats. The DRD1 gene is involved in the regulation of mothering in rats (Mileva-Seitz et al., 2012), while the HTR2A molecule binds to serotonin, a neurotransmitter upregulated during the transition to motherhood with strong effects on maternal care and pup survivability (Angoa-Pérez et al., 2014; Pawluski et al., 2019). However, it is worth emphasizing that interpreting the observed patterns of mRNA expression from a functional perspective is not completely straightforward. For instance, in the olfactory bulb of pregnant goats we have observed an upregulation of the histamine receptor H1 (HRH1, |log₂FC| 1.43, q-value < 0.05) and the prion protein (PRNP, |log₂FC| 0.60, q-value < 0.05) mRNAs. The inactivation of these genes leads to increased aggression (Tomaz et al., 2014; Naganuma et al., 2017), while in parallel we have detected an upregulation of histamine receptor H3 (HRH3, |log₂FC| 3, q-value < 0.05) mRNA, which in zebrafish has the opposite effect (Reichmann et al., 2020). Moreover, we have found that the neurotensin (NTS, |log₂FC| 0.86, q-value < 0.05) mRNA is upregulated in the olfactory bulb of pregnant females. Interestingly, the intracerebroventricular injection of neurotensin in mice decreases aggression associated with the maternal protection of the offspring (Gammie et al., 2009).

With regard to olfactory bulb DEGs linked to depression and anxiety, we have found an increased mRNA expression of both antidepressant (Son et al., 2012; Zhou and Zhou, 2014), e.g., lipid droplet biogenesis associated seipin (BSCL2, |log₂FC| 0.64, q-value < 0.05), neuritin (NRN1, |log₂FC| 1.21, q-value < 0.05) and glycine receptor β (GLRB, |log₂FC| 0.75, q-value < 0.05), and pro-depressive or pro-anxiety genes (Knoll and Carlezon, 2010; Su et al., 2018; Liu et al., 2021), e.g., salt inducible kinase 2 (SIK2, |log₂FC| 1.44, q-value < 0.05), prodynorphin (PDYN, |log₂FC| 0.95, q-value < 0.05) and histamine receptor H1 (HRH1, |log₂FC| 1.41, q-value < 0.05) in pregnant goats. Similarly, we have found increased mRNA levels of genes which promote feeding, e.g., neuropeptide Y, NPY(|log₂FC| 1.58, q-value < 0.05); C-X-C motif chemokine ligand 14, CXCL14 (|log₂FC| 0.99, q-value < 0.05); and neural EGFL like 2, NELL2 (|log₂FC| 3.13, q-value < 0.05); or inhibit it, e.g., G protein-coupled receptor 45, GPR45 (|log₂FC| 1.56, q-value < 0.05); and acyl-CoA binding domain containing 7, ACBD7 (|log₂FC| 0.63, q-value < 0.05), as reported in previous studies (Tanegashima et al., 2010; Cui et al., 2016; Jeong et al., 2017; Lanfray and Richard, 2017; Engström Ruud et al., 2020). These findings illustrate the complexity of ascertaining the functional implications of changes in olfactory bulb mRNA expression associated with pregnancy.

There is also a subset of genes with increased expression in the olfactory bulb of pregnant goats which is particularly interesting because they display dual functions related with both olfaction and behavior. Examples of this category are genes encoding hippocalcin (HPCA, |log₂FC| 1.60, q-value < 0.05), a mediator of second messenger signaling in the olfactory epithelium which also regulates spatial and associative memory (Mamman et al., 2004; Kobayashi et al., 2005), and the brain derived neurotrophic factor (BDNF, |log₂FC| 1.78, q-value < 0.05), which is important for olfactory bulb neurogenesis, learning and memory (Yuan, 2008; Miranda et al., 2019).

The results of the differential expression analysis should be interpreted with caution because sample size is small. Despite this caveat, we consider that the disproportionate number of DEGs detected in the olfactory bulb when compared to other brain regions is a robust result. Importantly, this anatomical structure is strongly involved in the development of maternal behavior in mice, i.e., its complete removal eliminates such capacity and, very often, implies that litters are cannibalized by their mothers soon after parturition (Gandelman et al., 1971). Moreover, intrabulbar infusions of an oxytocin antagonist delays maternal behavior in rats, while oxytocin itself has the opposite effect (Yu et al., 1996). In sheep, which are phylogenetically close to goats, disruption of the noradrenergic projections to the olfactory bulb suppresses the establishment of a maternal bond between mother and offspring, probably because such event strongly depends on odor cues (Pissonnier et al., 1985). These results support a critical role of the olfactory bulb in the induction of maternal behavior in mammals, although discrepant findings have also been published (Pissonnier et al., 1985).

It has been proposed that in sheep and goats the increase of estrogen levels in blood during prepartum, combined with the vaginocervical stimulation provoked by the expulsion of the fetus, trigger modifications in a neural network including the main olfactory system, the medial preoptic area and the paraventricular nucleus of the hypothalamus, in order to induce maternal behavior (Poindron et al., 2007). Our transcriptomic data demonstrate that extensive areas of the goat brain manifest changes in their expression profiles much before the peripartum period. These alterations might represent an adaptive response toward preparing the pregnant goat for maternity, and we hypothesize that the olfactory bulb might hold a very relevant role in such process.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA808876.

Ethics statement

Ethical review and approval was not required for the animal study because goats were slaughtered as a result of a routine culling procedure managed by the Service of Farms and Experimental Fields of the Universitat Autònoma de Barcelona (UAB). Such procedure was completely unrelated with our project, so permission from the Ethics Committee on Animal and Human Experimentation of the UAB was not required.

Author contributions

Conceptualization—MA and ML-B; formal analysis—ML-S; funding acquisition—MA; investigation—ML-S and DG; project administration—MA; resources—ML-B, MA, ML-S, DG, EC, SO-M, JG, MP, UIO, ULO, and EM-S; visualization—ML-S; writing—original draft—MA and ML-S; writing—review & editing—EM-S, MP, DG, ML-B, MA, ML-S, EC, SO-M, JG, UIO, ULO. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by the European Regional Development Fund (FEDER)/Ministerio de Ciencia e Innovación—Agencia Estatal de Investigación/Project Reference grants: AGL2016-76108-R and also by project PID2019-105805RB-I00 funded by MCIN/AEI/10.13039/501100011033. We also acknowledge the support of the CERCA Programme/Generalitat de Catalunya and of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo Ochoa 2020-2023 (CEX2019-000902-S) grant awarded to the Centre for Research in Agricultural Genomics (CRAG, Bellaterra, Spain). Maria Luigi-Sierra was funded with a PhD fellowship “Formación de Personal Investigador” (BES-2017-079709) awarded by the Spanish Ministry of Economy and Competitivity and by “ESF Investing in your future”. Dailu Guan was funded by a PhD fellowship from the China Scholarship Council (CSC). Emilio Mármol-Sánchez was funded with a FPU Ph.D. grant from the Spanish Ministry of Education (FPU15/01733).

Acknowledgments

Thanks to José Luis Ruiz de la Torre Casañas from the Service of Farms and Experimental Fields of the UAB and Dr Xavier Such from the Department of Animal Sciences of the UAB for providing the animal material used in this experiment.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1114749/full#supplementary-material

Abbreviations

C, combined score; cDNA, complementary desoxyribonucleic acid; CPM, counts per million; DEGs, differentially expressed genes; FC, fold change; FDR, false discovery rate; GTEx, genotype-tissue expression; KEGG, Kyoto encyclopedia of genes and genomes; mRNA, messenger ribonucleic acid; PC, principal component; PCA, principal component analysis; QC, quality control; RIN, RNA integrity number; RNA, ribonucleic acid.

References

Abdul, H. M., Sama, M. A., Furman, J. L., Mathis, D. M., Beckett, T. L., Weidner, A. M., et al. (2009). Cognitive decline in Alzheimer’s disease is associated with selective changes in calcineurin/NFAT signaling. J. Neurosci. 29, 12957–12969. doi:10.1523/jneurosci.1064-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Angoa-Pérez, M., Kane, M. J., Sykes, C. E., Perrine, S. A., Church, M. W., and Kuhn, D. M. (2014). Brain serotonin determines maternal behavior and offspring survival. Genes. Brain Behav. 13, 579–591. doi:10.1111/gbb.12159

PubMed Abstract | CrossRef Full Text | Google Scholar

Arbeitman, M. N. (2019). Maternal experience leads to lasting gene expression changes in some regions of the mouse brain. G3 Genes., Genomes, Genet. 9, 2623–2628. doi:10.1534/g3.119.400249

CrossRef Full Text | Google Scholar

Ardlie, K. G., DeLuca, D. S., Segrè, A. v., Sullivan, T. J., Young, T. R., Gelfand, E. T., et al. (2015). Human genomics. The genotype-tissue expression (GTEx) pilot analysis: Multitissue gene regulation in humans. Science 348, 648–660. doi:10.1126/science.1262110

PubMed Abstract | CrossRef Full Text | Google Scholar

Balmer, C. W., and Lamantia, A. S. (2004). Loss of Gli3 and Shh function disrupts olfactory axon trajectories. J. Comp. Neurology 472, 292–307. doi:10.1002/cne.20053

CrossRef Full Text | Google Scholar

Barrière, D. A., Ella, A., Szeremeta, F., Adriaensen, H., Même, W., Chaillou, E., et al. (2021). Brain orchestration of pregnancy and maternal behavior in mice: A longitudinal morphometric study. Neuroimage 230, 117776–117815. doi:10.1016/j.neuroimage.2021.117776

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate - a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Benton, M. J., Donoghue, P. C. J., Asher, R. J., Friedman, M., Near, T. J., and Vinther, J. (2015). Constraints on the timescale of animal evolutionary history. Palaeontol. Electron. 18, 1–106. doi:10.26879/424

CrossRef Full Text | Google Scholar

Bickhart, D. M., Rosen, B. D., Koren, S., Sayre, B. L., Hastie, A. R., Chan, S., et al. (2017). Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome. Nat. Genet. 49, 643–650. doi:10.1038/ng.3802

PubMed Abstract | CrossRef Full Text | Google Scholar

Camargo, A. P., Vasconcelos, A. A., Fiamenghi, M. B., Pereira, G. A. G., and Carazzolle, M. F. (2020). tspex: a tissue-specificity calculator for gene expression data. Res. Sq., 1–7. doi:10.21203/rs.3.rs-51998/v1

CrossRef Full Text | Google Scholar

Carlson, H. E. (2011). “The pituitary gland in pregnancy and the puerperium,” in The pituitary. Editor S. Melmed (Cambridge, MA: Academic Press), 397–412.

CrossRef Full Text | Google Scholar

Chen, E. Y., Tan, C. M., Kou, Y., Duan, Q., Wang, Z., Meirelles, G. V., et al. (2013). Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinforma. 14, 128. doi:10.1186/1471-2105-14-128

CrossRef Full Text | Google Scholar

Chhetri, K. P., and Das, J. M. (2021). Neuroanatomy, neural tube development and stages. Treasure Island (FL): StatPearls Publishing.

Google Scholar

Clark, E. L., Bush, S. J., McCulloch, M. E. B., Farquhar, I. L., Young, R., Lefevre, L., et al. (2017). A high resolution atlas of gene expression in the domestic sheep (Ovis aries). PLoS Genet. 13, 1006997. doi:10.1371/journal.pgen.1006997

CrossRef Full Text | Google Scholar

Cui, J., Ding, Y., Chen, S., Zhu, X., Wu, Y., Zhang, M., et al. (2016). Disruption of Gpr45 causes reduced hypothalamic POMC expression and obesity. J. Clin. Investigation 126, 3192–3206. doi:10.1172/JCI85676

CrossRef Full Text | Google Scholar

Dalmaz, C., Barth, B., Pokhvisneva, I., Wang, Z., Patel, S., Quillfeldt, J. A., et al. (2021). Prefrontal cortex VAMP1 gene network moderates the effect of the early environment on cognitive flexibility in children. Neurobiol. Learn Mem. 185, 107509. doi:10.1016/j.nlm.2021.107509

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, D. A., Griffis, K. T., and Tardivel, C. (1990). Olfactory bulb removal: Effects on sexual behavior and partner-preference in male rats. Physiol. Behav. 48, 447–450. doi:10.1016/0031-9384(90)90342-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Engström Ruud, L., Pereira, M. M. A., de Solis, A. J., Fenselau, H., and Brüning, J. C. (2020). NPY mediates the rapid feeding and glucose metabolism regulatory functions of AgRP neurons. Nat. Commun. 11, 442. doi:10.1038/s41467-020-14291-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Erö, C., Gewaltig, M. O., Keller, D., and Markram, H. (2018). A cell atlas for the mouse brain. Front. Neuroinform 12, 84. doi:10.3389/fninf.2018.00084

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferraz, A. L. J., Ojeda, A., López-Béjar, M., Fernandes, L. T., Castelló, A., Folch, J. M., et al. (2008). Transcriptome architecture across tissues in the pig. BMC Genomics 9, 173. doi:10.1186/1471-2164-9-173

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, A. O., Vasconcelos, B. G., Favaron, P. O., Santos, A. C., Leandro, R. M., Pereira, F. T. V., et al. (2018). Bovine central nervous system development. Pesqui. Veterinaria Bras. 38, 147–153. doi:10.1590/1678-5150-pvb-5020

CrossRef Full Text | Google Scholar

Freeman, T. C., Ivens, A., Baillie, J. K., Beraldi, D., Barnett, M. W., Dorward, D., et al. (2012). A gene expression atlas of the domestic pig. BMC Biol. 10, 90. doi:10.1186/1741-7007-10-90

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallego Romero, I., Pai, A. A., Tung, J., and Gilad, Y. (2014). RNA-seq: Impact of RNA degradation on transcript quantification. BMC Biol. 12, 42. doi:10.1186/1741-7007-12-42

PubMed Abstract | CrossRef Full Text | Google Scholar

Gammie, S. C., D’Anna, K. L., Gerstein, H., and Stevenson, S. A. (2009). Neurotensin inversely modulates maternal aggression. Neuroscience 158, 1215–1223. doi:10.1016/j.neuroscience.2008.11.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Gandelman, R., Zarrow, M. X., Denenberg, V. H., and Myers, M. (1971). Olfactory bulb removal eliminates maternal behavior in the mouse. Science 171, 210–211. doi:10.1126/science.171.3967.210

PubMed Abstract | CrossRef Full Text | Google Scholar

Glickstein, M., and Doron, K. (2008). Cerebellum: Connections and functions. Cerebellum 7, 589–594. doi:10.1007/s12311-008-0074-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Gratz, M. J., Stavrou, S., Kuhn, C., Hofmann, S., Hermelink, K., Heidegger, H., et al. (2018). Dopamine synthesis and dopamine receptor expression are disturbed in recurrent miscarriages. Endocr. Connect. 7, 727–738. doi:10.1530/EC-18-0126

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, D., Landi, V., Luigi-Sierra, M. G., Delgado, J. V., Such, X., Castelló, A., et al. (2020). Analyzing the genomic and transcriptomic architecture of milk traits in Murciano-Granadina goats. J. Anim. Sci. Biotechnol. 11, 35. doi:10.1186/s40104-020-00435-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, S. J., Parrish, M., and Monaghan, A. P. (2008). Sall3 is required for the terminal maturation of olfactory glomerular interneurons. J. Comp. Neurology 507, 1780–1794. doi:10.1002/cne.21650

CrossRef Full Text | Google Scholar

Harhay, G. P., Smith, T. P. L., Alexander, L. J., Haudenschild, C. D., Keele, J. W., Matukumalli, L. K., et al. (2010). An atlas of bovine gene expression reveals novel distinctive tissue characteristics and evidence for improving genome annotation. Genome Biol. 11, R102. doi:10.1186/gb-2010-11-10-r102

PubMed Abstract | CrossRef Full Text | Google Scholar

Hébert, J. M., Lin, M., Partanen, J., Rossant, J., and McConnell, S. K. (2003). FGF signaling through FGFR1 is required for olfactory bulb morphogenesis. Development 130, 1101–1111. doi:10.1242/dev.00334

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoekzema, E., Barba-Müller, E., Pozzobon, C., Picado, M., Lucco, F., García-García, D., et al. (2017). Pregnancy leads to long-lasting changes in human brain structure. Nat. Neurosci. 20, 287–296. doi:10.1038/nn.4458

PubMed Abstract | CrossRef Full Text | Google Scholar

Inestrosa, N. C., and Varela-Nallar, L. (2015). Wnt signalling in neuronal differentiation and development. Cell. Tissue Res. 359, 215–223. doi:10.1007/s00441-014-1996-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeong, J. K., Kim, J. G., Kim, H. R., Lee, T. H., Park, J. W., and Lee, B. J. (2017). A role of central NELL2 in the regulation of feeding behavior in rats. Mol. Cells 40, 186–194. doi:10.14348/molcells.2017.2278

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M. (2016). KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462. doi:10.1093/nar/gkv1070

PubMed Abstract | CrossRef Full Text | Google Scholar

Kavoi, B. M., and Jameela, H. (2011). Comparative morphometry of the olfactory bulb, tract and stria in the human, dog and goat. Int. J. Morphol. 29, 939–946. doi:10.4067/s0717-95022011000300047

CrossRef Full Text | Google Scholar

Keller, M., and Lévy, F. (2012). The main but not the accessory olfactory system is involved in the processing of socially relevant chemosignals in ungulates. Front. Neuroanat. 6, 39. doi:10.3389/fnana.2012.00039

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). Hisat: A fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Knoll, A. T., and Carlezon, W. A. (2010). Dynorphin, stress, and depression. Brain Res. 1314, 56–73. doi:10.1016/j.brainres.2009.09.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, M., Masaki, T., Hori, K., Masuo, Y., Miyamoto, M., Tsubokawa, H., et al. (2005). Hippocalcin-deficient mice display a defect in cAMP response element-binding protein activation associated with impaired spatial and associative memory. Neuroscience 133, 471–484. doi:10.1016/j.neuroscience.2005.02.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozareva, V., Martin, C., Osorno, T., Rudolph, S., Guo, C., Vanderburg, C., et al. (2020). A transcriptomic atlas of the mouse cerebellum reveals regional specializations and novel cell types. Nature 98, 214–219. doi:10.1101/2020.03.04.976407

CrossRef Full Text | Google Scholar

Lanfray, D., and Richard, D. (2017). Emerging signaling pathway in arcuate feeding-related neurons: Role of the Acbd7. Front. Neurosci. 11, 328. doi:10.3389/fnins.2017.00328

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Tissier, P. R., Hodson, D. J., Lafont, C., Fontanaud, P., Schaeffer, M., and Mollard, P. (2012). Anterior pituitary cell networks. Front. Neuroendocrinol. 33, 252–266. doi:10.1016/j.yfrne.2012.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Tang, W., Ji, C., Gu, J., Chen, Y., Huang, J., et al. (2021). The Selective SIK2 inhibitor ARN-3236 produces strong antidepressant-like efficacy in mice via the hippocampal CRTC1-CREB-BDNF pathway. Front. Pharmacol. 11, 624429. doi:10.3389/fphar.2020.624429

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mamman, A., Simpson, J. P., Nighorn, A., Imanishi, Y., Palczewski, K., Ronnett, G. V., et al. (2004). Hippocalcin in the olfactory epithelium: A mediator of second messenger signaling. Biochem. Biophys. Res. Commun. 322, 1131–1139. doi:10.1016/j.bbrc.2004.07.123

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsushima, S., Shimizu, A., Kondo, M., Asano, H., Ueno, N., Nakayama, H., et al. (2020). Anosmin-1 activates vascular endothelial growth factor receptor and its related signaling pathway for olfactory bulb angiogenesis. Sci. Rep. 10, 188. doi:10.1038/s41598-019-57040-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Melé, M., Ferreira, P. G., Reverter, F., DeLuca, D. S., Monlong, J., Sammeth, M., et al. (2015). Human genomics. The human transcriptome across tissues and individuals. Science 348, 660–665. doi:10.1126/science.aaa0355

PubMed Abstract | CrossRef Full Text | Google Scholar

Mileva-Seitz, V., Fleming, A. S., Meaney, M. J., Mastroianni, A., Sinnwell, J. P., Steiner, M., et al. (2012). Dopamine receptors D1 and D2 are related to observed maternal behavior. Genes. Brain Behav. 11, 684–694. doi:10.1111/j.1601-183X.2012.00804.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, M., Morici, J. F., Zanoni, M. B., and Bekinschtein, P. (2019). Brain-derived neurotrophic factor: A key molecule for memory in the healthy and the pathological brain. Front. Cell. Neurosci. 13, 363. doi:10.3389/fncel.2019.00363

PubMed Abstract | CrossRef Full Text | Google Scholar

Muriuki, C., Bush, S. J., Salavati, M., McCulloch, M. E. B., Lisowski, Z. M., Agaba, M., et al. (2019). A mini-atlas of gene expression for the domestic goat (Capra hircus). Front. Genet. 10, 1080. doi:10.3389/fgene.2019.01080

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, H. M., and Wildeman, C. H. (1986). Pup survival and development following hippocampal lesions in the female rat. Ohio J. Sci. 86, 89–94.

Google Scholar

Murtagh, F., and Legendre, P. (2014). Ward’s hierarchical agglomerative clustering method: Which algorithms implement ward’s criterion? J. Classif. 31, 274–295. doi:10.1007/s00357-014-9161-z

CrossRef Full Text | Google Scholar

Naganuma, F., Nakamura, T., Yoshikawa, T., Iida, T., Miura, Y., Kárpáti, A., et al. (2017). Histamine N-methyltransferase regulates aggression and the sleep-wake cycle. Sci. Rep. 7, 15899. doi:10.1038/s41598-017-16019-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagayama, S., Homma, R., and Imamura, F. (2014). Neuronal organization of olfactory bulb circuits. Front. Neural Circuits 8, 98. doi:10.3389/fncir.2014.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Nevard, R. P., Pant, S. D., Broster, J. C., Norman, S. T., and Stephen, C. P. (2022). Maternal behavior in beef cattle: The physiology, assessment and future directions-a review. Vet. Sci. 10, 10. doi:10.3390/vetsci10010010

PubMed Abstract | CrossRef Full Text | Google Scholar

Numan, M., and Numan, M. J. (1991). Preoptic-brainstem connections and maternal behavior in rats. Behav. Neurosci. 105, 1013–1029. doi:10.1037/0735-7044.105.6.1013

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortiz, C., Navarro, J. F., Jurek, A., Märtin, A., Lundeberg, J., and Meletis, K. (2020). Molecular atlas of the adult mouse brain. Sci. Adv. 6, abb3446. doi:10.1126/sciadv.abb3446

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, S., Rahmani, B., Gandhi, J., Seyam, O., Joshi, G., Reid, I., et al. (2020). Revisiting the pineal gland: A review of calcification, masses, precocious puberty, and melatonin functions. Int. J. Neurosci. 130, 464–475. doi:10.1080/00207454.2019.1692838

PubMed Abstract | CrossRef Full Text | Google Scholar

Pawluski, J. L., Li, M., and Lonstein, J. S. (2019). Serotonin and motherhood: from molecules to mood. Front. Neuroendocrinol. 53, 100742. doi:10.1016/j.yfrne.2019.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667. doi:10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Pissonnier, D., Thiery, J. C., Fabre-Nys, C., Poindron, P., and Keverne, E. B. (1985). The importance of olfactory bulb noradrenalin for maternal recognition in sheep. Physiol. Behav. 35, 361–363. doi:10.1016/0031-9384(85)90309-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Poindron, P., Lévy, F., and Keller, M. (2007). Maternal responsiveness and maternal selectivity in domestic sheep and goats: The two facets of maternal attachment. Dev. Psychobiol. 49, 54–70. doi:10.1002/dev.20192

PubMed Abstract | CrossRef Full Text | Google Scholar

Prosser, H. M., Bradley, A., and Caldwell, M. A. (2007). Olfactory bulb hypoplasia in Prokr2 null mice stems from defective neuronal progenitor migration and differentiation. Eur. J. Neurosci. 26, 3339–3344. doi:10.1111/j.1460-9568.2007.05958.x

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core team (2015). R: A language and environment for statistical computing. R foundation for statistical computing. Available at: http://www.R-project.org/.

Google Scholar

Ragancokova, D., Rocca, E., Oonk, A. M. M., Schulz, H., Rohde, E., Bednarsch, J., et al. (2014). TSHZ1-dependent gene regulation is essential for olfactory bulb development and olfaction. J. Clin. Investigation 132, 1214–1227. doi:10.1172/JCI72466

CrossRef Full Text | Google Scholar

Ray, S., Tzeng, R. Y., DiCarlo, L. M., Bundy, J. L., Vied, C., Tyson, G., et al. (2016). An examination of dynamic gene expression changes in the mouse brain during pregnancy and the postpartum period. G3 Genes., Genomes, Genet. 6, 221–233. doi:10.1534/g3.115.020982

CrossRef Full Text | Google Scholar

Reichmann, F., Rimmer, N., Tilley, C. A., Dalla Vecchia, E., Pinion, J., Al Oustah, A., et al. (2020). The zebrafish histamine H3 receptor modulates aggression, neural activity and forebrain functional connectivity. Acta Physiol. 230, e13543. doi:10.1111/apha.13543

CrossRef Full Text | Google Scholar

Reinert, S., Hübener, M., Bonhoeffer, T., and Goltstein, P. M. (2021). Mouse prefrontal cortex represents learned rules for categorization. Nature 593, 411–417. doi:10.1038/s41586-021-03452-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Renault, V. M., Rafalski, V. A., Morgan, A. A., Salih, D. A., Brett, J. O., Webb, A. E., et al. (2009). FoxO3 regulates neural stem cell homeostasis. Cell Stem Cell 5 (5), 527–539. doi:10.1016/j.stem.2009.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rolls, A., Schori, H., London, A., and Schwartz, M. (2008). Decrease in hippocampal neurogenesis during pregnancy: A link to immunity. Mol. Psychiatry 13, 468–469. doi:10.1038/sj.mp.4002126

PubMed Abstract | CrossRef Full Text | Google Scholar

Rowe, F. A., and Edwards, D. A. (1971). Olfactory bulb removal: Influences on the aggressive behaviors of male mice. Physiol. Behav. 7, 889–892. doi:10.1016/0031-9384(71)90059-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Santos Fontanez, S. E., and De Jesus, O. (2021). Neurohypophysis. Treasure Island, FL: StatPearls Publishing.

Google Scholar

Schurch, N. J., Schofield, P., Gierliński, M., Cole, C., Sherstnev, A., Singh, V., et al. (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22, 839–851. doi:10.1261/rna.053959.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Sjöstedt, E., Zhong, W., Fagerberg, L., Karlsson, M., Mitsios, N., Adori, C., et al. (2020). An atlas of the protein-coding genes in the human, pig, and mouse brain. Science 367, aay5947. doi:10.1126/science.aay5947

CrossRef Full Text | Google Scholar

Son, H., Banasr, M., Choi, M., Chae, S. Y., Licznerski, P., Lee, B., et al. (2012). Neuritin produces antidepressant actions and blocks the neuronal and behavioral deficits caused by chronic stress. Proc. Natl. Acad. Sci. U. S. A. 109, 11378–11383. doi:10.1073/pnas.1201191109

PubMed Abstract | CrossRef Full Text | Google Scholar

Stocker, A. M., and Chenn, A. (2015). The role of adherens junctions in the developing neocortex. Cell. Adh Migr. 9, 167–174. doi:10.1080/19336918.2015.1027478

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, W. J., Zhang, T., Jiang, C. L., and Wang, W. (2018). Clemastine alleviates depressive-like behavior through reversing the imbalance of microglia-related pro-inflammatory state in mouse hippocampus. Front. Cell. Neurosci. 12, 412. doi:10.3389/fncel.2018.00412

PubMed Abstract | CrossRef Full Text | Google Scholar

Sullivan, R. M., Wilson, D. A., Ravel, N., and Mouly, A. M. (2015). Olfactory memory networks: From emotional learning to social behaviors. Front. Behav. Neurosci. 9, 36. doi:10.3389/fnbeh.2015.00036

PubMed Abstract | CrossRef Full Text | Google Scholar

Sultana, A., Nakaya, N., Dong, L., Abu-Asab, M., Qian, H., and Tomarev, S. I. (2014). Deletion of olfactomedin 2 induces changes in the AMPA receptor complex and impairs visual, olfactory, and motor functions in mice. Exp. Neurol. 261, 802–811. doi:10.1016/j.expneurol.2014.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkin, S. M., Ng, L., Lau, C., Dolbeare, T., Gilbert, T. L., Thompson, C. L., et al. (2013). Allen brain atlas: An integrated spatio-temporal portal for exploring the central nervous system. Nucleic Acid. Res. 41, D996–D1008. doi:10.1093/nar/gks1042

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanegashima, K., Okamoto, S., Nakayama, Y., Taya, C., Shitara, H., Ishii, R., et al. (2010). CXCL14 deficiency in mice attenuates obesity and inhibits feeding behavior in a novel environment. PLoS One 5, e10321. doi:10.1371/journal.pone.0010321

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, S. A., and Palmiter, R. D. (1997). Impaired maternal behavior in mice lacking norepinephrine and epinephrine. Cell. 91, 583–592. doi:10.1016/S0092-8674(00)80446-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomaz, B., Aljaz, M., Ales, J., and Gregor, M. (2014). Deletion of the prion gene Prnp affects offensive aggression in mice. Behav. Brain Res. 266, 216–221. doi:10.1016/j.bbr.2014.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkatraman, A., Edlow, B. L., and Immordino-Yang, M. H. (2017). The brainstem in emotion: A review. Front. Neuroanat. 11, 15. doi:10.3389/fnana.2017.00015

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., Eraslan, B., Wieland, T., Hallström, B., Hopf, T., Zolg, D. P., et al. (2019). A deep proteome and transcriptome abundance atlas of 29 healthy human tissues. Mol. Syst. Biol. 15, e8503. doi:10.15252/MSB.20188503

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Michaelis, M. L., and Michaelis, E. K. (2010). Functional genomics of brain aging and alzheimer’s disease: Focus on selective neuronal vulnerability. Curr. Genomics 11, 618–633. doi:10.2174/138920210793360943

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanai, I., Benjamin, H., Shmoish, M., Chalifa-Caspi, V., Shklar, M., Ophir, R., et al. (2005). Genome-wide midrange transcription profiles reveal expression level relationships in human tissue specification. Bioinformatics 21 (5), 650–659. doi:10.1093/bioinformatics/bti042

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, S., Omichi, K., Yanazawa, M., Kitamura, K., and Yoshihara, Y. (2005). Arx homeobox gene is essential for development of mouse olfactory system. Development 132 (4), 751–762. doi:10.1242/dev.01619

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G. Z., Kaba, H., Okutani, F., Takahashi, S., and Higuchi, T. (1996). The olfactory bulb: a critical site of action for oxytocin in the induction of maternal behaviour in the rat. Neuroscience 72 (4), 1083–1088. doi:10.1016/0306-4522(95)00600-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, T. F. (2008). BDNF signaling during olfactory bulb neurogenesis. J. Neurosci. 28 (20), 5139–5140. doi:10.1523/JNEUROSCI.1327-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zapala, M. A., Hovatta, I., Ellison, J. A., Wodicka, L., Del Rio, J. A., Tennant, R., et al. (2005). Adult mouse brain gene expression patterns bear an embryologic imprint. Proc. Natl. Acad. Sci. U. S. A. 102 (29), 10357–10362. doi:10.1073/pnas.0503357102

PubMed Abstract | CrossRef Full Text | Google Scholar

Zemla, R., and Basu, J. (2017). Hippocampal function in rodents. Curr. Opin. Neurobiol. 43, 187–197. doi:10.1016/j.conb.2017.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, S., and Zhou, J. (2014). Neuritin, a neurotrophic factor in nervous system physiology. Curr. Med. Chem. 21 (10), 1212–1219. doi:10.2174/0929867321666131218093327

PubMed Abstract | CrossRef Full Text | Google Scholar