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-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. In some cases, tactics are fixed for life (fixed tactics) [6] and can represent distinct life histories. Fixed tactics can occur through either inherited genetic polymorphisms [7-10], condition-dependent switches that are triggered prior to sexual maturation [1,6,11], or a combination of genetic and environmental factors [12,13]. In other cases, individuals can exhibit different tactics throughout their reproductive life, either as they grow or in response to changing social or environmental context (plastic tactics or status-dependent tactics) [1,4,6,14]. Advances in sequencing technology, such as RNA sequencing (RNAseq), now allow behavioral ecologists to explore how variation in gene expression contributes to behavioral variation among mating tactics and examine if genes associated with these behaviors differ across species with ARTs.

Next-generation sequencing has led to more in-depth research into the molecular mechanisms driving ARTs [9,15-20]. For example, development of independent (territorial) males and two alternative tactics, satellite males and female-mimicking (faeder) males in a shorebird (the ruff, Philomachus pugnax) is driven by a supergene resulting from a chromosome inversion that contains 125 genes potentially influencing ART traits [9,10]. However, due to the lack of reference genomes for most teleosts, much of the work on ARTs in this group has focused on examining differential gene expression to identify genes associated with these tactics. These studies have found a large number of genes that vary among tactics in expression in the brain during mating. In the ocellated wrasse (Symphodus ocellatus), 1,048 genes were differentially expressed when comparing sneakers to two other male tactics (nesting and satellite) and to females [19]. In the black-faced blenny (Tripterygion delaisi) and peacock blenny (Salaria pavo), RNAseq identified approximately 600 transcripts differentially expressed within the brains of ‘sneaker’ versus other male tactics [18,20]. In another study, approximately 2,000 transcripts were differentially expressed between intermediate-sized sailfin molly (Poecilia latipinna) males performing courtship behaviors compared to small sneaker males [17].

With the increase in genomic studies examining differences among ARTs, there are a growing number of candidate genes associated with these tactics. Schunter et al. [18] proposed a list of potential candidate genes based on a number of studies (Table 1). Many of these genes are involved in hormone regulation and vertebrate mating behavior, and differences in expression levels have been observed among mating tactics in different fish species. For example, the product of the cyp19a1b gene is aromatase B, the key enzyme responsible for the conversion of androgens to estrogens within the brain of vertebrates [e.g., 24,31]. Higher levels of cyp19a1b brain expression have been observed in territorial males compared to sneaker males in the peacock blenny [23], black-faced blenny [18], and an African cichlid (Astatotilapia burtoni) [16] but higher levels have been observed in sneaker male plainfin midshipman (Porichthys notatus) [25]. As more data become available, the number of candidate genes in this list will likely increase and evaluating gene expression across teleosts will aid in determining whether similar molecular pathways drive ART behaviors across different species.

Table 1

Proposed candidate genes (from [18]) influencing teleost alternative reproductive tactics (ARTs).

POA = Pre-optic area

Proposed Candidate Genes	Function	Relationship to ARTs
Arginine vasotocin (avt)	Non-mammalian homolog of vasopressin. Activates some aspects of sexual behavior	↑ in posterior POA of territorial cichlid males, but ↑ anterior POA of non-territorial [21]; ↓ density of avt mRNA in POA in parental blenny males [22]
Gonadotropin releasing hormone (gnrh)	Regulates release of luteinizing hormone and follicle-stimulating hormone from the pituitary gland	↑ in territorial cichlid males [16]
Cytochrome P450 family 19, subfamily A, polypeptide 1 (cyp19a1)	Brain aromatase. Key enzyme in estrogen biosynthesis	↑ in territorial cichlid males [16]; ↑ territorial blenny males [23]; ↑ territorial black-faced blenny males [18]; ↓ in the sonic motor nucleus and ventromedial nucleus of nesting type I (territorial) male plainfin midshipman compared to type II (sneaker) males [24,25]
Ependymin (epd)	Glycoprotein associated with neuroplasticity and neuronal regeneration. Also affects aggression levels in zebrafish [26]; Associated with stress in trout [27]	↑ in territorial cichlid males [16]; ↓ in subordinate trout males [26]
Galanin/GMAP prepropeptide (gal)	Neuropeptide that influences neurotransmitters. Associated with male sexual behaviors [28] and parental care [29]	↑ in territorial cichlid males [16]
Somatostatin (sst)	Neuropeptide that regulates endocrine pathways. Also affects neurotransmitters	↑ in territorial blenny males [18]; ↑in territorial cichlid males [16]
Early growth response 1 (egr1)	Transcription factor that influences neural plasticity	↑ when subdominant cichlid males switch to dominant [30]

