Personal and transgenerational cues are nonadditive at the phenotypic and molecular level.

Organisms can gain information about their environment from their ancestors, their parents or their own personal experience. 'Cue integration' models often start with the simplifying assumption that information from different sources is additive. Here, we test key assumptions and predictions of cue integration theory at both the phenotypic and molecular level in threespined sticklebacks (Gasterosteus aculeatus). We show that regardless of whether cues about predation risk were provided by their father or acquired through personal experience, sticklebacks produced the same set of predator-adapted phenotypes. Moreover, there were nonadditive effects of personal and paternal experience: animals that received cues from both sources resembled animals that received cues from a single source. A similar pattern was detected at the molecular level: there was a core set of genes that were differentially expressed in the brains of offspring regardless of whether risk was experienced by their father, themselves or both. These results provide strong support for cue integration theory because they show that cues provided by parents and personal experience are comparable at both the phenotypic and molecular level, and draw attention to the importance of nonadditive responses to multiple cues.

Recent evolutionary theory seeks to understand how cues from ancestors, parents and personal experience are integrated together to produce adaptive phenotypes[1-6]. The central problem is that organisms in natural populations must decide how and whether to attend to cues from different sources, and those sources might not always agree with each other. For example, an animal might obtain cues from their father that the environment is safe, while personal experience suggests otherwise. Recent theory identifies the conditions that favour the evolution of reliance on some sources of cues over others[1-6], and highlights the importance of cue reliability during cue integration. The relative weight given to a cue depends on its accuracy as a predictor of selective conditions in the future[3]. For example, a cue might not give entirely reliable information on current conditions, and/or the cue might give information on current conditions but the environment might change during the interval between cue detection and when selection acts on the phenotype[7,8].

Evidence that different sources of cues (e.g. genetic and environmental) trigger similar adaptive phenotypic responses[9] provides support for information integration theory, but a key assumption of several cue integration models concerns the way that organisms respond to cues from different sources that are in agreement with each other. Several models start with the simplifying assumption that cues from different sources are additive[1-3]. Under this assumption, additional information increases an individual’s confidence in its assessment of the environment, which results in a linear relationship between the number of sources of consistent cues and the adaptive phenotype[3]. For example, assume a wide range of anti-predator phenotypes available to a developing individual and that greater elaboration of those phenotypes confers greater fitness benefits[10]. An individual receiving cues from its parent that the environment is dangerous might begin to develop anti-predator phenotypes. In an additive model, if personally-acquired cues confirm that the environment is dangerous, then the individual will further develop those phenotypes[1-3], but if personally acquired cues indicate that the environment is safe, the individual will stop developing those phenotypes.

However, there are also several reasons to expect that organisms receiving consistent cues from different sources will respond in a nonadditive manner. For example, there might be underlying constraints (epistasis, fundamental biochemical or biophysical constraints) that limit the most extreme phenotypes. Nonadditivity is also expected for threshold traits, i.e. when a single source of cues is sufficient to push a phenotype past a threshold[11,12]. Another possibility is that if organisms integrate cues in a Bayesian fashion, i.e. they update personal information by continuously sampling their environment[4,5], then they might not respond to a personally-acquired cue if it is consistent with their strong prior expectation, i.e., that was set by their evolutionary history or their parents[5]. Alternatively, additional cues might disproportionally increase the individual’s confidence in the state of the environment, causing a multiplicative effect on the phenotype. Finally, nonadditivity is expected when the absence of cues provides an unreliable assessment of the environment. For example, imagine two different sources that provide highly reliable cues about predation risk, and both sources indicate the same level of risk. If the absence of cues about predation risk is unreliable – perhaps because predators come and go – then organisms might be better off always strongly responding to cues of predation risk, even if they are only from a single source[13]. This scenario might be especially likely to occur when the costs of failing to respond to cues about risk is high, or even deadly (“smoke detector principle”[14]). In contrast, additive responses might be more likely to occur in response to environmental information that is not as immediately threatening, such as weather, food availability, etc.

Here, we investigate the independent and combined influence of personal and paternal experience with danger at both the phenotypic and molecular level in threespined sticklebacks. Specifically, using a 2×2 factorial experiment with a split-clutch design (Figure 1), we explore how juvenile sticklebacks combine personally- and paternally-acquired cues about predation risk. In this species, parental care is necessary for offspring survival, males are the sole providers of parental care for embryos and offspring for approximately two weeks, and the way fathers behave toward their offspring influences offspring phenotypic development[15-17].

