Predators exert a powerful selective force, however, predator avoidance can conflict with other important activities such as attracting mates. Decisions over whether to court mates versus avoid predators are vital to fitness, yet the mechanistic underpinnings of how animals manage such tradeoffs are poorly understood. Here, we investigate the flexibility of behavior and gene regulation in response to a tradeoff between avoiding predators (survival) and courting potential mates (reproduction) in three-spined stickleback (Gasterosteus aculeatus). We compared behavioral and transcriptomic responses of male sticklebacks faced with a courtship opportunity and cues of a predator simultaneously to males faced with a courtship opportunity or cues of a predator alone and found that males behaviorally compromised courtship in favor of predator avoidance when faced with a tradeoff between them. The need to manage this tradeoff elicited dynamic changes in brain gene expression, and sets of functionally connected genes were organized into discrete modules based on co-expression. Additionally, we found that behavioral flexibility in response to tradeoffs corresponded to flexibility in gene regulatory network structure. Combined, these results uncover the coordinated response by the brain to a fundamental ecological tradeoff, providing insight into the structure and function of genetic networks underpinning how animals make fitness-influencing decisions.

Study population

Adult male and female sticklebacks were collected using minnow traps from South Rolly Lake, AK (61.667N, 150.138W) on June 1, 2022, and shipped to the University of Illinois Urbana. Upon arrival, fish were placed in groups of 15 or 30 in 10- and 25-gallon tanks, respectively, and quarantined for one week. Fish were kept on a 16:8h light/dark photoperiod and water temperatures were maintained at 62-68 degrees C. Housing tanks were given 25% water changes every other day and all fish were fed ad libitum once per day.

Behavioral trials

Behavioral trials were conducted between July 12 – July 29, 2022. Males were placed in experimental tanks (53L x 33W x 24H cm) with sand and algae provided for nesting material. Once males built their nests, they were assigned one of four treatments: courtship opportunity alone, predation risk alone, courtship opportunity and predation risk together (hereafter ‘tradeoff’), and a control. Stimuli were confined to cages made from clear plastic (7.5x5x2.5”) with 0.5cm wire mesh at each end to allow olfactory cues to permeate. For each trial, a male was presented with two cages placed 10cm from the nest at either end of the tank. In the courtship opportunity treatment, a live gravid female was confined to one of the cages. In the predation risk treatment, trout visual and olfactory cues (see below) were presented in one of the two cages. In the tradeoff treatment, a live gravid female was confined to one cage and trout cues were presented in the other. Lastly, the control consisted of two empty cages. Visual and olfactory cues of a rainbow trout predator (Oncorhynchus mykiss) were used to simulate predation risk. Rainbow trout are a common predator of stickleback, and native to South Rolly Lake (Baker et al. 2010). Trout visual cues consisted of a realistic 6” rainbow trout fishing lure (Savage Gear, Denmark, EU) with the hook removed. To simulate predator movement, during the trial, the front of the trout was lifted up and down approximately 2cm once every 10 seconds using a clear fishing line tied around the head (as in Näslund et al. 2016). Trout odor cues consisted of water collected from a raceway housing several hundred rainbow trout at Jake Wolf Memorial Fish Hatchery. Injured conspecific cues were generated by macerating eight conspecifics and soaking them in 0.6 liters of water (as in Landeira-Dabarca 2019). Both the trout odor cues and the injured conspecific cues were stored at -20C in 100ml aliquots and thawed immediately before use. The 100ml olfactory cue consisted of 50ml trout odor cues and 50ml injured conspecific cues.

