Brain Transcriptional Profiles of Male Alternative Reproductive Tactics and Females in Bluegill Sunfish

Bluegill sunfish (Lepomis macrochirus) are one of the classic systems for studying male alternative reproductive tactics (ARTs) in teleost fishes. In this species, there are two distinct life histories: parental and cuckolder, encompassing three reproductive tactics, parental, satellite, and sneaker. The parental life history is fixed, whereas individuals who enter the cuckolder life history transition from sneaker to satellite tactic as they grow. For this study, we used RNAseq to characterize the brain transcriptome of the three male tactics and females during spawning to identify gene ontology (GO) categories and potential candidate genes associated with each tactic. We found that sneaker males had higher levels of gene expression differentiation compared to the other two male tactics. Sneaker males also had higher expression in ionotropic glutamate receptor genes, specifically AMPA receptors, compared to other males, which may be important for increased spatial working memory while attempting to cuckold parental males at their nests. Larger differences in gene expression also occurred among male tactics than between males and females. We found significant expression differences in several candidate genes that were previously identified in other species with ARTs and suggest a previously undescribed role for cAMP-responsive element modulator (crem) in influencing parental male behaviors during spawning.


Introduction
Understanding the molecular mechanisms that influence variation in behavior can provide insight into how different behavioral phenotypes within populations evolve and are maintained. One important area of research on behavioral phenotypes focuses on alternative reproductive tactics (ARTs), which are found in a wide array of taxa [1][2][3][4][5]. ARTs typically consist of larger males practicing a "territorial" tactic that maintain and protect breeding territories and smaller "sneaking" males that sneak fertilizations rather than compete with territorial males [6]. The mechanisms underlying the expression of ARTs can differ significantly across species. history, mating tactics are developmentally plastic, with males apparently transitioning from the sneaker tactic to the satellite tactic as they age [35].
While the spawning behavior, reproductive success, and hormone profiles of bluegill have been studied extensively [35,[37][38][39][40][41], the genes influencing behavioral differences during spawning are less clear [42]. Thus, for this study, we used RNAseq to characterize the brain transcriptome of the three spawning male tactics (parental, sneaker, and satellite), non-spawning parental males, and spawning females to examine how differences in gene expression may relate to behavioral variation among these groups. Specifically, we aim to (1) assess whether or not there is a greater difference in gene expression profiles between tactics in different life histories (parental versus the two cuckolder tactics) than between tactics within the same life history (sneaker versus satellite), (2) identify specific gene ontology categories that are expressed for each tactic, (3), examine the expression of potential candidate genes identified from other fish species to determine if they also differentiate ARTs in bluegill, and (4) compare expression differences between male and female bluegill.

Bluegill Sampling
In June 2013, bluegill sunfish were collected via dip net from Lake Opinion near Queen's University Biological Station (QUBS), Elgin, Ontario, Canada. A total of 12 parental males, 12 sneaker males, 13 satellite males, and 12 females were collected on the same day directly from the bluegill colony while in the act of spawning. All spawning fish used in this study were behaviorally verified as to tactic by snorkelers prior to collection. An additional 12 non-nesting parental males were collected off of the colony four days prior to spawning (as determined once spawning at these colonies began) and were used as our non-spawning parental males. These males were reproductively mature but were in between spawning bouts. Individuals were euthanized using clove oil, total body length was measured, and whole brains were immediately dissected out and stored in RNAlater (Life Technologies, Carlsbad, CA). The total amount of time required for euthanasia, brain dissection, and brain storage in RNAlater was Table 1. Proposed candidate genes (from [18]) influencing teleost alternative reproductive tactics (ARTs). POA = Pre-optic area

De novo Transcriptome Assembly and Reference Transcriptome
Prior to assembly, read quality was assessed using FastQC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc). Nucleotides whose quality score was below PHRED = 2 were trimmed using Trimmomatic version 0.32 [43], following recommendations from MacManes [44]. The reference transcriptome was assembled de novo using Trinity version 2.04 [45,46]. One representative of each of the five groups (spawning parental male, non-spawning parental male, sneaker male, satellite male, and female) was used to construct a combined reference transcriptome. The five representatives selected for the reference were the individuals with the highest number of reads within their group and a total of 85 million paired-end reads were assembled. The assembly was normalized using Trinity's (version 2.04) in silico normalization program. The fully assembled transcriptome consisted of 235,547 transcripts. To determine whether this was an appropriate representation of the bluegill brain transcriptome, reads from samples not used in the assembly were mapped back to the transcriptome using Burrows-Wheeler Aligner (bwa)-mem version 0.7.12 [47], and >90% of those reads aligned, which is comparable to the rate of mapping for the individuals used in the assembly (92%). TransDecoder [45] was used to identify protein-coding regions within the assembled transcriptome. Transcripts were blasted using Blastn to a custom database containing complete coding sequences (cds) and non-coding RNA (ncRNA) from spotted green puffer (Tetraodon nigroviridis), spotted gar (Lepisosteus oculatus), southern platyfish (Xiphophorus maculatus), medaka (Oryzias latipes), Japanese pufferfish (Takifugu rubripes), West Indian Ocean coelacanth (Latimeria chalumnae), Mexican tetra (Astyanax mexicanus), zebrafish (Danio rerio), and Amazon molly (Poecilia formosa) (downloaded from Ensembl). Transcripts that contained protein coding regions or those that blasted to the customized fish database with an e-value less that 1x10 -3 comprised the reference transcriptome and this was used for read alignment and to estimate transcript counts. This reference consisted of 72,189 transcripts, including isoforms, with a mean transcript length of 2,024 bp, a N50 = 3,106 bp and a N90 = 1,018 bp.

