Galanin and prolactin expression in relation to parental care in two sympatric cichlid species from Lake Tanganyika General and Comparative Endocrinology

Our understanding of the hormonal mechanisms underlying parental care mainly stems from research on species with uniparental care. Far less is known about the physiological changes underlying motherhood and fatherhood in biparental caring species. Here, using two biparental caring cichlid species ( Neolamprologus caudopunctatus and Neolamprologus pulcher ), we explored the relative gene-expression levels of two genes implicated in the control of parental care, galanin (gal) and prolactin (prl). We investigated whole brain gene expression levels in both, male and female caring parents, as well as in non-caring individuals of both species. Caring males had higher prl and gal mRNA levels compared to caring females in both fish species. Expression of gal was highest when young were mobile and the need for parental defense was greatest and gal was lowest during the more stationary egg tending phase in N. caudopunctatus . The onset of parenthood was associated with lower expression of prl and higher expression of gal in N. pulcher , but this pattern was not observed in N. caudopunctatus . Our study demonstrates that gal gene expression is correlated with changes in parental care in two biparental cichlid species and extends both knowledge and taxonomic coverage of the possible neurogenetic mechanisms underlying parental care.


Introduction
Across species, parental care varies widely, with great diversity in both, the type of care performed and in the sex of the care giver (Clutton-Brock, 2002;Reynolds et al., 2002;Royle et al., 2012;Sefc, 2011). Our current understanding of the neural substrates underlying parental care largely comes from research on the brain of female mammals (Choleris et al., 2013;Dulac et al., 2014;Lee et al., 2009). Although studies of the mechanisms underlying paternal, biparental and alloparental care (e.g., care by non-parents or helpers-at-the-nest) are increasing, we still know little about the neural substrates associated with these less common forms of care (Bendesky et al., 2017;Bukhari et al., 2019;Dulac et al., 2014;Feldman et al., 2019;. Research conducted to date on the hormonal mechanisms underlying care suggests considerable evolutionary conservation across vertebrates with specific hormones like prolactin playing key roles in care modulation within and across species (Angelier and Chastel, 2009;Schradin and Anzenberger, 1999;Whittington and Wilson, 2013). Recently, another neuropeptide, galanin (Gal) has been shown to have a role in parental care regulation in mice . Wu et al. (2014) found that during parenting, a subset of Gal expressing neurons in the medial preoptic area of the brain (MPOA) were activated, and that genetic ablation of these neurons resulted in impairment of parental care. Further research by Kohl and colleagues (Kohl et al., 2018) showed that if virgin males, which normally attack pups, had optogenetic activation of these Gal expressing neurons, they instead groomed pups Wu et al., 2014). Expression of galanin is correlated to parental care in other taxa as well, including some fish species (Bukhari et al., 2019;Tripp et al., 2020); however, nearly all studies conducted to date have focused on species with uniparental care (but see Fischer et al., 2019 for an exception). To improve understanding of the neurogenomic mechanisms underlying biparental care and to expand taxonomic coverage of this topic, we used two cichlid fish species to explore the relative gene-expression levels of prl and gal, two genes potentially implicated in the control of parental care.
Galanin is a small 29 amino-acid peptide first isolated from the porcine intestine (Tatemoto et al., 1983) and subsequently was identified in birds, reptiles, and fishes (Wynick et al., 2001). Among its many functions, galanin stimulates the release of prolactin , growth hormone (GH)  and luteinizing hormone (LH) (López et al., 1991) from the pituitary gland, all of which have been associated with parental care (Mammals: Ziegler, 2000; for review see Ziegler, 2000;Birds: Sharp et al., 1988, 1979 for review see Buntin, 1996). The anterior commissural nucleus (ACN) in the MPOA revealed a positive response by Gal neurons during suckling behavior, increasing maternal motivation in rats (Rattus norvegicus) (Cservenák et al., 2017). Among teleost fishes, gal is highly expressed in the diencephalon (including the preoptic area) during the nest, egg, and early hatching stages of care in paternal male three-spined sticklebacks (Bukhari et al., 2019). However, the gal neurons in the POA-AH of parental midshipman fish (Porichthys notatus) were not highly activated by direct brood care (Tripp et al., 2020(Tripp et al., , 2019. Aside from these correlative findings from uniparental species, increased neural activation of galanin positive POA neurons has also been reported in parents during tadpole transport in a biparental species of poison frog (Dendrobates imitator, Fischer et al., 2019).
Among fishes, uniparental male-only care is the most common type of care, followed by female-only care, while biparental care is less common (Balshine, 2012;Sefc, 2011). In one large family of fishes, the cichlids, there is a variety of care patterns found, with 2/3 of cichlid species showing female-only uniparental care and 1/3 showing biparental care (Goodwin et al., 1998;Sefc, 2011). Cichlids have rapidly diverged (Gante et al., 2016) making it possible to find many closely related species that differ in care but that live in similar habitats and under similar ecological conditions. These characteristics make the cichlids an excellent model system for comparative studies on the proximate causes of parental care.
In this manuscript, we investigated the possibility that prl1 and gal expression are linked to parenting behavior in both biparental and even in cooperative species, as has been previously shown for uniparental species. There are several forms or variants of prolactin in fishes; here we examined only prl1 as it is common to all vertebrates (Whittington and Wilson, 2013) while prl2 is highly expressed in the eye and is thought to be involved in retinal development (Huang et al., 2009). The gal gene structure is highly conserved across vertebrates (Mensah et al., 2010) and for this study, we used the long-form of the gal gene (Martins et al., 2014). We used two cichlid species from the same genus N. caudopunctatus and N. pulcher, that differ in caring strategies. In N. caudopunctatus, a male and a female raise the young together, and so this species is classified as a biparental breeding fish (Ochi and Yanagisawa, 1999). In contrast, in N. pulcher the breeding pair is joined by both related and unrelated non-reproducing group members in raising offspring and is classified as a cooperative breeder (Stiver et al., 2005;Wong and Balshine, 2011). In both of these cichlid species, the male and the female share the parental task of defending the brood and the territory (Cunha-Saraiva et al., 2019, 2018Wong and Balshine, 2011). However, the role of non-breeding fish differs significantly between these species. In an established social group, non-breeding subordinate N. pulcher usually help look after young, while non-breeding N. caudopunctatus are voracious egg consumers and are chased off by the breeding pairs (Balshine et al., 2001;Cunha-Saraiva et al., 2018;Wong and Balshine, 2011).
We predicted that caring individuals would have higher prl1 mRNA levels than non-caring individuals in both species. Given that the optogenetic activation of gal neurons is connected with the inhibition of infanticidal behaviour in male virgin rats Wu et al., 2014), we also predicted that gal would be down-regulated in non-caring individuals compared to caring individuals. Furthermore, we predicted that correlated prl1 and gal gene transcript expression levels would be linked with the degree of parental care behaviours performed.

Study animals and housing conditions
Fifty-six Neolamprologus caudopunctatus (28 males and 28 females) and 42 Neolamprologus pulcher (21 males and 21 females) were used in the study. Biparental N. caudopunctatus pairs excavate a breeding cavity in crevices or under stones where they spawn and guard their young together for up to six weeks (Ochi and Yanagisawa, 1999). N. pulcher is also biparental, both parents look after young by providing brood care and defense, but they live in social groups and breed cooperatively. Groups consist of a dominant breeding pair and 1-20 smaller subordinates that typically don't breed (Hellmann et al., 2015;Stiver et al., 2005). Both breeders and subordinates participate in territory defense, territory maintenance, direct brood care of the young (by cleaning and fanning the eggs and defending the young) and can remain in the same territory for several years (Balshine et al., 2001;Jungwirth et al., 2019;Stiver et al., 2005). All fish were wild caught from the most southern tip of Lake Tanganyika/Republic of Zambia in autumn 2015 with an unknown parental history.
Each fish was sexed and measured for standard length (SL), total length (TL), and mass (M) (see Supplementary Table 1 for morphological details). Fish were fed daily with either frozen food (a mixture of artemia, cyclops, and daphnia species plus red mosquito larvae) or with tropical fish flakes. Holding and experimental tanks were maintained at a constant water temperature of 26 ± 1 • C under a 12/12 h light/dark cycle. Prior to using the fish in experiments, the fish were housed in single-sex stock tanks (400 L) containing a maximum of 40 individuals for the smaller N. caudopunctatus and a maximum of 30 individuals for the larger N. pulcher. In the behavioral assays described below, the F1 sexually immature juvenile offspring from N. caudopunctatus wildcaught parents (fish that were between 30 and 45 mm, standard length SL) and Telmatochromis vittatus were used as territorial intruders. Two different intruder species were used because previous studies have identified sexually immature N. caudopunctatus juveniles as both territory intruder and egg/offspring consumers that elicit a strong defensive response from N. caudopunctatus parents (Cunha-Saraiva et al., 2019, 2018Ochi and Yanagisawa, 1999;Schaedelin et al., 2015). In the case of N. pulcher pilot studies in our lab showed that T. vittatus elicits a strong defensive response from N. pulcher parents and subordinate helpers of this species and a similar response to T. vittatus has been reported in a closely related and also cooperative species, N. savoryi (Josi et al., 2019).

General procedures and experimental protocols
A male and female pair of N. caudopunctatus were each placed in a 45 L experimental tank (n = 21) containing a breeding shelter (half a flowerpot), a sponge filter, a heater and a 3 cm sand layer. Breeding pairs in the experimental tanks were checked daily for signs of nest preparation and maintenance, spawning, pair stability and general health. Pairs were given a minimum of 7-days prior to any behavioral observations and were randomly selected for observation at different stages of the breeding cycle (Fig. 1a). Randomization was performed by giving each pair a unique ID, which was used to design a randomization table by using the cluster_ra function from the randomizr package in R. 1) Seven pairs were observed and sampled before egg laying began (classified as pre-spawning pairs); 2) seven pairs were observed and sampled within the first 24 hrs following egg laying (classified as spawning pairs) and 3) seven pairs were sampled 10 days after egg laying, during the free-swimming care phase of offspring development (classified as post spawning pairs). To compare breeding pairs to nonpaired non-breeding individuals we also sampled 4) seven non-paired males and seven non-paired females from same-sex 160 L holding stock tanks (classified as non-breeding individuals). All fish used in this experiment were kept under similar environmental conditions and water quality and fish general health were checked daily. A series of behavioral assays were performed with the fish in the different treatment or classification groups, and then immediately afterwards the fish were sampled. The sampling time points across the breeding cycle are described above and took into consideration tight connection between egg laying and the cessation of egg cannibalism in N. caudopunctatus (Cunha-Saraiva et al., 2018). These discrete sampling points (nonpaired, pre-spawning, spawning and post-spawning) allowed the genomic dynamics of prl1 and gal expression and possible connections with parental care to be explored.
We created 7 N. pulcher breeding pairs, by selecting males and females from the single sex holding tanks (400 L). One male and one female were placed together in an 80 L experimental tank . The experimental tanks contained 3 cm of sand as substrate, half a flowerpot as a breeding shelter, a heater and a sponge filter. In each pair the male breeder was always at least 5 mm larger than the female breeder as this is the size difference found in natural pairs (Balshine et al., 2001). N. pulcher were placed in larger tanks because they are larger and have larger territory sizes compared to the N. caudopunctatus (Balshine et al., 2001;Schädelin et al., 2012). The breeding pairs were checked daily, for territory maintenance, spawning, pair stability and general health for a period of 7-days before behavioral observations. Once spawning (within the first 24 hrs of egg laying) was observed, a series of behavioral assays were performed with the seven breeding pairs and then the fish were terminally sampled.
To compare breeding N. pulcher pairs to non-breeders ( Fig. 1b) and to disentangle the effect of social status on whole-brain mRNA levels of prl1 and gal, a well-established laboratory protocol was followed (Taborsky et al., 2012;von Siemens, 1990;Zöttl et al., 2013) where 7 male-male and 7 female-female social dyads of non-reproductive individuals were established. One individual of the dyad was always 5 mm larger than the other as this distinct size difference leads to a quick and consistent dominance hierarchy, with the larger fish of the same sex dyad emerging as dominant and the other one as the subordinate (Hamilton et al., 2005;Heg et al., 2008;Reddon et al., 2011;Wong and Balshine, 2011). These dyads were kept together (two fish in 80 L experimental tanks) for 7 days and then they were terminally sampled.

