Paternal effects in a wild‐type zebrafish implicate a role of sperm‐derived small RNAs

While the importance of maternal effects has long been appreciated, a growing body of evidence now points to the paternal environment having an important influence on offspring phenotype. Indeed, research on rodent models suggests that paternal stress leaves an imprint on the behaviour and physiology of offspring via nongenetic information carried in the spermatozoa; however, fish have been understudied with regard to these sperm‐mediated effects. Here, we investigated whether the zebrafish was subjected to heritable influences of paternal stress by exposing males to stressors (conspecific‐derived alarm cue, chasing and bright light) before mating and assessing the behavioural and endocrine responses of their offspring, including their behavioural response to conspecific‐derived alarm cue. We found that after males are exposed to stress, their larval offspring show weakened responses to stressors. Small RNA sequencing subsequently revealed that the levels of several small noncoding RNAs, including microRNAs, PIWI‐interacting RNAs and tRNA‐derived small RNAs, were altered in the spermatozoa of stressed fathers, suggesting that stress‐induced alterations to the spermatozoal RNA landscape may contribute to shaping offspring phenotype. The work demonstrates that paternal stress should not be overlooked as a source of phenotypic variation and that spermatozoal small RNAs may be important intergenerational messengers in fish.


| INTRODUC TI ON
The ancestral environment is already well-established as a factor which can shape the development of offspring phenotypes, and the existence of parental effects is now suspected to have a large bearing on the long-term responses of organisms to selective pressures (Räsänen & Kruuk, 2007). The contribution of the maternal environment has long been appreciated, for instance during gestation (Feil & Fraga, 2012;Lupien, McEwen, Gunnar, & Heim, 2009) or via the nutritional and endocrine composition of the oocyte (Best, Kurrasch, & Vijayan, 2017;Davis & Sandman, 2010;Redfern et al., 2017). Although relatively less studied, the paternal environment too can have important influences on offspring development in species with high levels of paternal investment postfertilization (Mcghee & Bell, 2014;Otero-Ferrer et al., 2020). Meanwhile, the sperm epigenome and its role in paternal nongenetic inheritance currently represent a highly active area of research (Champroux, Cocquet, Henry-Berger, Drevet, & Kocer, 2018;Chen, Yan, & Duan, 2016).
Stress can be defined as a state of threatened homeostasis induced by real or perceived danger, necessitating a suite of adaptive physiological and behavioural responses (stress responses) to restore it (Barton, 2002;Godoy, Rossignoli, Delfino-Pereira, Garcia-Cairasco, & de Lima Umeoka, 2018). Neuroendocrine pathways involving glucocorticoid signalling typically (but not exclusively) underlie these responses, although the functions of glucocorticoids are far from limited to mediation of stress responses (MacDougall-Shackleton, Bonier, Romero, & Moore, 2019). The perceived threat of predation can be considered a form of stress for prey species which rely on cues to detect nearby predators, as exposure to these cues alone is often sufficient to induce behavioural and physiological responses reflective of a neuroendocrine stress response. For instance, several species of teleost fish release a chondroitin-based alarm cue from the skin upon damage (Stensmyr & Maderspacher, 2012) which, when detected, induces behavioural responses and elevated cortisol in conspecifics (Mathuru, 2016;Stephenson, 2016).
Although stress responses are inherently adaptive in the short term, overactivation of these responses over a prolonged period, known as chronic stress, can lead to a maladaptive state through overproduction of glucocorticoids (McEwen, 2004). Aside from regulating the neuroendocrine stress response, glucocorticoids are also important in the regulation of foraging behaviour, energy mobilization and use and the circadian rhythm (Rao & Androulakis, 2019;Vitousek, Taff, Hallinger, Zimmer, & Winkler, 2018). Excessive glucocorticoid secretion can disrupt these functions and shift many physiological parameters away from their normal homeostatic ranges, leading to dysregulation (Maestripieri & Hoffman, 2011). Although previously regarded as a phenomenon inherent to industrialized society, it has become clear that the chronic stress concept is applicable in the context of wild animals and their fitness. For example, chronic stress induced by elevated predation threat can impact the reproductive success of prey species through physiological disruption and can thus influence population dynamics even in the absence of direct killing (Clinchy, Sheriff, & Zanette, 2013).
Stress, particularly in the context of predation, can be a prominent source of maternal and paternal effects, as documented in a range of species with high levels of parental investment postfertilization including fish (Bell, McGhee, & Stein, 2016;Ord, Holmes, Holt, Fazeli, & Watt, 2020), reptiles (Bestion, Teyssier, Aubret, Clobert, & Cote, 2014) and mammals (Chen et al., 2018). However, evidence for effects of paternal stress mediated by changes to the spermatozoa before fertilization has been limited to studies on laboratory rodents.
For instance, male mice stressed by maternal separation in early life sire offspring with reduced anxiety-like behaviours and glucose production (Gapp et al., 2014), while males exposed to chronic stress throughout spermatogenesis sire offspring with reduced corticosterone production (Rodgers, Morgan, Bronson, Revello, & Bale, 2013).
Meanwhile, evidence for sperm-mediated nongenetic inheritance in fish has been increasingly documented in recent years, albeit not explicitly in the context of stress. For instance, rearing environment induces specific changes to salmon sperm DNA methylation which appear to be inherited by the offspring (Rodriguez Barreto et al., 2019), while sperm-mediated paternal effects on phenotype have been demonstrated in zebrafish in response to increased sperm competition (S. Zajitschek, Hotzy, Zajitschek, & Immler, 2014) and paternal social status (Zajitschek, Herbert-Read, Abbasi, Zajitschek, & Immler, 2017).
The zebrafish has long justified its place in the pantheon of vertebrate models in translational research, particularly given its short life cycle, ease of manipulation and relative genetic similarity to humans. Core components of the neuroendocrine stress response pathways are highly conserved between mammals and fish, both of which, for instance, employ cortisol as the principle glucocorticoid regulator of these pathways (Barton, 2002). Readily measurable anxiety-like behavioural patterns influenced by glucocorticoid signalling are also similar between mammals and fish, such as edge preference behaviour, termed thigmotaxis (Steenbergen, Richardson, & Champagne, 2011). Furthermore, the spermatogenic cycle in adult zebrafish is remarkably rapid, lasting only six days from the onset of meiosis (Leal et al., 2009), allowing for experiments concerning sperm-mediated effects on the offspring to be carried out in a short timeframe. Therefore, zebrafish readily lend themselves to the study of paternal stress effects.
Although several possible epigenetic mechanisms to explain sperm-mediated paternal effects have been postulated, the transmission of noncoding RNA, such as microRNA (miRNA), in spermatozoa has garnered the most experimental evidence (Chen, Yan, & Duan, 2016). MiRNAs, which are 18-24 nucleotides in length, have the capacity to alter gene expression via the incorporation into RNA-induced silencing complexes (RISC) and subsequent cleavage of target transcripts. Remarkably, phenotypic effects of paternal stress have been shown to be recapitulated by injection of either total spermatozoal RNA from stressed males (Gapp et al., 2014) or specific stress-associated spermatozoal miRNAs into normal zygotes (Rodgers, Morgan, Leu, & Bale, 2015). MiRNAs are abundant in teleost fish (Babiak, 2014) and have important functions in early zebrafish development (Wienholds et al., 2005). It is therefore possible that paternal miRNA and/or other noncoding RNAs delivered in spermatozoa could also influence developmental outcomes in fish.
Although they could also be maladaptive, changes to the RNA consignment or other components of the spermatozoal epigenome could facilitate intergenetational adaptive responses to precondition offspring if necessitated by changes to environmental conditions during the parents' lifetime (Immler, 2018). Such soft inheritance has been suggested to provide a "buffer" which protects underlying genetic variation against environmental changes which occur too quickly for adaptation to occur via conventional mutations and in doing so could even facilitate later genetic adaptation (O'Dea, Noble, Johnson, Hesselson, & Nakagawa, 2016). However, intergenerational effects are predicted to be adaptive only when the parental environment is a good predictor of offspring environment (Guillaume, Monro, & Marshall, 2016).
Here, we hypothesized that paternal stress could induce measurable influences on stress response physiology and behaviour in the offspring of zebrafish and that such alterations would be linked to alterations in the small RNA composition of the paternal sperm.
We exposed males to stressors, which included alarm substance derived from the skin of conspecifics, and measured the responses to stressors in their offspring, specifically thigmotaxis (edge preference behaviour, which we took to be reflective of stimulus avoidance), and cortisol production. We further applied small RNA sequencing of paternal spermatozoa to investigate whether any intergenerational effects could be linked with differential expression of small RNAs putatively delivered to the offspring at fertilization. Although miR-NAs are the most well-known class of small RNA to be implicated in paternal effects, we also investigated whether the spermatozoal composition of other small RNA classes was subjected to stress-induced alterations, specifically PIWI-interacting RNAs (piRNAs) and tRNA-derived small RNAs (tsRNAs).

