Heat induces multiomic and phenotypic stress propagation in zebrafish embryos

Abstract

Heat alters biology from molecular to ecological levels, but may also have unknown indirect effects. This includes the concept that animals exposed to abiotic stress can induce stress in naive receivers. Here, we provide a comprehensive picture of the molecular signatures of this process, by integrating multiomic and phenotypic data. In individual zebrafish embryos, repeated heat peaks elicited both a molecular response and a burst of accelerated growth followed by a growth slowdown in concert with reduced responses to novel stimuli. Metabolomes of the media of heat treated vs. untreated embryos revealed candidate stress metabolites including sulfur-containing compounds and lipids. These stress metabolites elicited transcriptomic changes in naive receivers related to immune response, extracellular signaling, glycosaminoglycan/keratan sulfate, and lipid metabolism. Consequently, non-heat-exposed receivers (exposed to stress metabolites only) experienced accelerated catch-up growth in concert with reduced swimming performance. The combination of heat and stress metabolites accelerated development the most, mediated by apelin signaling. Our results prove the concept of indirect heat-induced stress propagation toward naive receivers, inducing phenotypes comparable with those resulting from direct heat exposure, but utilizing distinct molecular pathways. Group-exposing a nonlaboratory zebrafish line, we independently confirm that the glycosaminoglycan biosynthesis-related gene chs1 and the mucus glycoprotein gene prg4a, functionally connected to the candidate stress metabolite classes sugars and phosphocholine, are differentially expressed in receivers. This hints at the production of Schreckstoff-like cues in receivers, leading to further stress propagation within groups, which may have ecological and animal welfare implications for aquatic populations in a changing climate.

Keywords: multiomics; stress cues; stress propagation; stress response; thermal stress.

Fig. 1.

Scheme of experimental design. A) Hypotheses: abiotic stress causes (1) stress responses in aquatic animals, which (2) leads to the release of stress metabolites that (3) propagate stress responses to naive receivers. B) Zebrafish embryos (D. rerio) < 3 hpf were incubated according to a two-way factorial design represented by two predictors: thermal stress (0, control temperature of 27°C; 1, repeated thermal stress at a sublethal temperature of 32°C as shown in left graph) and SM (0, fresh medium free of metabolites; 1, SM released by embryos exposed to thermal stress). Treatments were CM, control metabolites at 27°C; C, control in fresh medium at 27°C; SM, stress metabolites at 27°C; TS, fresh medium in thermal stress; and TS + SM, stress metabolites in thermal stress. Arrows indicate medium transfer from metabolite donor (asterisk symbols) to metabolite receiver (plain black tube caps). C) Endpoints: individually exposed embryos were sampled for both molecular (pooled samples for RNA-seq and metabolomics of medium at 1 dpf and cortisol and HSP70 at 4 dpf) and phenotypic (1 dpf and 4 dpf) endpoints (circled numbers 1–5). D) Confirmatory experiment with LAMP data of candidate genes found in RNA-seq in group-exposed (20 embryos per petri dish) embryos incubated in stress metabolites (SM) or fresh medium (C) until 1 dpf. Endpoints 1/2/5 utilized laboratory inbred AB and endpoints 3/6 wilder PET embryos.

Fig. 2.

TS and SM alter the transcriptome and its functions in 1-dpf zebrafish embryos. Top row: volcano plots showing the DEGs in response to A) TS, B) SM, and C) their combination TS + SM compared with control C. Genes of interest are shown in red when they have significant raw P values (above horizontal line) and an absolute FC (representing the effect size) > 2 (|log 2 FC| > 1, vertical lines). DEGs left to the left vertical line and right to the right vertical line are respectively significantly underexpressed and overexpressed compared with the control C. Middle row: gene functional categories from KEGG, Reactome, and GO BP analysis for D) TS, E) SM, and F) those uniquely present in the combined treatment TS + SM. Top and middle rows show transcriptomic data from individually raised AB strain embryos. Bottom row shows LAMP data of candidate genes associated with stress metabolites showing the normalized time of detection (−ΔΔC_T) in C and SM in more realistic environmental conditions (genetically diversified outcrossing PET strain raised in groups). G) chs1, chitin synthase 1; H) ldha, lactate dehydrogenase A4; I) ora3, olfactory receptor class A related 3. ¹ND, mRNA amplification nondetected suggesting a depletion of ora3 in SM; J) otofa, otoferlin a; K) prg4a, proteoglycan 4a; L) tlr18, toll-like receptor 18. Student's t tests compared SM with C with significant comparisons shown by horizontal bars with *P ≤ 0.05 and **P ≤ 0.01. Treatments were SM, stress metabolites at 27°C; TS, fresh medium in thermal stress; and TS + SM, stress metabolites in thermal stress, compared with C, control in fresh medium at 27°C.

