High mirror-induced aggression phenotype is heritable over generations

To establish high mirror aggression zebrafish (HAZ) and low mirror aggression zebrafish (LAZ) lines, we randomly selected AB wild-type zebrafish from our local breeding colony and characterized their response to their own mirror image (the F0 generation). These zebrafish exhibited a large variation in aggression ranging from 0 to 178.06 s during a standard 5 min observation period with no differences between males and females (Fig. 1a, b). However, the aggression level of the 10 most aggressive and the 10 least aggressive F0 zebrafish was clearly separable (Fig. 1c). To investigate whether low and/or high mirror-induced aggression is heritable, we measured this behaviour in the offspring of the 3 most aggressive HAZ males and females and the 3 least aggressive LAZ males and females (F1 generation) respectively. Analysis revealed that F1 HAZ are on average more aggressive towards their mirror image than F1 LAZ at 1 month of age (i.e. juvenile zebrafish, Fig. 1d). Although F2 juvenile HAZ and LAZ did not show differences in mirror aggression (Fig. 1e), F3 and F4 juvenile fish did, again with HAZ being on average more aggressive than LAZ (Fig. 1f, g). This suggests that genetic factors are contributing to this behaviour. We also evaluated mirror aggression in young adult fish (3 months of age) across generations. For this, we selected the 20 most aggressive juvenile HAZ in each generation, raised them to adulthood and compared them to the 20 least aggressive juvenile LAZ of the corresponding generation, in order to better characterize the most and least aggressive individuals in each generation and maximize selection for and against mirror aggression over generations. As can be seen in Fig. 1h–k, these preselected adult HAZ were more aggressive than the preselected LAZ in all generations investigated suggesting that the propensity to be aggressive in the mirror test remains stable from adolescence to adulthood. Interestingly, when comparing data across generations, it appears that HAZ (especially preselected adult fish) became more aggressive over generations as suggested by increasing average times spent in an aggressive display. In contrast, LAZ remained stable (juvenile fish) or tended to increase in average aggression levels (adult fish) suggesting that high mirror-induced aggression phenotype, but not low mirror-induced aggression phenotype is heritable over generations. On the other hand, it is also clear that there is a considerable overlap in time spent in mirror aggression between LAZ and HAZ after prolonged selective breeding (F4 generation) with high and low mirror aggression individuals found in both groups (Fig. 1g, k). However, the fraction of non-aggressive or low aggressive individuals is considerably higher in LAZ and the fraction of highly aggressive individuals is markedly greater in HAZ. To analyse how long-term mirror exposure affects aggressive behaviour of F4 HAZ and LAZ, we next recorded 30 mixed sex and size-matched adult HAZ (13 males, 17 females, size: 2.70 +/− 0.05 cm) and LAZ (14 males, 16 females, size: 2.58 +/− 0.03 cm) during a prolonged MIA session (1h). Similar to our findings in the 5-min mirror-induced aggression (MIA) assay (Fig. 1g, k), HAZ were significantly more aggressive than LAZ without changes to locomotor activity (Fig. 1l–n).

Fig. 1
figure 1

Transgenerational effects of selective breeding on mirror fighting. a Representative heatmap images of an aggressive (top part) and non-aggressive (bottom part) zebrafish individual in the mirror-induced aggression (MIA) setup. b Variation in time spent interacting with the mirror of arbitrarily-selected male and female zebrafish from the local breeding colony (F0 fish) during a 5 min MIA assay. c The 10 most aggressive F0 zebrafish spent significantly more time in aggressive display than the 10 least aggressive F0 zebrafish. n = 10/group (dg). Increased mirror aggression levels of d 1-month-old (n = 85–96) F1 HAZ (high aggression zebrafish) derived from F0 HAZ, e 1-month-old (n = 100/group) F2 HAZ derived from F1 HAZ, f 1-month-old (n = 38–98) F3 HAZ derived from F2 HAZ and g 1-month-old (n = 34–90) F4 HAZ derived from F3 HAZ compared to respective LAZ (low aggression zebrafish). hk Mirror aggression levels of the most aggressive juvenile HAZ at adulthood (3 months of age) compared to the mirror aggression levels of the least aggressive juvenile LAZ at adulthood. Analysis revealed higher mirror aggression levels of h the most aggressive F1 HAZ at 3 months of age (n = 19–20), i the most aggressive F2 HAZ at 3 months of age (n = 14–20), j the most aggressive F3 HAZ at 3 months of age (n = 17–19) and k the most aggressive F4 HAZ at 3 months of age (n = 15–19) compared to the least aggressive LAZ of the respective generation at 3 months of age. ln Behaviour of adult F4 high aggression zebrafish (HAZ) and low aggression zebrafish (LAZ) during prolonged mirror exposure. l Time spent interacting with the mirror, m distance travelled and n time spent immobile of HAZ and LAZ exposed to a 1h mirror-induced aggression assay (n = 30/group). Mann-Whitney U test. ***, P < 0.001; **, P < 0.01 and *, P < 0.05 vs. respective LAZ. Data are presented as mean ± SEM. Source data and individual data values are available in Additional file 2