| Animals and housing
Mature male London wild-type (LWT) zebrafish (aged 6-10 months) were selected from healthy stocks, reared in tanks on a recirculatory system with water heated to 26°C and kept on a 12:12-hr light/dark cycle. Animals were fed brine shrimp or flake food twice daily.
Adult males were housed in individual compartments (30 × 15 × 25 cm) of the recirculatory system for at least one week prior to the onset of the 12-to 14-day treatment period and for the entire duration thereof. Five individuals in both control and stress groups comprised one experimental batch. A total of six experimental batches were carried out. Table S1 gives the number of adult male zebrafish used in each batch of the experiment.

| Alarm substance extraction
Alarm substance was derived from mature LWT zebrafish (indiscriminate of sex), using a previously described method .
Briefly, for every 3 ml of extract, five fish were euthanized, and 7-10 lacerations were made to the epidermis on both sides of each fish.
All five fish were then placed in a single 50-ml tube with 3 ml water and gently shaken to allow the alarm substance to seep out. The water containing the extract was then eluted, incubated at 95°C for 16 hr, centrifuged at 10,000 × g for 10 min to remove debris and filtered through a microfilter.

| Paternal chronic stress
Stress treatments were administered during a 12-to 14-day treatment period. To administer stress, each male was removed from its home tank into a separate exposure tank (17 × 11 × 12 cm) once per day and administered one of three randomly allocated stressors: 1: alarm substance (200 µl in 600 ml water) for 20 min, 2: chasing with a small net in 5× 1-min episodes over a 20-min period or 3: exposure to bright light via a 10w LED tube light suspended above the tank for 20 min. Individuals were not visible to each other during exposures. Following alarm substance exposure, animals were moved to separate "washing" tanks to remove traces of alarm substance before returning them to their home tanks. Males exposed to stressors were classed as paternal chronic stress (PCS) animals. Control animals were handled in an identical fashion but did not receive stress treatments.