Fig. 1

Experimental design. The effects of personal and paternal experience with risk on offspring phenotypes were compared in a 2×2 factorial design. Fathers either were (N = 5) or were not (N = 5) exposed to predation risk while they were caring for their offspring (paternal experience: unexposed versus exposed). Within each family, siblings either were or were not personally exposed to predation risk as juveniles (personal experience: unexposed versus exposed). Juveniles were then measured for either brain gene expression (no paternal cue/no personal cue: N = 7, paternal only: N = 10, personal only: N = 10, both: N = 9) or behaviour (no paternal cue/no personal cue: N = 6, paternal only: N = 10, personal only: N = 10, both: N = 9).

Adult males (one year of age) were randomly assigned to either a predator-exposed or control (unexposed) treatment. While males were providing care for their offspring, fathers in the “predator-exposed” group were chased by a model sculpin predator for two minutes (unexposed: not chased)[18]. Sculpin are a fish predator that primarily prey on stickleback nests and juveniles[19]. At two months of age, half of the offspring within each family were chased by a model sculpin predator for one minute a day for seven days (personal experience: exposed), while the other half of the family was undisturbed (personal experience: unexposed).

This design resulted in four different conditions: offspring that were not exposed to risk and whose fathers were also unexposed, offspring that were not exposed to risk but whose fathers were exposed, offspring that were exposed to risk but whose fathers were unexposed and offspring that were exposed to risk and whose fathers were also exposed (Figure 1). At three months of age, offspring were measured for size, weight, latency to emerge from a refuge (timidity), and brain gene expression using RNA-Seq. We infer that differences between offspring of predator-exposed versus unexposed fathers reflect transgenerational plasticity, while differences between predator-exposed and unexposed offspring reflect developmental plasticity. We investigated additivity by comparing offspring with both personal and paternal experience with risk to the three other conditions.

RESULTS AND DISCUSSION

When fathers were exposed to predation risk while they were caring for offspring, they decreased parental behaviour (Supplementary Figure 1, consistent with [17, 18, 20]). This behavioural shift suggests that fathers can provide cues to their offspring through their behavioural interaction with them, similar to the way that mothering influences the behavioural development of offspring in mammals[21]. Juvenile offspring of predator-exposed fathers were relatively small, had lower mass for a given length, and took more time to emerge from a refuge compared to juvenile offspring of unexposed fathers (Supplementary Table 1; Figure 2), consistent with a previous study on sticklebacks[17], and with other studies on both evolved and developmental response to risk in small fishes[22-25]. It is possible that offspring of predator-exposed fathers had these phenotypes because they received less fanning (oxygen) from their fathers, which caused altered growth patterns during embryonic development. As these phenotypes align with anti-predator phenotypes arising from selection and from developmental plasticity[17,22-25], it is unlikely that they are due to poor parenting from fathers, and instead might reflect adaptive anti-predator phenotypes.

Fig. 2

The effect of personal and paternal experience with predation risk on offspring phenotypes was nonadditive. Box plots indicate median, interquartile range (IQR), and 1.5*IQR at both the upper and lower ranges (whiskers). Dots indicate raw data points. There was a significant interaction between personal and paternal experience on (A) standard length (linear mixed model; F1,20.86 = 5.08, p = 0.035); (B) body mass relative to length (linear mixed model; F1,23.60 = 10.23, p = 0.004); and (C) latency to emerge from a refuge (linear mixed model; F1,25.83 = 11.79, p = 0.002).

In general, the phenotypes of offspring with personal experience with predation risk resembled the phenotypes of offspring whose fathers had been exposed to predation risk (Figure 2). It is possible that the personal experience of being chased by the model sculpin caused offspring to hide more and forage less, again resulting in smaller, more timid phenotypes[13]. These results support the hypothesis that regardless of its source, cues about risk cause sticklebacks to produce a similar set of predator-adapted phenotypes. Moreover, the combined influence of personal and paternal experience on body size and timidity was nonadditive: offspring that received cues about risk from two sources were statistically indistinguishable from offspring that received cues about risk from a single source (Figure 2). In general, offspring of predator-exposed fathers had lower body mass relative to length compared to the control group. Personal experience with risk by itself strongly decreased body mass relative to length. Interestingly, personal experience with risk combined with paternal experience with risk appeared to attenuate the negative effects of personal experience with risk by itself on body mass relative to length.