Read Alignment and Transcript Counts
Reads from each individual were separately aligned to the reference transcriptome using bwamem 0.7.10 [47]. At least 85% of all reads from each individual mapped back to the reference, with the majority aligning 90% of reads or higher. The sequence alignment/map (sam) files were then converted to a binary format (bam) using Samtools 0.1.19 [48]. Transcript counts for each individual were obtained using the program eXpress 1.5.1 [49]. Unless otherwise indicated, all programs were run using the default options. Differential gene expression was determined using the R statistical package edgeR 3.6.8 [50]. Transcripts with cpm values of <10 for at least 4 individuals were filtered out prior to analysis, leaving 19,084 transcripts. While this filtering process is conservative, we are less likely to observe false positives due to outliers with highly variable expression, which is common for transcripts with low counts. Transcript counts were normalized to account for differences in cDNA library size among individuals and dispersion parameters were estimated using Tagwise dispersion estimates. Differences in gene expression between groups were calculated using an Exact-test for binomial distribution. Genes with p-values lower than 0.05 after false discovery rate (FDR) correction were determined to be statistically significant. All fold changes are reported as log2 fold change. Hierarchical cluster analysis to visualize overall group differences was performed on only those transcripts with FDR values below 0.05 and with log2 fold changes greater than 1.5 (equaling 1,400 transcripts) using the R package ggplot2 (2.1.0) [51].
For the transcripts that were differentially expressed among behavioral groups, enrichment analysis was conducted using a Fisher Exact test in the R Stats package (v 3.3.1) to examine whether the proportion of genes within each GO category was significantly higher than expected based upon the proportion of expressed genes assigned to that GO term within the reference transcriptome. To ensure adequate statistical power, only GO terms with at least 10 transcripts within each category were included in the statistical analysis. A FDR correction was applied to control for multiple testing and GO terms with p-values < 0.05 were considered to be significant. Visual representations of enriched GO terms were generated using REVIGO [52].

Differential Gene Expression across All Groups
Hierarchical cluster analysis of the top differentially expressed transcripts showed sneaker males grouped separately from the other male tactics (Fig 2). Satellite males tended to have expression profiles intermediate between sneakers and the other groups.
When comparing across all groups, five transcripts consistently displayed higher expression in spawning parental males compared to all other groups ( Table 2). Fourteen transcripts were differentially expressed in satellite males compared to all other groups. Expression for these transcripts in satellite males was higher compared to parental males (spawning and nonspawning) and females, but lower compared to sneaker males ( Table 2). There were 2,253 transcripts differentially expressed between sneaker males and all other groups (S2 Table). The majority of these transcripts with higher expression in sneakers were related to ion transport, ionotropic glutamate signaling pathway, and mRNA processing (Fig 3). Two transcripts were differentially expressed in females compared to the other groups and both of these were expressed at lower levels than in the other groups (Table 2).