Selective breeding for mirror aggression co-selects for anxiety

To evaluate whether selective breeding for aggression also affects anxiety, the preselected adult HAZ and LAZ that had been tested in the 5-min MIA assay were also tested in the novel tank diving test (NTD [25];) in each generation. Although these fish did not show differences in anxiety-like behaviour in the F1 and F2 generations (Fig. 2a, b), further selective breeding led to changes in anxiety-like behaviour between the groups. Specifically, in the F3 and F4 generations HAZ spent on average more time in the top compartment of the NTD than LAZ indicating reduced anxiety (Fig. 2c, d). To evaluate whether higher aggression levels and lower anxiety levels are correlated, we used Pearson correlation analysis. As expected, there was no correlation between aggression and anxiety in the F1 and F2 generations given that HAZ did not differ from LAZ in anxiety phenotype in these generations (Fig. 2e–j). However, analysis of the F3 and F4 generation revealed a significant negative correlation between aggression duration and time in the top zone for HAZ (Fig. 2k, l). This suggests that aggressive HAZ are also less anxious.

Fig. 2
figure 2

Transgenerational effects of selective breeding on anxiety-like behaviour. ad Time spent in the top zone of the novel tank diving test across generations (F1–F4) of adult high aggression zebrafish (HAZ) and adult low aggression zebrafish (LAZ), which had been previously tested in the mirror-induced aggression assay. Mann-Whitney U test. n = 13–21. ***, P < 0.001 and **, P < 0.01 vs. LAZ. Data are presented as mean ± SEM. eh Scatter plots visualizing Pearson correlation analysis between time spent attacking the mirror and time spent in the top zone of the novel tank diving test of LAZ across generations (F1–F4). n = 13–21. il Scatter plots visualizing Pearson correlation analysis between time spent attacking the mirror and time spent in the top zone of the novel tank diving test of HAZ across generations (F1–F4). n = 18–20. Significant correlations are shown by displaying the correlation coefficient (rp) and p value in the respective Figure panels. Source data and individual data values are available in Additional file 2

To evaluate whether F4 HAZ and LAZ differ in other behavioural domains as well, these fish were tested for locomotion, social behaviour and boldness. We did not observe any other differences in behaviour including locomotor activity in the open field test (distance moved, velocity, time spent immobile, and angular velocity), social behaviour (shoaling and social interaction) and novel object boldness (Additional file 1: Fig. S1).

Selective breeding and mirror exposure strongly alter the brain transcriptome

Next, we investigated how selective breeding and mirror fighting affect the brain transcriptome. We collected the brains of HAZ and LAZ founder fish used in this study (F0 generation), brains from F4 HAZ and LAZ individuals not exposed to the mirror and brains from F4 HAZ and LAZ exposed to the mirror for 1h. Analysis of the F0 fish brains collected at baseline revealed only minor transcriptomic differences between HAZ and LAZ founders. Although the two groups were very different during the mirror assay (Fig. 1c), differential expression analysis showed that only 6 protein-coding genes (prostaglandin reductase 1 (ptgr1), si:ch211-213a13.2, secretogranin V (scg5), si:dkey-238k10.1, LON peptidase N-terminal domain and ring finger 1 (lonrf1) and si:ch211-76m11.3) were differentially expressed (padj < 0.05 and LFC > |2|). In addition, no clustering of the HAZ and LAZ samples was evident by principal component analysis (PCA; Additional file 1: Fig. S2). This suggests that the F0 fish used to generate high aggressive and low aggressive lines were similar in their baseline genetic makeup and that their different behaviour during the MIA assay was most likely the result of non-genetic factors (e.g. tank hierarchy) in the founders.

When analysing the brain transcriptome of F4 HAZ and LAZ fish, the situation was different. Already at baseline (F4 fish not exposed to a mirror), there were considerable gene expression differences between the groups. PCA revealed a clear separation of HAZ and LAZ along PC1 of the plot explaining 46% of the total variance (Fig. 3a) and differential expression analysis showed 85 differentially expressed genes (DEGs) between the two groups. Twenty-two of these DEGs were upregulated in HAZ compared to LAZ, while 63 DEGs were downregulated (Fig. 3b). Hierarchical clustering of samples also revealed a clear separation between HAZ and LAZ indicating similar transcriptional patterns within the 2 lines (Fig. 3c).

Fig. 3
figure 3