One of the best-studied vertebrate species with male ARTs is the bluegill sunfish (Lepomis macrochirus). Male bluegill have two distinct life histories: parental and cuckolder. In Lake Opinicon (Ontario, Canada), all bluegill tactics spawn within large breeding colonies. Parental males are part of the parental life history and mature at around seven years of age (Fig 1). These males construct nests, court females, and provide care to young [32]. Males in the cuckolder life history become reproductively mature around two years of age [32]. Initially these males use a “sneaking” tactic (i.e., sneakers) to dart in and out of nests within the colony to cuckold fertilizations while parental males and females are spawning. As they grow, sneakers transition into a “satellite” tactic and take on female-like coloration and behaviors [32, 33]. Satellite males use this female mimicry to enter a parental male nest and cuckold fertilizations [34]. The parental and cuckolder life histories are fixed–once a male adopts the parental or cuckolder life history, he remains in that life history [35]. However, within the cuckolder life history, mating tactics are developmentally plastic, with males apparently transitioning from the sneaker tactic to the satellite tactic as they age [35].

Fig 1

Alternative reproductive tactics of bluegill sunfish.

This figure has been modified from Gross and Charnov [35] and Neff and Knapp [36]. Ages are based on data from Gross and Charnov [35] but may differ among populations [33].

While the spawning behavior, reproductive success, and hormone profiles of bluegill have been studied extensively [35, 37–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.

Materials and Methods

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 under 5 minutes. Brains remained in RNAlater at 4°C for 24 hours and were then transferred to fresh cryovials, flash frozen, and kept in liquid nitrogen until they were transported on dry ice to the University of Western Ontario. Samples were then stored at -80°C until total RNA extraction. The Animal Care Committee at Western University (UCC) approved all procedures performed in this study (AUP # 2010–214).

Total RNA Extraction

Total RNA was extracted from whole brains using a standard Trizol (Life Technologies, Carlsbad, CA) extraction protocol (https://tools.thermofisher.com/content/sfs/manuals/trizol_reagent.pdf). Total RNA was submitted to the London Genomics Center at the University of Western Ontario and quality was assessed using a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Four individuals from each group (spawning parental males, non-spawning parental males, sneaker males, satellite males, and females), for a total of 20 individuals, were submitted to the Michigan State University Research Technology Support Facility—Genomics Center for cDNA library construction and sequencing. Individuals used for this study had RIN (RNA Integrity Number) values ranging from 9.2–9.9.

cDNA Library Construction and Sequencing

The cDNA libraries were constructed for each individual using Illumina TrueSeq Stranded mRNA Library Preparation Kits LT (Illumina, San Diego, CA), with each individual receiving a uniquely identifiable index tag. The quality of each library was evaluated and the 20 individuals were multiplexed into a single sample that was subsequently run on two lanes of an Illumina HiSeq2500 Rapid Run flow cell (v1). Sequencing was performed on paired end 2 x 150 bp format reads and bases were called using Illumina Real Time Analysis software (v1.17.21.3). Reads from each individual were identified based on their unique index tag, separated, and converted to fastq files using Illumina Bcl2fastq v1.8.4. Sequencing produced an average of 14.5 million reads per individual, with over 90% of the reads having a Q-score >30.

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 bwa-mem 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].

Gene Annotation and Enrichment Analysis

For gene annotation, all transcripts were blasted using the program Blastx against a custom-assembled fish protein database. This database consisted of Ensembl protein databases of 13 different fish species: Amazon molly (Poecilia formosa), zebrafish (Danio rerio), Mexican tetra (Astyanax mexicanus), Atlantic cod (Gadus morhua), West Indian Ocean coelancanth (Latimeria chalumnae), Japanese pufferfish (Takifugu rubripes), sea lamprey (Petromyzon marinus), medaka (Oryzias latipes), southern platyfish (Xiphophorus maculatus), spotted gar (Lepisosteus oculatus), three-spined stickleback (Gasterosteus aculeatus), green spotted pufferfish (Tetradon nigroviridis), and Nile tilapia (Oreochromis niloticus). Blast hits with e-values less than 1x10-10 were considered significant. All annotated transcripts used for differential expression analysis are listed in S1 Table. Ensembl IDs from the blast hits were then converted into GO term identifiers using Biology Database Network (bioDBnet) (http://biodbnet.abcc.ncifcrf.gov/db/dbFind.php).

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].

Results

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.

Fig 2

Heatmap of transcripts differentially expressed in at least one group comparison.

Only transcripts with a log2 fold change of 1.5 or greater are included in the heatmap, representing 1,400 transcripts. Count values were averaged within each group and are scaled by row. Sneak = sneaker males, Sat = satellite males, Parent = parental males, Fem = females, NS_P = non-spawning parental males.

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 non-spawning) 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).