Between Life History Comparisons
Spawning Parental Males versus Sneaker Males. A total of 9,279 transcripts were differentially expressed between spawning parental males and sneaker males (Fig 4). Of these, 4,537 transcripts showed higher expression in parental males (S3 Table) and 4,742 transcripts showed higher expression in sneaker males (S4 Table).
Enrichment analysis of GO terms associated with differentially expressed genes showed that the biological functions most enriched in parental males included translation initiation, translation elongation, proteolysis involved in cellular protein catabolism, and oxidationreduction processes (S5 Table). The 27 molecular processes most enriched in parental males Satellite Males-Fold changes are higher compared to parental males (spawning and nonspawning) and females but lower compared to sneaker males. Sneaker Males-• 2,253 differentially expressed transcripts (S2 Table) • Transcripts with high expression related to ion transport, ionotropic glutamate receptor signaling pathway, and mRNA processing (Fig 3) Females-Expression levels are lower compared to other groups compared to sneaker males included ribosomal structure, oxidoreductase activity and catalytic activity (S5 Table). Biological processes enriched with genes displaying higher expression in sneaker males included ion transport, homophilic cell adhesion, protein phosphorylation, ionotropic glutamate receptor signaling pathway, and synaptic transmission (S5 Table). The 10 molecular processes enriched in sneaker males included ion channel activity, protein binding, and ionotropic glutamate receptor activity (S5 Table).
Spawning Parental Males versus Satellite Males. A total of 1,141 transcripts were differentially expressed between spawning parental males and satellite males (Fig 4). Of these, 676 transcripts had higher expression in parental males (S6 Table) and 465 transcripts showed higher expression in satellite males (S7 Table).
One GO term related to biological function, oxidation-reduction processes, was enriched in parental males compared to satellite males (S5 Table). Six GO terms related to molecular processes were enriched in parental males (S5 Table). These were iron ion binding, two types of oxidoreductase activity, heme binding, acylCoA dehydrogenase activity and catalytic activity (S5 Table).
Only one GO term related to biological function, ion transport, was enriched in satellite males compared to spawning parental males (S5 Table). Three GO terms related to molecular processes were enriched in satellite males relative to spawning parental males. These were nucleic acid binding, ion channel activity, and GTP binding (S5 Table).

Differential Expression within Life Histories
Satellite Males verses Sneaker Males. There were 2,590 transcripts differentially expressed between satellite males and sneaker males (Fig 4). Of these, 2,480 transcripts were also differentially expressed between spawning parental and sneaker males (Fig 4) and all showed expression to be in the same direction for parental and satellite males compared to sneakers (i.e. those with higher expression in parental males compared to sneaker males were also higher in satellite males compared to sneakers). Only 110 transcripts were differentially expressed in satellite males compared to sneaker males that were not also differentially expressed between parental and sneaker males. Seventy-six transcripts had higher expression levels in satellite males (S8 Table) and 34 transcripts had higher expression in sneaker males (S9 Table). The number of transcripts differentially expressed was too low to have adequate statistical power to perform enrichment analysis for GO terms. However, many of the transcripts with higher expression in satellite males are associated with GTP catabolism, while transcripts with higher expression in sneaker males are involved in signal transduction, neural crest cell migration, and DNA integration.
Spawning Parental Males verses Non-Spawning Parental Males. A total of 137 transcripts were differentially expressed between spawning and non-spawning parental males. The majority of these transcripts (132 transcripts) showed higher expression in spawning males (S10 Table). Genes with the highest expression in spawning parental males compared to nonspawning males were MHC II beta antigen, cytosolic 5'-nucleotidase II (nt5c2), cAMP responsive element modulator a (crem), cysteine dioxygenase type 1 (cdo1), and an uncharacterized protein. Only 8 transcripts showed higher expression in non-spawning parental males. These were nuclear receptor subfamily 1 group D member 4b (nr1d4b), neuronal tyrosine-phosphoinositide-3-kinase adaptor 2 (nyap2), sphingosine-1-phosphate receptor 4 (s1pr4), gamma-aminobutyric acid A receptor beta 3 (gabrb3), and four uncharacterized proteins (S11 Table). Again, the number of transcripts assigned to GO terms was too small to have adequate statistical power to perform an enrichment analysis for this comparison.