| Mate pairing, embryo collection and rearing
Two days after the final stress treatment, male zebrafish were paired with mature naïve females randomly selected from healthy stock. Dishes of marbles (12 cm diameter × 5 cm height) were placed in the compartments to provide a substrate for egg laying and to protect embryos from cannibalism. Eggs were collected the following morning, and dead eggs were removed. Twenty-four hours later, 1 day postfertilization (DPF) embryos from each male were placed in 5 cm petri dishes with 15 ml aquarium water in densities of 13-30 for assessment of larval stress response phenotypes. All embryo dishes were incubated at approx. 26°C under a 12:12 light cycle.

| Larval exposure to alarm substance and physical stress
5-DPF offspring were tested for behavioural and physiological responses to stress. Two sets of offspring from each male were housed in petri dishes, and one was exposed to alarm substance (50 µl), while the other was given the same volume of aquarium water as a blank. Petri dishes were placed on a platform of frosted glass which was illuminated from below with LED lights to maximize contrast between fish and background. Dishes were video-recorded from above using a Panasonic HC-X920 digital camcorder beginning from 5 min prior to alarm or control treatment. A smaller number of sets of control and PCS larvae were exposed to physical stress, which was established by stirring the petri dish water by hand with a Pasteur pipette for 1 min (~2 rotations per second).

| Larval cortisol measurements
Fifteen minutes after exposure to stressors, larvae were immediately euthanized by immersion in ice water. For cortisol extraction, euthanized larvae were collected into 1.5-ml microfuge tubes and stored at −20°C until cortisol extraction. Pools of larvae were homogenized in 2-ml tubes in 120 µl distilled water using a pellet mixer. Cortisol was ethyl acetate-extracted from homogenates and quantified by ELISA, both following a previously published protocol (Yeh, Glöck, & Ryu, 2013).

