Transcriptomic changes associated with maternal care in the brain of mouthbrooding cichlid Astatotilapia burtoni reflect adaptation to self-induced metabolic stress

Parental behavior has evolved multiple times in animals. Although it is a costly reproductive strategy for parents, it increases fitness through higher offspring survival rate. These parental costs come in the form of time and energy, as time spent caring for and protecting offspring must be balanced against competing behavioral interests such as feeding and mating (Numan and Woodside, 2010; Olazábal and Young, 2006). For most species, the diverse set of activities that compose ‘feeding’ and ‘parental care’ must be co-regulated through partially overlapping neuroendocrine circuits, and evidence indicates substantial cross-talk between these circuits in at least some species (e.g. Olszewski et al., 2010; Crossin et al., 2012; Shahjahan et al., 2014). While some species increase feeding during parental care to compensate for the added demands, others inhibit feeding behavior to favor parental care, a tradeoff that comes at the expense of parental condition (Fischer and O'Connell, 2017; O'Rourke and Renn, 2015).

The neural circuits and molecular mechanisms regulating parental care and feeding behaviors are well studied, primarily in mammalian systems. Ongoing work provides inference regarding homology across vertebrates at both the circuit (Goodson, 2005; Goodson and Kingsbury, 2013; O'Connell and Hofmann, 2011) and molecular signaling (Volkoff, 2016) levels. The close connection between reproduction and energy balance is well documented in mammals with several nuclei of the ‘maternal brain’ receiving either direct or indirect input from the arcuate nucleus feeding regulatory neurons (Schneider, 2004). Within the arcuate nucleus of the hypothalamus, two specific populations of neurons, one producing orexigenic neuropeptide Y (NPY) and agouti related peptide (AgRP), and the other producing anorexigenic pro-opiomelanotcortin (POMC) and cocaine- and amphetamine-regulated transcript (CART), are directly responsive to the peripheral signals of satiety, leptin (Cone, 2005) and ghrelin (Cowley et al., 2003).

Mouthbrooding, a form of parental care seen in several species of fish in which the parent holds the developing eggs and fry in the buccal cavity, represents an extreme example of food restriction and parental investment. In an ultimate sense, mouthbrooding is an adaptation in response to predation pressure (Kuwamura, 1997) and has evolved independently in cardinalfishes, gobies, betas (Rüber et al., 2004) and several clades of cichlid fishes (Barlow, 2000; Fryer and Iles, 1972), where it is thought to have promoted rapid speciation (Salzburger et al., 2005). In a proximate sense, mouthbrooding severely limits food intake, such that most species completely abstain from food intake while others obtain food only for provisioning the fry (Yanagisawa et al., 1996). This brooding phase is plastically adjusted in response to predation threat (Taborsky and Foerster, 2004), noise-induced perturbations (Butler and Maruska, 2021) and parental energy reserve as well as the quality and availability of mates (Iwao et al., 2004).

In the African cichlid fish Astatotilapia burtoni, it is the female alone that provides parental care in the form of mouthbrooding, abstaining from food intake for 2–3 weeks. This leads to a dramatic reduction in body mass (Grone et al., 2012; Renn et al., 2009), and reduced gonadosomatic index associated with lowered estradiol and 11-KT (a fish-specific androgen) (Maruska et al., 2020b). Unlike a state of involuntary starvation, feeding motivation is decreased during mouthbrooding even when the developing eggs are removed from the buccal cavity (Mrowka, 1984), and thus the female does not anticipate caloric intake. The capacity to withstand and recover from prolonged periods of food deprivation requires behavioral, physiological, morphological and molecular adaptations to conserve energy and maintain immune function (Secor and Carey, 2016). We have previously identified an enhanced reduction in gut epithelial turnover during mouthbrooding as one such potentially adaptive mechanism for energy conservation (Faber-Hammond et al., 2019). While the full extent of metabolic changes during mouthbrooding is not known, it may share similarities with changes seen in animals that engage in hibernation or other related extreme adjustments when caloric intake is not anticipated. During hibernation, shifts in mitochondrial ATP production and calcium uptake may improve energy efficiency (Barger et al., 2003), neuroprotective mechanisms may mitigate oxidative damage (Yin et al., 2016), and changes to brain water homeostasis can also impact neuronal signal transduction (Osborne and Hashimoto, 2008).

The expression of genes encoding some neuropeptides and their receptors is known to differ between brooding, non-brooding and starved females, suggesting that different regulatory mechanisms control this special case of feeding inhibition (Grone et al., 2012). Furthermore, recent work examining neuronal activation across the brain during mouthbrooding in A. burtoni provides evidence for a large degree of shared circuitry involved in the regulation of maternal care, food intake and energy balance with nuanced differences in specific brain regions (Maruska et al., 2020a). Such work predicts global changes in gene expression across the brain with specifically regulated genes in some regions.

In the current study, we analyzed whole-brain transcriptomes of female A. burtoni from multiple time points across the reproductive cycle. The time points selected represent dramatically different stages in terms of nutrition state, energy requirements and parental care behaviors. We conducted differential expression analyses and parallel enrichment analyses for functional gene sets to identify patterns of upregulation or downregulation at either of two specific mouthbrooding stages relative to each other and to two post-release time points, one of which includes intensive maternal care. Genes that are differentially regulated across the maternal stages (both mouthbrooding and post-release) may underlie the specific behaviors performed at these time points, while genes showing trends that correlate with decreasing body condition may reflect energy reserves. In addition to providing valuable insight into the evolution of parental behavior, our results reveal previously unknown changes in the maternal cichlid brain that could have clinical impacts in the study of brain hypoxic ischemia, inflammation and neuroprotection. Our holistic approach using whole-brain samples aims to identify novel mechanisms beyond a candidate gene approach and capture global changes of a metabolically active tissue.

Astatotilapia burtoni (Günther 1894) used in this experiment were derived from the same lineage as was sequenced for the reference genome assembly (GCA_000239415). Females of reproductive age were maintained in 30 gallon (∼114 l) mixed-sex stock tanks with 2–4 males and 8–12 females. Tanks were kept on a 12 h:12 h day:night cycle with 0.5 h dawn and dusk settings. The salinity and pH levels of tanks were matched to the average levels of Lake Tanganyika. Tanks were checked daily, and when brooding females were captured and fry development stage was confirmed. Females were assigned to one of four treatment groups: 2 days mouthbrooding (B02), 14 days mouthbrooding (B14), 2 days post-release of fry (R02) and 14 days post-release (R14) (Fig. 1). To avoid confounding date with treatment groups, females were sequentially assigned one to each group, with the exceptions that only females confirmed to be holding stage 0 fry were assigned to the B02 group. Brooding females were placed into a 5 gallon (∼19 l) tank with a terracotta pot for shelter and gravel. A single male stimulus fish (∼3 g±1 g) was housed in the same tank, retained by a cylindrical mesh barrier (10 cm diameter) to allow visual and olfactory contact. The male stimulus fish was fed standard cichlid pellets daily within the confined cylinder. Because brooding females are not expected to eat, the addition of food would have fouled the tank. Therefore, females had visual and olfactory cues regarding food availability without actual access to it. Once females released the fry, they were fed standard fish flakes daily. All protocols complied with Renn lab IACUC approval #01-2017.

Females were monitored daily for brooding or release of fry. All fry were removed 2 days after the initial release date. Any female that cannibalized her fry prior to 2 days post-release was reassigned to the R14 group, so that all R02 fish were actively caring for free-swimming fry. At B02, B14, R02 and R14 time points, females were anesthetized in MS-222, weighed, measured for standard length and rapidly decapitated. Whole brains were removed and preserved in 1 ml RNAlater (Ambion) and stored at −80°C. All dissections were done between 16:00 and 18:00 h. While whole-brain gene expression profiling avoids the technical variation introduced by brain punch techniques, this approach is not without caveats and limitations. Primarily, region-specific regulation and genes expressed in a limited number of cells are less likely to be detected. Certain genes are known to be differentially regulated in brain region-specific patterns that are lost with whole-brain expression profiling. For example, Butler et al. (2020) studied the role of galanin in maternal care and infanticide in A. burtoni, and found that expression in preoptic area (POA) neurons is associated with maternal care while expression in nucleus lateralis tuberis (NLT) neurons is associated with feeding. While this finer scale resolution is lost in analysis of whole brains, our method does capture a snapshot of the broad changes in the brain throughout the reproductive cycle associated with this metabolically demanding behavior.

Brains were removed from RNAlater and placed in 200 μl of 1-Thioglycerol/Homogenization Solution from the MAxwell 16 LEV Simply RNA Tissue Kit (Promega AS1280). Samples were homogenized using a hand-held pestle motor and plastic pestle in a microcentrifuge tube. Homogenized samples were transferred to the MAxwell 16 LEV Simply RNA Tissue Kit cartridge and processed in the Maxwell 16 Instrument (Promega AS1000) using the High Strength LEV Magnetic Rod and Plunger Bar Adaptor (Promega SP1070). Standard DNAse treatment was included and samples were eluted in 50 μl H₂O including RNAseOUT (Invitrogen). RNA quality was initially checked by electrophoresis and quantity by NanoDrop 1000 (Thermo Scientific). RNA (1–2 μg per sample) was sent to Novogene (Sacramento, CA, USA) for library preparation and paired-end 150 bp sequencing (Illumina Platform PE150).

Paired-end RNA-seq reads were preprocessed with TrimGalore v0.6.4 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) for adaptor trimming and FastQC v0.11.3 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality filtering. Remaining paired reads were aligned to the A. burtoni reference genome (GCF_000239415.1) with STAR v2.5.3 (Dobin et al., 2013), using the associated RefSeq GTF file to guide the alignment. Aligned reads were assembled and quantified using StringTie v1.3.4 (Pertea et al., 2015), and count tables were constructed as input for differential expression analysis. All methods, results and discussion of a WGCNA analysis (v1.69: Langfelder and Horvath, 2008) that was conducted in parallel are available from GitHub (github.com/jfaberha/cichlid_transcriptome/tree/main/lab_stock_analysis).