Selective breeding-induced and mirror-induced neurotranscriptomic differences of male and female high aggression zebrafish (HAZ) and low aggression zebrafish (LAZ). a Principal component analysis plot of the top 200 most variable genes after differential expression analysis between F4 HAZ and LAZ without mirror exposure. n = 3/group. b Volcano plot of differentially expressed genes (DEGs; padj < 0.05 and LFC > |2|) between F4 HAZ and LAZ at baseline. DEGs with the lowest adjusted p values (padj) are highlighted. n = 3/group. c Heat map displaying DEGs between F4 HAZ and LAZ at baseline. Hierarchical clustering of samples and genes reveals large differences between HAZ and LAZ, but similar transcriptional patterns within the two lines. n = 3/group. d Principal component analysis plot of the top 200 most variable genes after differential expression analysis between F4 HAZ and LAZ displaying similar aggression levels after mirror exposure. n = 6/group. e Heatmap of DEGs between female HAZ (HAZf) and female LAZ (LAZf) after mirror exposure. Hierarchical clustering of samples and genes reveals large differences between HAZ and LAZ, but similar transcriptional patterns within aggression subgroups. n = 3-4/group. f Heatmap of DEGs between male HAZ (HAZm) and male LAZ (LAZm) after mirror exposure. Hierarchical clustering of samples and genes reveals similar differences like in corresponding male animals. n=2–3/group. g Volcano plot displaying DEGs between HAZf and LAZf after mirror exposure. DEGs with lowest adjusted p values (padj) are highlighted. n = 3–4/group. h Volcano plot displaying DEGs between HAZm and LAZm after mirror exposure. DEGs with the lowest adjusted p values (padj) are highlighted. n = 2–3/group. Golden dots in Volcano plots indicate genes upregulated in HAZ more than log fold change 2, blue dots represent genes downregulated in HAZ more than LFC − 2 and black dots represent genes not passing these thresholds. Source data and individual data values are available at the ebrains data repository, DOI: 10.25493/VTP5-8J9 and in Additional file 2

To analyse how baseline gene expression differences between HAZ and LAZ relate to mirror aggression behaviour, we also sequenced F4 LAZ and HAZ brain samples after 1h mirror exposure. As can be seen in Fig. 1l, these fish showed a large spectrum of behaviour and there was a partial overlap of time spent interacting with the mirror between the groups. We thus decided to analyse two cohorts of fish from this experiment. We first evaluated transcriptomic differences between HAZ and LAZ that showed similar mirror aggression levels (cohort 1). We then evaluated transcriptomic differences between the behavioural extremes of the two groups, the most aggressive HAZ versus the least aggressive LAZ of the experiment (cohort 2).

Fish from the first cohort did not show any significant difference in behaviour during MIA (Additional file 1: Fig. S3), but did show strong transcriptomic differences. PCA again revealed a clear separation between HAZ and LAZ along PC1 of the plot, but no separate clusters for male and female fish from the respective lines (Fig. 3d). Like in F4 fish not exposed to the mirror, hierarchical clustering revealed a clear separation of HAZ and LAZ samples for both male and female fish (Fig. 3e, f), but differential expression analysis revealed considerably more DEGs. Specifically, female high aggression zebrafish (HAZf) had 243 protein-coding DEGs (110 upregulated and 133 downregulated) compared to female low aggression zebrafish (LAZf) from this cohort (Fig. 3g) and male high aggression zebrafish (HAZm) had 209 protein-coding DEGs (97 upregulated and 112 downregulated) compared to male low aggression zebrafish (LAZm) (Fig. 3h). In contrast, only 4 genes were differentially expressed between HAZf and HAZm and 12 genes between LAZf and LAZm indicating only minor sex-specific effects within the groups (Additional file 1: Fig. S4).

When analysing the behaviour of the most aggressive HAZ compared to the least aggressive LAZ during prolonged MIA (cohort 2), we detected extreme differences in mirror interaction between HAZ and LAZ (main effect: F(1,20) = 110.487; P < 0.001) but no sex differences. HAZ spent on average 1117.95 s in the mirror zone during the testing period, whereas LAZ interacted with the mirror for an average of only 3.93 s with half of the LAZ fish not interacting with the mirror at all (Additional file 1: Fig. S5). These stark differences in behaviour between these HAZ and LAZ were accompanied by major transcriptomic differences in the brains of these fish. PCA revealed a clear separation of HAZ and LAZ along PC1 of the plot explaining 48% of the total variance. In contrast, male and female zebrafish did not form separate clusters (Fig. 4a). In line with this finding, differential expression analysis showed a large number of DEGs between aggression subgroups, but not between males and females. Specifically, female mirror fighters (HAZf) of this cohort had 549 protein-coding DEGs (266 upregulated and 283 downregulated) compared to the least aggressive female LAZ (LAZf), which exhibited no or very little mirror fighting behaviour (Fig. 4b, c). Hierarchical clustering of samples revealed a clear separation between HAZ and LAZ indicating similar transcriptional patterns within aggression subgroups (Fig. 4b). Similar to female animals, male mirror fighters (HAZm) of this cohort had 494 protein-coding DEGs compared to LAZm (251 upregulated and 243 downregulated; Fig. 4d, e) and hierarchical clustering of samples also revealed a clear separation between HAZ and LAZ (Fig. 4d).

Fig. 4
figure 4