Potential Candidate Genes Associated with ART Spawning Behavior
We observed differential expression of a number of transcripts previously identified as candidate genes associated with differences in ART behaviors (described in Table 1) ( Table 3). In our data set, the candidate genes cyp19a1b, epd, and gal showed higher expression in spawning parental males compared to sneaker males. Epd also had higher expression in satellite males compared to sneakers. Egr1 showed higher expression in both satellite and sneaker males relative to spawning parental males. Sst1 showed higher expression in satellite males compared to sneaker males, but no differences in other comparisons between tactics. No differences in expression related to gnrh, avt, or sst3 were observed between any of our groups.
Another transcript that displayed large differences in expression between spawning parental males and all other groups (including non-spawning males) was cAMP-responsive element modulator (crem) ( Table 2). Multiple isoforms of the transcript were expressed, with log2 fold changes ranging from 1.3-2.6 times higher in spawning parental males compared to other groups. Consistent with the findings for GO term enrichment, transcripts that showed the highest levels of expression in sneaker males compared to other groups were related to glutamate receptor genes, particularly AMPA ionotropic glutamate receptors (S2 Table).
In addition to the candidate genes listed in Table 1, a number of endocrine genes were differentially expressed among two or more male tactics. Among these are a number of genes that we consider bluegill candidate genes based on documented male tactic differences in circulating steroid hormone levels on the day of spawning [37,38,41]. Some of these include oxytocin, pro-melanin concentrating hormone-like, prostaglandin, and corticotropin releasing hormone receptor 2 (S2 Table, S3 Table). Further investigation of these specific hormone-associated genes is currently in progress and will be reported elsewhere. Table 3. Gene expression differences (log2 fold change) among male tactics for proposed candidate genes (see Table 1).

Sex Differences
Two transcripts were differentially expressed between females and all of the male groups (sneaker, satellite, spawning parental male, and non-spawning parental male) ( Table 2). These corresponded to galanin/GMAP prepropeptide (gal) and protachykinin (tac) and both were expressed at lower levels in females. The number of differentially expressed genes between females and satellite males was higher than between females and parental males (Fig 4), despite females and satellites exhibiting some similarity in spawning behavior. The relatively low number of differentially expressed genes observed between males and females may be due to higher variation in gene expression among females compared to males (S1 Fig). The datasets supporting the conclusions of this article are available on the Sequence Read Archive (SRA) through BioProject ID: PRJNA287763. Environmental data, RNA quality information, the assembled transcriptome, the transcript count matrix, and R code for differential gene analysis are available on Dryad (doi: 10.5061/dryad.82fd8).

Discussion
Bluegill sunfish are a classic system for examining behavioral differences in ARTs. In this study, we generated and assembled the first bluegill brain transcriptome and identified candidate genes associated with different male spawning tactics. The main differences in gene expression were found between sneaker males when compared to the two other male tactics and females. Generally, sneaker males showed higher expression in transcripts influencing neural activity, whereas parental and satellite males exhibited higher expression in genes related to translation and oxidoreductase activity. There were larger differences in transcript expression among different male tactics than between males and females.

Overall Expression Differences among ARTs
One of our main findings is that shared life history is not a driving factor influencing similarity in brain gene expression among tactics. In bluegill, parental and cuckolder life histories are fixed, but within the cuckolder life history, males transition from the sneaking to the satellite tactic as they age [32,35]. Our data showed that, regardless of whether comparisons were made across fixed (parental versus sneaker or parental versus satellite) or plastic (sneaker versus satellite) tactics, sneaker males showed the highest level of differentiation in gene expression. The expression differences in sneakers may be partially due to age and size considering sneaker males are both younger and smaller than satellite and parental males. Genes associated with increased age in other fish species, such as translation elongation and ribosomal proteins [53], had higher levels of expression in parental and satellite males compared to sneaker males in our dataset. However, age and size are not likely the only factors contributing to these differences. The differences in expression observed in this study are also likely to be reflective of behavioral differences exhibited by these tactics. For example, sneaker males in Atlantic salmon, Salmo salar, show higher expression of genes related to neural activity [15] compared to immature males of similar age and size. While we were not able to separate the effects of age, size, or behavioral tactic for our data, many of the genes with higher expression in bluegill sneakers are related to similar gene pathways (synaptic transmission) that were observed in the Atlantic salmon study. Thus, while age and size are likely responsible for some of these expression differences, they are also a reflection of different behavioral tactics and not just exclusively the result of different life histories.