One possible explanation for these nonadditive patterns is that they reflect constraints on the maximum phenotype that can be produced in response to cues about risk. For example, it might not be possible to be much smaller or have lower weight relative to body size and still function. There might also be a constraint imposed by the tradeoff between foraging and predation risk that limits timidity: an animal can only hide in the refuge for so long before eventually venturing out to feed[13]. The results could also be consistent with a threshold model: once a certain threshold of information about the environment is reached, one of only a few alternate states is induced[11,12], perhaps because there are few benefits to having an intermediate phenotype.

Another potential explanation for the nonadditive patterns is that sticklebacks combined cues from their fathers and their personal experience in a Bayesian fashion[5]. In this population, fathers are likely to have highly reliable information about the extent to which sculpin are likely to be a threat to their offspring. Fathers have opportunities to perceive visual and/or olfactory cues of sculpin without being threatened themselves because sculpin tend to specialize on juveniles soon after they emerge[19], before juveniles have had time to sample their environment. Under this Bayesian scenario, after receiving highly reliable cues from their fathers, offspring in this experiment maximally produced anti-predator phenotypes, but additional cues (based on personal experience) that also indicated that sculpin were present did not provide any additional information about predation risk to those subjects. Similarly, when offspring were chased by a model sculpin for several days, this provided highly reliable cues that sculpin were present, and this information over-rode the effects of unreliable cues from their father indicating that predation risk was low, and those offspring also maximally produced anti-predator phenotypes. Indeed, because sticklebacks are a prey species highly vulnerable to predation[26], they might be better off responding to a false alarm than not responding at all (the “smoke detector principle”[14]). Our results suggest that once a response is triggered in response to paternal information indicating that the environment is dangerous, it remained “on”, perhaps because the costs of reversal were higher than the costs of failing to respond to an unpredictable, but potentially deadly threat.

Other studies that have examined how organisms integrate information from their parents and personal experience have also found that responses to multiple cues tend to be nonadditive[27-34], but the precise nonadditive pattern is variable across studies. For example, personal and parental responses to cues of predation risk are synergistic in snails, such that snails only mounted a phenotypic response when they received cues of predation risk from both sources[27]. Another recent study found that phenotypic responses to both personal and maternal experience with food availability was highly variable among clones of Daphnia[28]. An important consideration is that different types of patterns are likely to be expected in studies where the environment simply acts as a cue, e.g. cues of predation risk, versus in studies where the environment also influences state, e.g. food availability. A challenge for theory is to incorporate experiences that not only act as cues but also affect state.

Personal and parental experiences also produced similar responses[9] at the molecular level: there was a core set of genes that were differentially expressed in the brain in response to risk, regardless of whether the risk was experienced by fathers, their offspring, or both (Figure 3B), and the number of shared genes between the three pair-wise contrasts is greater than expected due to chance (Shared genes across all treatments: 208; hypergeometric test: p < 1e–10). Moreover, the brain gene expression pattern of the core set of genes was remarkably concordant (Figure 3C). The brain gene expression profile of offspring with both personal and paternal experience with predation risk resembled the brain gene expression profile of offspring that independently received either source of information on its own. These results suggest that for this core set of genes, both sources of information trigger the same response at the molecular level, and that personally and paternally-acquired information share some “equivalence” at the molecular level. While West-Eberhard[9] discussed “equivalence” in the context of the exchangeability of genetic and environmental effects, our findings suggest that the same concept applies to different environmental effects acting over different timescales (transgenerational versus developmental). This is in contrast to a study in Daphnia, which found few similarities between personal experience and maternal experience at either the phenotypic or molecular level[35], highlighting the need for future work to examine patterns of information integration across organisms with differing life histories, sensory inputs and development.

Fig. 3