Neurotranscriptomic differences of the most aggressive male and female high aggression zebrafish (HAZ) compared to the least aggressive low aggression zebrafish (LAZ). a Principal component analysis plot of the top 200 most variable genes after differential expression analysis. n = 6/group. b Heatmap of differentially expressed genes (DEGs; padj < 0.05 and LFC > |2|) between female HAZ (HAZf) and female LAZ (LAZf) of cohort 2. Hierarchical clustering of samples and genes reveals large differences between HAZ and LAZ, but similar transcriptional patterns within aggression subgroups. n = 6/group. c Volcano plot displaying DEGs between HAZf and LAZf. DEGs with the lowest adjusted p value (padj) are highlighted. n = 6/group. d Heatmap of differentially expressed genes (DEGs; padj < 0.05 and LFC > |2|) between male HAZ (HAZm) and male LAZ (LAZm) of cohort 2. Hierarchical clustering of samples and genes reveals large differences between HAZ and LAZ, but similar transcriptional patterns within aggression subgroups. n = 6/group. e Volcano plot displaying DEGs between HAZm and LAZm. DEGs with the lowest adjusted p value (padj) are highlighted. n = 6/group. f Venn diagram showing the overlap of DEGs between the HAZf vs. LAZf comparison and the HAZm vs. LAZm comparison. g Volcano plot displaying DEGs between HAZf and HAZm. DEGs are highlighted. n = 6/group. h Volcano plot displaying DEGs between LAZf and LAZm. DEGs are highlighted. n = 6/group. i Venn diagram showing the overlap of DEGs between the HAZf vs. HAZm comparison and the LAZf vs. LAZm comparison. Golden dots in Volcano plots indicate genes upregulated in HAZ more than log fold change 2, blue dots represent genes downregulated in HAZ more than LFC -2 and black dots represent genes not passing these thresholds. Source data and individual data values are available at the ebrains data repository, DOI: 10.25493/VTP5-8J9 and in Additional file 2

Next, we analysed which genes were differentially expressed between HAZ and LAZ with and without mirror exposure and the degree of DEG overlap between conditions. When looking at the DEGs of cohort 2 animals that showed the largest expression changes between the groups, we found that the 10 most significant DEGs of female animals include both up- and downregulated genes in the HAZf group compared to LAZf (Fig. 4c, Additional file 1: Fig. S6). These DEGs are conserved in mammals and have been associated with diverse important biological functions such as arsenic metabolism (arsenite methyltransferase; as3mt), cell migration (ADP-ribosylation factor-like 4ab; arl4ab), immune system activity (ptgr1), actin cytoskeletal remodelling (WD repeat domain 1; wdr1), corticogenesis (DiGeorge syndrome critical region gene 2; dgcr2), protein dephosphorylation (ubiquitin-like domain containing CTD phosphatase 1; ublcp1), sialic acid metabolism (ST6 (alpha-N-acetyl-neuraminyl-2,3-beta-galactosyl-1,3)-N-acetylgalactosaminide alpha-2,6-sialyltransferase 3; st6galnac3) and ketone body metabolism (acetoacetyl-CoA synthetase; aacs) [26,27,28,29,30,31,32] as well as with unknown biological functions (slc4a1ap, sccpdha.1). The 10 most significant DEGs for male animals in this cohort are highlighted in Fig. 4e and shown in Additional file 1: Fig. S7. Strikingly, the most significant DEGs between conditions in both males and females is as3mt indicating an important role during mirror fights independent of sex. We also found a large degree of overlap in highly significantly changed DEGs in HAZm vs LAZm compared to HAZf vs LAZf, with slc4a1ap, arl4ab, ptgr1, sccpdha.1, dgcr2, ublcp1, aacs and st6galnac3 appearing in both lists. In fact, most DEGs between HAZf and LAZf (313 genes) of this cohort are also differentially expressed between HAZm and LAZm and vice versa (Fig. 4f). This overlap is also very high when analysing, for example, the 25 most significant (21 differentially expressed in HAZf vs. LAZf and HAZm vs. LAZm) or the 100 most significant DEGs in both comparisons (71 differentially expressed in HAZf vs. LAZf and HAZm vs. LAZm) (Additional file 2). This suggests a strong common neurotranscriptional basis of mirror fighting in males and females in agreement with our behavioural data. In contrast, considerably fewer genes were differentially expressed between male and female zebrafish within the HAZ and LAZ groups of cohort 2. HAZf zebrafish had 7 DEGs (5 upregulated and 2 downregulated) compared to HAZm (Fig. 4g) and LAZf showed 5 DEGs (3 upregulated and 2 downregulated) compared to LAZm (Fig. 4h). Three of these DEGs (coagulation factor XIII, A1 polypeptide a, tandem duplicate 1 (f13a1a.1), iodothyronine deiodinase 2 and phosphate cytidylyltransferase 1B, choline b) were differentially expressed between both HAZf vs HAZm and LAZf vs LAZm (Fig. 4i). Taken together, this analysis suggests that there are neurotranscriptional sex differences, but these are minor compared to the differences detected between HAZ and LAZ and they are not detectable at the behavioural level.