Gene Categories Associated with ARTs
Identifying distinct gene categories expressed by ART types provides information regarding which functional gene categories may be associated with behavioral differences during spawning. As mentioned above, previous studies in Atlantic salmon and sailfin mollies, Poecilia latipinna, indicate that sneaker males have increased expression of genes related to neurotransmission and learning [15,17]. We found that the GO terms enriched in bluegill sneaker males compared to all other groups were the ionotropic glutamate signaling pathway and ionotropic glutamate receptor activity. Ionotropic glutamate receptors are primarily excitatory neurotransmitter receptors and play an important role in fast synaptic transmission (reviewed in [54]). Two of these receptors, NMDA and AMPA, play important roles in memory function and spatial learning (reviewed in [55]). Blocking NMDA receptors impairs learning new spatial locations in goldfish [56] and mice with impaired AMPA receptors show normal spatial learning but have impaired spatial working memory (i.e. their ability to alter their spatial choice in response to changing environments is impaired) [57]. We propose that increased expression of genes related to spatial memory, particularly spatial working memory, could be important for bluegill sneakers during spawning as they attempt to enter nests while avoiding detection by parental males and common predators around the colony [58]. Bluegill sneakers must also position themselves in close proximity to females to time sperm release to coincide with female egg release [59]. Similarly, sailfin molly sneakers, who also show enrichment in ionotropic glutamate related genes [17], probably benefit from increased spatial working memory as they position themselves by the female for quick and successful copulations. In this context, increased expression in gene pathways that improve neural function related to spatial working memory would be especially beneficial for sneaking tactics to increase their reproductive success.
While ARTs with fixed tactics maintain the same mating tactic over their lifetime, ARTs with plastic tactics can alter their behavior and, in some cases, their phenotype when switching from one tactic to another. Different phenotypes can be accomplished without altering the underlying genomic sequence through a number of mechanisms including epigenetic regulation, alternative gene splicing, and post-translational modification of proteins [60,61]. A number of genes involved in these processes showed higher expression in the plastic tactics (satellite and sneaker) compared to the fixed parental tactic (Table 2). For example, ogt plays a key role in chromatin restructuring and post-translational modification of proteins [62]. It has been also implicated in a number of different processes including nutrient and insulin signaling [63,64], sex-specific prenatal stress [65], and behavioral plasticity [66]. Genes associated with alternative splicing that were expressed at higher levels in plastic tactics included isoforms of serine/arginine-rich proteins (SR proteins), a family of proteins involved in RNA splicing [67], and CLK-4 like proteins, which are kinases that function in regulating SR protein activity [68]. Similarly, differential expression of RNA splicing genes has also been observed in two other teleost species with plastic tactics, the black-faced blenny and intermediate-sized sailfin mollies [17,18]. While the mechanisms influencing how ART males switch between tactics is currently unresolved, epigenetic regulation, alternative gene splicing, and post-transcriptional modifications could be important for plastic tactics in altering their phenotype in response to environmental or developmental cues.