Behavioral assays and scoring
To control for diurnal variation in behavior and physiology (Desjardins et al., 2011;Lema et al., 2010;Werner et al., 2003), behavioral observations always took place between 10 h and 13 h. Parental behavior was observed and scored ~ 1 h before the whole-brain sampling by conducting 10-min behavioral observation followed by a 2-min nest intrusion assay. To ensure that the fish were not disturbed by observer presence, the observer sat still 1 m from the tank and the observation period only began after a two-minute habituation period to the presence of the observer. Following the habituation period each pair was observed for 10 min and all behaviors including any aggressive (e. g., head down or opercula spread) and submissive behaviors (e.g., tail quiver or flee), as well as any parental care/brood care (e.g., egg cleaning and fanning) and territory maintenance (e.g., nest cleaning and Fig. 1. a) A diagram illustrating the different reproductive stages of Neolamprologus caudopunctatus's breeding cycle. Non-breeding individuals and pre-spawning breeding pairs are egg cannibals, while spawning and post-spawning pairs are noncannibals (Cunha-Saraiva et al., 2018); b) The diagram illustrates the different social and reproductive stages of Neolamprologus pulcher's breeding cycle. Non-breeding individuals live in a social group with a breeding pair that is always dominant to the non-breeders . building) behaviors were recorded. A full description of all the behaviors recorded for N. caudopunctatus can be found in Cunha-Saraiva et al., 2018, andfor N. pulcher in Sopinka et al., 2009. Thereafter, a nest defense assay was performed by placing three juvenile N. caudopuncatus or two adults Telmatochromis vittatus in a transparent plexiglas cylinder (diameter × height: 18 × 40 cm) placed at ~ 10 cm from the entrance of the breeding shelter in N. caudopunctatus or N. pulcher pair tanks, respectively. The nest defense assay lasted 2 min and began as soon as one of the fishes performed an aggressive display or act towards the egg predators in the cylinder (mean latency to attack the cylinder: ~30 s, range: 0 -120 s). All aggressive behaviors, such as head-down display, approach or opercula spread, towards the N. caudopunctatus juveniles and towards the T. vittatus were recorded and classified as defense behavior.