To disentangle which of the DEGs between HAZ and LAZ are related to selective breeding and which DEGs are related to mirror fighting, we analysed the DEG overlap between mirror-exposed individuals (cohorts 1 and 2) and fish not exposed to a mirror. These comparisons revealed an overlap of 32-42 DEGs (Fig. 5a–d) with a core set of 28 genes present in all comparisons (Fig. 5e). When subtracting genes that were already differentially expressed in the F0 generation (ptgr1 and lonrf1), 26 genes remained, reflecting stable gene expression changes induced by selective breeding independent of mirror exposure (Additional file 3: Table 1). However, when considering the large amount of DEGs between cohort 1 HAZ and LAZ fish (males: 209 DEGs, females: 243 DEGs), which spent a similar amount of time interacting with the mirror (LAZ: 608.9 ± 108.2 s vs HAZ: 429.5 ± 91.58 s), it seems likely that selective breeding led to additional gene expression changes unmasked by mirror exposure. Finally, we investigated, which genes are induced by mirror exposure. For this, we subtracted the stably expressed genes induced by selective breeding mentioned above (Additional file 3: Table 1) and also the overlap of DEGs between mirror-exposed individuals (cohort 1 versus cohort 2 fish; Fig. 5f, g), reflecting additional genes induced by selective breeding, from the DEG lists of cohort 2. This results in two lists of genes of 395 and 372 DEGs (Additional file 4: Table 2, Additional file 5: Table 3) for females and males respectively, which are most likely induced by mirror fighting.

Fig. 5
figure 5

Overlap of differentially expressed genes (DEGs) between high aggression zebrafish (HAZ) and low aggression zebrafish (LAZ) at baseline and after mirror exposure. a Venn diagram showing the overlap of DEGs between the HAZ vs. LAZ comparison at baseline and the comparison of female cohort 2 fish. b Venn diagram showing the overlap of DEGs between the HAZ vs. LAZ comparison at baseline and the comparison of male cohort 2 fish. c Venn diagram showing the overlap of DEGs between the HAZ vs. LAZ comparison at baseline and the comparison of female cohort 1 fish. d Venn diagram showing the overlap of DEGs between the HAZ vs. LAZ comparison at baseline and the comparison of male cohort 1 fish. e Venn diagram depicting common DEGs between the comparisons of male and female cohort 1 and cohort 2 fish. f Venn diagram showing the overlap of DEGs between the comparison of female cohort 1 fish and the comparison of female cohort 2 fish. g Venn diagram showing the overlap of DEGs between the comparison of male cohort 1 fish and the comparison of male cohort 2 fish. For all comparisons, cohort 1 denotes HAZ and LAZ exposed to the mirror for 1h and displaying similar aggression levels, whereas cohort 2 denotes the comparison between the most aggressive HAZ during prolonged mirror exposure and the least aggressive LAZ. Source data and individual data values are available in Additional file 2

To validate the RNAseq findings, we performed RT-qPCR analysis of 6 DEGs between HAZm and LAZm across levels of significance, fold change and direction of change using an independent group of animals taken from the prolonged MIA experiment. There was a significant positive correlation between the RNAseq and qPCR expression changes (Pearson correlation coefficient of 0.869, p = 0.02) suggesting that our RNAseq data is reliable (Additional file 1: Fig. S8).

Pathway analysis reveals an enrichment of immune-related pathways

To investigate biological pathways underlying the observed gene expression changes, we used DAVID pathway analysis and weighted gene expression analysis (WGCNA). For this, we focused on cohort 2 of the prolonged MIA experiment, which showed the largest differences in behaviour and also gene expression. We used DAVID pathway analysis [33] to classify DEGs into functional clusters revealing 8 significantly enriched annotation clusters (enrichment score ≥1.3) in HAZf animals compared to LAZf from this cohort (Fig. 6). The most enriched annotation cluster for the HAZf vs LAZf comparison (cluster 1, enrichment score: 4.35) contains the GO term “immune response” and the INTERPRO term “Lymphocyte function associated antigen 3” suggesting differences in immunity-related pathways between the two groups. This notion is supported by the significant enrichment of two other clusters (Fig. 6a) containing the immunoglobulin-related terms “Immunoglobulin-like domain”, “Immunoglobulin-like fold”, “Immunoglobulin subtype”, “Immunoglobulin V-set” and “IG” (cluster 3, enrichment score: 2.56) and the chemotaxis-related terms “Chemokine receptor family”, “chemokine receptor activity”, “chemotaxis” and “Cytokine-cytokine receptor interaction” (cluster 5, enrichment score: 1.89). DEGs with important immunomodulatory functions present in one or more of these enriched clusters include chemokine (C-X-C motif) ligand 12b (cxcl12b), tumour necrosis factor receptor superfamily member 1B (tnfrsf1b) and various chemokine receptors (chemokine (C motif) receptor 1a, duplicate 1 (xcr1a.1), chemokine (C motif) receptor 1b, duplicate 1 (xcr1b.1), chemokine (C-C motif) receptor 11.1(ccr11.1), chemokine (C-C motif) receptor 8.1 (ccr8.1)). (Fig. 6b–g, Additional file 1: Table 4).

Fig. 6
figure 6