Table 2

Differentially expressed transcripts associated with each male mating tactic and females.

Sunfish Focal Group	Differentially Expressed Genes
Parental Males (Spawning)—Expression levels are higher in parental males compared to other groups, except for MHC class 1 antigen.	• Pancreatic progenitor cell differentiation and proliferation (ppdpf): FC = 1.8• Potassium voltage-gated channel, Isk-related family, member 4 (kcne4): FC = 1.4• Cysteine dioxygenase type 1 (cdo1), 3 isoforms: FC = 1.7• cAMP-responsive element modulator (crem), 2 isoforms: FC = 2.1• MHC class 1 antigen: FC = -5.5
Satellite Males—Fold changes are higher compared to parental males (spawning and non-spawning) and females but lower compared to sneaker males.	• Serine/arginine-rich splicing factor 4-like (srsf4): FC = 1.4 (other groups)/-0.7 (sneaker)• Arginine/serine-rich protein 1 (rsrp1), 2 isoforms: FC = 0.8/-0.6• CLK4-associating serine/arginine rich protein (clasrp): FC = 0.8/-0.8• RNA binding motif protein, X-linked (rbmx): FC = 0.9/-0.9• Dual specificity protein kinase CLK4-like (clk4), 2 isoforms: FC = 0.9/-0.8• Ultraconserved element locus (57322): FC = 1.0/-0.9• Cat eye syndrome chromosome region, candidate 2 (cecr2): FC = 1.6/-1.0• Luc7-like protein 3-like (luc7l3): FC = 0.8/-0.8• SET domain, bifurcated 1 (setdb1): FC = 0.8/-0.8• SUZ12 polycomb repressive complex 2 subunit (suz12): FC = 1.4/-1.2• O-linked N-acetylglucosamine (GLCNAc) transferase (ogt): FC = 0.7/-0.8• Serine/arginine-rich splicing factor 3-like (srsf3): FC = 1.1/-0.9• RNA-binding protein 25-like (rbm25): FC = 0.7/ -0.7• Uncharacterized protein: FC = 1.5/-0.9
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	• Protachykinin-like (tac): FC = -1.6• Galanin/GMAP prepropeptide (gal): FC = -1.7