Tissue sampling and processing
One hour after the observations the individuals were captured and euthanized with an overdose of anesthetic (MS222, Sigma, 1000 mg/L) and then the spinal cord was severed. Fish were then measured, dissected, and their sex was confirmed. The whole brains were preserved in vials of RNAlater (EURx, Gdansk, Poland) that were refrigerated at 4 • C for 12 hrs before storing at − 80 • C. Brains were defrosted and homogenized and total RNA was isolated using a Maxwell ® 16 total RNA Purification kit. RNA integrity and purity were assessed by 1% agarose gel electrophoresis and was quantified using a NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, USA). Prior to cDNA synthesis, total RNA was treated with DNase (DNA-free kit, Ambion, UK) to eliminate genomic DNA contamination, according to the manufacturer's instructions. First strand cDNA was synthesized from 500 ng DNase-treated total RNA in a 20 µl reaction containing 500 ng of DNasetreated RNA, 200 ng of random hexamers (Jena Biosciences, Germany), 100 U of RevertAid (Fermentas, Thermo Fischer Scientific, USA) reverse transcriptase and 8 U of RiboLockRNase Inhibitor (Fermentas). Reactions were incubated for 10 min at 25 • C and 60 min at 42 • C, followed by enzyme inactivation for 10 min at 70 • C, and then stored at -20 • C until use.

