TYPE Original Research PUBLISHED 14 July 2023 DOI 10.3389/fgene.2023.1114749 A protein-coding gene expression OPEN ACCESS atlas from the brain of pregnant EDITED BY Monica Correa Ledur, Embrapa Suínos e Aves, Brazil and non-pregnant goats REVIEWED BY Jan Mulder, María Gracia Luigi-Sierra 1, Dailu Guan 1, Manel López-Béjar 2, Karolinska Institutet (KI), Sweden Stephen J. Bush, Encarna Casas 2, Sergi Olvera-Maneu 2, Jaume Gardela 2, University of Oxford, United Kingdom María Jesús Palomo 3, Uchebuchi Ike Osuagwuh 3, *CORRESPONDENCE Uchechi Linda Ohaneje 3, Emilio Mármol-Sánchez 1 and Marcel Amills, marcel.amills@uab.cat Marcel Amills 1,4* RECEIVED 02 December 2022 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Bellaterra, Spain, ACCEPTED 27 June 2023 2 Department of Animal Health and Anatomy, Universitat Autònoma de Barcelona, Bellaterra, Spain, PUBLISHED 14 July 2023 3 Department of Animal Medicine and Surgery, Universitat Autònoma de Barcelona, Bellaterra, Spain, 4 Departament de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, Bellaterra, Spain CITATION Luigi-Sierra MG, Guan D, López-Béjar M, Casas E, Olvera-Maneu S, Gardela J, Palomo MJ, Osuagwuh UI, Ohaneje UL, Mármol-Sánchez E and Amills M (2023), A protein-coding gene expression atlas Background: The brain is an extraordinarily complex organ with multiple from the brain of pregnant and non- anatomical structures involved in highly specialized functions related with pregnant goats. behavior and physiological homeostasis. Our goal was to build an atlas of Front. Genet. 14:1114749. doi: 10.3389/fgene.2023.1114749 protein-coding gene expression in the goat brain by sequencing the transcriptomes of 12 brain regions in seven female Murciano-Granadina goats, COPYRIGHT © 2023 Luigi-Sierra, Guan, López-Béjar, from which three of them were 1-month pregnant. Casas, Olvera-Maneu, Gardela, Palomo, Osuagwuh, Ohaneje, Mármol-Sánchez Results: Between 14,889 (cerebellar hemisphere) and 15,592 (pineal gland) and Amills. This is an open-access article protein-coding genes were expressed in goat brain regions, and most of them distributed under the terms of the displayed ubiquitous or broad patterns of expression across tissues. Principal Creative Commons Attribution License (CC BY). The use, distribution or component analysis and hierarchical clustering based on the patterns of mRNA reproduction in other forums is expression revealed that samples from certain brain regions tend to group permitted, provided the original author(s) according to their position in the anterior-posterior axis of the neural tube, and the copyright owner(s) are credited and that the original publication in this i.e., hindbrain (pons and medulla oblongata), midbrain (rostral colliculus) and journal is cited, in accordance with forebrain (frontal neocortex, olfactory bulb, hypothalamus, and hippocampus). accepted academic practice. No use, Exceptions to this observation were cerebellum and glandular tissues (pineal gland distribution or reproduction is permitted which does not comply with these terms. 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 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. Frontiers in Genetics 01 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 olfactory bulb, a finding consistent with the key role of this anatomical structure on the development of maternal behavior. KEYWORDS goat, RNA-Seq, gestation, embryonic vesicle, differential gene expression (DGE), encephalon 1 Introduction systematically. In goats, mRNA expression of two neural tissues, frontal lobe cortex and cerebellum, has also been characterized by The mammalian brain is an extraordinarily complex organ RNA-Seq (Muriuki et al., 2019). The main goal of the current work integrating multiple highly specialized structures involved in the was to establish an atlas of protein-coding gene expression of the regulation of memory, behavior, learning, sensory function, motor caprine brain by sequencing the transcriptomes of 12 encephalic skills and body homeostasis, amongst others. During embryo regions in seven Murciano-Granadina goats. Given that three of development, three vesicles emerge from the neural tube (Chhetri these goats were pregnant at the time of slaughter, we have also and Das, 2021), which roughly correspond to the forebrain investigated whether gestation affects the transcriptomic profiles of (prosencephalon), midbrain (mesencephalon) and hindbrain the 12 brain regions under study. (rhombencephalon). Further subdivisions take place in the forebrain (telencephalon and diencephalon) and hindbrain (metencephalon and myelencephalon), resulting in the formation 2 Materials and methods of five encephalic vesicles (Chhetri and Das, 2021). In the bovine embryo, these five encephalic vesicles are visible at 24 days after 2.1 Sample collection conception, but they do not become fully developed until 110 days of gestation (Ferreira et al., 2018). Multiple anatomical structures with Seven adult non-lactating Murciano-Granadina goats raised in specialized functions (Chhetri and Das, 2021) are subsequently the experimental farm of the Faculty of Veterinary Sciences at the derived from the telencephalon (cerebral hemisphere, basal Universitat Autònoma de Barcelona (UAB) were slaughtered ganglia, hippocampus, olfactory bulb, lateral ventricles, etc.), because of reasons unrelated with this project. Sampled goats diencephalon (thalamus, hypothalamus, pineal body, were kept under the same management and environmental neurohypophysis, infundibulum, third ventricle, etc.), conditions and had similar ages (6.28 ± 1.38 years). Besides, three mesencephalon (rostral colliculus, tegmentum, crus cerebri, of these goats were 1-month pregnant at the time of slaughtering. To cerebral aqueduct, etc.), metencephalon (pons, cerebellum, upper minimize pain, goats were administered pentobarbital (150 mg/kg) part of the fourth ventricle, etc.) and myelencephalon (medulla in the jugular vein. A lethal intravenous dose of sodium oblongata, spinal cord lower part of the fourth ventricle, etc.). pentobarbital causes a mammal to lose consciousness within Transcriptomic analyses may hold the key to significantly seconds and results in clinical death within just minutes, but we advance our knowledge about the biological functions of brain cannot completely rule out the possibility that the euthanasia regions. In mice, patterns of gene expression have been used to procedure used in the current work might have altered the brain establish a molecular atlas of the adult brain (Ortiz et al., 2020), and mRNA expression of certain genes. Since animals were killed due to as much as 737 brain structures have been identified with the routine culling process implemented in this experimental farm, microscopy techniques complemented with other approaches, no permission from the Ethics Committee on Animal and Human providing a comprehensive view about the high functional Experimentation at UAB was required. complexity of this organ (Erö et al., 2018). An important feature After slaughtering, goats were transported to the Necropsy room of the brain is its high plasticity, so any atlas of gene expression is of the Faculty of Veterinary Sciences at UAB. The cranial vault of necessarily dynamic, not only in space but also in time. For instance, each goat was opened with a bone saw and twelve brain anatomical pregnancy in mice is associated with extensive changes in gene structures were carefully dissected and biopsied by an expert expression in the neocortex, cerebellum, hippocampus and anatomist. Samples were drawn from the adenohypophysis, hypothalamus, and there is evidence that several of such cerebellar hemisphere, cerebellar trunk or peduncle, frontal modifications might be long lasting (Ray et al., 2016; Arbeitman, neocortex, hippocampus, hypothalamus, medulla oblongata, 2019). neurohypophysis, pineal gland, pons, olfactory bulb and rostral Very few atlases of brain gene expression have been generated in colliculus. Biopsies were immediately submerged in RNAlater domestic animals. Recently, Sjöstedt et al. (2020) investigated the (Thermofisher Scientific, Barcelona, Spain) to be stored at -80°C gene expression profiles of 10 major mammalian brain regions in until processing. humans, mice, and pigs. They found that global transcriptomic To purify total RNA, tissue samples were mixed with 1 ml profiles are, in general, well conserved in these three species, with QIAzol (QIAGEN Inc., Barcelona, Spain) and homogenized using cerebrum and brainstem regions clustering apart in hierarchical the Lysing Matrix D reagent (MP Biomedicals, Santa Ana, CA) in a trees, and the cerebellum showing a highly divergent profile of gene Precellys 24 tissue homogenizer (Bertin Instruments, Rockville, expression (Sjöstedt et al., 2020). The mRNA expression of several MD). The extraction of total RNA was performed using the brain regions has also been reported in sheep (Clark et al., 2017) and RNeasy lipid tissue mini kit (QIAGEN Inc., Barcelona, Spain) cattle (Harhay et al., 2010), although not comprehensively or following the protocol described by the manufacturer. The Frontiers in Genetics 02 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 concentration and purity of extracted RNA molecules were analyzed 2.4 Construction of a brain atlas based on using the Nanodrop ND-1000 spectrophotometer (Thermofisher the expression of protein-coding genes Scientific, Barcelona, Spain) and RNA integrity was assessed with a Bioanalyzer-2100 equipment (Agilent Technologies, Santa Clara, Genes expressed in each tissue were cataloged according to their CA) using the RNA 6000 Nano Kit 4.2 (Agilent Technologies, Santa Ensembl biotype. To remove systematic technical effects (e.g., Clara, CA). library size), count data from protein-coding genes were normalized to CPM with the edgeR package (Robinson et al., 2009). Finally, data were log2 transformed adding a pseudo-count 2.2 Sequencing of total RNA of 1. Each tissue was normalized separately. Normalization factors were estimated by jointly considering pregnant and non-pregnant Paired-end sequencing (2 x 50 bp) of total RNA was carried out goats because library size features were quite similar in both groups. at the Centre Nacional de Anàlisi Genòmica (CNAG). Sequencing To perform downstream analyses, only genes with a CPM over 0.5 in methods have been reported by Guan et al. (2020). Briefly, the RNA- at least two samples were considered. Seq library was prepared with the KAPA Stranded mRNA-Seq In order to detect errors during the sampling procedure, which Illumina Platforms Kit (Roche, Sant Cugat, Spain) by using might alter our results (accidental cross-contaminations when 500 ng total RNA as template. Oligo-dT magnetic beads were retrieving samples, errors when annotating the anatomical origin used to enrich the poly-A fraction and subsequently, RNA was of samples etc), the quality of the regional sampling of the brain was fragmented. Strand cDNA synthesis was performed in the presence evaluated in two ways. First, Euclidean distances were estimated in a of dUTP to enforce strand-specificity. The blunt-ended double pairwise manner based on the corrected level of gene expression stranded cDNA was 3′-adenylated and ligated to Illumina (log2(CPM + 1)) considering all protein-coding genes (17,054 genes) adaptors with unique dual indexes and unique molecular to identify outlier samples. The distribution of distances between identifiers (Integrated DNA Technologies, Coralville, IA). samples from the same tissue was analyzed using the boxplot.stats Enrichment of the ligation product was ensured by performing function from the grDevices package in R (https://r-universe.dev/ 15 cycles of polymerase chain reaction amplification. An Agilent manuals/grDevices.html), and samples deviating more than 1.5 2100 Bioanalyzer equipment was employed to verify the quality of times from the interquartile range were considered as outliers the final library by using the DNA 7500 assay (Agilent Technologies, and, in consequence, eliminated. Additionally, we built a Inc., Santa Clara, CA). Library sequencing was carried out with a heatmap based on the gene expression level of 61 genes selected HiSeq 4000 equipment (Illumina, San Diego, CA) in accordance by their high tissue-specificity in a set of 12 human tissues (https:// with the protocol for dual indexing advised by the manufacturer. www.proteinatlas.org/, version 22) equivalent to those investigated Image analysis, base calling and quality scoring of the sequencing in goats herewith. This analysis was performed to make sure that the run were checked with the Real-Time Analysis (RTA 2.7.7) tool expression patterns of the goat samples match what is expected (Illumina, San Diego, CA) and FASTQ sequence files were based on human transcriptomic data. subsequently generated. 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 2.3 Quality control, alignment and applying the Ward D clustering method with the hclust function quantification from the R stats package (Murtagh and Legendre, 2014). For each brain region, the top 1,000 genes with the highest expression levels The quality of the sequences was assessed with the FastQC were selected to perform pathway enrichment analyses with the software v.0.11.9 (http://www.bioinformatics.babraham.ac.uk/ Enrichr R package (Chen et al., 2013). The KEGG release projects/fastqc/) and adapters were trimmed using the 99.0 human database (Kanehisa et al., 2016) was used as TrimGalore v.0.6.6 tool (https://www.bioinformatics.babraham.ac. reference to annotate pathways enriched in the 12 sets of genes uk/projects/trim_galore/). In addition, reads with more than five (one for each brain region). In order to identify such enriched ambiguous bases (Ns) were removed. Filtered sequences were pathways, a combined score (c) was estimated from the output of a aligned to the goat ARS1 reference genome (Bickhart et al., 2017) traditional Fisher exact test and a Z-score was estimated from with the HISAT2 v.2.2.1 aligner using paired end option (-1 and -2) expected ranks and variances of the set of genes. In accordance and default parameters (Kim et al., 2015). Gene assembly and with Chen et al. (2013) the following formula was applied: expression quantification were performed with StringTie v.2.1. c  logp × z 0 by extracting, with the -e option, those genes annotated in the reference GTF file as reported in a previous protocol (Pertea et al., where c is the combined score, p corresponds to p-values estimated 2016). The estimated count matrix was obtained from coverage with the Fisher exact test, and z is the Z-score from the values of each feature (genes in our case) using the dedicated script expected rank. “prepDE.py” from the StringTie pipeline and based on the following The tissue-specificity expression of protein-coding genes was formula (https://github.com/gpertea/stringtie/blob/master/ assessed by calculating the tissue specificity index tau (τ) defined by prepDE.py): Yanai et al. (2005): coverage × gene length ni1 (1 − xi ) Reads per gene  τ read length n−1 Frontiers in Genetics 03 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 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. Brain region RIN Average alignment Average number of reads Number of expressed protein-coding Mean ± SD rate (%) (millions) Mean ± SD genes (CPM > 0.5) Adenohypophysis 7.900 ± 0.428 84.44 37.786 ± 2.348 15,220 Cerebellar hemisphere 8.100 ± 0.486 94.11 41.461 ± 4.571 14,889 Cerebellar trunk/ 8.143 ± 0.550 94.17 42.117 ± 9.052 14,898 peduncle Frontal neocortex 7.400 ± 0.569 94.26 41.361 ± 7.438 15,098 Hippocampus 7.629 ± 0.411 93.98 41.717 ± 7.825 15,437 Hypothalamus 6.943 ± 0.326 93.92 38.618 ± 4.817 15,557 Medulla oblongata 6.843 ± 0.181 94.27 40.141 ± 5.510 15,420 Neurohypophysis 7.771 ± 0.482 87.8 39.443 ± 5.931 15,366 Olfactory bulb 7.214 ± 0.422 94.16 40.139 ± 3.938 15,581 Pineal gland 8.029 ± 0.594 94.23 39.228 ± 7.142 15,592 Pons 7.186 ± 0.654 94.35 39.505 ± 4.233 15,273 Rostral colliculus 7.243 ± 0.412 93.63 38.832 ± 3.179 15,153 where n is the number of brain regions, and xi corresponds to the non-pregnant goats were considered as the baseline reference, expression profile component normalized by the maximal meaning that any upregulation or downregulation of gene component value (Yanai et al., 2005). This calculation was expression is referred to changes in the pregnant goats with performed using the tspex v.0.6.1 program (Camargo et al., regard to the non-pregnant. The annotation of differentially 2020). A τ-value of 0 would correspond to a housekeeping gene expressed genes (DEGs) was performed as explained in section with ubiquitous tissue expression. In contrast, genes with a restricted 2.4, while the enrichment pathway analyses were carried out regional expression show τ-values above 0.85, while genes with τ = independently in the sets of upregulated and downregulated 1 can be assumed to be expressed in a single brain region. genes. The level of enrichment was estimated by calculating a c score (Chen et al., 2013), as previously defined. 2.5 Differential expression and pathway enrichment analyses 3 Results Before doing the differential expression analysis, the level of 3.1 Sequencing of total RNA from 12 goat variability across replicates from the same condition (pregnant or brain regions non-pregnant) was assessed by calculating the coefficient of variation (CV) of the expressed genes and genes differentially The average RNA integrity number (RIN) of the 84 RNA expressed between conditions: extractions was 7.53 ± 0.62, ranging from 6.2 to 9 (Table 1). σ These values are in agreement with the conservative threshold of CV  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 where μ is the average expression of each gene and σ is the Table 1, the average sequencing depth ranged from 37.78 to corresponding standard deviation. 41.71 million paired-end reads per sample across all 12 brain The software DEseq2 (Love et al., 2014) was used to perform a regions, and the average alignment rate ranged from 84.44% differential expression analysis comparing the transcriptomic (adenohypophysis) to 94.35% (pons). The read counts per gene profiles of pregnant vs. non-pregnant goats. Genes with a per sample, including pregnant and non-pregnant goats are number of estimated counts below 10 were removed. Correction available in Supplementary Table S1. 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) 3.2 Clustering of brain regions according to absolute logarithm of the fold change |log2FC| > 0.58, and 2) their mRNA expression profile Statistical significance threshold: q-value < 0.05. Differentially expressed genes were used to perform pathway enrichment The total number of protein-coding loci expressed (counts per analyses using the Enrichr R package (Chen et al., 2013). The million, CPM > 0.5 in at least two samples) in each tissue ranged Frontiers in Genetics 04 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 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. from 14,889 genes for the cerebellar hemisphere to 15,592 genes for marker genes in close correspondence to expectations based on the pineal gland. Considering all brain regions, 17,054 protein- human data, with the only exception of pons and medulla oblongata. coding genes passed the quality control (QC) filters and were Samples from the medulla oblongata display high mRNA levels of expressed in at least two samples. One pineal gland sample genes PHOX2A and SLC6A2, which are considered markers for the deviated more than 1.5 times from the interquartile range, so it pons based on human data. In congruence with our results, the was identified as an outlier and excluded from any further analysis. expression of these genes is reported to be high in the medulla As a complementary measure to detect sampling errors, we built a oblongata of mice (Sunkin et al., 2013). Furthermore, the gene EXD1 heatmap based on the expression levels of 61 marker genes which in selected as a marker for pons was highly expressed in the pineal humans show high tissue specificity for each one of the brain regions gland of goats. This gene is also expressed at high levels in the considered in this research. As shown in Supplementary Figure S1, pituitary gland of pigs, and at moderate levels in multiple brain the expression patterns of the samples based on the selected 61 genes structures of pigs and mice (Sjöstedt et al., 2020). We also evaluated is concordant across replicates from the same brain region, thus the variability of the mRNA expression patterns across replicates by evidencing that dissection and sampling were correctly performed. calculating the corresponding coefficient of variation per group. In In general, all brain regions display high expression levels of the general, the coefficient of variation was lower within the pregnant Frontiers in Genetics 05 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 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. and not pregnant groups than when all samples were pooled compared to the remaining encephalic regions (Figure 1A; Figure 2). together in a single group (Supplementary Figure S2). This result The dendrogram indicated that pons and medulla oblongata, which was even more evident when analyzing the sets of differentially are hindbrain structures, tended to group together. Moreover, expressed genes in brain regions moderately or highly affected by frontal neocortex, olfactory bulb, and hippocampus, which derive pregnancy (Supplementary Figure S2B). from the forebrain telencephalon vesicle, also tended to group In close agreement with data reported above, the PCAs, with together. In strong contrast, cerebellum trunk and hemisphere (Figure 1A) and without (Figure 1B) cerebellum and glandular clustered far apart from any other brain region (including pons samples, and the dendrogram (Figure 2) showed that most despite its anatomical proximity), and adenohypophysis (derived samples clustered according to their region of origin. We also from the oral ectoderm) and neurohypophysis (derived from the observed that four regions, namely adenohypophysis, diencephalon), which have completely different embryonic origins, neurohypophysis, pineal gland, and cerebellum (hemisphere and also grouped jointly and clearly separated from the remaining brain trunk), displayed highly divergent gene expression patterns when structures. Frontiers in Genetics 06 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 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. 3.3 Regional specificity of mRNA expression tissue-specific genes (τ ≥ 0.85) identified for each tissue is shown in Supplementary Figure S3. The distribution of the tissue-specificity index τ-values was As shown in Supplementary Table S2B, the gene pathway highly skewed to the left (Figure 3A). Indeed, 5,655 (33.16%) enrichment analysis for loci with an ubiquitous pattern of expression protein-coding genes showed τ-values below 0.15, indicative of a (τ = 0–0.15) revealed that they are mainly involved in general biological highly ubiquitous pattern of expression across all brain regions processes such as RNA transport (combined score, c = 110.230, q-value = (Supplementary Table S2A). In contrast, 10,036 (58.84%) 1.16 E-13), autophagy (c = 105.214, q-value = 5.66 E-12), oxidative protein-coding genes displayed intermediate τ-values (τ = phosphorylation (c = 110.941, q-value = 3.34 E-12), spliceosome (c = 0.15–0.85) and 1,363 protein-coding genes (7.99%) had a 124.950, q-value = 1.16 E-13), proteasome (c = 185.773, q-value = 5.52 E- highly region-specific profile (τ > 0.85) of mRNA expression 10) and ubiquitin mediated proteolysis (c = 127.938, q-value = 1.56 E- (Figures 3A, B). A heatmap based on the level of expression of 13), as well as in the development of neurological diseases such as Frontiers in Genetics 07 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 amyotrophic lateral sclerosis (c = 213.777, q-value = 5.57 E-26), sample decreased the number of replicates from the pregnant Parkinson disease (c = 195.845, q-value = 1.86 E-21), and group to two, which is insufficient to make a sound statistical Huntington disease (c = 142.384, q-value = 3.91 E-19), to mention a few. analysis. When comparing pregnant vs. non-pregnant goats, As indicated in Supplementary Table S3 and Supplementary Figures brain regions showed strong differences in terms of the number S4–S15, we have also investigated whether genes with highly region- of differentially expressed genes. Six regions displayed little changes specific gene expression patterns (τ = 0.85–1, listed in Supplementary in their expression levels in response to 1 month-pregnancy: in the Table S3A) are enriched in specific functional pathways in all analyzed cerebellar hemisphere, cerebellar trunk, hypothalamus, medulla regions (Supplementary Table S3B–S3M). The neuroactive ligand- oblongata, neurohypophysis and rostral colliculus, only 2, 1, 1, receptor interaction pathway was highly enriched in genes with τ = 12, 4 and 1 DEGs were identified, respectively (Figure 4B, 0.85–1 in several of the analyzed regions, i.e., frontal neocortex, Supplementary Table S5). In strong contrast, we observed hippocampus, hypothalamus, and olfactory bulb. Other pathways remarkable changes in the expression profiles of five brain associated with neuronal transduction were also enriched in genes regions (Figure 4A, Supplementary Table S5): adenohypophysis with a tissue-restricted expression, e.g., serotonergic synapse (201 DEGs, 13 downregulated and 188 upregulated in pregnant (cerebellar hemisphere and cerebellar trunk with q-values < 0.05 and goats, Supplementary Table S5A), frontal neocortex (82 DEGs, c of 74.868, and 65.095, respectively) and GABAergic synapse (cerebellar 37 downregulated and 45 upregulated in pregnant goats, hemisphere and trunk, with q-values < 0.05 and c of 69.725 and 60.823, Supplementary Table S5D), hippocampus (70 DEGs, respectively). Besides, the cholinergic synapse pathway showed 15 downregulated and 55 upregulated in pregnant goats, enrichment in the rostral colliculus (c = 43.169, q-value = 0.05). Supplementary Table S5E), pons (190 DEGs, 57 downregulated Several brain regions showed an enrichment pattern that agreed well and 133 upregulated in pregnant goats, Supplementary Table S5J) with their biological functions (Supplementary Tables S3B–S3M). For and, most remarkably, olfactory bulb (1207 DEGs, instance, in the pineal gland the phototransduction pathway showed a c 381 downregulated and 826 upregulated in pregnant goats, of 1101.90 (q-value = 1.06 E-13). 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 3.4 KEGG pathway enrichment analysis of (Supplementary Table S6D), the majority of biological functions highly expressed genes enriched in the set of upregulated DEGs in pregnant goats were associated with axon guidance (c = 89.246; q-value = 7.29 E-07), The results of the Kyoto encyclopedia of genes and genomes dopaminergic synapse (c = 79.526; q-value = 6.07 E-06), neuroactive (KEGG) pathway enrichment analysis for the 1,000 mRNA genes ligand-receptor interaction (c = 53.217; q-value = 5.27 E-06), with the highest expression in each region (Supplementary Tables cholinergic synapse (c = 55.712; q-value = 1.32 E-04) and S4A–S4K) revealed several pathways broadly shared across brain phenylalanine, tyrosine and tryptophan biosynthesis (c = 53.087; regions and also reaching strong statistical significance (q-value < q-value = 8.51 E-02). In contrast, downregulated DEGs in the 0.05), e.g., synaptic vesicle cycle, endocrine and other factor- olfactory bulb of pregnant goats were enriched in pathways regulated calcium reabsorption, endocytosis, phagosome, long-term related to nervous development and control of cellular activities, potentiation, dopaminergic synapse, glutamatergic synapse, circadian including migration, differentiation, and proliferation entrainment, gap junction, gastric acid secretion and adrenergic (Supplementary Table S6D). A heatmap based on the expression signaling in cardiomyocytes. We also detected many pathways levels of genes differentially expressed in five tissues moderately or related to neurological conditions that attained high statistical highly affected by pregnancy (i.e., adenohypophysis, frontal significance (q-value < 0.05) and were represented in a broad array neocortex, hippocampus, olfactory bulb, and pons) is shown in of encephalic regions, e.g., Alzheimer disease, Huntington disease, Supplementary Figure S21. It can be seen that the main determinant Parkinson disease, prion disease, amyotrophic lateral sclerosis, and, of gene expression patterns is the region-of-origin of samples rather to a much lesser extent, spinocerebellar ataxia. Pathways related with than pregnancy. Vibrio cholerae and Salmonella infection were also commonly found. By making a literature search for each one of the genes showing Protein processing in endoplasmic reticulum (q-value < 0.05), RNA differential expression in the olfactory bulb of pregnant vs. non- transport (q-value < 0.05) and ribosome (q-value < 0.05) were amongst pregnant goats, we observed that many of them are related to human the most significant pathways in the adenohypophysis (Supplementary behavioral traits (Supplementary Table S7), including, but not Table S4A) and neurohypophysis (Supplementary Table S4H), but not limited to: maternal behavior (e.g., downregulated: DDC; in the other regions. As said before, in the pineal gland (Supplementary upregulated: DBH, DRD1, NTS and HTR2A), affective behavior, Table S4J) phototransduction was also a highly significant pathway (q- sociability and exploration (e.g., downregulated: FGFR4, GRIP2, value = 0.006) while in the remaining brain regions it did not reach NTRK2 and RYR3; upregulated: HTR2C, FBX O 45, GPR3, statistical significance. 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., 3.5 Analysis of differential gene expression in downregulated: MOCOS, CACNA1D, MBD6 and AUTS2; pregnant and non-pregnant goats upregulated: RAB39B, BTBD11, KCNQ3, CDH9, CADPS2, DOCK4 and EXT1), aggression (e.g., downregulated: NOS1; Samples from the pineal gland were not included in the upregulated: HRH3, PRNP, HNMT and GRIA3), cognition, differential expression analysis, since the removal of an outlier memory and learning (e.g., downregulated: EPHA10, STAT5 and Frontiers in Genetics 08 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 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. MMP28; upregulated: EPHA6, SYP, SORBS2, ARHGEF4, HCN1, As shown in Supplementary Table S8, our literature search also CAMK2N2, MMP17, BTBD9, CLSTN3, STAU2RIMKLA, PAK6, revealed several DEGs associated with the migration and maturation SLC22A4 and NEURL1), response to stress (e.g., downregulated: of olfactory bulb interneurons (e.g., PROKR2 and SALL3), olfactory IFIT1; upregulated: PDYN, EPOP, DPYSL2 and HCN2), feeding bulb angiogenesis (e.g., ANOS1), olfactory bulb morphogenesis and behavior (e.g., upregulated: NELL2, CXCL14, GPR45, GPR162, NPY establishment of a functional olfactory neural circuitry (e.g., ARX, and ACBD7) and diverse neuropsychiatric disorders (e.g., TSHZ1, OLFM2 and FGFR1), axonal growth and guidance (e.g., downregulated: SLITRK6, IL1RAPL1, BAHCC1, WDR62 and RAP1GAP2, SHH, GLI3 and DSCAM) and sensitivity of olfactory SLC6A1; upregulated: SCN1A, GABBR2, RTN4R, ACOT7, KIF5A. sensory neurons to food cues (e.g., NPY). A heatmap displaying the HECW2, NRG1, SNCA, ATP8A2, WFS1, HTR5A, STX1A and JPH3). gene expression patterns of 207 genes related to behavior Frontiers in Genetics 09 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 (Supplementary Table S7) and 53 genes involved in olfaction while pons and medulla oblongata come from the rhombencephalon. (Supplementary Table S8) that show differential expression in the In the comprehensive analysis of porcine tissue expression carried out olfactory bulb of pregnant vs. non-pregnant goats is shown in by Ferraz et al. (2008), tissues clustered according to the germ layer Supplementary Figure S22. (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 4 Discussion anterior-posterior axis of the neural tube and gene expression in different regions from the adult brain (Zapala et al., 2005), as also 4.1 The clustering of samples from certain published by Ortiz et al. (2020). These published results and those brain regions follows an anterior-posterior obtained by us agree with the hypothesis that embryogenesis might pattern consistent with embryonic origin leave a durable footprint in the profile of mRNA expression of mammalian brain regions. 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 4.2 Evolutionary conservation of the 469 genes (2.75%) were expressed in just one of the 12 sampled brain genome-wide patterns of expression across regions (τ = 1). Results obtained in the Genotype-Tissue Expression brain regions (GTEx) project were consistent with this finding (Melé et al., 2015), since only 200 genes showed tissue-specific expression (95% of these Sjöstedt et al. (2020) described the patterns of genome-wide genes were exclusively expressed in the testis). expression of protein-coding genes in a number of brain regions Brain region of origin, rather than pregnancy, was the major from humans, mice and pigs. They found that, in general, factor explaining the clustering of samples. By using microarrays, transcriptomic patterns are evolutionarily conserved in these three Ferraz et al. (2008) investigated the profiles of expression of species, with the three forebrain structures (cerebral cortex, 16 porcine tissues and observed that the factor “tissue of origin” hippocampus, and amygdala) grouping together, and the midbrain, accounted for ~11 times more variability than sex or breed. In thalamus, and pons and medulla samples forming another cluster close humans, the GTEx Consortium retrieved 1641 post-mortem to the hypothalamus (Sjöstedt et al., 2020). They also found that the samples covering 54 body sites from 175 individuals and cerebellum behaves as an outlier, a finding fully coherent with ours. reported that tissue type was the primary factor explaining Although transcriptomic profiles of brain regions seem to be well differences in gene expression (Ardlie et al., 2015), a result that is conserved across species, relevant differences also exist. For instance, fully consistent with ours. Similar findings have been obtained when Sjöstedt et al. (2020) indicated that in humans the olfactory bulb clusters building atlases of gene expression in cattle (Harhay et al., 2010), with other forebrain structures (cerebral cortex, hippocampus, pigs (Freeman et al., 2012) and sheep (Clark et al., 2017). We have amygdala, etc.), while in pigs and mice this brain structure almost observed that a few samples do not cluster with the remaining behaves as an outlier. They reasoned that this might be due to the fact samples from the same region of origin (see for instance the two that these two latter species have more evolved olfactory systems than hippocampus samples in the upper left quadrant of the PCA shown humans (Sjöstedt et al., 2020). In goats, the olfactory bulb is also much in Figure 1B). This lack of consistent clustering has been observed in more developed than in humans (Kavoi and Jameela, 2011), a feature other atlases of gene expression (Harhay et al., 2010; Ardlie et al., consistent with the fundamental role of the olfactory system in the 2015; Clark et al., 2017). Rather than mix-ups or contaminations, maternal and social behaviors of ungulates (Keller and Lévy, 2012). unexpected patterns of clustering could be due to the expression of However, our data indicate that in goats the pattern of mRNA thousands of housekeeping genes, that are transcribed at similar expression of the olfactory bulb is closely aligned with that of other rates in different tissues, and, additionally, to the existence of forebrain regions (Figure 1B; Figure 2). This result agrees much more “communication pathways” between tissues which might with the PCA reported by Sjöstedt et al. (2020) for human brain contribute to generate a high correlation of expression profiles structures than for those corresponding to pigs or mice. A higher between them, as evidenced in previous studies (Wang et al., 2019). resemblance of goats to humans, rather than to pigs, is unexpected We have observed that, in goats, frontal neocortex, olfactory bulb because goats and pigs are ungulates (Benton et al., 2015). Sjöstedt et al. and hippocampus samples tend to group together, and the same trend (2020) stated that the functional importance of olfaction might be a key is observed for rostral colliculus, medulla oblongata and pons samples, factor explaining the differential patterns of olfactory bulb mRNA while hypothalamus samples would be placed in the middle of these expression observed in pigs/mice vs. humans. Such hypothesis is not two clusters. In strong contrast, cerebellum, pineal gland and adeno/ fully consistent with our results and should be interpreted cautiously neurohypophysis show highly differentiated patterns of gene until further data are available in other mammalian species. 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 4.3 Cerebellum, pineal gland and hypophysis located closer to hypothalamus than to the hippocampus, or that the have highly differentiated patterns of gene distance between the rostral colliculus and the hippocampus is shorter expression than that between the olfactory bulb and the hippocampus. The common denominator between frontal neocortex, olfactory bulb and Cerebellum displayed a highly differentiated pattern of mRNA hippocampus is that these three organs derive from the telencephalon, expression when compared to other brain structures (Figure 1A; Frontiers in Genetics 10 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 Figure 2), a finding that, as mentioned before, is consistent with Further studies are needed to assess whether such interpretation is previous reports (Ortiz et al., 2020; Sjöstedt et al., 2020). Moreover, correct. 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 4.4 Heterogeneous effects of pregnancy on it has been suggested that it is scarcely connected with cognitive the expression profiles of eleven goat brain areas of the brain cortex (Glickstein and Doron, 2008). One of the regions most distinctive features of the cerebellum is its extraordinarily high cell density: while this organ represents 10% of brain volume, it Amongst the seven individuals sampled in our study, there were encompasses 42% of all brain cells and 59.8% of all excitatory three pregnant goats providing the opportunity to investigate the neurons, mainly due to the conspicuous abundance of tightly effect of pregnancy on the brain transcriptome. Since the number of packed granular cells (Erö et al., 2018). Another relevant feature replicates for each pregnant (N = 3) and non-pregnant (N = 4) of the cerebellum is that the same circuit, composed by mossy fibers category is low, sensitivity to detect DEGs is expected to be quite that excite granule cells that, in turn, excite Purkinje cells, forms a limited (Schurch et al., 2016). To circumvent this difficulty, we have fundamental unit that is replicated thousands of times (Kozareva used the total number of DEGs as the main criterion to identify et al., 2020). These biological particularities might have contributed which brain regions are mostly affected by 1-month pregnancy. In to transform the cerebellum into an organ with a very specific principle, low number of replicates should affect the number of transcriptomic profile, not only in goats but also in other mammals. detected DEGs to a similar extent in all 11 brain regions (pineal We also observed highly divergent patterns of gene expression in gland was excluded from the differential expression analysis), so the three glandular tissues: the pineal gland and the adeno/ number of DEGs seems to be an appropriate indicator of which neurohypophysis (Figure 1A), a finding consistent with results regions are more affected by pregnancy. reported by Harhay et al. (2010). The strong functional In comparisons involving cerebellar hemisphere, cerebellar specialization of these three anatomical structures in hormonal trunk, hypothalamus, medulla oblongata, neurohypophysis and secretion might explain these findings. In response to light, the rostral colliculus, the total number of DEGs when comparing pineal gland, which develops from an evagination of pregnant vs. non-pregnant goats was very low or inexistent neuroepithelium in the dorsal midline of the diencephalon, (Figure 4B). In contrast, between 70 and 201 DEGs were synthesizes and releases melatonin, which is a key regulator of detected in the frontal neocortex, hippocampus, adenohypophysis the circadian sleep-wake cycle and seasonal rhythms (Patel et al., and pons (Figure 4A). By far, the organ which displayed the largest 2020). Remarkably, in our study the pineal gland was the second number of DEGs was the olfactory bulb (826 upregulated and organ showing the highest number of genes with region-specific 381 downregulated genes in pregnant goats). From these data, we expression (τ > 0.85, Figure 3B and Supplementary Table S3A). conclude that 1-month pregnancy does not have the same effect on Neurohypophysis and adenohypophysis are also specialized in the all goat brain regions. Ray et al. (2016) investigated changes in the secretion of molecules with key physiological roles, and data mRNA expression of four brain structures (hypothalamus, collected in sheep and cattle indicate that their profiles of mRNA neocortex, hippocampus and cerebellum) in virgin, pregnant, and expression are highly differentiated from those of other brain postpartum mice, and they found that in the virgin vs. pregnant structures (Harhay et al., 2010; Clark et al., 2017). The main comparison the number of DEGs was much higher in the hormones produced by the adenohypophysis are prolactin, hippocampus than in the cerebellum, while hypothalamus and adrenocorticotropic hormone, luteinizing hormone, follicle- neocortex showed intermediate values. These results support the stimulating hormone, growth hormone and thyroid-stimulating notion that pregnancy does not affect all brain regions to the same hormone, which regulate a very diverse set of biological processes extent. including growth, metabolism, lactation, stress and reproduction (Le Brain gene expression during pregnancy is dynamic and it might Tissier et al., 2012). One particular feature of the adenohypophysis is change depending on the time point under consideration. By using that it does not develop from any of the five neural vesicles but from magnetic resonance imaging, it has been shown that the brain of the oral ectoderm (Le Tissier et al., 2012). In contrast, the primiparous women experiences substantial morphological changes neurohypophysis has a diencephalic origin and stores vasopressin during gestation which mostly affect the right middle temporal and oxytocin (both are synthesized in the hypothalamus), two gyrus, inferior frontal gyrus and posterior cingulate cortex hormones regulating diuresis and a broad array of reproduction (Hoekzema et al., 2017). A longitudinal morphometric study in and behavioral processes, respectively (Santos Fontanez and De mice also provided evidence of transient hypertrophy associated Jesus, 2021). Despite having completely different embryological with gestation and/or lactation in the medial preoptic area, bed origins, histological structure and biological functions, our data nucleus of the stria terminalis, amygdala, caudate nucleus, and indicate that neurohypophysis and adenohypophysis share similar hippocampus (Barrière et al., 2021). These findings indicate that profiles of mRNA expression (Figure 1A; Figure 2). Our certain brain structures undergo significant alterations in their interpretation is that high functional specialization in hormonal profiles of expression and morphology because of gestation, and secretion might erase, at least partially, the transcriptomic footprint that biological changes are highly dynamic. Thus, the lists of DEGs associated with embryogenesis, a hypothesis supported by the detected for 11 goat brain regions at 1-month of gestation might not distinctive patterns of gene expression observed in both pituitary be representative of the whole gestation, but just of this specific time structures when compared with the remainder encephalic regions. point. Frontiers in Genetics 11 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 In the goat frontal neocortex and the hippocampus, 82 DEGs (Inestrosa and Varela-Nallar, 2015; Renault et al., 2009; Stocker (45 upregulated and 37 downregulated in pregnant goats) and and Chenn, 2015). 70 DEGs (55 upregulated and 15 downregulated in pregnant Amongst this list of 94 loci, one gene that is particularly goats) were detected, respectively (Supplementary Tables S5D interesting is that encoding the dopa decarboxylase (DDC) mRNA, and S5E). The hippocampus and the frontal neocortex play which was consistently downregulated in the pons and olfactory bulb important roles in modulating memory and learning as well as (log2FC -2.08 and -1.3, respectively). This gene is responsible for the social behavior (Zemla and Basu, 2017; Reinert et al., 2021), and synthesis of dopamine, essential for a healthy pregnancy (Gratz et al., there is evidence that hippocampal neurogenesis is affected by 2018). Dopamine is a key modulator of maternal behavior in rodents, pregnancy (Rolls et al., 2008). In this way, genes included in partly because it inhibits the pituitary release of prolactin which pathways involved in the immune response, including NFATC3 modulates lactation, corpora lutea function (in certain species), and CCL2, and apoptosis (CTSH) were upregulated in these tissues and maternal behavior (Nevard et al., 2022). in pregnant goats (Supplementary Table S6C). Interestingly, these Notably, the direction of differential gene expression was not genes have been associated with neurodegenerative processes in the always consistent across tissues. For instance, the vesicular hippocampus and memory loss (Abdul et al., 2009; Wang et al., associated membrane protein 1 (VAMP1) gene was significantly 2010). Besides, hippocampal lesion in the female rat results in upregulated in the olfactory bulb of pregnant goats (log2FC = 1.15), deficient nest construction and reduced pup survival, suggesting while in the frontal neocortex was downregulated in pregnant goats that this organ has a key influence on the development of maternal (log2FC = -0.67). The VAMP1 gene encodes a protein that plays an behavior (Murphy and Wildeman, 1986). important role in synapsis and synapsis plasticity and forms part of In the adenohypophysis 201 DEGs (188 upregulated and the soluble N-ethylmaleimide-sensitive factor attachment protein 13 downregulated in pregnant goats) were identified receptor complex (Dalmaz et al., 2021). These proteins are involved (Supplementary Table S5A). Adenohypophysis provides in neurocognitive processes and associated with the Wnt pathways hormonal signals that are fundamental for the maintenance of that, among other things, are important for neurogenesis and cell gestation. As listed in Supplementary Table S6A, genes included proliferation (Dalmaz et al., 2021). 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 4.6 Dramatic changes in the mRNA upregulated in pregnant goats. The anterior pituitary is strongly expression of the olfactory bulb in response enlarged during gestation in response to hormones produced by the to early gestation placenta and ovaries, and the secretion of prolactin increases gradually to prepare the female for lactation (Carlson, 2011). The olfactory bulb was the brain region that displayed, by far, the Differential mRNA expression in the pons, with 190 DEGs highest number of DEGs (Supplementary Table S5I), implying that (133 upregulated and 57 downregulated in pregnant goats) is this organ is strongly affected by 1 month-pregnancy. The olfactory harder to interpret (Supplementary Table S5J). This brain region bulb processes smell information transmitted from olfactory sensory is mostly involved in the regulation of breathing and sleep, and in neurons, which express odorant receptors in the olfactory relaying information to or from the cerebellum to other brain epithelium, and relays it to the olfactory cortex (Nagayama et al., regions (Venkatraman et al., 2017). However, there is evidence 2014). Importantly, the sense of olfaction is strongly connected to a that the pons has a basic role in the generation and experience of broad array of behavioral responses related with aversion to food, emotions through the integration of arousal, autonomic function, avoidance of predators, sex arousal, partner preference, and motor control, and somatosensory signals (Venkatraman et al., aggression, amongst others (Sullivan et al., 2015). For instance, 2017). Moreover, it has been reported that preoptic area bilateral bulbectomy of male rats is associated with indifference to projections to lower brainstem regions affect maternal behavior the sexual status of females and suppression of mating (Edwards in postpartum rats (Numan and Numan, 1991). 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 4.5 Several genes are differentially expressed the olfactory bulb observed in pregnant goats, when compared to in two or more brain regions in pregnant vs. their non-pregnant counterparts, are not triggered by odorant non-pregnant goats signals delivered by the offspring. More likely, this remarkable change of gene expression might be promoted by other brain We have detected 94 genes that are differentially expressed in regions and/or endocrine glands delivering chemical signals to more than one tissue when comparing pregnant and non-pregnant the olfactory bulb. As expected, several DEGs are related with the goats (90 genes are differentially expressed in two tissues and four sense of olfaction (Supplementary Table S8). For instance, the genes are differentially expressed in three tissues). The most inactivation of the zinc finger homeodomain factor teashirt zinc significant pathways detected in this list of 94 genes are Wnt finger family member 1 (TSHZ1, |log2FC| 1.36, q-value < 0.05) or of signalling pathway, adherens junctions signalling pathways, and the prokineticin receptor 2 (PROKR2, |log2FC| 1.28, q-value < 0.05) FoxO signalling pathway (Supplementary Table S9). These genes results in olfactory bulb hypoplasia (Prosser et al., 2007; pathways are related to cellular growth and proliferation and Ragancokova et al., 2014); and the loss of the Spalt like they are involved in the formation and modulation of neuronal transcription factor 3 (SALL3, |log2FC| 1.02, q-value < 0.05) circuits, neurogenesis, neuronal migration and axon guidance inhibits the formation of the glomerular layer in the olfactory Frontiers in Genetics 12 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 bulb (Harrison et al., 2008). The aristaless related homeobox (ARX, pregnant goats. Similarly, we have found increased mRNA levels of |log2FC| 0.97, q-value < 0.05) gene is fundamental for the proper genes which promote feeding, e.g., neuropeptide Y, NPY(|log2FC| development of mouse olfactory system (Yoshihara et al., 2005), and 1.58, q-value < 0.05); C-X-C motif chemokine ligand 14, CXCL14 loss of sonic hedgehog signaling molecule (SHH, |log2FC| 1.02, q- (|log2FC| 0.99, q-value < 0.05); and neural EGFL like 2, NELL2 value < 0.05) or olfactomedin 2 (OLFM2, |log2FC| 1.09, q-value < (|log2FC| 3.13, q-value < 0.05); or inhibit it, e.g., G protein-coupled 0.05) impairs olfactory function (Balmer and Lamantia, 2004; receptor 45, GPR45 (|log2FC| 1.56, q-value < 0.05); and acyl-CoA Sultana et al., 2014). Moreover, the formation of the olfactory binding domain containing 7, ACBD7 (|log2FC| 0.63, q-value < bulb requires fibroblast growth factor (FGF) signaling mediated 0.05), as reported in previous studies (Tanegashima et al., 2010; Cui by the FGF receptor one FGFR1 (|log2FC| 0.59, q-value < 0.05) et al., 2016; Jeong et al., 2017; Lanfray and Richard, 2017; Engström (Hébert et al., 2003), and anosmin 1 (ANOS1, |log2FC| 1.68, q- Ruud et al., 2020). These findings illustrate the complexity of value < 0.05), which is essential for blood vessel formation in the ascertaining the functional implications of changes in olfactory olfactory bulb (Matsushima et al., 2020). So, modifications in the bulb mRNA expression associated with pregnancy. expression of these and other genes might result in changes of smell There is also a subset of genes with increased expression in the perception and/or of behavioral responses associated with it as a olfactory bulb of pregnant goats which is particularly interesting transition to maternity. because they display dual functions related with both olfaction and More enigmatic are changes in the mRNA expression of a large behavior. Examples of this category are genes encoding hippocalcin number of genes known to influence behavioral traits in humans (HPCA, |log2FC| 1.60, q-value < 0.05), a mediator of second messenger (Supplementary Table S7). For instance, we have observed an signaling in the olfactory epithelium which also regulates spatial and upregulation, in the olfactory bulb of pregnant goats, of the associative memory (Mamman et al., 2004; Kobayashi et al., 2005), and dopamine β-hydroxylase gene (DBH, |log2FC| 0.93, q-value < 0.05), the brain derived neurotrophic factor (BDNF, |log2FC| 1.78, q-value < which encodes an oxidoreductase that catalyzes the conversion of 0.05), which is important for olfactory bulb neurogenesis, learning and dopamine to norepinephrine. In mice, the inactivation of the DBH memory (Yuan, 2008; Miranda et al., 2019). gene causes a profound impairment of maternal behavior The results of the differential expression analysis should be characterized by the abandonment of pups which die a few days interpreted with caution because sample size is small. Despite this after birth (Thomas and Palmiter, 1997). The dopamine D1 receptor caveat, we consider that the disproportionate number of DEGs detected (DRD1, |log2FC| 1.91, q-value < 0.05) and the 5-hydroxytryptamine in the olfactory bulb when compared to other brain regions is a robust receptor 2A (HTR2A, |log2FC| 0.91, q-value < 0.05) show an increased result. Importantly, this anatomical structure is strongly involved in the expression in the olfactory bulb of pregnant goats. The DRD1 gene is development of maternal behavior in mice, i.e., its complete removal involved in the regulation of mothering in rats (Mileva-Seitz et al., eliminates such capacity and, very often, implies that litters are 2012), while the HTR2A molecule binds to serotonin, a cannibalized by their mothers soon after parturition (Gandelman neurotransmitter upregulated during the transition to motherhood et al., 1971). Moreover, intrabulbar infusions of an oxytocin with strong effects on maternal care and pup survivability (Angoa- antagonist delays maternal behavior in rats, while oxytocin itself has Pérez et al., 2014; Pawluski et al., 2019). However, it is worth the opposite effect (Yu et al., 1996). In sheep, which are phylogenetically emphasizing that interpreting the observed patterns of mRNA close to goats, disruption of the noradrenergic projections to the expression from a functional perspective is not completely olfactory bulb suppresses the establishment of a maternal bond straightforward. For instance, in the olfactory bulb of pregnant between mother and offspring, probably because such event strongly goats we have observed an upregulation of the histamine receptor depends on odor cues (Pissonnier et al., 1985). These results support a H1 (HRH1, |log2FC| 1.43, q-value < 0.05) and the prion protein (PRNP, critical role of the olfactory bulb in the induction of maternal behavior |log2FC| 0.60, q-value < 0.05) mRNAs. The inactivation of these genes in mammals, although discrepant findings have also been published leads to increased aggression (Tomaz et al., 2014; Naganuma et al., (Pissonnier et al., 1985). 2017), while in parallel we have detected an upregulation of histamine It has been proposed that in sheep and goats the increase of receptor H3 (HRH3, |log2FC| 3, q-value < 0.05) mRNA, which in estrogen levels in blood during prepartum, combined with the zebrafish has the opposite effect (Reichmann et al., 2020). Moreover, vaginocervical stimulation provoked by the expulsion of the fetus, we have found that the neurotensin (NTS, |log2FC| 0.86, q-value < 0.05) trigger modifications in a neural network including the main mRNA is upregulated in the olfactory bulb of pregnant females. olfactory system, the medial preoptic area and the paraventricular Interestingly, the intracerebroventricular injection of neurotensin in nucleus of the hypothalamus, in order to induce maternal behavior mice decreases aggression associated with the maternal protection of (Poindron et al., 2007). Our transcriptomic data demonstrate that the offspring (Gammie et al., 2009). extensive areas of the goat brain manifest changes in their expression With regard to olfactory bulb DEGs linked to depression and profiles much before the peripartum period. These alterations might anxiety, we have found an increased mRNA expression of both represent an adaptive response toward preparing the pregnant goat antidepressant (Son et al., 2012; Zhou and Zhou, 2014), e.g., lipid for maternity, and we hypothesize that the olfactory bulb might hold droplet biogenesis associated seipin (BSCL2, |log2FC| 0.64, q-value < a very relevant role in such process. 0.05), neuritin (NRN1, |log2FC| 1.21, q-value < 0.05) and glycine receptor β (GLRB, |log2FC| 0.75, q-value < 0.05), and pro-depressive or pro-anxiety genes (Knoll and Carlezon, 2010; Su et al., 2018; Liu Data availability statement et al., 2021), e.g., salt inducible kinase 2 (SIK2, |log2FC| 1.44, q- value < 0.05), prodynorphin (PDYN, |log2FC| 0.95, q-value < 0.05) The datasets presented in this study can be found in online and histamine receptor H1 (HRH1, |log2FC| 1.41, q-value < 0.05) in repositories. The names of the repository/repositories and accession Frontiers in Genetics 13 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 number(s) can be found below: https://www.ncbi.nlm.nih.gov/, Maria Luigi-Sierra was funded with a PhD fellowship “Formación de PRJNA808876. 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 Ethics statement China Scholarship Council (CSC). Emilio Mármol-Sánchez was funded with a FPU Ph.D. grant from the Spanish Ministry of Ethical review and approval was not required for the animal Education (FPU15/01733). 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 Acknowledgments was completely unrelated with our project, so permission from the Ethics Committee on Animal and Human Experimentation of the Thanks to José Luis Ruiz de la Torre Casañas from the Service of UAB was not required. 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. Author contributions Conceptualization—MA and ML-B; formal analysis—ML-S; Conflict of interest funding acquisition—MA; investigation—ML-S and DG; project administration—MA; resources—ML-B, MA, ML-S, DG, EC, SO- The authors declare that the research was conducted in the M, JG, MP, UIO, ULO, and EM-S; visualization—ML-S; absence of any commercial or financial relationships that could be writing—original draft—MA and ML-S; writing—review & construed as a potential conflict of interest. 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. Publisher’s note All claims expressed in this article are solely those of the authors Funding and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the This research was funded by the European Regional reviewers. Any product that may be evaluated in this article, or Development Fund (FEDER)/Ministerio de Ciencia e claim that may be made by its manufacturer, is not guaranteed or Innovación—Agencia Estatal de Investigación/Project Reference endorsed by the publisher. 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/ Supplementary material Generalitat de Catalunya and of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo The Supplementary Material for this article can be found online Ochoa 2020-2023 (CEX2019-000902-S) grant awarded to the Centre at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1114749/ for Research in Agricultural Genomics (CRAG, Bellaterra, Spain). full#supplementary-material References Abdul, H. M., Sama, M. A., Furman, J. L., Mathis, D. M., Beckett, T. L., Weidner, A. Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate - a M., et al. (2009). Cognitive decline in Alzheimer’s disease is associated with selective practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300. changes in calcineurin/NFAT signaling. J. Neurosci. 29, 12957–12969. doi:10.1523/ doi:10.1111/j.2517-6161.1995.tb02031.x jneurosci.1064-09.2009 Benton, M. J., Donoghue, P. C. J., Asher, R. J., Friedman, M., Near, T. J., and Vinther, Angoa-Pérez, M., Kane, M. J., Sykes, C. E., Perrine, S. A., Church, M. W., and Kuhn, J. (2015). Constraints on the timescale of animal evolutionary history. Palaeontol. D. M. (2014). Brain serotonin determines maternal behavior and offspring survival. Electron. 18, 1–106. doi:10.26879/424 Genes. Brain Behav. 13, 579–591. doi:10.1111/gbb.12159 Bickhart, D. M., Rosen, B. D., Koren, S., Sayre, B. L., Hastie, A. R., Chan, S., et al. (2017). Arbeitman, M. N. (2019). Maternal experience leads to lasting gene expression Single-molecule sequencing and chromatin conformation capture enable de novo reference changes in some regions of the mouse brain. G3 Genes., Genomes, Genet. 9, assembly of the domestic goat genome. Nat. Genet. 49, 643–650. doi:10.1038/ng.3802 2623–2628. doi:10.1534/g3.119.400249 Camargo, A. P., Vasconcelos, A. A., Fiamenghi, M. B., Pereira, G. A. G., and Ardlie, K. G., DeLuca, D. S., Segrè, A. v., Sullivan, T. J., Young, T. R., Gelfand, E. T., Carazzolle, M. F. (2020). tspex: a tissue-specificity calculator for gene expression et al. (2015). Human genomics. The genotype-tissue expression (GTEx) pilot analysis: data. Res. Sq., 1–7. doi:10.21203/rs.3.rs-51998/v1 Multitissue gene regulation in humans. Science 348, 648–660. doi:10.1126/science. Carlson, H. E. (2011). “The pituitary gland in pregnancy and the puerperium,” in The 1262110 pituitary. Editor S. Melmed (Cambridge, MA: Academic Press), 397–412. Balmer, C. W., and Lamantia, A. S. (2004). Loss of Gli3 and Shh function disrupts Chen, E. Y., Tan, C. M., Kou, Y., Duan, Q., Wang, Z., Meirelles, G. V., et al. (2013). olfactory axon trajectories. J. Comp. Neurology 472, 292–307. doi:10.1002/cne.20053 Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Barrière, D. A., Ella, A., Szeremeta, F., Adriaensen, H., Même, W., Chaillou, E., et al. Bioinforma. 14, 128. doi:10.1186/1471-2105-14-128 (2021). Brain orchestration of pregnancy and maternal behavior in mice: A longitudinal Chhetri, K. P., and Das, J. M. (2021). Neuroanatomy, neural tube development and morphometric study. Neuroimage 230, 117776–117815. doi:10.1016/j.neuroimage. stages. Treasure Island (FL): StatPearls Publishing. 2021.117776 Frontiers in Genetics 14 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 Clark, E. L., Bush, S. J., McCulloch, M. E. B., Farquhar, I. L., Young, R., Lefevre, L., Knoll, A. T., and Carlezon, W. A. (2010). Dynorphin, stress, and depression. Brain et al. (2017). A high resolution atlas of gene expression in the domestic sheep (Ovis Res. 1314, 56–73. doi:10.1016/j.brainres.2009.09.074 aries). PLoS Genet. 13, 1006997. doi:10.1371/journal.pgen.1006997 Kobayashi, M., Masaki, T., Hori, K., Masuo, Y., Miyamoto, M., Tsubokawa, H., et al. Cui, J., Ding, Y., Chen, S., Zhu, X., Wu, Y., Zhang, M., et al. (2016). Disruption of (2005). Hippocalcin-deficient mice display a defect in cAMP response element-binding Gpr45 causes reduced hypothalamic POMC expression and obesity. J. Clin. protein activation associated with impaired spatial and associative memory. Investigation 126, 3192–3206. doi:10.1172/JCI85676 Neuroscience 133, 471–484. doi:10.1016/j.neuroscience.2005.02.034 Dalmaz, C., Barth, B., Pokhvisneva, I., Wang, Z., Patel, S., Quillfeldt, J. A., et al. (2021). Kozareva, V., Martin, C., Osorno, T., Rudolph, S., Guo, C., Vanderburg, C., et al. Prefrontal cortex VAMP1 gene network moderates the effect of the early environment (2020). A transcriptomic atlas of the mouse cerebellum reveals regional specializations on cognitive flexibility in children. Neurobiol. Learn Mem. 185, 107509. doi:10.1016/j. and novel cell types. Nature 98, 214–219. doi:10.1101/2020.03.04.976407 nlm.2021.107509 Lanfray, D., and Richard, D. (2017). Emerging signaling pathway in arcuate feeding- Edwards, D. A., Griffis, K. T., and Tardivel, C. (1990). Olfactory bulb removal: Effects related neurons: Role of the Acbd7. Front. Neurosci. 11, 328. doi:10.3389/fnins.2017. on sexual behavior and partner-preference in male rats. Physiol. Behav. 48, 447–450. 00328 doi:10.1016/0031-9384(90)90342-2 Le Tissier, P. R., Hodson, D. J., Lafont, C., Fontanaud, P., Schaeffer, M., and Mollard, Engström Ruud, L., Pereira, M. M. A., de Solis, A. J., Fenselau, H., and Brüning, J. C. P. (2012). Anterior pituitary cell networks. Front. Neuroendocrinol. 33, 252–266. doi:10. (2020). NPY mediates the rapid feeding and glucose metabolism regulatory functions of 1016/j.yfrne.2012.08.002 AgRP neurons. Nat. Commun. 11, 442. doi:10.1038/s41467-020-14291-3 Liu, Y., Tang, W., Ji, C., Gu, J., Chen, Y., Huang, J., et al. (2021). The Selective Erö, C., Gewaltig, M. O., Keller, D., and Markram, H. (2018). A cell atlas for the mouse SIK2 inhibitor ARN-3236 produces strong antidepressant-like efficacy in mice via the brain. Front. Neuroinform 12, 84. doi:10.3389/fninf.2018.00084 hippocampal CRTC1-CREB-BDNF pathway. Front. Pharmacol. 11, 624429. doi:10. 3389/fphar.2020.624429 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. Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change doi:10.1186/1471-2164-9-173 and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/ s13059-014-0550-8 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. Mamman, A., Simpson, J. P., Nighorn, A., Imanishi, Y., Palczewski, K., Ronnett, G. V., et al. Veterinaria Bras. 38, 147–153. doi:10.1590/1678-5150-pvb-5020 (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 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- Matsushima, S., Shimizu, A., Kondo, M., Asano, H., Ueno, N., Nakayama, H., et al. 7007-10-90 (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- Gallego Romero, I., Pai, A. A., Tung, J., and Gilad, Y. (2014). RNA-seq: Impact of 019-57040-3 RNA degradation on transcript quantification. BMC Biol. 12, 42. doi:10.1186/1741- 7007-12-42 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. Gammie, S. C., D’Anna, K. L., Gerstein, H., and Stevenson, S. A. (2009). Neurotensin Science 348, 660–665. doi:10.1126/science.aaa0355 inversely modulates maternal aggression. Neuroscience 158, 1215–1223. doi:10.1016/j. neuroscience.2008.11.045 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 Gandelman, R., Zarrow, M. X., Denenberg, V. H., and Myers, M. (1971). Olfactory behavior. Genes. Brain Behav. 11, 684–694. doi:10.1111/j.1601-183X.2012.00804.x bulb removal eliminates maternal behavior in the mouse. Science 171, 210–211. doi:10. 1126/science.171.3967.210 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 Glickstein, M., and Doron, K. (2008). Cerebellum: Connections and functions. brain. Front. Cell. Neurosci. 13, 363. doi:10.3389/fncel.2019.00363 Cerebellum 7, 589–594. doi:10.1007/s12311-008-0074-4 Muriuki, C., Bush, S. J., Salavati, M., McCulloch, M. E. B., Lisowski, Z. M., Agaba, M., Gratz, M. J., Stavrou, S., Kuhn, C., Hofmann, S., Hermelink, K., Heidegger, H., et al. et al. (2019). A mini-atlas of gene expression for the domestic goat (Capra hircus). Front. (2018). Dopamine synthesis and dopamine receptor expression are disturbed in Genet. 10, 1080. doi:10.3389/fgene.2019.01080 recurrent miscarriages. Endocr. Connect. 7, 727–738. doi:10.1530/EC-18-0126 Murphy, H. M., and Wildeman, C. H. (1986). Pup survival and development Guan, D., Landi, V., Luigi-Sierra, M. G., Delgado, J. V., Such, X., Castelló, A., et al. (2020). following hippocampal lesions in the female rat. Ohio J. Sci. 86, 89–94. 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 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. Harrison, S. J., Parrish, M., and Monaghan, A. P. (2008). Sall3 is required for the 1007/s00357-014-9161-z terminal maturation of olfactory glomerular interneurons. J. Comp. Neurology 507, 1780–1794. doi:10.1002/cne.21650 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. Harhay, G. P., Smith, T. P. L., Alexander, L. J., Haudenschild, C. D., Keele, J. W., Sci. Rep. 7, 15899. doi:10.1038/s41598-017-16019-8 Matukumalli, L. K., et al. (2010). An atlas of bovine gene expression reveals novel distinctive tissue characteristics and evidence for improving genome annotation. Nagayama, S., Homma, R., and Imamura, F. (2014). Neuronal organization of Genome Biol. 11, R102. doi:10.1186/gb-2010-11-10-r102 olfactory bulb circuits. Front. Neural Circuits 8, 98. doi:10.3389/fncir.2014.00098 Hébert, J. M., Lin, M., Partanen, J., Rossant, J., and McConnell, S. K. (2003). FGF Nevard, R. P., Pant, S. D., Broster, J. C., Norman, S. T., and Stephen, C. P. (2022). signaling through FGFR1 is required for olfactory bulb morphogenesis. Development Maternal behavior in beef cattle: The physiology, assessment and future directions-a 130, 1101–1111. doi:10.1242/dev.00334 review. Vet. Sci. 10, 10. doi:10.3390/vetsci10010010 Hoekzema, E., Barba-Müller, E., Pozzobon, C., Picado, M., Lucco, F., García-García, Numan, M., and Numan, M. J. (1991). Preoptic-brainstem connections and maternal D., et al. (2017). Pregnancy leads to long-lasting changes in human brain structure. Nat. behavior in rats. Behav. Neurosci. 105, 1013–1029. doi:10.1037/0735-7044.105.6.1013 Neurosci. 20, 287–296. doi:10.1038/nn.4458 Ortiz, C., Navarro, J. F., Jurek, A., Märtin, A., Lundeberg, J., and Meletis, K. (2020). Inestrosa, N. C., and Varela-Nallar, L. (2015). Wnt signalling in neuronal Molecular atlas of the adult mouse brain. Sci. Adv. 6, abb3446. doi:10.1126/sciadv. differentiation and development. Cell. Tissue Res. 359, 215–223. doi:10.1007/s00441- abb3446 014-1996-4 Patel, S., Rahmani, B., Gandhi, J., Seyam, O., Joshi, G., Reid, I., et al. (2020). Revisiting Jeong, J. K., Kim, J. G., Kim, H. R., Lee, T. H., Park, J. W., and Lee, B. J. (2017). A role the pineal gland: A review of calcification, masses, precocious puberty, and melatonin of central NELL2 in the regulation of feeding behavior in rats. Mol. Cells 40, 186–194. functions. Int. J. Neurosci. 130, 464–475. doi:10.1080/00207454.2019.1692838 doi:10.14348/molcells.2017.2278 Pawluski, J. L., Li, M., and Lonstein, J. S. (2019). Serotonin and motherhood: from Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M. (2016). KEGG molecules to mood. Front. Neuroendocrinol. 53, 100742. doi:10.1016/j.yfrne.2019. as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, 03.001 D457–D462. doi:10.1093/nar/gkv1070 Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript- Kavoi, B. M., and Jameela, H. (2011). Comparative morphometry of the olfactory level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. bulb, tract and stria in the human, dog and goat. Int. J. Morphol. 29, 939–946. doi:10. Nat. Protoc. 11, 1650–1667. doi:10.1038/nprot.2016.095 4067/s0717-95022011000300047 Pissonnier, D., Thiery, J. C., Fabre-Nys, C., Poindron, P., and Keverne, E. B. (1985). Keller, M., and Lévy, F. (2012). The main but not the accessory olfactory system is The importance of olfactory bulb noradrenalin for maternal recognition in sheep. involved in the processing of socially relevant chemosignals in ungulates. Front. Physiol. Behav. 35, 361–363. doi:10.1016/0031-9384(85)90309-9 Neuroanat. 6, 39. doi:10.3389/fnana.2012.00039 Poindron, P., Lévy, F., and Keller, M. (2007). Maternal responsiveness and maternal Kim, D., Langmead, B., and Salzberg, S. L. (2015). Hisat: A fast spliced aligner with selectivity in domestic sheep and goats: The two facets of maternal attachment. Dev. low memory requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317 Psychobiol. 49, 54–70. doi:10.1002/dev.20192 Frontiers in Genetics 15 frontiersin.org Luigi-Sierra et al. 10.3389/fgene.2023.1114749 Prosser, H. M., Bradley, A., and Caldwell, M. A. (2007). Olfactory bulb hypoplasia in inflammatory state in mouse hippocampus. Front. Cell. Neurosci. 12, 412. doi:10.3389/ Prokr2 null mice stems from defective neuronal progenitor migration and fncel.2018.00412 differentiation. Eur. J. Neurosci. 26, 3339–3344. doi:10.1111/j.1460-9568.2007.05958.x Sullivan, R. M., Wilson, D. A., Ravel, N., and Mouly, A. M. (2015). Olfactory memory R Core team (2015). R: A language and environment for statistical computing. R networks: From emotional learning to social behaviors. Front. Behav. Neurosci. 9, 36. foundation for statistical computing. Available at: http://www.R-project.org/. doi:10.3389/fnbeh.2015.00036 Ragancokova, D., Rocca, E., Oonk, A. M. M., Schulz, H., Rohde, E., Bednarsch, J., et al. Sultana, A., Nakaya, N., Dong, L., Abu-Asab, M., Qian, H., and Tomarev, S. I. (2014). (2014). TSHZ1-dependent gene regulation is essential for olfactory bulb development Deletion of olfactomedin 2 induces changes in the AMPA receptor complex and impairs and olfaction. J. Clin. Investigation 132, 1214–1227. doi:10.1172/JCI72466 visual, olfactory, and motor functions in mice. Exp. Neurol. 261, 802–811. doi:10.1016/j. expneurol.2014.09.002 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 Sunkin, S. M., Ng, L., Lau, C., Dolbeare, T., Gilbert, T. L., Thompson, C. L., et al. and the postpartum period. G3 Genes., Genomes, Genet. 6, 221–233. doi:10.1534/g3.115. (2013). Allen brain atlas: An integrated spatio-temporal portal for exploring the central 020982 nervous system. Nucleic Acid. Res. 41, D996–D1008. doi:10.1093/nar/gks1042 Reichmann, F., Rimmer, N., Tilley, C. A., Dalla Vecchia, E., Pinion, J., Al Oustah, A., Tanegashima, K., Okamoto, S., Nakayama, Y., Taya, C., Shitara, H., Ishii, R., et al. et al. (2020). The zebrafish histamine H3 receptor modulates aggression, neural activity (2010). CXCL14 deficiency in mice attenuates obesity and inhibits feeding behavior in a and forebrain functional connectivity. Acta Physiol. 230, e13543. doi:10.1111/apha. novel environment. PLoS One 5, e10321. doi:10.1371/journal.pone.0010321 13543 Thomas, S. A., and Palmiter, R. D. (1997). Impaired maternal behavior in mice Reinert, S., Hübener, M., Bonhoeffer, T., and Goltstein, P. M. (2021). Mouse lacking norepinephrine and epinephrine. Cell. 91, 583–592. doi:10.1016/S0092- prefrontal cortex represents learned rules for categorization. Nature 593, 411–417. 8674(00)80446-8 doi:10.1038/s41586-021-03452-z Tomaz, B., Aljaz, M., Ales, J., and Gregor, M. (2014). Deletion of the prion gene Prnp Renault, V. M., Rafalski, V. A., Morgan, A. A., Salih, D. A., Brett, J. O., Webb, A. E., affects offensive aggression in mice. Behav. Brain Res. 266, 216–221. doi:10.1016/j.bbr. et al. (2009). FoxO3 regulates neural stem cell homeostasis. Cell Stem Cell 5 (5), 527–539. 2014.03.006 doi:10.1016/j.stem.2009.09.014 Venkatraman, A., Edlow, B. L., and Immordino-Yang, M. H. (2017). The brainstem in Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: A Bioconductor emotion: A review. Front. Neuroanat. 11, 15. doi:10.3389/fnana.2017.00015 package for differential expression analysis of digital gene expression data. Wang, D., Eraslan, B., Wieland, T., Hallström, B., Hopf, T., Zolg, D. P., et al. (2019). A Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616 deep proteome and transcriptome abundance atlas of 29 healthy human tissues. Mol. Rolls, A., Schori, H., London, A., and Schwartz, M. (2008). Decrease in hippocampal Syst. Biol. 15, e8503. doi:10.15252/MSB.20188503 neurogenesis during pregnancy: A link to immunity. Mol. Psychiatry 13, 468–469. Wang, X., Michaelis, M. L., and Michaelis, E. K. (2010). Functional genomics of brain doi:10.1038/sj.mp.4002126 aging and alzheimer’s disease: Focus on selective neuronal vulnerability. Curr. Genomics Rowe, F. A., and Edwards, D. A. (1971). Olfactory bulb removal: Influences on the 11, 618–633. doi:10.2174/138920210793360943 aggressive behaviors of male mice. Physiol. Behav. 7, 889–892. doi:10.1016/0031- Yanai, I., Benjamin, H., Shmoish, M., Chalifa-Caspi, V., Shklar, M., Ophir, R., et al. 9384(71)90059-X (2005). Genome-wide midrange transcription profiles reveal expression level Santos Fontanez, S. E., and De Jesus, O. (2021). Neurohypophysis. Treasure Island, FL: relationships in human tissue specification. Bioinformatics 21 (5), 650–659. doi:10. StatPearls Publishing. 1093/bioinformatics/bti042 Schurch, N. J., Schofield, P., Gierliński, M., Cole, C., Sherstnev, A., Singh, V., et al. Yoshihara, S., Omichi, K., Yanazawa, M., Kitamura, K., and Yoshihara, Y. (2005). Arx (2016). How many biological replicates are needed in an RNA-seq experiment and homeobox gene is essential for development of mouse olfactory system. Development which differential expression tool should you use? RNA 22, 839–851. doi:10.1261/rna. 132 (4), 751–762. doi:10.1242/dev.01619 053959.115 Yu, G. Z., Kaba, H., Okutani, F., Takahashi, S., and Higuchi, T. (1996). The olfactory Sjöstedt, E., Zhong, W., Fagerberg, L., Karlsson, M., Mitsios, N., Adori, C., et al. bulb: a critical site of action for oxytocin in the induction of maternal behaviour in the (2020). An atlas of the protein-coding genes in the human, pig, and mouse brain. Science rat. Neuroscience 72 (4), 1083–1088. doi:10.1016/0306-4522(95)00600-1 367, aay5947. doi:10.1126/science.aay5947 Yuan, T. F. (2008). BDNF signaling during olfactory bulb neurogenesis. J. Neurosci. 28 Son, H., Banasr, M., Choi, M., Chae, S. Y., Licznerski, P., Lee, B., et al. (2012). Neuritin (20), 5139–5140. doi:10.1523/JNEUROSCI.1327-08.2008 produces antidepressant actions and blocks the neuronal and behavioral deficits caused Zapala, M. A., Hovatta, I., Ellison, J. A., Wodicka, L., Del Rio, J. A., Tennant, R., et al. by chronic stress. Proc. Natl. Acad. Sci. U. S. A. 109, 11378–11383. doi:10.1073/pnas. (2005). Adult mouse brain gene expression patterns bear an embryologic imprint. Proc. 1201191109 Natl. Acad. Sci. U. S. A. 102 (29), 10357–10362. doi:10.1073/pnas.0503357102 Stocker, A. M., and Chenn, A. (2015). The role of adherens junctions in the Zemla, R., and Basu, J. (2017). Hippocampal function in rodents. Curr. Opin. developing neocortex. Cell. Adh Migr. 9, 167–174. doi:10.1080/19336918.2015. Neurobiol. 43, 187–197. doi:10.1016/j.conb.2017.04.005 1027478 Zhou, S., and Zhou, J. (2014). Neuritin, a neurotrophic factor in nervous system Su, W. J., Zhang, T., Jiang, C. L., and Wang, W. (2018). Clemastine alleviates physiology. Curr. Med. Chem. 21 (10), 1212–1219. doi:10.2174/ depressive-like behavior through reversing the imbalance of microglia-related pro- 0929867321666131218093327 Frontiers in Genetics 16 frontiersin.org