Each trial consisted of a 5-minute acclimation period, after which the cages with their corresponding stimuli were introduced and the behavior of the focal male was recorded. In the predation risk and tradeoff treatments, olfactory cues were poured gently above the predator cage. To control for disturbance caused by pouring water, 100ml of tank water was added simultaneously above the female cage in the courtship opportunity and tradeoff treatments and the two empty cages in the control treatment. Predator avoidance was quantified as the number of retreats by the focal male, where the male would swim rapidly backward with dorsal spines raised (Näslund et al. 2016). Male courtship was quantified by recording the number of ‘zigzag’ displays performed by the male, a common courtship display in sticklebacks (van Iersel 1953), as well as the number of ‘pecks’ at the female cage (or a randomly assigned empty cage in the control). Because females were confined to cages, we assumed the number of pecks at the female cage reflected courtship effort by the males. We focus on the number of pecks because zigzags occurred at a low frequency relative to the number of pecks at the female cage, and rates of zig-zags and pecks were positively correlated (correlation estimate = 0.40, t = 2.318, df = 29, p-value = 0.028), suggesting that the number of pecks is a reasonable proxy for courtship behavior. There was no indication that variation in female behavior impacted male behavior; females were confined to the cage and therefore had little opportunity to behaviorally respond to male courtship.

Behavioral analysis

All analyses were conducted in R version 4.2.2 (R Core Team 2022). To investigate behavioral responses to the tradeoff between predation risk and courtship two linear models were fit using the ‘lme4’ package: 1) the number of pecks at the female (or empty control) cage and 2) the number of retreats from the predator cage, with treatment as a fixed effect. Response variables were +1 ln-transformed to fit the assumptions of linear models. Tukey post-hoc tests were performed on significant treatment effects using the ‘emmeans’ package.

RNA extraction and sequencing

One hour after the trial, males were sacrificed (n = 8 after a courtship opportunity, n = 8 after predation risk, n = 8 after a tradeoff, n = 7 after a control) and their telencephalon and diencephalon were rapidly dissected and deep frozen at -80^°C until extraction (as in Barbasch et al. 2023). RNA was extracted using Invitrogen PureLink RNA extraction kits (Invitrogen Corporation, Carlsbad, CA) and quality was checked using Agilent Bioanalyzer chips. Library prep and sequencing were performed by the Functional Genomics Unit of the W.M. Keck Center (University of Illinois Urbana Champaign). Libraries were prepared using the ‘Kapa Hyper Stranded mRNA library kit’ (Roche) and pooled, quantitated by qPCR, and sequenced on one S4 lane on a NovaSeq6000 to a depth of >30M reads. Reads were aligned using STAR (Dobin et al. 2013) to the G. aculeatus reference genome (Ensembl release 95; Jones et al. 2012; Nath et al. 2021; Peichel et al. 2020). Gene transcript counts were generated using HTSeq (Anders et al. 2015). Counts from the diencephalon and telencephalon samples were analyzed separately. Out of 62 samples, we removed two samples from the telencephalon (one ‘Control’ and one ‘Predation Risk’) due to overlap with the diencephalon samples on a PCA plot, potentially because not all of the diencephalon tissue was cut away while dissecting the smaller telencephalon.

 Differential gene expression

Genes were filtered to include those with >0.5 CPM in four or more samples, resulting in 18302 and 18588 genes in the telencephalon and diencephalon, respectively. Filtered read counts were TMM (timed mean of M-values) normalized using a tagwise dispersion estimate and used to perform a GLM to identify the set of differentially expressed genes in the three experimental treatments relative to the control, with an empirical FDR P-value correction (as in Barbasch et al. 2023). Briefly, sample names were permuted 500 times to generate a null distribution of p-values and compared to the actual p-values. A cut-off of eFDR < 0.01 was used. Differentially expressed genes (DEGs) were annotated with ‘biomaRt’ and ‘TopGo’ was used to identify genes and biological processes recruited in response to each treatment relative to the control. Here, we included only the DEGs unique to a specific treatment, to test whether these distinct genes represent shared or distinct biological processes.

Next, to investigate the overlap between the unique DEGs associated with managing the tradeoff between courtship and predation risk (this study) and the unique DEGs associated with managing the tradeoff between courtship and territorial defense (Barbasch et al. 2023), we used a hypergeometric overlap test comparing gene lists between the two experiments in the diencephalon and telencephalon. Significantly enriched GO terms (p-value < 0.05) for the DEGs unique to each experimental treatment (Courtship Opportunity, Predation Risk, and Tradeoff) were identified using Fischer’s exact tests and the ‘parentchild’ algorithm (Grossmann et al. 2007), which takes into account parent-child relationships among terms.