Enriched clusters and functional categories between the most aggressive female high aggression zebrafish (HAZf) and the least aggressive female low aggression zebrafish (LAZf). a Functional annotation clustering using DAVID pathway analysis revealed 8 significantly enriched clusters (enrichment score ≥ 1.3). Annotation terms related to each cluster are displayed in pink boxes. bg DEGs with known important immunomodulatory functions present in one or more of these enriched clusters include b chemokine (C-X-C motif) ligand 12b (cxcl12b), c tumour necrosis factor receptor superfamily member 1B (tnfrsf1b), d chemokine (C motif) receptor 1a, duplicate 1 (xcr1a.1), e chemokine (C motif) receptor 1b, duplicate 1 (xcr1b.1), f chemokine (C-C motif) receptor 11.1(ccr11.1) and g chemokine (C-C motif) receptor 8.1 (ccr8.1). n = 6/group. Source data and individual data values are available in Additional file 2

Similarly, we detected an enriched cluster containing the immunoglobulin-related terms “immune response”, “Immunoglobulin-like fold”, “IG”, “Immunoglobulin subtype”, “Immunoglobulin-like domain”, “Immunoglobulin V-set” and “Lymphocyte function associated antigen 3” in male animals when comparing HAZm with LAZm (cluster 4, enrichment score: 1.66, Fig. 7a). In addition, this comparison showed another enriched cluster (cluster 5, enrichment score: 1.48) related to cytokine action containing the terms “TNFR” and “TNFR/NGFR cysteine-rich region”, “response to lipopolysaccharide” and “tumor necrosis factor-activated receptor activity”. DEGs in these enriched clusters with known important immunomodulatory functions include chemokine (C-C motif) ligand 35 duplicate 2 (ccl35.2), interleukin 12B, c (il12bc), leukotriene B4 receptor 2b (ltb4r2b) and chemokine receptors (xcr1b.1, ccr11.1, ccr8.1) (Fig. 7b–g, Additional file 1: Table 5)

Fig. 7
figure 7

Enriched clusters and functional categories between the most aggressive male high aggression zebrafish (HAZm) and the least aggressive male low aggression zebrafish (LAZm). a Functional annotation clustering using DAVID pathway analysis revealed 5 significantly enriched clusters (enrichment score ≥ 1.3). Annotation terms related to each cluster are displayed in blue boxes. bg DEGs with known important immunomodulatory functions present in one or more of these enriched clusters include b chemokine (C-C motif) ligand 35 duplicate 2 (ccl35.2), c interleukin 12B, c (il12bc), d leukotriene B4 receptor 2b (ltb4r2b), e chemokine (C motif) receptor 1b, duplicate 1 (xcr1b.1), f chemokine (C-C motif) receptor 11.1 (ccr11.1) and g chemokine (C-C motif) receptor 8.1 (ccr8.1). n = 6/group. Source data and individual data values are available in Additional file 2

As well as enrichment of immune system-related processes, both comparisons (HAZf vs LAZf and HAZm vs LAZm) revealed an enriched cluster (enrichment score: 2.71 and 2.57, respectively) containing the GO terms “oxidoreductase activity”, “oxidation-reduction process” and the UniProt term “Oxidoreductase” (Figs. 6a and 7a). Most DEGs related to these terms were found when comparing both HAZf vs LAZf and HAZm vs LAZm (Additional file 1: Table 6). DEGs within these clusters encompass enzymes with a number of diverse biological functions such as nitric oxide synthesis (nitric oxide synthase 2a, inducible, nitric oxide synthase 2b, inducible), prostaglandin metabolism (ptgr1, prostaglandin reductase 2, LOC566996, prostaglandin-endoperoxide synthase 2b, 15-hydroxyprostaglandin dehydrogenase) or steroid hormone metabolism (cytochrome P450, family 11, subfamily A, polypeptide 1, cytochrome P450, family 2, subfamily AA, polypeptide 8, cytochrome P450, family 2, subfamily K, polypeptide16, cytochrome P450, family 3, subfamily C, polypeptide 4, hydroxysteroid (17-beta) dehydrogenase 2) (Additional file 1: Table 6).

Finally, we detected the enrichment of a cluster containing the membrane-associated GO terms “integral component of membrane” and “membrane” as well as the UniProt terms “Transmembrane helix”, “Transmembrane” and “Membrane” (enrichment score: 1.99 and 2.13, respectively) in both HAZf vs LAZf and HAZm vs LAZm; Figs. 6a and 7a). More than 100 DEGs are part of this cluster in both HAZf vs LAZf and HAZm vs LAZm, with a large overlap again (Additional file 1: Table 7). Among the annotated DEGs within this cluster, we found a number of interesting genes encoding, for example, receptor proteins (e.g. adrenoceptor alpha 1Bb (adra1bb) and corticotropin releasing hormone receptor 2 (crhr2)), tight junction proteins (e.g. claudin 11b (cldn11b)) and cell adhesion proteins (e.g. CD164 molecule, sialomucin, selectin E (cd164)). Moreover, in line with the findings described above, many of the DEGs from these clusters have reported immunological functions including leukotriene C4 synthase (ltc4s), toll-like receptor 4b, duplicate a (tlr4ba) or tnfrsf1b (Additional file 1: Table 7).