FC = Mean log2 fold change across comparisons. Positive numbers indicate expression levels that were higher in focal group compared to other groups, negative numbers indicate expression is lower in focal group. For satellite males, the first FC value is the mean log2 fold change of satellite males compared to spawning, non-spawning parental males, and females. The second FC value is satellite males compared to sneaker males. When transcripts had multiple isoforms, FC values were averaged across isoforms.

Fig 3

Biological process GO terms enriched by genes with higher expression in sneaker males compared to all other groups.

Boxes of similar color are grouped into the same GO term hierarchy. Box size reflects the–log10 p-value of the GO term.

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).

Fig 4

UpSet graph displaying the number of differentially expressed transcripts among groups.

For each comparison, the number on top of each bar represents the number of differentially expressed transcripts (intersection size). The comparison in question is indicated by the dots or connected dots below its respective bar. Set size indicates the total number of transcripts for each comparison.

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 oxidation-reduction processes (S5 Table). The 27 molecular processes most enriched in parental males 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 non-spawning 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.

Table 3

Gene expression differences (log2 fold change) among male tactics for proposed candidate genes (see Table 1).

		Comparison between Male Tactics (Log2 Fold Change)
Proposed Candidate Gene	Isoform ID	Parent vs Sneak	Parent vs Sat	Sat vs Sneak	Spawn Parent vs NonSpawn
Arginine vasotocin (avt)	c34708_g2_i1	0.45 (0.32)	-0.98 (0.09)	0.54 (0.33)	0.74 (0.50)
Gonadotropin releasing hormone (gnrh)	c63124_g1_i1	0.76 (0.50)	0.32 (0.87)	0.44 (0.77)	0.77 (0.88)
Cytochrome P450 19a1b (cyp19a1b)	c48084_g2_i1	0.93 (0.0002)	0.64 (0.06)	0.28 (0.40)	0.39 (0.58)
Ependymin (epd)	c44195_g1_i5	1.54 (1.4 x 10–8)	0.66 (0.07)	0.89 (0.007)	-0.51 (0.45)
Galanin/GMAP prepropeptide (gal)	c41071_g5_i2	1.12 (0.0001)	0.53 (0.91)	-0.59 (0.10)	0.09 (0.97)
Somatostatin 1 (sst1)	c30013_g1_i1	0.53 (0.15)	-0.39 (0.49)	0.93 (0.03)	0.27 (0.88)
Somatostatin 3 (sst3)	c46547_g6_i1	0.001 (1.00)	-0.25 (0.54)	0.25 (0.48)	0.15 (0.90)
Early growth response 1 (egr1)	c37907_g1_i1	-0.74 (0.02)	-0.91 (0.03)	0.16 (0.72)	-0.63 (0.42)

Values in brackets represent p-values after false discovery rate correction. Values in bold are significant at p < 0.05. Parent = parental male, Sneak = sneaker male, Sat = satellite male, NonSpawn = non-spawning parental male.

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.

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 (11-KT) 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–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-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-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–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.

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-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.

Multi-dimensional space (MDS) plot based on the biological coefficient of variation (bcv) among bluegill male ARTs and females.

(PDF)

Click here for additional data file.

Annotated transcripts used for differential gene expression.

(TXT)

Click here for additional data file.

Transcripts differentially expressed in sneaker males compared to all other groups.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill parental males compared to sneaker males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill sneaker males compared to parental males.

(TXT)

Click here for additional data file.

Biological process and molecular function GO terms that are significantly enriched with genes differentially expressed between tactics.

(XLSX)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill parental males compared to satellite males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill satellite males compared to parental males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill satellite males compared to sneaker males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in bluegill sneaker males compared to satellite males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in spawning parental males compared to non-spawning parental males.

(TXT)

Click here for additional data file.

Transcripts with significantly higher expression in non-spawning parental males compared to spawning parental males.

(TXT)

Click here for additional data file.