Fig. 3.

Heat-induced DEGs are functionally connected to stress metabolites, which interact with the transcriptome of receivers. Hypothetical CPI network analysis of significant DEGs (circles) of TS (top row, A) or genes of SM (bottom row, B) and stress metabolite compounds (in red without circle). Line width represents the increasing confidence score with a minimum threshold of 0.7 (TS) or 0.4 (SM). Computed and drawn independently for A) and B) in STITCH using Cytoscape and merged in Inkscape. A) Compound-protein interaction (CPI) enrichment P < 0.0001, 64 edges observed vs. 5 expected edges. B) CPI enrichment P < 0.0001, 14 CPI edges observed vs. 3 expected. GPC, glycerophosphocholine. Treatments were SM, stress metabolites at 27°C, and TS, fresh medium in thermal stress compared with C, control in fresh medium at 27°C. Genes with red bold circles were used for LAMP measurements. Genes with red bold circles were selected for LAMP measurements based on evidence from multiomic data of their possible involvement in stress propagation and functional links with stress metabolites.

Fig. 4.

Multiomic data evidence that stress metabolites (SM) differ in classes and functions from control metabolites (CM). A) Correlation plot showing filtered (n = 89) masstags that are possible biomarkers of CM (blue circles, bottom right) and SM (green squares, top left) groups. Data are presented as blank-corrected total percent intensity (BC%), which represents the relative concentrations of the compounds in the medium sample. Masstag inside the black dotted ellipse are unlikely to be biomarkers of SM or CM since they are present in both conditions. Masstags inside the dot-dashed blue and dashed green ellipses represent CM- and SM-specific biomarkers, respectively. FC, fold-change between SM and CM where the numerator is the medium with the highest concentration for each compound. Labels show the masstags with possible hits among the top 23 compounds of SM, top 14 compounds of CM based on fold-change (see Table S3 for associated chemicals), and top 10 compounds based on difference between SM and CM (delta Δ identifiers). Superclasses of representative hits for the biomarkers of B) CM and C) SM. Bubble sizes are proportional to counts of superclasses per medium type. Candidate cause-and-effect pathways derived from multiomic analysis of KEGG pathways integrating the representative hits of the SM with the DEGs of D) TS or E) SM. The top five most evident pathways are given within plots, respectively. Significant terms are shown by **** (P-adj ≤ 0.0001). PPP, pentose phosphate pathway. See Dataset S2 for class-sorted compounds.

Fig. 5.

TS and SM alter the development and behavior of zebrafish embryos. A) Principal component analysis (PCA) of 1-dpf morphology. HPF, final stage in hours post fertilization; Stage, % in segmentation or pharyngula; EA, eye area; EL, eye length; DVL, dorsal–ventral length; HTA, head–trunk angle; LEL, longest embryo length; OVL, otic vesicle length; SEL, shortest embryo length; WBA, whole-body area; YBA/L, yolk ball area/length; YEA/L, yolk extension area/length; YE/YB, yolk extension/yolk ball ratio. Individual variables significantly altered by TS are shown by black arrows. B) Burst count per minute in active 1-dpf embryos. C) SEL (mm) increment from 1 to 4 dpf. D) Mean acceleration (m/s²) of the first complete burst following the first and third touch stimuli. E) Density distribution of total burst count and total distance in cm F) after three touch stimuli. Effects of TS and SM across the two-way factorial design (C, SM, TS, and TS + SM) were computed using PERMANOVA A), ANOVA B–D, F), or negative binomial generalized linear model (E, confirmed by a density test) with significant predictors and covariates (“batch”) reported on top-right corners. Pairwise tests with FDR P value correction compared CM to C and SM, or C to SM, TS, and TS + SM, with significant comparisons shown by horizontal bars. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, and ****P ≤ 0.0001. Treatments were CM, control metabolites at 27°C; C, control in fresh medium at 27°C; SM, stress metabolites at 27°C; TS, fresh medium in thermal stress; and TS + SM, stress metabolites in thermal stress.

Fig. 6.

Conceptual summary of stress-induced communication mediated through SM in zebrafish embryos.