DOI: 10.5061/dryad.cvdncjtc6
title: Managing the tradeoff between reproduction and survival requires flexibility in behavior and gene regulation in three-spined stickleback

description: 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.
description: <em>Study population</em> 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 <em>ad
 libitum </em>once per day. <em>Behavioral
 trials</em> 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 (<em>Oncorhynchus mykiss</em>) 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. <em>Behavioral
 analysis</em> 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. <em>RNA extraction and
 sequencing</em> 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<sup>°</sup>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 &gt;30M reads.
 Reads were aligned using STAR (Dobin et al. 2013) to the <em>G.
 aculeatus</em> 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.
  <em>Differential gene expression</em>
 Genes were filtered to include those with &gt;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 &lt; 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 &lt; 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.
 <em>WGCNA</em> 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.
 <em>Gene regulatory network analysis</em>
 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 &lt; 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.
description: # Managing the tradeoff between reproduction and survival requires
 flexibility in behavior and gene regulation in three-spined stickleback
 --- [Access this dataset on
 Dryad](https://doi.org/10.5061/dryad.cvdncjtc6) ## Description of the data
 and file structure ### File List: **A)**
 Courtship_Predation_Tradeoffs_Final.csv **B)** PredCourt_sampleinfo.csv
 **C)** counts:
 01D_AAGCTCGTGG-TCTTCAGGTC_L001_R1_001.fastq.gz_htseq_counts.txt,
 04D_TGTGCACCAA-TACCTGCCTT_L001_R1_001.fastq.gz_htseq_counts.txt,
 06D_TCGCGAAGCT-TGTTGACCGG_L001_R1_001.fastq.gz_htseq_counts.txt,
 07D_CTAGACTTCG-CCAACTCCGA_L001_R1_001.fastq.gz_htseq_counts.txt,
 11D_ACAGACCACG-CGTAACCTGG_L001_R1_001.fastq.gz_htseq_counts.txt,
 12D_CGATCTCAGG-TGGAGGTTCG_L001_R1_001.fastq.gz_htseq_counts.txt,
 13D_GGTTCCTATT-TAATTCCAGC_L001_R1_001.fastq.gz_htseq_counts.txt,
 15D_CACAAGATCC-TGTCCGGTTG_L001_R1_001.fastq.gz_htseq_counts.txt,
 17D_GGCACAACCT-CAGGAGTCTA_L001_R1_001.fastq.gz_htseq_counts.txt,
 18D_CCGGTCATAC-TTGTCTATGG_L001_R1_001.fastq.gz_htseq_counts.txt,
 20D_CTTGGTAGCA-CCAACCTATA_L001_R1_001.fastq.gz_htseq_counts.txt,
 21D_ACCGGTCGAA-CTATGTTCCA_L001_R1_001.fastq.gz_htseq_counts.txt,
 23D_GAGCTGTATA-GTACTCCTGG_L001_R1_001.fastq.gz_htseq_counts.tx,
 24D_ACCGCGGATA-CAACACCGTA_L001_R1_001.fastq.gz_htseq_counts.txt,
 25D_AAGGAAGGAA-CACGATTCCG_L001_R1_001.fastq.gz_htseq_counts.txt,
 26D_CCATCCTGTG-CGCCTTGAGA_L001_R1_001.fastq.gz_htseq_counts.txt,
 28D_GCGCGTAGTT-TAGTGACCAG_L001_R1_001.fastq.gz_htseq_counts.txt,
 29D_CTTAGCGCTG-TCCATGATCT_L001_R1_001.fastq.gz_htseq_counts.txt,
 30D_TGCAGTGCTC-AGGTGTTCTA_L001_R1_001.fastq.gz_htseq_counts.txt,
 32D_TCGAAGTCAA-CATTCCTCCG_L001_R1_001.fastq.gz_htseq_counts.txt,
 33D_CAGAAGATGG-TTCCGGAGAA_L001_R1_001.fastq.gz_htseq_counts.txt,
 34D_TCCTGAAGTG-GCCTAATTCC_L001_R1_001.fastq.gz_htseq_counts.txt,
 35D_TCACATGAGA-GGCTTGTCGA_L001_R1_001.fastq.gz_htseq_counts.txt,
 37D_ACCTGGACAA-CAGATAAGCG_L001_R1_001.fastq.gz_htseq_counts.txt,
 38D_GGCTCTACTG-AGAAGCGGAA_L001_R1_001.fastq.gz_htseq_counts.txt,
 39D_CTCCGCAGTT-TGTTCTACGG_L001_R1_001.fastq.gz_htseq_counts.txt,
 40D_AGCGGTCTTC-AAGTACCACT_L001_R1_001.fastq.gz_htseq_counts.txt,
 41D_GATCTCGTCC-CTGCGAAGAC_L001_R1_001.fastq.gz_htseq_counts.txt,
 42D_TTGTTGCAGA-GTAAGCGAGG_L001_R1_001.fastq.gz_htseq_counts.txt,
 44D_GCTCTTATTG-TGATACCAGA_L001_R1_001.fastq.gz_htseq_counts.txt,
 46D_CTTGTCGTCG-CTCGAAGGAA_L001_R1_001.fastq.gz_htseq_counts.txt **D)**
 eFDR null p-values:
 DcephalonSouthRolly_nullCourtshipDControlBROAD2023-04-09.txt,
 DcephalonSouthRolly_nullPredationCourtshipDControlBROAD2023-04-09.txt,
 DcephalonSouthRolly_nullPredationDControlBROAD2023-04-09.txt,
 T2cephalonSouthRolly_nullCourtshipT2ControlBROAD2023-10-05.txt,
 T2cephalonSouthRolly_nullPredationCourtshipT2ControlBROAD2023-10-05.txt,
 T2cephalonSouthRolly_nullPredationT2ControlBROAD2023-10-05.txt **E)** eFDR
 empirical p-values:
 DiencephalonSouthRolly_contrast_Courtship_eFDRDControlBROAD2023-04-09.txt,
 DiencephalonSouthRolly_contrast_Predation_eFDRDControlBROAD2023-04-09.txt,
 DiencephalonSouthRolly_contrast_PredationCourtship_eFDRDControlBROAD2023-04-09.txt, T2cephalonSouthRolly_contrast_Courtship_eFDRT2ControlBROAD2023-10-05.txt, T2cephalonSouthRolly_contrast_Predation_eFDRT2ControlBROAD2023-10-05.txt, T2cephalonSouthRolly_contrast_PredationCourtship_eFDRT2ControlBROAD2023-10-05.txt **F)** Code_PredCourt.R #### A) Behavioral Data: Courtship\_Predation\_Tradeoffs\_Final.csv Description of Variables * Male_ID: Unique identifier for each focal male * Treatment: Treatment the male was given. F = Courtship Opportunity, P = Predation Risk, B = Tradeoff, C = Control * date: Calendar date that data was collected * Female_bites: The number of pecks at the female cage (or randomly assigned empty cage) performed by the focal male * Predator_bites: Number of bites at the predator cage (or a randomly assigned empty cage) performed by the focal male * Female_zigzags: Number of zigzags performed by the focal male * Retreat: the number of retreats performed by the focal male #### B) Behavioral Data: PredCourt\_sampleinfo.csv * Sample Info: PredCourt_sampleinfo.csv * Sample: Unique identifier for each brain sample * Experiment: Identifier for which experiment the samples belong * Treatment: Treatment the male was given (as with Behavioral Data) * Brain_Region: The brain region the sample is from. T = Telencephalon, D = Diencephalon #### C) counts Count data for all samples #### D) eFDR\_null\_pvalues null distribution of p-values for DEG analysis #### E) eFDR\_contrast\_pvalues eFDR corrected p-values for DEG analysis #### F) Code/Software Code_PredCourt.R Code for the statistical analysis of behavioral and gene expression data. All analyses were conducted using R version 4.2.2.