| Larval group thigmotaxis
Thigmotaxis (edge preference behaviour) in response to alarm substance was recorded directly from videos using automated tracking software (Viewpoint® ZebraLab) which calculated the average number of animals present in a defined central zone (covering most of the dish area apart from the outer 6mm or approx. 1.5 larval standard lengths; Figure 1, inset) for every 10-s time interval.
These values were subtracted from the total number of animals to calculate the percentage of animals in the peripheral zone (% thigmotaxis). These values were then summarized into 1-min time bins (i.e. the mean of six values). As the average thigmotaxis response appeared to reach a peak at around 5 min after exposure and subsequently began a gradual decrease, we considered only the first 10 min after exposure in the statistical analyses (i.e. 5 min either side of the observed peak).

| Statistical analyses of phenotype data
All statistical analyses were carried out in "R" 3.6.1 (R Core Team, 2018). We used linear mixed-effects models fit by restricted maximum likelihood (REML) with the lme4 package version 1.1.21 (Bates, Maechler, Bolker, & Walker, 2015) because of the hierarchical structure of the experimental design and the fact that experiments were carried out in batches. In all models, "Batch" was included as a random term. "Father ID" was included as a random effect (nested within batch) as the analyses included multiple samples from the same parent. "Dish ID" was included as a random effect (nested within Father ID) in the case of thigmotaxis as observations were made at multiple time points. "Assay" was included as a separate random effect in the case of cortisol analyses due to samples being run over three plates. Significance of fixed effect terms was evaluated using F tests and t tests with Kenward-Roger approximation of degrees of freedom, derived from summary tables obtained using the ANOVA() and summary() functions from the stats package in conjunction with lmertest version 3.1.1 (Kuznetsova, Brockhoff, & Christensen, 2017) and pbkrtest version 0.4.7 (Hakekoh & Hojsgaard, 2014) packages.
Post hoc pairwise t tests of model terms were computed using the emmeans() function from the emmeans package version 1.4.3 F I G U R E 1 Visual summary of larval behavioural and physiological stress responses. (a) Thigmotaxis (% of animals within 1.5 standard lengths of the dish edge) in 1-min time periods for control (circles, solid lines) and paternal chronic stress (PCS) offspring (triangles, dashed lines) exposed to alarm substance (black) or a "blank" (water; white) at 5 days postfertilization. The dashed vertical lines demarcate the period during which the alarm substance (or blank) was added to the dishes (between minutes 5 and 6). Inset: schematic representation of the measurement of thigmotaxis. (b) Cortisol levels in control and PCS offspring exposed to alarm substance (grey) or blank (white). (c) Cortisol levels from control and PCS offspring exposed to a stirring stimulus (grey) or blank (white). Each point represents a pool of 13-30 animals. Stars represent p < .05, derived from post hoc pairwise comparisons of mixed-effects models by estimated marginal means [Colour figure can be viewed at wileyonlinelibrary.com] (Lenth, 2019). To approximate normal distributions, the thigmotaxis response variable was subjected to a Gaussianization transformation using the Gaussianize() function from the lambertw package version 0.6.4 (Goerg,2016), while cortisol response variables were log-transformed. Details of the models, including model structure and sample sizes can be found in Table S2. The R code for the models is also provided in the Supporting information. Plots of phenotype data were generated using ggplot2 version 3.2.1 (Wickham, 2011).  Table S3 for the sperm RNA sampling scheme.

| RNA-seq library preparation and sequencing
Small RNA sequencing libraries were prepared using the SomaGenics (Santa Cruz) Real-Seq Biofluids kit following the manufacturers' instructions, using a minimum of 2.5 ng of input total RNA. Twenty PCR cycles were performed to amplify the reverse transcription product. The quality of the prepared libraries was examined using the Agilent Hi-sensitivity DNA chip prior to the individual libraries being pooled to a concentration of 4 nM for Illumina sequencing.
Sequencing was carried out on an Illumina HiSeq 2,500 using sbs v3 chemistry and single end 65 base reads, producing on average 7.59 ± 2.22 million (mean ± SD) reads per sample.

| miRNA alignment and quantification
For miRNAs, reads in the length range of 18-24 nt were retrieved and aligned to miRbase mature zebrafish miRNA sequences using bowtie2 version 2.3.4.1 (Langmead & Salzberg, 2012) allowing no mismatches in the seed region but allowing the default maximum of two mismatches in the remainder of the alignment. Primary alignments were counted using samtools version 1.7 (Li et al., 2009) with a mapping quality (MAPQ) filter of 10. To account for possible alternative miRNA isoforms derived from differential clipping of hairpin precursors (isomiRs), reads not aligning to mature miRNAs were recovered and aligned to known hairpin sequences retrieved from miRbase using the same Bowtie2 parameters. Primary alignments to hairpin sequences were also counted in the same manner, as we did not seek to identify specific isomiRs.