As expected, due to the low number of DEGs, pathway analysis did not show any enriched pathways for HAZf vs HAZm and LAZf vs LAZm.

Weighted gene co-expression network analysis reveals unique gene modules associated with high and low mirror aggression

Weighted gene co-expression network analysis (WGCNA) was used to identify modules of highly correlated genes that are significantly associated with aggression duration or experimental group. Similar to DAVID pathway analysis, we decided to use the differentially expressed genes of cohort 2 for this analysis and were able to identify 18 co-expression gene modules containing between 41 and 5093 genes (Fig. 8a). To evaluate which of these gene modules are significantly associated with experimental groups or mirror interaction, associations between the eigengene values of each module (weighted average module expression profile) and experimental group or mirror interaction were determined yielding 12 gene co-expression modules that were significantly associated with aggression duration and/or the experimental groups (HAZm, HAZf, LAZm and LAZf). Only one of these gene modules (magenta) was significantly associated with aggression duration and all of the experimental groups. This module is enriched in genes related to a number of processes including cell projections, regulation of neurotransmitter levels, oxidative stress, oxidoreductase activity and lipid metabolism (Additional file 6: Table 8) and was positively correlated with aggression duration, HAZf and HAZm, but negatively correlated with LAZf and LAZm (Fig. 8a, b). Two other gene modules showed a similar association pattern to module magenta; that is in an opposite manner between HAZ and LAZ subgroups. Module light yellow (Fig. 8c), which is enriched in genes related to multicellular signalling, transmembrane transporter activity and the MAPK signalling pathway was positively correlated with HAZf, but negatively correlated with LAZm. In contrast, gene module dark red (Fig. 8d; enriched in genes related to cation transmembrane transporter activity, cell differentiation and the hedgehog signalling pathway) was positively correlated with LAZf, but negatively correlated with HAZm (Fig. 8a–d; Additional file 6: Table 8). Most experimental groups were also characterized by a significant association with at least one unique gene module that was not associated with the other experimental groups, suggesting that each experimental condition has a unique neurotranscriptomic state. Specifically, HAZf was characterized by a positive correlation with the violet module (Fig. 8e; enriched in genes related to regulation of gene expression, transcription factor binding and arginine and proline metabolism) and a negative correlation with the dark turquoise module (Fig. 8f; enriched in genes related to RNA processing, mitochondrial function, oxidoreductase activity and protein modification). Moreover, LAZm and LAZf were characterized by a significant association with three unique gene modules each. Module dark slate blue (enriched in genes related to ribonucleoprotein complex biogenesis, calmodulin binding and amino sugar and nucleotide sugar metabolism), module dark green (enriched in genes related to response to external stimulus, photo receptor activity, tryptophan and purine metabolism) and module light cyan1 (enriched in genes related to protein alkylation, oxidoreductase activity and beta-Alanine metabolism) were positively correlated with LAZm (Fig. 8g–i). Module light green (enriched in genes related to peripheral nervous system development, protein folding, lipid transporter activity and histidine metabolism) was negatively correlated with LAZf (Fig. 8j) and the modules dark orange2 (enriched in genes related to vesicle organization, endocytosis and autophagy) and plum2 (enriched in genes related to second-messenger-mediated signalling) were positively correlated with LAZf (Fig. 8k, l; Additional file 7: Table 9).

Fig. 8
figure 8