Data normalization and subsequent differential expression analysis were performed with DESeq2 v1.20.0 (Love et al., 2014) in R-studio v1.1.383 (http://www.rstudio.com/). A likelihood ratio test (LRT) was performed to generate a list of genes that showed significant expression differences between any of the four sampled time points, followed by pairwise binomial Wald tests to discover which pair of time points differ. As a supplementary analysis, brooding (B02, B14) and release (R02, R14) stages were pooled and a pairwise binomial Wald test was performed to discover candidate genes that differentiate by associated parental behaviors and response to food intake. Genes showing significant differential expression from both tests are reported at a false discovery rate (FDR)<0.1 and for exploratory purposes we report genes trending towards differential expression at an unadjusted P<0.05 (results reported in Table S1). This less stringent threshold provided a sufficient number of genes for downstream functional analyses. LRT differentially expressed genes (DE genes) were tested for gene ontology (GO) term enrichment with BiNGO v3.0.4 (Maere et al., 2005), implemented in Cytoscape v3.7.2 (Shannon et al., 2003). While A. burtoni transcripts were assigned names based on their NCBI RefSeq GTF annotation, they were reannotated by predicted protein sequence to obtain a comprehensive list of GO terms (Table S2). In order of priority, GO annotations were preferentially used from the Swiss-Prot, TrEMBL (The UniProt Consortium, 2019) and EggNOG (Huerta-Cepas et al., 2019) protein databases (accessed February 2018).

We performed hierarchical clustering analysis on subsets of genes showing significant differential expression and trending towards differential expression that address original and novel exploratory hypotheses with the goal of discovering co-expression and possible interaction of candidate genes. For all clustering analyses, we used variance stabilizing transformation (VST)-normalized count tables generated by DESeq2, and constructed heatmaps using the made4 v1.56.0 R-package (Culhane et al., 2005). This package utilizes the correlation similarity metric and average linkage clustering for expression data, which allows for clustering of expression patterns regardless of the overall magnitude of expression. Clustering was performed on individuals to generate expression-based dendrograms, while the associated heatmaps display time-point medians for easier interpretation of expression data. Bootstrap values for dendrograms were generated with the PVclust v2.2-0 R-package (Suzuki and Shimodaira, 2006) following 10k clustering iterations.

In addition to the detection of individual DE genes, we aimed to identify changes in functional gene categories across time points; therefore, we ran Gene Set Enrichment Analysis (GSEA v4.0.3) for sets of genes representing GO categories (Subramanian et al., 2005). Pairwise comparisons were done with gene-level VST-normalized counts between each of the four time points in both directions, resulting in 12 possible comparisons. Only six comparisons are possible between four time points; however, reciprocal comparisons (e.g. B02 versus B14 and B14 versus B02) were performed because GSEA results slightly differ based on the direction of the test, due to significance being determined by ranking of sets, and the union is reported. For each GSEA comparison, 1000 permutations of gene set randomization were performed because there were fewer than 10 samples per group. Sets with FDR<0.05 were considered significantly enriched, and the GO categories had to meet this threshold in either of two reciprocal tests per time-point comparison. GSEA documentation recommends an FDR cutoff of 0.25 for determining significance in exploratory analyses because of false negatives being less desirable than false positives in such analyses; however, we chose the more restrictive alpha threshold of 0.05 to limit the number of significant GO sets. GSEA tests were followed by leading-edge analysis to determine both the direction of change and the overlap between enriched gene sets.

Significantly enriched GO gene sets are not fully independent in the gene accessions they contain, as single genes often have a number of associated GO terms. To account for this, we took steps to reduce complexity in the list of significant GO terms by merging terms with minimum 50% overlap of member genes into GO families. Overlap coefficients (CO) were calculated as the proportion of shared genes between pairs of terms divided by the smaller of the two gene sets. The inverse values of these coefficients (1−CO) were converted into a distance matrix for all significant GO terms from GSEA, then used to build a neighbor joining tree with PHYLIP v3.698 (http://evolution.genetics.washington.edu/phylip.html). TreeClust v1.00 (Balaban et al., 2019) was used to cluster nodes that had a maximum distance of 0.5 from each other, equivalent to a minimum 50% overlap coefficient between GO terms. Merged GO families were named using the most significant (lowest FDR) term from all GSEA time-point pairwise comparisons.

For hierarchical clustering of significant GO terms and families, GSEA was re-run to obtain normalized enrichment scores (NES) of each of four time points versus all other samples, with positive values representing upregulation and negative values representing downregulation of entire gene sets in an analog to log₂-fold change of individual genes. Hierarchical clustering was performed using the Ward.D2 clustering and Euclidean distance methods implemented by the Heatmap tool from the ComplexHeatmaps R package v2.6.2 (Gu et al., 2016; Ward, 1963).

The time points selected for the current study represent dramatically different states with regard to behavioral, feeding and metabolic activity. We used the R14 stage as a baseline in which there is active feeding behavior, adequate energy reserves and normal metabolic activity but no maternal care behavior. The early brooding stage, B02, is characterized by modest buccal churning as maternal care, no feeding behavior, adequate energy reserves and potentially reduced metabolic activity to conserve energy; and while the later brooding stage, B14, shares the lack of feeding behavior, it is also characterized by more intense buccal-based maternal care and dramatically reduced energy reserves, which likely necessitates dramatic regulation of metabolic activity. The early release stage, R02, shares this increased maternal care but of a very different nature in that it includes aggressive territory defense, and while feeding behavior has resumed, energy stores remain depleted, again potentially requiring substantial metabolic regulation.

Twenty-two samples passed QC filters for analysis, with 5 samples representing B02 and R02 time points and 6 samples representing B14 and R14 time points. From ∼610 million total high-quality reads (mean per sample ∼27.7 million±7.8 million), 84.6% reads were aligned to the A. burtoni genome. These A. burtoni brain transcriptomes yielded hits to 25,572 out of the 44,653 RefSeq annotated genes (57.3%), which accounted for 92.9% of final gene counts. Conversely, 7.1% of gene counts represented novel unannotated transcripts within the genome assembly, and upon GO annotation for the entire transcriptome we found only a small portion (17.7%) of these novel transcripts to be putative protein coding genes, many of which appeared to code for transposable element machinery. Therefore, only annotated NCBI RefSeq genes in the assembly were used for downstream analyses.

We found 128 genes to be differentially expressed across the four reproductive time-points through DESeq2 LRT analysis with an FDR<0.1, making up 0.5% of RefSeq annotated genes expressed in the brain. For exploratory purposes, we also examined a secondary list of 1643 genes (1515 additional genes) trending toward differential expression at an unadjusted P<0.05 (Table S1). With post-hoc pairwise comparisons, we explored patterns among the less stringent set of genes trending toward differential expression (Table 1). The highest number of DE genes, 105, occurred between late brooding B14 and early release R02 stages in the brooding cycle (44 B14 biased; 61 R14 biased) when maternal behavior and feeding are dramatically changing. Conversely, the fewest number of DE genes, 20, was detected in the contrast between early brooding B02 and early release R02 (11 B02 biased and 9 R02 biased). Among the full set of DE genes identified for a given brooding stage (both upregulated and downregulated), only a small percentage were biased toward that brooding stage relative to all three other stages (B02 2.56%, B14 0.56%, R02 1.36%, R14 0%), meaning that there were few stage-specific gene expression levels. The same dynamic pattern was seen among the larger list of genes trending toward differential expression, with a greater percentage of differential expression being stage specific (2.5–6.8%). While the two brooding stages are very different in terms of duration of starvation and other physiological factors, they share a mouthbrooding maternal component that differs from either of the release stages; therefore, we also combined the two brooding stages and the two release stages and performed a pairwise pooled brooding versus pooled release binomial Wald test. This analysis revealed 31 DE genes at FDR<0.1, 21 of which were significant DE genes, with the remaining 10 trending towards significance based on LRT analysis. An additional 1124 genes trended towards differential expression in the brooding versus release analysis at an unadjusted P<0.05 (Table S1). Functionally interesting candidate genes appearing in both LRT and brooding versus release differential expression lists are highlighted in expression heatmaps while the following discussion primarily reports on the LRT DE gene list across all four time points.

The DE genes discussed below include several that are classically associated with parental care such as the teleost ortholog to oxytocin, isotocin neurophysin (itnp), and others that have more recently emerged as potential candidates in parental care such as neurotensin (nts). We found strong correlation of expression that clustered genes involved in parental care with others classically associated with feeding, such as galanin (gal). As expected, given the voluntary starvation phase, many of the known orexigenic genes and anorexigenic genes varied dramatically across the time points sampled, and we saw compensatory response genes that constitute metabolic regulatory mechanisms, suggesting hypoxic conditions and activation of neuroprotective genes.

GO enrichment for the 128 genes determined to be differentially expressed across stages by LRT analysis (regardless of comparison or directionality) yielded only six significant terms (Table 2). However, none of the enriched GO terms were typically associated with maternal or voluntary fasting behavior in our searches of current literature. Instead, five were related to oxygen transport, and the sixth was ‘negative regulation of protein phosphatase type 2B activity’. The oxygen transport GO categories were detected based on differential expression of four globin genes, and ‘negative regulation of protein phosphatase type 2B activity’ was represented by two DE genes, repressor of calcipressin1-like (rcan1l) and peptidyl-prolyl cis-trans isomerase fkbp1-like. Although the number of significantly differentially expressed genes in both of these groups was low and significance levels may be inflated as a result, there were only 25 and three non-significant genes in these GO categories, respectively, out of the 22,924 annotated genes in the brain transcriptome. In addition, these GO groups are represented by genes that show dramatic changes across time points, with globin gene hemoglobin subunit beta-A-like as the third-most significant DE gene (FDR=4.68E−6) and rcan1l as the second-most significant gene (FDR=8.75E−14) in the dataset. When running GO enrichment analysis on the 1643 genes that trend toward differential expression for an unadjusted P<0.05, we found a list of six different enriched terms (Table 2), including ‘hormone activity’ as the top GO term, represented by 36 accessions. These findings suggest that observed behavioral changes may not require vast changes in expression of entire gene pathways, but rather they could result from dramatic changes in a small number of genes or minor expression changes in larger gene networks.

We used GSEA as a parallel approach to discover enriched pathways that differ between time points, which allowed us to detect differences in aggregate signal for groups of genes related in function even when individual member genes failed to meet significance thresholds in differential expression analysis. In fact, there were a number of enriched sets that failed to contain any identified DE genes. Because there was partial overlap of gene accessions between sets, we merged significant GO sets into families, and after this data reduction step we found a total of 23 enriched GO families composed of 107 GO sets and an additional 32 unmerged GO terms for a total of 55 independent sets (Fig. 2; Table S3). These enriched GO terms include several that are potentially involved in modifying behavior, including several related to neuronal growth/death and neuropeptide receptor activity. Other notable enriched families and sets are related to ion channel regulation, oxygen transport, energy metabolism and inflammatory response. Following NES clustering for these sets, the two dominant patterns that emerged were both highly dynamic, alternating between upregulation and downregulation across adjacent time points. The majority of independent GO sets and families (38/55) showed these patterns with 19 sets upregulated in B02 and R02 relative to B14 and R14, and another 19 sets showing the opposite pattern. Of the remaining 17/55 sets, 14 showed more gradual expression changes across the approximately cyclical time points, with adjacent pairs of time points having either higher or lower expression than the two distal time points; the other three sets were families composed of GO terms with inconsistent NES patterns between time points.

Among cichlids, mouthbrooding is thought to have arisen 10–14 times (Goodwin et al., 1998) and once at the base of the remarkably diverse haplochromine clade (DeVos and Seegers, 1998). Thus, the mouthbrooding haplochromine cichlid species A. burtoni, an important model system for the physiological and genomics studies of complex social behavior (Maruska and Fernald, 2013; Brawand et al., 2014), offers a useful perspective to identify mechanisms underlying mouthbrooding that may have contributed to diversification in this cichlid clade and possibly also in others.

Across the mouth-brooding reproductive cycle, female A. burtoni undergo two important transitions: an initial transition to starvation during the obligatory mouthbrooding phase and a second transition to increased aggression and maternal defense following release of fry (Renn et al., 2009). The care given to the developing offspring during the mouthbrooding phase includes metabolically demanding buccal churning to aerate the fry. As the obligatory phase of mouthbrooding continues, this metabolic demand is coupled with a dramatic reduction in body condition as a result of prolonged voluntary starvation. Following the release of fry from the buccal cavity, the metabolically demanding maternal care is coupled with the resumption of feeding behavior. From the perspective of parental care, while considering fundamental distinctions, these two transitions can be compared to initiation of mammalian pregnancy and parturition, or avian incubation of eggs, which involves an early decrease in activity for incubation followed by increased activity for provisioning of chicks. From the perspective of feeding behavior, these transitions are dramatically different from mammalian transitions of pregnancy and parturition given the complete lack of food intake and dramatic reduction in body condition during brooding.

The time points selected for the current study highlight these important transitions that encompass behavioral, feeding and metabolic changes. The late release time point, R14, provides a rough baseline during which extensive maternal care and aggression have subsided. The early brooding phase, B02, describes the initial transition to the maternal brooding stage, while late brooding stage, B14, describes an individual that has undergone substantial loss of body condition and yet retains brood care motivation while sustaining suppression of feeding behavior. The early release stage, R02, describes the social stages of maternal care as the female interacts with offspring, shows heightened aggression and has resumed feeding behavior. That we found the greatest number of differentially regulated genes between B14 and R02 likely reflects the dramatic differences with regard to both maternal behavior and body condition.

The results discussed below include both known and novel candidates for the regulation of maternal care, feeding behavior and metabolism. In addition, oxygen transport pathways were found to be dramatically downregulated throughout mouthbrooding while neuroprotective pathways and immune response were upregulated, likely to mitigate the negative consequences of metabolic depression brought about by fasting. These results offer new insights into the evolution of mouthbrooding behavior. The individual genes identified by differential expression analysis are corroborated by GO term and GSEA analysis implicating full pathways in these transitions.

Some of the GO and GSEA analysis results (e.g. ‘equilibrio-perception’ and ‘sensory perception of light’) indicate biological processes known to be controlled by specialized brain regions (Kasumyan, 2004; Rosa Salva et al., 2014), thus pointing to the need for future studies focused on more localized gene expression in the brain using micropunch or single-cell approaches. With regard to maternal mouthbrooding, it is known that activation of galanin neuropeptide-expressing neurons in the POA plays a role in maternal care and infanticide while activation of neurons expressing the same peptide in the NLT brain region is associated with feeding regulation (Butler et al., 2020). To the extent that neuronal activation correlates with gene regulation, when not studied in isolation, gene expression profiles from these regions controlling different aspects of behavior will be lost.

Numerous studies show a positive correlation between parental care and expression of the teleost ortholog of oxytocin, isotocin neurophysin (itnp) (DeAngelis et al., 2017; Keverne and Kendrick, 1992; O'Connell et al., 2012; Pedersen et al., 1982), making it a candidate gene for regulation in the current study. When compared with the R14 baseline, late stage brooding females in B14 showed significantly lower expression of itnp (FDR=7.32E−2) (Fig. 3). This negative, rather than positive, correlation between itnp and care may be analogous to the reduced expression seen in female stickleback following spawning (Kulczykowska and Kleszczyńska, 2014), associated with oocyte maturation and ovulation, or may simply reflect the low social engagement during this phase and suggest it is not a primary regulator of maternal mouthbrooding behavior in A. burtoni.

The neuropeptide neurotensin (nts) has emerged as a candidate for involvement in parental care, showing significantly elevated expression (FDR=4.1E−2) throughout the mouthbrooding stages, with a sharp decrease after release of fry (Figs 3 and 4). Expression of nts is positively and negatively associated with different forms of parental behavior in different taxa. In mice, injected nts inversely modulates aggressive behavior, potentially through interaction with its receptor neurotensin receptor 1 (ntsr1; Fig. 4) (Gammie et al., 2009), and causes increased neuronal activity in several brain regions including the lateral septal nucleus (LS) and bed nucleus of the stria terminalis (BNST), which are involved in social behavior and anxiety (Clauss et al., 2019; Deng et al., 2019). Localized expression studies in mice suggest a complex relationship with elevated nts expression in postpartum females in the medial preoptic area (MPOA) and lower expression in the LS (Driessen et al., 2014). Meanwhile, data from birds (domesticated Japanese quails) show higher nts expression in the female hypothalamus associated with parental behavior including reduced aggression towards juveniles (Lopes and de Bruijn, 2021).

The involvement of the nts pathway, beyond differential expression of the neuropeptide itself, is corroborated by our GSEA enrichment results that reveal the GO set ‘G protein-coupled neurotensin receptor activity’ to be increased in B14 compared with B02 fish (FDR=1.18E−2), suggesting receptor expression may be subject to feedback from early expression of the neuropeptide (Fig. 2B; Table S3). This GO set represented many genes of low expression level that failed to meet significance thresholds independently, though three of the 17 nts receptors or receptor-interacting ligands did reach an unadjusted threshold (P<0.05) (Fig. 4). Two of these genes, sortilin 1/neurotensin receptor 3 (sort1) and sortilin-related VPS10 domain containing receptor 3-like (sorcs3l) follow a pattern similar to that of nts itself, although their highest expression was in B14 instead of B02. Spatial co-expression of nts and sort1 is also well documented in mammals, with both genes upregulated in the MPOA while being downregulated in the LS brain region of postpartum female mice relative to virgins (Driessen et al., 2014). These region-specific changes in mice further highlight the need for more targeted analysis in the future in A. burtoni in order to disentangle the interactive role of nts and its receptors in the regulation of maternal behavior. In our study, a third nts receptor, neuromedin U receptor 1 (nmur1) followed an inverted pattern relative to nts with its highest expression at release, R02. This gene's functional relationship to nts signaling is not well studied as its primary affinity is for the pleiotropic peptide neuromedin U (nmu), which has a brain expression profile that parallels that of nmur1, with highest expression in R02, suggesting they are co-expressed in the maternal brain (Martinez and O'Driscoll, 2015).

Among the differentially expressed neuropeptides, hormones and receptors shown in Fig. 4, hierarchical clustering revealed that expression of nts was most strongly correlated with gal, with both expressed at significantly higher levels in brooding time points relative to release time points (Fig. 4; Table S1). Studies of fish, mammals and amphibians have found gal to be associated with parental behavior (Butler et al., 2020; Fischer et al., 2019; Wu et al., 2014). Specific MPOA galanin neurons are activated during parental behavior in male and female mice, and optogenetic activation induces parental pup-grooming while genetic ablation suppresses parental behavior (Wu et al., 2014). In A. burtoni, galanin neuron activity is associated with maternal care in the POA and feeding in the arcuate nucleus (Butler et al., 2020). Correlated increased expression of gal and nts is also reported in POA of Dendrobates tinctorius males involved in tadpole transport (Fischer et al., 2019). That expression of nts and gal wanes in A. burtoni as maternal aggression increases, just prior to and following the release of fry (Renn et al., 2009), suggests a conserved role in recognition and protection of offspring rather than overt defense. Additional parallels with paternal care in frogs include the upregulation of myelin basic protein-like (mbpl) (unadjusted P=1.51E−2) and downregulation of progesterone receptor (pgr) (unadjusted P=8.83E−3) and prolactin receptor b (prlrb) (unadjusted P=8.80E−3) at the two mouthbrooding time points in A. burtoni and in parental D. tinctorius males, yet prolactin (prl) itself was not differentially expressed in either study (Renn et al., 2009; Fischer et al., 2019). Several genes associated with parental behavior in both studies show opposite regulation, including neuroligin-3 (nlgn3) and secretogranin II (scg2), which are downregulated in D. tinctorius parental males and upregulated in A. burtoni mouthbrooding females, while androgen receptor (ar) shows the opposite pattern in each species. Assuming these normalized read counts correlate with translated peptide levels, the opposing patterns of differential expression between A. burtoni and D. tinctorius may reflect the species-specific demands of parental behavior or result from parental behavior being expressed in the opposite sex.

While inconsistent patterns of regulation between D. tinctorius and A. burtoni may result from the different techniques (brain punch versus whole brain) or taxa-specific parental strategies, the similarity in suites of genes at play is of note. In the frog POA, cocaine- and amphetamine-regulated transcript (cart) is upregulated in parental males while corticotropin-releasing hormone binding protein (crhbp) is downregulated and in the medial pallium (MP), corticotropin-releasing factor receptor 2 (crhr2) as well as cart are upregulated while pro-thyrotropin-releasing hormone and thyrotropin-releasing hormone (trh) receptor are downregulated (Fischer et al., 2019). In A. burtoni, homologs of these ligands and receptors (cartl, trh, crhbpl and crh) are strongly correlated in expression with increased expression in B02 and R02 surrounding the transition from mouth-brooding care to defense of fry. In mice, Crhbp inhibits endogenous Crh to stimulate maternal aggression and deletion of the gene inhibits maternal aggression (Gammie et al., 2008), suggesting it may play a similar role in the A. burtoni maternal transition.

Under normal conditions, the arcuate nucleus contains two specific populations of neurons, one producing orexigenic NPY and AgRP to stimulate feeding behavior, and the other producing POMC and CART to inhibit feeding, which are directly responsive to the peripheral signals of satiety – leptin (Cone, 2005) and ghrelin (Cowley et al., 2003) – and indicate energy reserves and body condition. However, in brooding females, the reduction in body mass (Grone et al., 2012; Renn et al., 2009) that should signal to increase orexigenic signals does not result in feeding behavior during the mouthbrooding time points. Among the many hormone/neuropeptide genes related to feeding behavior, we found several to be among those trending toward differential expression across the reproductive time points, including relaxin-3-like (rln3l), gal, cartl, cholecystokinin-like (cckl), pro-opiomelanocortin B (pomcb) and peptide YY (pyy) (Fig. 3). While many of these showed large variance within time points, they appeared to be responsive to the resumption of feeding post-release.

The anorexigenic genes pomcb and cart showed highest expression after resumption of feeding, at R02, and were low during brooding, again reflecting energetic state. These genes are regulated by the peripheral signals leptin and insulin (Lau and Herzog, 2014; Valen et al., 2011). The orexigenic genes gal and rln3l clustered with nts (Fig. 3) and showed increased expression during mouthbrooding when females are not eating. As these genes are expected to promote food intake (Ganella et al., 2013; Volkoff and Peter, 2001), their upregulation during brooding suggests that the brains may be receiving the typical signals of hunger from the periphery, although the females suppress or do not respond to the signals with increased feeding. Alternatively, these genes may be playing a role directly promoting maternal mouthbrooding care.

While anorexigenic neuropeptide Y (npy) did not show differential expression across time points in the current study, a related peptide, pyy, was elevated in mouthbrooding stages (unadjusted P=0.01). In the digestive system, pyy expression increases with food intake to signal satiety (Karra et al., 2009); however, in the brain, its response to fasting is variable across species (Volkoff, 2016). In a study of the mouse enteroendocrine system, pyy was shown to be co-expressed and co-secreted with both nts and glucagon-like peptide-1 (glp1) (Grunddal et al., 2016). The genes pyy and glp1 were each found to work synergistically with nts in these mice to inhibit gastric emptying and decrease food intake. In A. burtoni maternal brains, pyy expression correlates with nts expression (Fig. 3) across samples through a clustering of candidate differentially expressed neurohormones and neuropeptides. Its role in this system and potential interaction with nts in the regulation of fasting behavior warrants further study.

In addition to the potential role of nts in parental behavior, this neuropeptide has a pleiotropic role and also can regulate fasting. In a number of studies in mice and rats, researchers have shown that increases in nts in the brain suppressed feeding behavior through interaction with its associated receptor, ntsr1 (Kim et al., 2008; Schroeder and Leinninger, 2018). In A. burtoni, we saw the highest expression of nts in mouthbrooding fish and a sizable drop in non-brooders but no change in low-expressed ntsr1 (Fig. 4). This pattern is consistent with the hypothesis that nts plays a role in fasting behavior in mouthbrooding females by eliciting an anorectic effect despite unchanged expression of ntsr1. It is possible that mouthbrooding A. burtoni may actually be suppressing signals that would normally increase hunger and feeding behavior. Such suppression has been suggested by studies addressing the expression of feeding behavior-associated hormones, neuropeptides and their receptors in A. burtoni and Oreochromis niloticus brains (Das et al., 2019; Grone et al., 2012).

Many of our DE genes point to compensatory metabolic mechanisms to allow fish to withstand long periods of fasting. For example, the pyruvate dehydrogenase kinases pdk2 (FDR=9.71E−5) and pdk4 (FDR=0.024) show significantly higher expression in the brain in the two mouthbrooding time points relative to the two non-brooding time points. These genes are associated with food deprivation in fish, mammals and other taxa, and have been shown to increase during periods of starvation then decrease after resumption of feeding (Wu et al., 2000; Yang et al., 2019). This energy intensive pathway mediates the interaction of glycolysis and the tricarboxylic acid cycle (Sugden and Holness, 2003). The functional conservation of pyruvate dehydrogenase kinase activity across taxa may have been a necessary mechanism allowing for the evolution of mouthbrooding behavior in cichlids and other independent fish lineages as a result of the associated fasting behavior.

GSEA results also corroborate this change in mitochondrial function, with GO sets related to mitochondrial ATP production and electron transport showing upregulation in B02 fish relative to R14 fish (Fig. 2; Table S3). These GO terms signal a shift in energy metabolic processes at the start of a long period of prolonged fasting, and include gene sets such as ‘ATP biosynthetic process’, ‘mitochondrial fusion’, ‘mitochondrial electron transport, NADH to ubiquinone’, ‘proton-transporting ATP synthase activity, rotational mechanism’ and ‘mitochondrial proton-transporting ATP synthase complex’. Metabolic downshifts are common among starving or hibernating animals and are often characterized by a reduction in mitochondrial function in liver and skeletal muscle tissue for energy conservation (Barger et al., 2003). That we see these changes anticipated in the periphery occurring in brain samples may reflect the high metabolic demand of neural tissue.

Late-stage (post-hatching) mouthbrooding cichlids show reduced tolerance to hypoxia (P. multicolor: Corrie et al., 2008) possibly as a result of the costs of maternal behavior including buccal churning to aerate the eggs. However, we found a dramatic reduction in gene expression for several globin genes, which could also contribute to the observed sensitivity to hypoxia. Four genes reached, and one approached, a statistically significant reduction in expression at the late-brooding B14 and R02 time points (Fig. 5). These contribute to GO term enrichment (Table 2) and GSEA results (Fig. 2) for oxygen transport-related terms. These are among the most dramatic gene expression changes observed, likely reflecting consistent regulation across multiple brain regions. The reduction in globin expression may be a response to nutrient state as seen under starvation in seabream (Ntantali et al., 2020). The sustained expression could reflect ‘carryover’, the persistence of acquired neurogenomic states into subsequent stages (Burkhari et al., 2019). However, we also found hypoxia-inducible factor 1-alpha-like (FDR=0.12, unadjusted P=7.72E−4) elevated in both brooding time points along with two target genes, pdk4, known to be also induced by starvation (Lee et al., 2012), and round spermatid basic protein 1-like (rsbn1l), a less understood gene known to be expressed in testes as well as breast cancer tissue (Abu-Jamous et al., 2017). The correlated expression of rsbn1l and pdk4 supports the hypothesis that mouthbrooding cichlids are responding to hypoxia-related signals possibly due to reduced globin gene expression.

The reduction in nutrient intake and apparent reduction in neural oxygen transport capacity confer a large metabolic cost that could be mitigated by neuroprotective pathways as suggested by GSEA (Figs 2 and 4). We saw varied and dynamic patterns among the neuroprotection-related genes, some being upregulated in late brooding and release such as two butyrophilin-related DE genes (butyrophilin subfamily 2 member A1-like and V-set domain-containing T-cell activation inhibitor 1-like) that co-localize on a genome scaffold and are known immune system regulators that may have neuroprotective function in response to hypoxia, reduced nutrient intake or the action of other DE genes associated with foreign DNA/RNA (Abeler-Dörner et al., 2012). Conversely, other genes involved in neuroprotection were downregulated at late brooding time points and upregulated at early brooding time points. Downregulation of ECM/collagen-related DE genes including matrix metalloproteinases such as matrix metalloproteinase 9 (mmp9) could impact recovery via delayed promotion of neuroblast cell migration (Chen et al., 2009), while regulation of platelet-derived growth factor receptor-like (pdgfrl), similar to other PDGF proteins seen in cerebral hypoxia of stroke victims (Krupinski et al., 1997), could aid in neural regeneration in response to neuronal damage.

Nutrient and oxygen deficit impact the regulation of ion gradients and ion transport pathways in the brain that require energy and aerobic ATP synthesis. Both ion transport and cell junction genes are upregulated in whole-brain samples under starvation in seabream (Ntantali et al., 2020) and hypoxia-tolerant species are known to have evolved means to reduce ion leakage in the CNS (Boutilier, 2001), thus reducing energy demands in such conditions. Similarly, increased myelination mitigates the loss of action potentials under hypoxia (Waxman et al., 1990). These facts are consistent with the observed upregulation at B14 and R02 time points of cell adhesion molecule 4 (cadm4) (FDR=2.21E−14 and 3.2-fold change), which regulates myelination of the CNS (Elazar et al., 2019; Golan et al., 2013) and could contribute to the maintenance of action potentials under starvation and hypoxic conditions during mouthbrooding. These and other DE genes, including those coding for myosin, actin and troponin complex proteins, contribute both to molecular function GO terms, such as ion transport-related and action potential gradient-related GO sets, and biological process GO terms related to the structure and function of skeletal muscle and heart tissue (Fig. 2; Table S3), supporting the hypothesis of added demand to regulate ion gradients for neural function during mouthbrooding. Here again, the patterns of regulation among genes contributing to ion transport-related GO terms is highly variable as a result of the individual genes being negative versus positive regulators of the same or related processes.

In addition to oxygen transport and ion regulation, our GSEA analysis detected the inflammatory pathway as a large family of 16 GO sets under the umbrella term ‘interleukin-1 beta production’ (Table S3) that showed decreased expression at the early brooding B02 time point and increased expression at the late brooding B14 time point. While IL1B, a cytokine secreted by microglia in the brain in response to hypoxia or injury to the brain (Hewett et al., 2012), is not itself differentially expressed in our dataset, GSEA identifies genes involved in IL1B suppression via (NF-κB) signaling (e.g. protein nlrc3-like, NACHT, LRR and PYD domains-containing protein 12-like) (Allen, 2014). The upregulation of these pathways and static expression of cytokines may be adaptations that help avoid negative physiological consequences of prolonged inflammation associated with hypoxic conditions and starvation experienced during each repeated reproductive cycle (Martin and Król, 2017; Schäfer et al., 2021). To the best of our knowledge, these immune function pathways have not been implicated in other whole-brain studies of fish under starvation (e.g. Ntantali et al., 2020; Yang et al., 2019) but have been associated with paternal care in stickleback (Burkhari et al., 2019) suggesting their regulation may represent an adaptation in cichlids that reduce the costs of mouthbrooding.

Using female A. burtoni, we carried out a comprehensive analysis of whole-brain gene expression and how it changes across distinct brooding stages. The gene expression patterns reveal how the female brain may resolve trade-offs between metabolic demands and parental care. Similar conclusions can be drawn from lists of DE genes, GO analysis and the GSEA results that rely on differing levels of statistical significance for individual genes. In addition to identifying itnp, classically associated with parental care, our data support the involvement of nts, a gene that has more recently emerged as a candidate in the regulation of social phenotypes. These, and other identified DE genes, could be subjects for exciting functional analysis with gene editing techniques that are being made available in A. burtoni and other emerging model organisms (Juntti, 2019). Both GO enrichment and GSEA analysis further support the involvement of receptors and downstream activity of this pathway. Furthermore, correlation of gene expression ties these pathways to important neuropeptides involved in the regulation of feeding behavior and metabolism, highlighting the known tradeoffs between care and feeding (O'Rourke and Renn, 2015).

Interestingly, we found both orexigenic and anorexigenic genes that showed increased expression during mouthbrooding stages as well as anorexigenic genes that showed decreased expression during mouthbrooding, which underscores the complexity of feeding regulation under voluntary starvation. Future transcriptomic studies should include food-deprived individuals as an additional control group for comparison with mouthbrooding fish. Such a comparison would help elucidate the degree to which our observed increase in neuroprotective pathways and increased expression of negative regulators of inflammation represent evolved mechanisms to counter the hypoxic and metabolic demands of mouthbrooding. Our results underscore that cichlid maternal behavior and fasting are potentially regulated by pleiotropic activity of novel peptides such as neurotensin as well as the expected classical neurohormones and neuropeptides. It is noteworthy that these hormonal signals were dwarfed by dramatic changes in oxygen transport capacity, neuroprotection and ionic regulation pathways in the brain that mirror metabolic shifts from animals regularly subjected to starvation or hypoxic conditions. Here, additional time points through these transitions could be used to further test the degree to which gene expression changes persist to reflect previous neurogenomic states versus priming the brain for future neurogenomic states (Burkhari et al., 2019). While the observed changes in neurohormones are logically brain specific, many of the other gene expression changes are likely not unique to brain tissue and rather reflect systemic responses to the stress of starvation and energy demands.

Samples and preliminary data were collected by Kimberly Engeln with fish husbandry support from several Reed lab undergraduates.