| piRNA alignment and quantification
For piRNAs, reads in the 25-32 nt range were aligned to predicted zebrafish piRNA cluster sequences obtained from the piRNA cluster database (https://www.smallrnagroup.uni-mainz.de/piCdb/; Rosenkranz, 2016). Alignment was carried out using Bowtie2 as above, and up to three alignments for each read were reported. As the MAPQ score produced by Bowtie2 is influenced by multiple alignments but we were interested only in cluster-level expression (and therefore not concerned with within-cluster multimapping), we did not consider MAPQ when counting alignments. Instead, alignments were filtered to retain those with alignment scores (AS)>/= −11 (permitting one mismatch or one read gap of up to two bases), and for each read, the alignment with the highest AS was retained. Where reads had multiple alignments tied for the highest AS, these reads were kept if each alignment was within the same cluster and discarded if they aligned to multiple clusters. Multiple alignments were then assigned fractional values prior to counting the number of alignments to each cluster, such that the sum of all values counted equalled the total number of reads retained.

| tsRNA reference generation, alignment and quantification
Transfer RNA-derived small RNAs (tsRNAs) are typically derived from the 5' or 3' ends of mature tRNAs. Therefore, to obtain a reference database of putative spermatozoal tsRNAs, we pooled together 18-40 nt reads (which did not align to miRNAs/hairpins or piRNAs in previous steps) from eight biological samples (approx. 36 million reads) and aligned them (Bowtie2) to the first (5'ends) and last (3'ends) 40nt of mature tRNA sequences in the genomic tRNA database (Chan & Lowe, 2009). Perfect match primary alignments (i.e. no substitutions or indels; AS = 0) were extracted and further filtered to include only those that aligned to the edges of the reference sequences using custom R scripting. Unique sequences were considered as putative tsRNAs if they were represented by at least 4 reads. Due to many putative tsRNAs being redundant (i.e. could be completely encapsulated within another qualifying sequence) and so could not be pinpointed to a specific tsRNA, a clustering approach was adopted. Redundant sequences were removed, and nonredundant sequences were subjected to clustering using cd-hit-est from the cd-hit software package version 4.8.1 (Li & Godzik, 2006) with a similarity threshold of 90%.
Putative tsRNA sequences belonging to the same cluster were concatenated into single FASTA entries with each constitutive sequence separated by four Ns. Subsequently, the final reference FASTA contained both tsRNA sequence clusters and individual, nonredundant tsRNA sequences (i.e. "single-member" clusters).
18-40 nt reads from individual samples were then aligned to the putative tsRNA cluster reference using Bowtie2 as described for miRNAs. As a non-negligible portion of reads were expected to have multiple alignments within clusters, up to three alignments were reported per read which were subsequently filtered and counted using the same method described for piRNA clusters. were derived based on the trended dispersions, and statistical comparisons were performed using generalized linear models followed by likelihood ratio tests. We allowed results with an FDR-corrected p-value of < 0.1 and a log2-fold change > 1 or < −1 to be considered differentially expressed. The TMM-normalized CPM values of all differentially expressed small RNAs were converted to a single set of Z-scores which were plotted as a heatmap with hierarchical clustering using the pheatmap() function from the pheatmap package version 1.0.12 with default parameters (https://cran.r-proje ct.org/web/ packa ges/pheat map/index.html).

| MiRNA target analysis
Lists of predicted targets of differentially expressed mature miR-NAs were obtained using the miRmap web interface (Vejnar, Blum, & Zdobnov, 2013)

| Offspring of males exposed to stressors have compromised stress responses
Both control and PCS offspring exhibited increased thigmotaxis in response to alarm cue (Figure 1a), but the response in the PCS larvae was markedly weaker, the overall mean dropping by 3.6 ± 2.5 (SEM) percentage points. A mixed-effects model (Table S4) revealed an overall significant effect of alarm substance on larval thigmotaxis (F test of linear mixed model with Kenward-Roger approximation of degrees of freedom, F 1,36 = 28, p < .001), but no significant overall effect of PCS (F 1,32 = 1.61, p = .21) or larval alarm x PCS interaction, that is no effect of PCS on the larval response to alarm substance (F 1,36 = 2.1, p = .16). However, post hoc pairwise comparisons revealed that in the presence of alarm substance, thigmotaxis was significantly lower for PCS larvae (T 60.6 = 2.14, p = .036).
As expected, both alarm substance and the physical stressor, stirring, induced a significant increase in cortisol in control offspring F 1,10 = 3.1, p = .11). Cortisol in response to stirring was markedly lower in PCS compared to control larvae, although the difference was not significant (T 13 = 1.19, p = .08). Table S5 shows the results of the mixed-effects model.
PCS did not affect hatching success ( Figure S3), suggesting that the observed effects on stress response were not a consequence of general developmental impairments.
Approximately 34 ± 2.3% of reads in the 25-32nt range (or 9.4 ± 2.3% of total cleaned sRNA reads) aligned to predicted piRNA cluster sequences derived from the piRNA cluster database, although 31.6 ± 1.8% of piRNA-aligning reads with alignment score(s) of>/= −11 had equally valid alignments to multiple clusters (these reads were not counted).
Following alignment of the remaining 18-40 nt reads to the 5' and 3' ends of mature tRNA sequences in the genomic tRNA database, we identified 610 and 142 unique sequences which we considered putative 5' and 3'tsRNAs, respectively. In terms of total read abundance, most 5'tsRNAs were 32-34 nt in length, while almost all 3'tsRNAs were 36 nt (Figure 2c). The distribution of unique sequences followed the same general pattern, albeit sequence lengths were more evenly distributed than for total reads ( Figure 2d). When subjected to clustering by sequence similarity, 214 nonredundant tsRNA sequences (i.e. those which were not completely encapsulated by another sequence) formed 66 clusters. Re-alignment of 18-40nt reads (excluding miRNAs and piRNAs) from individual samples to clustered tsRNA sequences resulted in a 5.2 ± 1.2% alignment rate (or 2.8 ± 0.7% of total cleaned sRNA reads).

F I G U R E 3
Male stress alters small RNA expression in zebrafish spermatozoa. (a) Principal components analysis plots based on multidimensional scaling, showing Euclidean distances between samples according to differences in log2FC of mature miRNAs, hairpinderived sequences (putative isomiRs), piRNA clusters and tsRNA clusters in spermatozoa from paternal control (white points) and stress samples (black points), with experimental batch denoted by point shape. (b) Heatmap derived from hierarchical clustering of normalized expression values (Z-scores of counts per million) of 12 differentially expressed miRNAs, four hairpins exhibiting putative isomiR expression, six piRNA clusters and 12 differentially expressed putative tsRNAs or tsRNA clusters (FDR < 0.1). Putative tsRNAs are named after the mature tRNA sequence as well as the end (5' or 3') from which they are putatively derived. For clusters, the name of a representative tsRNA is shown, followed by the number of nonredundant tsRNAs in the cluster, in brackets. (c) Network plot of differentially expressed spermatozoal miRNAs, their predicted target genes (derived using miRmap) in the egg transcriptome and KEGG pathways enriched with miRNA targets (identified via pathway overrepresentation analysis). The plot is restricted to putative target genes that are included in enriched KEGG pathways (mitophagy and notch signalling) or those that are the predicted targets of at least three miRNAs (socs5b, phyhipla and rnf144aa) and associated miRNAs. Target genes that were not included in any enriched pathways or which were targeted by fewer than three miRNAs are not shown. The two miRNAs in the dre-miR-723 family (dre-miR-723-5p and dre-miR-723-3p) were combined as one for this analysis Together, miRNAs (including possible isomiRs), piRNAs and putative tsRNAs comprised approximately 17.3 ± 2.3% of cleaned sRNA reads.

| Differential expression of miRNAs, piRNAs and tsRNAs in paternal spermatozoa
Principal components analyses based on multidimensional scaling revealed that samples consistently segregated according to paternal treatment over the second principal component axis (Figure 3a), suggesting that male stress explained an appreciable proportion of variation in the populations of all three small RNA classes.
Out of 213 mature miRNAs that met the counts threshold for differential expression analysis, we detected 12 that were differentially expressed in response to stress (Figure 3b, Table 1 Among counts of 18-24nt reads which did not align to known mature miRNA sequences but which aligned to hairpin sequences, 79 such hairpins met the required counts threshold, of which four were found to be downregulated in response to stress (Figure 3b, Table 1). The topmost of these were dre-mir-16a (logFC = −2.3, FDR = 0.005) and dre-let-7b (logFC = −2.84, FDR = 0.035). Visual inspection of reads aligning to these hairpins using integrative genomics viewer (Robinson et al., 2011) indicated that most overlapped with the location of the known mature miRNA sequence, suggestive of isomiRs, as opposed to degraded fragments of hairpin precursors.
Out of 55 putative tsRNAs or tsRNA clusters that met the required counts threshold, 12 were differentially expressed in response to male stress (Figure 3b, Table S7). Most of these were downregulated, including a single-member cluster comprising a nonredundant

| Mitophagy and Notch signalling pathways are enriched among predicted targets of differentially expressed spermatozoal miRNAs
Out of 643 total predicted miRNA targets present in the unfertilized egg transcriptome data set (CPM>/= 5), 301 genes were targeted exclusively by upregulated miRNAs, while 289 were targeted exclusively by downregulated miRNAs. Among genes targeted exclusively by upregulated miRNAs, none were targeted by more than two. Among genes targeted exclusively by downregulated miRNAs, we identified three-socs5b, phyhipla and rnf144aa-which were predicted targets of all three of dre-miR-17a-5p, dre-miR-132-3p and dre-miR-200a (Figure 3c). A pathway overrepresentation analysis of the putative targets of differentially expressed miRNAs revealed two KEGG pathways which were enriched with target genes: mitophagy (dre04137; 8 genes; FDR = 0.08) and notch signalling (dre04330; 11 genes; FDR = 0.08), the latter including notch2 and notch3 ( Figure 3c).

| D ISCUSS I ON
Although the role of paternal effects has been increasingly studied in the context of paternal investment postfertilization, such as brood pouch influences in seahorses and pipefish (Cunha, Berglund, Mendes, & Monteiro, 2018;Otero-Ferrer et al., 2020), and nest guarding in sticklebacks (Bell et al., 2016;Mcghee & Bell, 2014), paternal effects mediated by alterations to the sperm epigenome have hitherto not been well characterized in nonmammalian vertebrates such as fish. Here, we show that in a fish with no parental care, heightened paternal stress before fertilization results in offspring with a reduced tendency to exhibit a behavioural response to a natural danger signal (alarm cue). Paternal chronic stress (PCS) offspring furthermore did not exhibit an endocrine stress response in terms of increased cortisol levels following exposure to either alarm cue TA B L E 1 List of mature miRNAs differentially expressed in zebrafish spermatozoa in response to male stress, as well as hairpin sequences exhibiting additional differential expression arising from alternative miRNA isoforms (putative isomiRs) In terms of experimental framework and the context of paternal stress, the previous studies most comparable to the present study were carried out in rodent models of paternal stress, in which the production of corticosterone (Rodgers et al., 2013) and glucose (Gapp et al., 2014) was ameliorated in offspring from stressed fathers. Not only do our results mirror those of these studies on laboratory rodents, but they also add to a growing body of evidence from a range of organisms, including invertebrates (Crean, Dwyer, & Marshall, 2013;Guillaume et al., 2016;Tauffenberger & Parker, 2014;Zajitschek, Zajitschek, & Manier, 2017), fish (Evans, Lymbery, Wiid, Rahman, & Gasparini, 2017;Rodriguez Barreto et al., 2019;Zajitschek et al., 2014) and wild mammals (Weyrich et al., 2016) that the prefertilization paternal environment can be a crucial determinant of offspring epigenome and/or phenotype.
It is unclear whether the suppressed responses to stressors exhibited by PCS offspring has any adaptive significance, as suppression of behavioural response to a natural danger signal would more immediately appear to be maladaptive. However, in the context of neuroendocrine stress, such a suppression could also be seen as adaptive within higher-stress environments with the propensity to induce chronic stress, as it would ameliorate the risk of excessive glucocorticoid production causing dysregulation to other essential processes that are mediated by these hormones (Maestripieri & Hoffman, 2011;Rao & Androulakis, 2019). Furthermore, activation of neuroendocrine stress responses incurs an energetic cost while simultaneously suppressing foraging behaviour-a combination which can be highly detrimental to fitness (Clinchy et al., 2013). Indeed, wild songbirds have lower reproductive success under heightened predation threat, in part due to starvation of nestlings due to reduced parental foraging (Zanette, White, Allen, & Clinchy, 2011).
Therefore, although weaker responses may be a disadvantage in the context of an isolated predator encounter, regular initiation of a strong response in the face of heightened predation risk may be ultimately unsustainable and detrimental to long-term fitness.
Alternatively, it is possible that other mechanisms which rely less extensively on cortisol are employed to respond to challenges most adaptively (MacDougall-Shackleton et al., 2019). Regardless, intergenerational mechanisms to facilitate "fine-tuning" of the intensity or mode of stress response may be ultimately beneficial for longterm offspring fitness.
The differential expression of spermatozoal small RNAs in spermatozoa in response to stress shows that the molecular composition of spermatozoa is sensitive to environmental disruption and is further concordant with patterns in rodent models, suggesting that spermatozoal small RNAs comprise a form of heritable environmental information in teleosts. Two downregulated mature miRNAs found in the spermatozoa of stressed males (dre-let-7b and dre-let-7d-5p, the former of which also exhibited putative isomiR downregulation) belonged to the let-7 family, which is among the most evolutionarily conserved (Lee, Han, Kwon, & Lee, 2016). Interestingly, downregulation of spermatozoal let-7 miRNAs has been implicated in intergenerational effects of nutritional stress (protein restriction) in a rodent model (Sharma et al., 2016), suggesting that the functions of intergenerational transmission of this miRNA family may be similarly conserved.
In addition to miRNAs, we detected increased relative expression of six piRNA clusters in the spermatozoa of stressed males. This particular class of small RNA has also been shown to play a role in intergenerational effects in other organisms, including maternal inheritance of an acquired trait in Drosophila (Grentzinger et al., 2012).
Altered piRNA expression levels were also detected in the spermatozoa of male mice that had been stressed by maternal separation (Gapp et al., 2014). Altered levels of piRNA cluster expression in spermatozoa may reflect stress-induced changes to the expression of transposable elements (TEs) during spermatogenesis, in which piRNAs are employed extensively in the suppression of "opportunistic" TE expression during chromatin remodelling (Ernst, Odom, & Kutter, 2017 to be important in early tissue specification in Xenopus embryos (Heasman, 2006).
Noncoding RNAs are believed play a role in chromatin remodelling during embryogenesis (Pauli, Rinn, & Schier, 2011), and mounting evidence suggests that piRNAs are capable of mediating DNA methylation (Calcagno et al., 2019). Therefore, paternally transmitted noncoding RNAs may influence the methylome of the early embryo. Indeed, zebrafish may be especially sensitive to paternal effects on methylation given that the methylation profile established in the early embryo recapitulates that of the sperm Potok, Nix, Parnell, & Cairns, 2013) Nevertheless, given the fast spermatogenic cycle in zebrafish (Leal et al., 2009) and assuming frequent spawning opportunities, a male could potentially produce several different spermatozoal epigenotypes throughout the course of a single spawning season, thus transmitting a large degree of cumulative epigenetic variation to the next generation. Although some traits inherited epigenetically may be maladaptive to offspring, the range of possible phenotypic variation generated might be to the overall advantage of the father in increasing the overall likelihood of transmitting his genes to subsequent generations. In other words, epigenetic inheritance may constitute a form of genetically encoded plasticity by which the genome itself can enjoy protection from selective pressures.
To conclude, stress-related behavioural and physiological responses are sensitive to intergenerational perturbation via the paternal germ line in wild-type zebrafish. Our data support the idea that spermatozoal small RNAs are important in this process, as found in rodent models, suggesting that RNA-mediated paternal inheritance is a shared between mammals and fish. Equivalent studies on a broader range of species are needed, however, to determine whether this process is widespread among vertebrates and whether it can facilitate adaptation. Moreover, the mechanistic role of spermatozoal small RNAs requires validation in functional experiments, and the role of other epigenetic factors, particularly DNA methylation and histone modification, should be explored. Finally, the work demonstrates that the zebrafish could serve great utility as a model for the further exploration of paternal epigenetic effects.

E TH I C A L A PPROVA L
All animal work was carried out under a UK Home Office licence for the use of animals in scientific procedures (licence number 40/3704).

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and scripts used for the phenotypic analyses (thigmotaxis and cortisol) as well as the complete workflow for the small RNA