Candidate gene primer validation
Whole-brain gene expression levels of prl1 and gal were studied, using quantitative real-time PCR (RT-qPCR). Primers for PCRamplification of a partial cDNA sequence for each of the genes of interest were designed based on the sequences of prl1 and gal of four cichlid species (Oreochromis niloticus, Neolamprologus brichardi, Astatotilapia burtoni and Pundamilia nyererei; Supplementary Material Table 2). Multiple sequence alignments were performed using ClustalX (v2.0.11) (Larkin et al., 2007), the alignments were edited using GeneDoc version 2.7.0 and the conserved 5' and 3 ′ flanking regions were selected for primer design. The selected primers were then used to amplify the ORF of prl1 and gal from a pool of brain cDNA. RT-PCR amplification was done using 1 µl of cDNA, 50 µM dNTPs, 10 pmol of each primer and 0.5 U DreamTaq DNA Polymerase (Fermentas). Cycling conditions were 5 min at 95 • C followed by 35 cycles of denaturation for 20 s at 95 • C, 20 s of annealing at the optimized temperature for primer pairs and 1 min of extension at 72 • C, followed by a final extension of 5 min at 72 • C. The amplified fragments were purified using a GFX DNA purification kit (GE healthcare) and sequenced to confirm their identity. The percentage of identity of the N. caudopunctatus and N. pulcher prl1 and gal nucleotide sequence, respectively was confirmed using the BLASTN and BLASTX tools to search against the Genebank database. Multiple-sequence alignments with other cichlid prl1 and gal were executed as described above to determine sequence identity and confirm gene identity.