Brain gene expression responses to personal experience with risk, paternal experience with risk. (A) Comparisons. The brain gene expression profiles (RNA-seq) of offspring in response to personally and paternally-acquired information about risk was compared relative to a control group of offspring that did not receive information about risk from either source (“double control”). The brain gene expression pattern of offspring of unexposed fathers was compared between offspring with and without personal experience with risk. This pair-wise contrast represents developmental plasticity genes (purple). The brain gene expression profile of offspring without personal experience risk, but whose fathers did experience risk, was compared to the double control. This pair-wise contrast represents transgenerational plasticity genes (blue). The brain gene expression profile of offspring with both personal and paternal experience with risk was compared to the double control. This pair-wise contrast includes both developmental and transgenerational plasticity, as well as their interaction (green). (B) Number of differentially expressed genes in each pairwise contrast, along with the number of overlapping genes between contrasts. The size of each circle is proportional to the number of genes. (C) Heat map showing the differential expression pattern of the 208 genes that were common to all three contrasts. Red=upregulated, purple=downregulated. Columns represent pairwise contrasts, rows represent genes. Note that genes that were upregulated in the brain in response to paternal information were upregulated in response to personal information and were also upregulated in animals that received information from both sources, and vice versa. The direction of regulation is more congruent than expected by chance (χ2 = 60.84, n=208, p<0.00001). The full gene lists and their functional enrichments are in Supplementary Tables 3 and 4, respectively.

Although the overlap between developmental and transgenerational plasticity at the molecular level was much greater than expected due to chance, there were also sets of genes that were unique to the different forms of plasticity. There were, for example, 322 genes that were differentially expressed in response to paternal experience with risk, but were not differentially expressed in response to personal experience with risk. Given the common response to personal and paternal cues about risk at the phenotypic level, it is tempting to speculate that the shared genes reflect the similar “output” in response to cues about risk from different sources, while the unique genes reflect differences in the “input” between developmental and transgenerational plasticity, i.e. whether the cue was acquired via paternal behavior versus from the experience of being personally chased by the model predator. The large number (n=425) of genes that were unique to the “both” comparison could reflect a variety of different mechanisms involved in weighing, processing and integrating cues. A previous study in stickleback found that maternal stress altered offspring brain gene expression in a sex-specific fashion[36]. An informal comparison of the gene lists suggests that there is little commonality between the genes associated with paternally-mediated transgenerational plasticity in our study and maternally-mediated transgenerational plasticity in [36], suggesting different molecular mechanisms responsive to cues from fathers versus mothers. As mothers do not provide care or interact with their offspring after fertilization, mothers and fathers provide different cues about environmental conditions to offspring. Future studies explicitly comparing cues from both parents may help resolve whether and how stickleback integrate cues from mothers and fathers differently. Given the effects of personally- and paternally-acquired information on nonbehavioral traits (e.g. body size), it would also be interesting for future studies to examine how cues from different sources are “read” by the genome in peripheral tissues and at different developmental timepoints.

Cue integration models offer a fresh framework for understanding why developing organisms sometimes pay more attention to their genes, their parents or their own personal experience to produce adaptive phenotypes. Key assumptions and predictions of these models are beginning to be empirically tested by studies that simultaneously manipulate cues from different sources[27-35]. Our study provides strong empirical support at both the phenotypic and molecular level for this theory[1-6] and suggests that future models should explore the consequences of relaxing the assumption of additivity.

METHODS

Study population and breeding