WGCNA

To identify modules of co-expressed genes and how they relate to individual behavior and treatment, weighted co-expression networks were built using WGCNA (Langfelder and Horvath 2008) with the genes that passed filtering in the two brain regions. Briefly, signed networks were built using a bi-weight mid-correlation function with a soft thresholding power chosen according to where the scale-free topology index curve flattened upon reaching 0.9 (power = 11 for the telencephalon, 5 for the diencephalon). The minimum block size was set to 30 and similar modules were merged using the dynamic tree cut method, which resulted in 23 modules in the telencephalon and 40 modules in the diencephalon.

First, to identify general patterns organizing functionally connected genes into co-expression modules, hypergeometric overlap tests were performed using the sets of DEGs from each experimental treatment (including overlapping genes) and the sets of genes composing each module. Overlapping genes were included to identify modules potentially enriched for DEGs associated with multiple treatments, as well as modules enriched for DEGs associated with a single treatment. A Bonferroni-corrected p-value threshold of 0.1 was used to assess significance. To investigate the biological processes associated with the genes composing these modules, we performed GO analysis on the modules significantly enriched for the DEG lists as described above. Next, to find modules whose expression patterns differed between treatments, we fit a linear model with module eigengene expression as the response variable and treatment as the independent variable. Empirical p-values were calculated as above, with sample names permuted 10000 times. We used an eFDR threshold of 0.1 and Tukey post-hoc tests were used to make pairwise comparisons for models with a significant treatment effect.

Lastly, to understand the relationship between module expression and individual behavior, linear models were fit with module eigengene expression as the response variable. To investigate female-directed behavior, we used the number of pecks at the female cage in the Courtship Opportunity and Tradeoff treatments as the independent variable. To improve power, we included treatment as a co-variate only if it significantly improved model fit, as determined using likelihood ratio tests. Similarly, to investigate predator avoidance behavior, we used the number of retreats in the Predation Risk and Tradeoff treatments, controlling for treatment where treatment significantly improved model fit. We used an eFDR threshold of 0.1 to control for multiple comparisons and Tukey post-hoc tests were used to make pairwise comparisons.

Gene regulatory network analysis

To investigate patterns of transcriptional regulation across our experimental treatments, gene regulatory networks were inferred from expression data using GENIE3 (Huynh-Thu et al. 2010). GENIE3 quantifies the strength of regulatory interactions using random forests to predict the expression pattern of each ‘target’ gene based on the expression of a set of ‘regulators.’ As regulators, we used the putative transcription factors for stickleback from the Animal Transcription Factor Database (Shen et al. 2023). A hard threshold of 0.017 was used to select the connections with the strongest importance values, resulting in a network of 1036 TFs with 27132 connections in the telencephalon and 1071 TFs with 13390 connections in the diencephalon. Transcription factors (TFs) whose targets were significantly enriched for the unique DEGs from each experimental treatment were identified using a Bonferroni-corrected hypergeometric overlap test (FDR < 0.1), including TFs with at least five targets to improve power to detect regulators involved in decision-making by omitting TFs with very few target genes. Here, as with the GO analysis, we used the sets of DEGs unique to each treatment, allowing us to identify distinct regulators involved in the response to tradeoffs, as well as test whether there are differences in the structure of gene regulation when animals are faced with multiple demands compared to only one. To test the hypothesis that the need to manage tradeoffs is associated with greater flexibility in GRN structure, characterized by fewer TF-TF edges, for each treatment we calculated the TF-TF edge ratio as the number of target genes for each TF that also encode TFs divided by the total number of target genes for that TF (as in Sinha et al. 2020). To understand network-level differences between our treatments, we used a p-value cutoff of 0.05, which was not corrected because we were comparing TFs across treatments to make inferences about TFs that were enriched in one treatment but not others. Thus a positive false discovery rate is not appropriate. We compared the ratios across treatments using Kruskal-Wallis tests and post hoc comparisons were made using Dunn’s tests.