Analysis of gene expression by quantitative real-time PCR
Transcript levels of prl1 and gal (candidate genes) and β-Actin and 18S (reference genes, (Ahi et al., 2017) were measured by real time quantitative PCR (RT-qPCR) in a CFX Connect qPCR thermocycler (BioRad). A sample of 14 fish were analyzed at each point in the reproductive cycle: for N. caudopunctatus (non-breeding, pre-spawning, spawning and post-spawning pairs) and for N. pulcher (breeding pairs, non-breeding dominant and subordinate). Quantification of transcripts was carried out in duplicate using the relative standard curve method and EvaGreen chemistry, following the procedure described in Martins et al. (2014). Briefly, prl1, gal, β-actin (actb) and 18S mRNA levels were measured in each sample in duplicate using 11 μl reactions that contained 2 μl of each cDNA (diluted 1:2), 300 nM of each specific primer (Table 3 Table 3), followed by a final melting curve between 65 and 95 • C. Specificity of qPCR reactions for prl1 and gal was confirmed by the presence of a single peak in melt curves, amplification of a product of the expected size (prl1: 129 bp, gal: 120 bp, actb: 189 bp, 18S: 103 bp) and the sequence of the amplified products. A cDNA control sample (a standard cDNA synthesis reaction without the addition of reverse transcriptase, − RT) was included for both genes to confirm there was no genomic DNA contamination and no product was detected. Standard curves relating amplification cycle with initial template quantity were generated using serial dilutions of specific RT-PCR products for each gene under quantification. The efficiency of qPCR reactions ranged from 91% to 98% for both prl1 and gal, with an R 2 > 0.98. Relative expression levels were calculated using the comparative C T method (2 − ΔΔCT , Livak and Schmittgen, 2001) and the expression of prl1 and gal were relative to actb for N. caudopunctatus and relative to 18S for N. pulcher. Different reference genes were used due to the lack of primer efficiency of actb for N. pulcher, notwithstanding, both reference genes did not vary in expression between samples. All samples for each of the species were extracted and reverse transcribed as a group, and each gene was measured on a single PCR plate for each species. The relative values for each gene can therefore be directly compared between each sampling point for each species separately.

Statistical analysis
We used R statistical software 3.6.1 to perform the analyses. Prior to any statistical test, all dependent variables were categorized according to their distribution and transformed when necessary to reach normal distribution of residuals and equal variance. prl1 mRNA levels were logtransformed to meet normality for both N. caudopunctatus and N. pulcher, whereas untransformed gal mRNA levels followed a gaussian distribution for both species (Shapiro-Wilk test, P > 0.05). To determine how prl1 and gal mRNA levels change across the defined sampling points (N. caudopunctatus: non-breeding, pre-spawning, spawning and postspawning; N. pulcher: non-breeder subordinate, non-breeder dominant, and breeding) for each species (see Section 2.2 of the methods), an ANOVA (aov function, from the stats package) was run for each gene with sampling point, sex (male and female) and their interactions as factors. Each model was validated by assessing the distribution of its residuals and variance (Johnson and Omland, 2004). Effect sizes estimates were computed as follows: Eta-square for analysis of variance (etaSquared function, from the lsr package; Navarro, 2015) and Cohen's d for pair-wise comparisons (cohen.d function, from the effsize package; Torchiano, 2018).
To test for differences between the different stages of the breeding cycle, we used glht as a general linear hypothesis testing function in R (Hothorn et al., 2008) and employed user-defined contrasts (N. caudopunctatus: non-breeding vs spawning; pre-spawning vs spawning; spawning vs post-spawning; non-breeding vs post-spawning; N. pulcher: breeding vs non-breeding dominant; breeding vs nonbreeding subordinate; non-breeding dominant vs non-breeding subordinate). This is a single-step method test, which adjusts the p-values to the number of comparisons using Tukey contrasts. User-defined contrasts were established as such to provide a comparison between the different sequential phases of the species' breeding cycle as described in Fig. 1a, b.
To test if males and females differed in the amount of parental care behavior provided, a brood defense score was calculated for the defense assay (the sum of all observed aggressive behaviors, such as head down, approach and opercula spread, towards the nest intruders) and a brood care score (the sum of all observed care behaviors including egg cleaning, fanning and cavity visits) for each individual in the study. Since egg cleaning/fanning and cavity visits were strongly correlated (Spearman correlation, rho = 0.65, P < 0.001, corr.test function was performed using the psych package, Revelle, 2019), we used cavity visits as our brood care measure. To analyze behaviors, all the behavioral variables were SQRT transformed and then used as the response variable in a generalized linear model. No statistical difference was found between the spawning and post-spawning pairs in terms of brood care or defense (Brood care: ANOVA, F-test = 0.01, P = 0.89; Brood defense: ANOVA, F-test = 0.40, P = 0.52). Thus, we explore the link between gene expression levels and the behavioral sex differences in N. caudopunctatus. We focused on the most active and long phase of care, Fig. 2. a) prl1 and b) gal relative mRNA levels for both male and female Neolamprologus caudopunctatus at each reproductive stage. Black boxes represent males (M) and white boxes represent females (F). Box-plots show the interquartile range (IQR) of each group analyzed with whiskers extending to 1.5 × the IQR. Horizontal lines represent medians. Outliers were included in the analysis and are depicted on the graphs as white circles. Significant differences as revealed by Tukey Contrasts are depicted on the graph as bold stars, all other contrasts were not significant. *P < 0.05, **P < 0.01. the post-spawning phase which last roughly 40-days (Cunha-Saraiva et al., 2019;Ochi and Yanagisawa, 1999). Each model was validated by assessing the distribution of its residuals and variance (Johnson and Omland, 2004).
Finally, to explore how the relative expression of each candidate gene might be linked to specific behaviors, a spearman rank correlation was performed. Because gene expression data was compared to several different behaviors (e.g. brood care and defense), correlations were corrected for multiple comparisons using the "holm" P-value adjustment method (Aickin and Gensler, 1996). More specifically, the relative mRNA levels of prl1 and gal were correlated with 1) brood care and 2) defense in both N. caudopunctatus and N. pulcher.