Adult threespined stickleback (approximately 1 year of age) were collected from Putah Creek, a dammed, regulated freshwater stream in northern California, in April 2013. Sculpin (Cottus spp), a fish predator known to prey on stickleback eggs, fry, and adults[19] are present at this site. Fish were shipped to the University of Illinois at Urbana-Champaign, and males were introduced into separate 9.5L (36 × 21 × 18 cm) tanks with a refuge (plastic “plant”), an open plastic box (13 × 13 × 3 cm) filled with fine sand, and filamentous algae for nest building. Following nest completion, males were presented with a gravid female and allowed to spawn. A previous study showed that there was no effect of previous breeding experience or previous experience with predation risk while breeding on subsequent paternal behavior[18]. Each male spawned with a unique female. After spawning, the female was removed. Fish were kept at 20 degrees Celsius on a summer (16L:8D) photoperiod in freshwater. Water was cleaned via a recirculating flow-through system that consists of a series of particulate, biological, and UV filters (Aquaneering, San Diego, USA). 10% of the water volume in the tanks was replaced each day. Fish were fed a mixed diet consisting of frozen bloodworm, brine shrimp and Mysis shrimp in excess each day. Experiments were carried out in accordance with institutional guidelines (University of Illinois IACUC protocol #15077). Animals were collected under a California Fish and Game Collecting permit #SC-3310 to AMB.

Exposing fathers to predation risk and recording paternal behaviour

A total of 20 males were randomly assigned to either the “unexposed” or “predator-exposed” treatment (N = 10 unexposed, N = 10 predator-exposed). The first five males from each treatment group to complete clutches were used in this experiment (N = 5 unexposed, N = 5 predator-exposed). Predator exposure did not increase the likelihood of a male’s nest failing. On the third day after males spawned (when the embryos were three days old), males in the “predator-exposed” treatment were chased with a 10 cm rubber model sculpin (Jewel Bait Company) for two minutes to simulate a nest predation attempt, as in [18]. Model predator exposure occurred at 11AM CST. A predator of this size is a threat to the eggs and fry, but not to the adult males[19]. Previous research has shown that male stickleback adjust their parenting behaviour in response to this predator model[17,18]. At this developmental stage, the optic cups of the embryo are still developing[37] and the eggs were covered by nesting material, thereby reducing the possibility of direct embryonic exposure to predation risk. For males in the “unexposed” treatment, we removed the top of the tank and gently splashed the water when the eggs were three days old to simulate the water disturbance caused when the model predator entered the tank. This splashing did not cause males to alter their paternal behavior[18].

After spawning, paternal behaviour was observed every day for ten minutes between 1000 and 1200 CST from one day after spawning through five days after the eggs hatched (when fry naturally disperse in this population). Eggs hatched on day 5 following fertilization (Supplementary Figure 1). We measured the total amount of time the male spent fanning his eggs. Fanning is a paternal behaviour that oxygenates the eggs[38], is important for offspring development[38], and consistently varies among fathers[18,20]. The simulated predation threat (or water splashing in the unexposed treatment) occurred after the daily observation of paternal behaviour. Five days after the eggs hatched, males were removed from the tank.

Exposing offspring to predation risk

Once fry were approximately one cm in length (at around one month of age), each full sibling family was evenly divided into two separate tanks and randomly assigned to either “unexposed” or “predator-exposed” treatments. Offspring were fed newly-hatched Artemia nauplii shrimp in excess each day until they reached three cm in length, at which time they were fed the adult slurry of frozen food.

At two months of age, juveniles in the predator-exposed treatment were briefly exposed to risk once a day for seven days. Specifically, they were chased with a 10 cm model sculpin for one minute at a random time each day (between 1000 and 1400 CST), once a day for seven days, to minimize the potential for habituation. For juveniles in the unexposed treatment, we removed the top of the tank and gently splashed the water once per day for seven days.

Offspring measurements and behaviour

At three months of age, we collected a subset of juveniles (N = 2 per treatment per family) and quickly measured standard length and body weight. Due to differences in clutch size and offspring mortality, the final sample sizes were no paternal cue/no personal cue: N = 7, paternal only: N = 10, personal only: N = 10, both: N = 9. We regressed length on weight and analyzed the residuals to obtain a measure of weight relative to length. We then euthanized the juveniles via rapid decapitation and flash froze the heads and bodies in supercooled ethanol (-110°C) for later RNA extraction. At this time we also removed the caudal fin and stored it in 70% ethanol for later determination of genetic sex using a male-specific genetic marker[39].

For behavioural testing of predator responses, another subset of juveniles (N = 2 per treatment per family) were measured at five months of age. Juveniles were transferred individually to an observation tank in an opaque cylinder (10 cm height, 10 cm diameter) plugged with a cork. After a 15-minute acclimation period, we removed the cork remotely and recorded latency to emerge from the refuge. Juveniles were returned to their home tanks following behaviour assays. Due to differences in clutch size and offspring mortality, the final sample sizes were no paternal cue/no personal cue: N = 6, paternal only: N = 10, personal only: N = 10, both: N = 9.

Phenotypic data analysis

We analyzed phenotypic data (length, mass relative to length, and latency to emerge from a refuge) using linear mixed models (LMMs). All models included paternal treatment (predator-exposed, unexposed), offspring treatment (predator-exposed, unexposed) and offspring sex as fixed effects, and father ID as a random effect. Analyses were conducted with R version 3.2.2[40]. LMMs were performed using the lmer function from the “lme4”[41] and “lmerTest”[42] packages. We used REML estimation and a diagonal covariance structure for our models, with Satterthwaite approximation for degrees of freedom. We determined whether levels of fixed factors differed from one another using Tukey’s HSD test.

RNA extraction and RNA-seq

Individuals for brain gene expression profiling were gently netted directly from their home tanks and rapidly decapitated with sharp scissors. Heads were flash frozen and stored at −80°C until dissection. We first scraped the skull with rongeurs to expose brain tissue. Heads were placed in RNALater for 24 hours at 4°C. We then dissected whole brains in RNAlater (Thermo Fisher Scientific) on dry ice and extracted RNA using the PicoPure RNA Isolation Kit with optional DNase treatment (Thermo Fisher Scientific).

Library Preparation

Poly-A RNA was enriched from 1–2 μg of total RNA by using Dynabeads Oligo(dT)25 (Life Technologies), following the manufacturer’s protocol. Two rounds of poly(A) enrichment were performed with a final elution in 14μL of water. The poly-A–enriched RNA was used to prepare RNAseq libraries, using the Illumina TruSeq kit (Illumina). Manufacturer’s instructions were followed and 13–15 cycles of PCR amplification were performed depending on the starting input of total RNA. All samples were barcoded, libraries were quantified on a Qubit fluorometer using the dsDNA High Sensitivity Assay Kit (Life Technologies), and library size was assessed on a Bioanalyzer High Sensitivity DNA chip (Agilent). Libraries were pooled and diluted to a final concentration of 10 nM. Final library pools were quantified using real-time PCR, using the Illumina compatible kit and standards (KAPA) by the W. M. Keck Centre for Comparative and Functional Genomics at the Roy J. Carver Biotechnology Centre (University of Illinois). Single-end sequencing was performed on an Illumina HiSeq 2500 instrument by the W. M. Keck Centre for Comparative and Functional Genomics at the Roy J. Carver Biotechnology Centre (University of Illinois). Not all individuals yielded enough RNA, therefore the final sample sizes for RNA-seq informatics were no paternal cue/no personal cue: N = 5, paternal only: N = 9, personal only: N = 9, both: N = 7.

RNA-seq Informatics

FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess the quality of the reads. Adaptor sequences and low quality bases were clipped from 100 bp single-end sequences using Trimmomatic. RNA-seq produced an average of 34 million reads per sample. We aligned reads to the Gasterosteus aculeatus reference genome (the repeat masked reference genome, Ensembl release 75), using TopHat (2.0.8)[43] and Bowtie (2.1.0)[44]. On average 26 million reads aligned to the genome that translate to ~76% alignment rate (Supplementary Table 2). Reads were assigned to features according to the Ensembl release 75 gene annotation file (http://ftp.ensembl.org/pub/release-75/gtf/gasterosteus_aculeatus/).

Defining differentially expressed genes (DEGs)

HTSeq-Count[45] was used to count reads mapped to gene features using stickleback genome annotation. Any reads that fell in multiple genes were excluded from the analysis. One sample from the transgenerational plasticity treatment group was excluded from the analysis based on high variability on MDS plot (Supplementary Figure 2), resulting in a final sample size for the transgenerational plasticity treatment group of N = 8. We included genes with at least 0.5 count per million (cpm) in at least five samples. Cpm values were log transformed and were analyzed using limma voom[46], a program which allowed us to control for the effect of Father). To assess differential expression, we fit a linear model ~Sex + Treatment and performed pairwise comparisons among Treatment levels to find differentially expressed genes (DEGs) due to fathers, offspring, or both exposure to predators relative to individuals who themselves nor their fathers had seen the predators. We also controlled for family by including father identity as a random factor. For false discovery rate (FDR) correction we used the “global” method in limma decideTests functionality (Limma user guide section 13.3), which adjusts p-values from all contrasts at once. An FDR cutoff of < 0.05 was used to call for differentially expressed genes (Supplementary Table 3).

To test for reproducibility of the results, we randomly permutated our sample labels 250 times and generated an empirical null distribution of coefficients by fitting a same model using limma voom. A permutation-based p-value was generated for each gene by comparing the observed model coefficient with the permutated ones (Supplementary Figure 3). A statistically significant overlap was observed between DE identified by limma voom alone and permutation tests, which suggests that our results were not biased by comparing the three experimental conditions to the same “double” control condition.

The significance of the pattern of congruent gene expression of the core set of genes was assessed with χ2 tests in each sex, where 25% of DEGs within each sex are expected to show a congruent pattern by chance alone.