Candidate Genes Associated with ARTs
A number of candidate genes have been proposed to influence the expression of ARTs in teleosts [18] (Table 1). In our study of bluegill, we corroborate some of these candidates. Similar to many other species, cyp19a1b, epd, and gal had higher expression levels in spawning parental males compared to sneaker males. Expression levels of cyp19a1b (brain aromatase) on the day of spawning initially seem contrary to what would be expected based on observed differences in circulating androgen and estrogen levels in male bluegill morphs. Estradiol (E2) and testosterone (T) levels have been shown to increase cyp19a1b expression in a number of teleosts [69,70], however 11-ketotestosterone  shows little to no effect [70]. In bluegill, sneaker males have higher circulating levels of E2 and T compared to parental males on the day of spawning, while 11-KT levels are higher in parental males during this time [41]. However, testosterone levels of parental males can peak just prior to or on the day spawning [37,71] possibly influencing the higher expression in cyp19a1b we observed.
The one candidate gene that was expressed opposite to expectations was egr1. Egr1 expression was lower in bluegill spawning parental males compared to sneaker or satellite males although previous work in cichlids found that expression of this gene increases when subdominant males transition into dominant males [30]. Egr1 is an important transcription factor involved in neural plasticity [72], so it may be one of a group of genes involved in regulating the switch from one tactic to another. Taken together, our results corroborate roles for cyp19a1b, epd, gal, and egr1 as candidate genes contributing to behavioral differences in ARTs across multiple species. Future work will explore how candidate genes are expressed across different brain regions, as some studies have found regional differences associated with genes, such as avt, in other species with ARTs [21,[73][74][75][76].
We also identified one transcript with a previously unrecognized function in influencing male spawning behavior for any teleost. Transcripts corresponding to isoforms of crem were expressed at significantly higher levels in spawning parental males compared to all other male groups, including non-spawning parental males. Crem plays a key role in modulating the hypothalamic-pituitary-gonadal (HPG) axis by regulating transcriptional responses to cAMP in neuroendocrine cells and also serves as an important activator of spermatogenesis in Sertoli cells of mice [77][78][79]. This gene can act as both transcriptional activator and inhibitor depending on the splice variant produced [77]. One splice variant is inducible cAMP early repressor (ICER), a powerful repressor of cAMP-regulated transcription [80]. ICER plays a key role in circadian melatonin synthesis by repressing the key enzyme that converts serotonin to melatonin [81]. High levels of these neurotransmitters have been associated with increased mating and cooperative behavior and decreased aggressive behavior [82][83][84]. ICER has not yet been well characterized in teleosts but one of our differentially expressed crem transcripts had a significant blast hit to an ICER variant from Epinephelus brunes (longtooth grouper). The relationship among crem, melatonin, and aggression is opposite to what would be expected if ICER is playing a role since parental males have darker pigmentation and are more aggressive than other groups [58,[85][86][87]. However, increased expression of crem, whether through ICER or another crem transcript variant, could be a candidate gene influencing behaviors associated with parental male spawning given its role in transcriptional regulation and its involvement in the HPG axis.

Sex Differences
Neural differences between the sexes are common and found in many taxa (reviewed in [88,89]). However, within ARTs, differences in neural expression profiles can often be larger among male tactics than between males and females [18][19][20]. In bluegill, only two transcripts were consistently differentially expressed in females when compared to all male groups and these corresponded to gal and tac. Gal and tac are neuropeptides and neurons expressing these genes have been associated with male sexual behavior and aggression [28,90]. Injections of gal into the preoptic area (MPOA) of the brain increase sexual behaviors in male rats [28] and stimulate both male-typical and female-typical sexual behaviors in females [91]. In male rats, testosterone can enhance the pituitary's response to galanin (endoded for by gal), which heightens gonadotropin releasing hormone's (GnRH) stimulation of luteinizing hormone. If gal is directly involved in regulating gnrh expression in bluegill, this neuropeptide may play an important role in behavioral differences between the sexes. In sequential hermaphroditic fish, surges in GnRH drive the switch from female to male [92]. Although bluegill are gonochoristic, gonadal sex is not evident until 30-60 days post hatch [93] and changes in sex can be hormonally induced [94]. Thus, gal expression, through its influence on gnrh expression, may play an important role in sex differences for this species.
The role of tac in influencing sexual behaviors in teleosts has not been addressed, but tac expression significantly increases in the brain of male eels (Anguilla anguilla) during sexual maturation [95] and leads to increased male aggression in Drosophila [90]. In bluegill, the primary role of tac expression may not be male-male aggression, considering higher expression levels of this gene are also observed in the non-aggressive satellite and sneaker males when compared to females. Although the ways in which gal and tac expression specifically influence sex-specific behaviors in bluegill is currently undefined, the fact that lower expression is consistently observed in females compared to all male groups suggests that these are important sex-specific neural genes.
In summary, our work describes differences in gene expression profiles in the brains of bluegill sunfish during spawning. The largest differences in expression levels were observed when comparing sneakers to parental males, satellite males, and females, suggesting that differences in gene expression are more related to male reproductive tactic than to life history. Consistent with other studies, our work demonstrates that sneaker males have greater expression of genes involved in neural function relative to more territorial-type males, particularly in relation to spatial working memory, as mediated by ionotropic glutamate receptors. We also found support for the previously identified candidate genes cyp19a1b, epd, gal, and egr1 contributing to behavioral differences in ARTs and identified a potential new candidate gene, crem, for regulating parental males' behavior during spawning.