Ethical note
The procedures described in this study were approved by the Uni-
Female and male N. caudopunctatus were not overall different in gal mRNA expression levels (Table 1), however, similar to prl1, parental males showed higher gal mRNA levels compared to females at the stage of care when young were mobile and free-swimming (post-spawning) (Fig. 2b, Table 1 Contrary to the patterns of male biased territorial defense reported previously in the field and in the lab for this species (Ochi and Yanagisawa, 1999;Cunha-Saraiva et al., 2018), in this study, male N. caudopunctatus did not defend the brood more than females (Table 2). However, as expected based on previous laboratory and field results, females did care for eggs and tended the nest more than males (Table 2;  see also Supplemental Table 4 for behavioral descriptive statistics).
prl1 appeared to be negatively correlated with brood care in postspawning males, but not females (Fig. 4, Table 3), while gal mRNA levels were not correlated with brood care or offspring defense during the post-spawning phase of care for both males and females (Table 3). Moreover, prl1 and gal mRNA levels were not correlated in postspawning pairs (Spearman correlation, N = 13, rho = 0.22, P = 0.45).

prl1 and gal mRNA expression in Neolamprologus pulcher
Non-breeding subordinate N. pulcher had the highest prl1 mRNA levels, and they were considerably higher than both non-breeding dominant fish and breeding N. pulcher, whilst no clear difference was found between breeding pairs and non-breeding dominant fish (Fig. 3a, b, c; Table 1, breeding pairs vs non-breeding subordinates, t = 4.46, P < 0.001, Cohen's D = 0.68; non-breeding dominants vs non-breeding subordinates, t = 6.14, P < 0.001, Cohen's D = 1.19; breeding pairs vs non-breeding dominants, t = − 1.62, P = 0.25, Cohen's D = 0.001). N. pulcher males and females differed in prl1 mRNA levels (Table 1), however, the effect of sex on prl1 mRNA levels was dependent on the time point of sampling. Male N. pulcher breeders had higher prl1 mRNA levels compared to breeder females whereas non-breeding subordinate females had higher prl1 mRNA levels compared to non-breeding subordinate males (Fig. 3a and c; Table 1; breeding: M vs F, t = 2.67, P = 0.03, Cohen's D = 0.52; non-breeding dominants: M vs F, t = − 1.02, P = 0.66; non-breeding subordinates: M vs F, t = − 3.08, P = 0.01).
gal mRNA levels were also dependent on the social or breeding status as spawning individuals had higher gal mRNA expression levels compared to non-breeders, both dominants and subordinate's females and males (Fig. 3d, e, f; Table 1; breeders vs non-breeding subordinates,   . Moreover, N. pulcher male breeders had higher gal mRNA levels compared to female breeders while no such sex differences were found between male and female non-breeding dominants or subordinate individuals (Fig. 3d, e, f;  Table 4 for behavioral descriptive statistics). Brood care (e.g., fanning, cleaning and cavity visits combined) was negatively correlated with gal mRNA levels in males but not females (Table 3, Fig. 4), but no relationship was observed with defense (Table 3). No relationship was detected between prl1 mRNA levels and parental behaviors in either males or females (Table 3). Notwithstanding, prl1 and gal mRNA levels were positively correlated in breeding pairs (Spearman correlation, N = 14, rho = 0.69, P = 0.01).

Discussion
This study revealed high gal mRNA levels, but not high prl1 mRNA levels, in the breeders of two cichlid species the biparental N. caudopunctatus and the cooperative breeder, N. pulcher. The caring males of both species expressed significantly higher prl1 and gal mRNA levels than did caring females. prl1 was negatively correlated with the frequency of brood care in male N. caudopunctatus, while the expression of gal was negatively correlated with the frequency of brood care in male N. pulcher.
Several studies have shown that increased gal neuron activation is associated with reproductive behavior, such as the onset of caring in mice (Bukhari et al., 2019;Wu et al., 2014), and the switch from a subordinate non-reproductive status to a reproductively active dominance status in male convict cichlids (Astatotilapia burtoni) (Renn et al., 2008). The increase in the relative expression of gal mRNA in breeder N. pulcher and caring N. caudopunctatus, might indicate that gal is also associated with reproductive behaviors in biparental species (Cunha-Saraiva et al., 2018).
Prolactin has long been associated with parental care (Bachelot and Binart, 2007). In uniparental male bluegill (Lepomis macrochirus), prolactin implants resulted in an increase in caring behavior while 11-ketotestosterone (11-KT) implants increased aggression towards a brood predator . Cunha and colleagues' study  revealed a trade-off between aggression and caring; a parent cannot be both caring and highly aggressive. prl1 was negatively correlated with brood care in male N. caudopunctatus, in the phase of post-spawning care, when free-swimming offspring are mobile and therefore extremely vulnerable to predators, with non-breeding conspecifics representing one of the biggest threats (Ochi and Yanagisawa, 1999). Therefore, caring males and females should be allocating most of their effort and energy to defending their offspring from predators. The reasons for higher prl1 expression in N. caudopunctatus post-spawning males, but not in females, is not known and does not seem to be linked to increased brood care as described in the bluegill . Moreover, the endocrine correlates we have identified in this study might be related to the regulation of specific gal or prl1 receptors (Bouret et al., 2000;Power, 2005). Although we don't have information about gal or prl1 receptors for the two cichlid species studied (but see Martins et al., 2014), these are obvious candidates for future research and additional studies are still needed to infer causality between the neuroendocrine system and parental behavior.
In the wild, dominant breeding N. pulcher males do not typically participate in direct caring behaviors (i.e., egg cleaning or fanning), their time is mostly spent protecting the territory (Balshine et al., 2001;Desjardins et al., 2008). Egg tending is performed by both dominant breeding N. pulcher females and non-breeding subordinate individuals . However, as our laboratory based behavioral paradigm did not fully replicate a natural N. pulcher social group, the lack of subordinate helpers might have resulted in males sharing the load, performing more direct care behaviors compared to what is typically observed in the wild. This may explain the absence of typical sex differences in brood care behavior observed in the N. pulcher of our study. In spawning parental males gal mRNA levels were significantly higher than in females, and the reasons for the negative relationship between care and gal mRNA are currently unknown. Future studies are needed to explore the function of gal in sex-specific parental care behaviors and how the dynamics of gal expression might be synchronized across the different phases of the breeding cycle.
The observation that males, but not females, in N. caudopunctatus and N. pulcher, expressed higher levels of prl1 within the first 24 hrs. after egg-laying is intriguing and conflicts with previous results that reported lower levels of prl1 in N. pulcher males compared to females (Bender et al., 2008). Several factors may explain the divergence between our results and those from the previous study. For example, Bender et al., 2008 sampled breeding N. pulcher pairs within social groups, while we sampled breeding pairs without a social group. There are also technical considerations linked to the PCR approach taken to measure gene expression in the two studies. The almost 200 times more prl1 measured in non-breeding subdominant individuals compared to the levels observed in caring parental fish may be associated with social position and stress. In support of this idea is the finding of similarity in prl1 levels in the dominant non-breeders and breeders. Similar, high prl1 mRNA levels in non-breeding subordinate individuals in comparison with breeding parents have previously been reported in stable social groups of N. pulcher (Bender et al., 2008) and other cooperatively breeding species (e.g. Meerkats, Suricata suricatta, Carlson et al., 2006;Florida scrub-jay, Aphelocoma coerulescens,Schoech, 1998;Schoech et al., 1996; Mexican Jay, Aphelocoma wollweberi, Brown and Vleck, 1998; Common marmosets, Callithrix jacchus, Roberts et al., 2001;and Harris' Hawks, Parabuteo unicinctus, Vleck and Goldsmith, 1991; but see Schradin and Pillay, 2004 for a counter example). As we did not measure the cortisol levels of the fish, we cannot rule out that the differences that we Table 3 Sex, sample sizes, Spearman rank correlation coefficients and P-values (Holm correction applied) for the behavioral responses and whole brain mRNA levels. observed in prl1 expression in the present study are due to stress. No prl1 gene expression difference existed between non-breeding and breeding biparental N. caudopunctatus. However, gal mRNA levels were higher in breeding than non-breeding N. caudopunctatus, suggesting that an increase in gal might be associated with the transition from a non-reproductive non-caring individual into a care giver (Cunha-Saraiva et al., 2018). Notwithstanding, we did find that, at the most intense and longest phase of care (at 7-days after egg laying), males had higher levels of both prl1 and gal when compared to females.
The results obtained in the present study together with our recent study (Cunha-Saraiva et al., 2019) in which post-spawning females of N. caudopunctatus had higher whole-brain bioactive arginine-vasotocin (AVT) levels compared to post-spawning males, contribute to increase knowledge of endocrine correlates of parental care (Aubin-Horth et al., 2007). AVT appears to be more important for maternal behavior and Gal appears to be important for paternal behavior in this biparental cichlid species (Butler et al., 2021), which suggests that although both parents perform brood care and defense the mechanism modulating care for Fig. 3. The prl1 mRNA expression levels for a) non-breeding subordinate Neolamprologus pulcher. Note, the y-axis is depicted in bold to highlight the difference in scale for prl1 relative expression in comparison with the other groups, b) non-breeding dominants and c) breeding individuals. The gal mRNA levels for d) nonbreeding subordinates, e) non-breeding dominants and f) breeding pairs of Neolamprologus pulcher. Black boxes represent males and white boxes represent females. Box-plots show the interquartile range (IQR) of each group analyzed with whiskers extending to 1.5 × the IQR. Horizontal lines represent medians. Outliers were included in the analysis and are depicted on the graphs as white circles. Sample sizes are depicted in the graphs. Significant differences between male and female of each sampling are depicted on the graph as bold stars. *P < 0.05. each sex may be different.

Conclusion
In conclusion, gal expression is associated with parental care and thus might be involved in the regulation of parental care in both N. caudopunctatus and N. pulcher, two common cichlid fishes endemic to Lake Tanganyika. As these species show somewhat different breeding systems, our results emphasize the importance of gal signaling in caring parents. Our results also underscore that gal expression might not only be important for the onset of caring in uni-parental species but also for biparental or cooperative species. Here the investigation of two behaviorally distinct species, has facilitated a more comprehensive understanding of the mechanisms underlying parental care.

Declaration of Competing Interests
We have no competing interests.

Fig. 4.
The negative relationship between a) prl1 and brood care in N. caudopunctatus males (black circles) but not females (white circles), and b) gal and brood care in N. pulcher males (black circles) but not females (white circles). The significant male regression line is depicted in the graph as a bold line. Regression coefficient and p-value are also depicted in the graph.