Weighted Gene Coexpression Network Analysis (WGCNA) of the most aggressive high aggression zebrafish (HAZ) compared to the least aggressive low aggression zebrafish (LAZ). a Correlations between gene coexpression modules identified by WGCNA and experimental group (male and female high aggression zebrafish (HAZm and HAZf), male and female low aggression zebrafish (LAZm and LAZf) as well as aggression duration. The colours of the boxes are scaled with the value of the correlation coefficient ranging from − 1 (green) to 1 (red). The p value of significant correlations and the respective correlation coefficient are shown in the Figure. bl Eigengene values of samples separated by group (HAZf, HAZm, LAZf and LAZm) for gene modules significantly associated to one or more of the experimental groups or aggression duration. n = 6/group. b Gene coexpression module magenta, c gene coexpression module light yellow, d gene coexpression module dark red, e gene coexpression module violet, f gene coexpression module dark turquoise, g gene coexpression module dark slate blue, h gene coexpression module dark green, i gene coexpression module light cyan1, j gene coexpression module light green, k gene coexpression module dark orange2 and l gene coexpression module plum2. Source data and individual data values are available in Additional file 2

Aggression duration was significantly associated with seven gene modules. Gene modules ivory (enriched in genes related to lipid binding, serotonin receptor activity and tyrosine metabolism), magenta and light yellow were positively correlated with mirror interaction, whereas gene modules dark slate blue, dark red, dark green and light cyan1 were negatively correlated with this behaviour (Fig. 8a). The direction of correlations (positive or negative) shows high similarity for HAZf and HAZm, but not with the LAZ groups. This is most obvious for module magenta that shows a very strong positive correlation with aggression duration (rp = 0.93), moderate positive correlations with HAZm (rp = 0.57) and HAZf (rp = 0.57), but also moderate negative correlations with LAZm (rp = − 0.59) and LAZf (rp = − 0.56). Given that HAZ were very aggressive against their mirror image and LAZ were not, this suggests that the genes within this module are induced by mirror aggression or are particularly important to determine the mirror aggression phenotype of zebrafish. A comparison of these genes with the differential expression analysis result reveals that 42.6% of the DEGs (234 genes) between HAZf and LAZf and 49.8% of the DEGs (246 genes) between HAZm and LAZm are part of this module (Additional file 7: Table 9, Additional file 8: Table 10) showing good overlap between the two different analysis methods used (DESeq2 and WGCNA). Strikingly, all of the 10 most significant DEGs from DESeq2 differential expression analysis are part of module magenta. In addition, 85 of the 100 most significant DEGs between HAZm and LAZm and 90 of the 100 most significant DEGs between HAZf and LAZf are part of module magenta as well. Two other modules (dark red and light yellow) were significantly correlated with aggression duration as well as two of the experimental groups. Module dark red, which is negatively correlated with aggression duration (rp = − 0.69) and HAZm (rp = − 0.52), but positively correlated with LAZf (rp = 0.65) might reflect genes repressed by mirror aggression and module light yellow, which positively correlated with aggression duration (rp = 0.57) and HAZf (rp = 0.61), but negatively correlated with LAZm (rp = − 0.45) might reflect another set of genes induced by mirror aggression. However, for these two modules, the overlap between the genes in each module and the DEGs that differ between HAZf and LAZf or HAZm and LAZm is much lower. Seventeen DEGs between HAZm and LAZm were also found in module dark red and 8 DEGs were also found in module light yellow. A similar overlap was found for DEGs between HAZf and LAZf, with 28 DEGs also appearing in module dark red and 15 DEGs also appearing in module light yellow.

HAZ and LAZ display morphometric differences

Some species of fish display dominance or aggression by changing their external appearance or displaying certain body features [34]. Recent evidence suggests a transient darkening of zebrafish during dyadic fights [35], but it is currently unclear if mirror fighting and selection for aggressive traits leads to morphometric changes in zebrafish. Thus, we investigated whether the observed changes in mirror aggression, anxiety-like behaviour and gene expression are accompanied by changes in body appearance. Stereomicroscopic images of the fish used for RNAseq revealed that aggressive HAZ are slightly bigger than non-aggressive LAZ independent of sex as assessed by standard length (main effect: F(1,20) = 9.269; p = 0.006; Fig. 9a) and height at nape (main effect: F(1,20) = 8.371; p = 0.009; Fig. 9b). Body coloration as well as stripe and interstripe width was largely similar for HAZ and LAZ (Fig. 9c; Additional file 1: Fig. S9), but HAZ had darker 2D stripes (main effect: F(1,20) = 20.508; p < 0.001; Fig. 9d), larger X1V interstripes independently of sex (Fig. 9e; main effect: F(1,20) = 10.195; p = 0.005) and female HAZ had smaller X1D interstripes than male HAZ and female LAZ (Fig. 9f; sex X genotype interaction: F(1,20) = 10.674; p = 0.004). In contrast, sex determined the width of several stripes and interstripes with females displaying larger 1V (sex main effect: F(1,20) = 4.956; p 0.038) and 2V (sex main effect: F(1,20) = 28.687; p <0.001) stripes than males independently of mirror aggression phenotype, but smaller X0 (sex main effect: F(1,20) = 43.685; p < 0.001) and X1V (sex main effect: F(1,20) = 10.827; p = 0.004) interstripes (Additional file 1: Fig. S9a–c; Fig. 9e). The size of the tail- and anal fin did not differ between HAZ and LAZ (Additional file 1: Fig. S9d, e), but female zebrafish had smaller anal fins than males independent of mirror interaction (sex main effect: F(1,20) = 15.313; p < 0.001). Whole body coloration and 1D darkness was similar in all groups (Additional file 1: Fig. S9f, g), but we observed sex main effects in the colour of the 1V and 2V stripes (1V: F(1,20) = 5.002; p = 0.037 and 2V: F(1,20) = 6.806; p = 0.017), likely due to the clear differences in body appearance between males and females (Additional file 1: Fig. S9h, i).

Fig. 9
figure 9

Morphological differences between high aggression zebrafish (HAZ) and low aggression zebrafish (LAZ). a Standard length and b height at nape measurements of male and female HAZ and LAZ exposed to the mirror-induced aggression setup for 1h and used for RNAseq (cohort 2). c Representative stripe pattern images of used HAZ and LAZ. d Stripe colouration of 2D stripes as assessed by grayscale measurements of stereomicroscopic images. (ef Width of X1V and X1D interstripes. Two-way ANOVA followed by Tukey post hoc test in the case of a significant interaction term. n = 6/group. ***, P<0.001; **, P<0.01 main effect HAZ vs. LAZ; bb, P<0.01 vs. LAZ male; ccc, P<0.001 vs. HAZ male. Data are presented as mean ± SEM. Source data and individual data values are available in Additional file 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Disclaimer:

This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.biomedcentral.com/)