Integrated cellular response of the zebrafish (Danio rerio) heart to temperature change

Ectothermic fish live in environments where temperatures can vary both daily and seasonally. These variations in temperature are due in part to changes in sun exposure, stream flow and air temperature (Turko et al., 2021; Fullerton et al., 2015). An acute change in physiological temperature has significant consequences for physiological performance. For example, a decrease in temperature leads to increased blood viscosity, decreased biochemical reaction rates, and increased cellular and tissue stiffness (Burton-Opitz, 1906; Graham and Farrell, 1989). These changes impair the capacity of the heart to contract while also requiring greater force generation to compensate for an increase in blood viscosity and capillary stiffness (Gillis et al., 2000; Graham and Farrell, 1989; Shiels et al., 2003). Such changes can lead to the heart becoming unable to effectively support the aerobic requirements of the animal (Keen et al., 2016; Gillis and Tibbits, 2002). Despite these challenges, several fish species remain active at environmental temperatures that can seasonally cycle by more than 10°C (Klaiman et al., 2011; Spence et al., 2006). To compensate for the effects of low temperature on the function of the heart, some fish undergo significant cardiac remodelling.

Compensation for the effects of a decrease in temperature on the heart requires integrated regulation across multiple levels of biological organization. For example, cold-acclimated rainbow trout (Oncorhynchus mykiss) increase the activity and Ca²⁺ sensitivity of the cardiac myofilament (Alderman et al., 2012; Klaiman et al., 2014) and the abundance of collagen in the myocardium (Keen et al., 2016, 2018, 2021; Klaiman et al., 2011). These changes are associated with increased thickness of myofibrillar bundles and phosphorylation of regulatory proteins, such as troponin T (TnT) (Keen et al., 2016; Klaiman et al., 2011; Kirkpatrick et al., 2011). It is thought that these changes offset the influence of low temperature on cardiac force generation. Cold acclimation also increases reliance on lipid utilization for energy production, as indicated by an increase in the activity of lipid catabolic enzymes without a change in protein concentration (Patey and Driedzic, 1997). This suggests that post-translational modifications (PTMs) of cellular proteins, such as phosphorylation, are used to modify their function. Together, these findings indicate that compensating for the effects of cold exposure on heart function involves comparatively expensive changes to the heart, such as tissue remodelling, coupled with comparatively lower cost modifications, such as PTMs of existing proteins.

In their natural habitat, zebrafish (Danio rerio), live in waters that can range in temperature from 6°C to 38°C (Spence et al., 2008) with the average seasonal variation being closer to 10°C (Das et al., 2022). There is also significant spatial variation where temperature can be 6°C cooler in the centre of flow versus along the bank of a river (Das et al., 2022). Previous work demonstrates that cold acclimation of zebrafish can decrease ventricular collagen content (Johnson et al., 2014; Shaftoe et al., 2023) and can cause changes to the expression of transcripts that encode proteins for collagen regulation (Johnson et al., 2014). This decrease in collagen content in response to cold acclimation is opposite to that which occurs in the trout heart with cold acclimation, mentioned above. It has been suggested that this difference in response is because the blood pressure of zebrafish is less than one-tenth that of trout (Keen et al., 2017). As a result, any increase in tissue stiffness, caused by a decrease in temperature, would impair the ability of the zebrafish heart to fill with blood during diastole thereby necessitating a decrease in tissue stiffness with cold acclimation (Johnson et al., 2014; Keen et al., 2016; Johnston and Gillis, 2022).

Previous work also demonstrates that cold acclimation (from 26–28°C to 18–20°C) causes zebrafish to alter the expression of gene transcripts encoding proteins that regulate the active properties of the myofilament. For example, cold acclimation alters the expression of troponin C (TnC), phospholamban (PLN) and sarcoendoplasmic reticulum Ca²⁺ ATPase (SERCA) (Genge et al., 2013; Little and Seebacher, 2014). Additionally, the use of cardiac ultrasound demonstrates that cold acclimation alters diastolic function in zebrafish (Lee et al., 2016) and that cardiac output is maintained at low temperature (Shaftoe et al., 2023). It is not known whether these functional changes are supported at the cellular level. Other work using cultured fibroblasts and larval zebrafish suggests that cold acclimation activates cytoprotection and metabolic reorganization (Long et al., 2013; Yan et al., 2020). However, to our knowledge, there has not been a comprehensive study examining the integrated cellular changes within the heart of zebrafish during either cold exposure or cold acclimation.

In the current study, zebrafish were exposed to either an acute decrease in temperature or a 6-week acclimation protocol and hearts were sampled at multiple time points in each experiment. Mass spectrometry-based quantitative proteomics followed by gene ontology (GO) analysis was then used to examine changes to the components of functional and structural pathways. We hypothesized that the initial response to a 7°C decrease in temperature would be via phosphorylation of critical regulatory proteins and that 6 weeks of cold acclimation would induce more permanent changes to cellular pathways involved in Ca²⁺-sensitive processes, collagen regulation, energy metabolism and cytoprotection.

Zebrafish [Danio rerio (Hamilton 1822); n=10] of both sexes were acquired from a local supplier (Big Al's Aquarium Services, Kitchener, Ontario, Canada) and raised to at least 6 months of age in a 39.8 l flow-through recirculating tank maintained at 27±1°C using submersible aquarium heaters (EHEIM GmbH and Co. KG., Deizisau, Germany). Water parameters were maintained throughout the experiment, within recommended limits for zebrafish, ammonia <0.1 ppm, nitrite <0.02 ppm and nitrate <10 ppm (Aleström et al., 2020). Zebrafish were fed Gemma 300 (Skretting, ME, USA) twice daily to satiation. At the beginning of the experiment, zebrafish were randomly assigned to 1 of 3 groups: controls (n=4); 24 h cold-exposed (n=3) or 1 week cold-exposed (n=3). All fish were transferred to 2 l flow-through aquaria maintained at the same holding conditions, with the exception of temperature.

Control fish were transferred to a new tank and held at 27±1°C while cold groups were transferred directly to 20±1°C. All fish remained active during and after the transfer. It has been shown that the whole-body cortisol of zebrafish recovers from net stress within an hour of handling (Ramsay et al., 2009). Even so, control fish were treated the same as cold-exposed fish other than the change in temperature. Cold-exposed fish were sampled after 24 h and 1 week. All procedures were approved by the University of Guelph Animal Care Committee (AUP: 4891) under the auspices of the Canadian Council for Animal Care.

On the day of sampling, fish were killed by immersion in an ice water bath followed by decapitation. The thoracic cavity was then opened and the ventricle was extracted with forceps, washed in physiological saline (in mmol l⁻¹: 94 NaCl, 24 NaHCO₃, 5 KCl, 1 MgSO₄, 1 Na₂HPO₄, 0.7 CaCl₂; buffered to 7.6 pH with HCl at 20°C) containing proteinase (cOmplete^TM Protease Inhibitor Cocktail, Roche Austria GmbH, Vienna, Austria) and phosphatase (PhosStop, Roche Austria GmbH, Vienna, Austria) inhibitors, then flash frozen on dry ice. Samples were stored at −80°C until processing. The mass and length of each fish without the heart was recorded to calculate Fulton's condition factor (K): (mass/length³)×100.

The methods used for preparing, managing and analysing samples for mass spectrometry are modified from Geddes-McAlister and Gadjeva (2019) and Sukumaran et al. (2022). Frozen, single hearts; n=4 for controls; n=3 for each of 24 h and 1-week cold-exposed treatments) were removed from dry ice to 300 µl ice-cold Tris-HCl (pH 8.5) containing proteinase and phosphatase inhibitors. Sodium dodecyl sulphate (SDS) was added to 2% to assist in cell lysis. These samples were then homogenized with a probe sonicator (Fisherbrand Model 120 Sonic Dismembrator), for 5 cycles of 30 s each on ice. Dithiothreitol was added to a final concentration of 10 mmol l⁻¹ and incubated for 10 min at 95°C on a Thermal Shake Touch Thermoshaker (VWR International, Radnor, PA, USA) at 800 rpm. The samples were then allowed to cool to room temperature before iodoacetamide was added to a final concentration of 55 mmol l⁻¹. These samples were then incubated in the dark for 20 min at room temperature. To precipitate the protein, 100% ice cold acetone was added to 80% and this mixture was incubated overnight at −20°C.

Following overnight incubation, samples were centrifuged at 13,500 rpm, 4°C for 10 min (rotor head diameter=16.8 cm). The resulting pellet was then washed twice with 500 µl of 80% acetone under the same conditions. The final pellet was dried for 15–20 min at 37°C in the Thermoshaker. The pellets were then resuspended in 100 µl of 8 mol l⁻¹ urea/40 mmol l⁻¹ HEPES and subsequently sonicated in a bath sonicator (VWR) on ice for 15 cycles of 30 s on/30 s off. The protein concentration of each sample was determined using bovine serum albumin (BSA) assay according to manufacturer's protocol (Pierce BSA Assay, Thermo Fisher Scientific). In brief, BSA was added to 8 mol l⁻¹ urea then serially diluted in duplicate to construct the standard curve. Sample absorbance was then measured at 280 nm in duplicate and the protein content was quantified from the standard curve. Each sample was then standardized to 400 µl volume by adding 300 µl ammonium bicarbonate. Then 170 µg of protein was removed and trypsin/Lys-C enzyme mix (0.5 µg µl⁻¹) was added (1:50) (Thermo Fisher Scientific). The samples were then incubated overnight at room temperature. One sample from each of the treatment groups did not yield sufficient protein to proceed and were removed from the study.

Stopping solution [20% acetonitrile (CAN)/6% trifluoroacetic acid (TFA)] was added in a 1:10 ratio to arrest digestion. The samples were divided into 20 µg peptide for total proteome analysis and 150 µg aliquot for phosphoproteome analysis. The total proteome samples were immediately desalted using three layers of C18 paper packed into a pipette tip and centrifuged at 1000 g (Sonation STC-V2, STAGE-tipping microcentrifuge). C18 filters were equilibrated by sequentially washing buffers through the filter. Filters were centrifuged with buffer at 1000 g until nearly drained through. The next buffer was added and likewise centrifuged until nearly drained and so on. Therefore, 100 µl 100% ACN was first added to the C18 filter column and centrifuged. Then 50 µl buffer B (2% ACN, 0.1% TFA, 0.5% acetic acid in ddH₂O) was added and centrifuged. Finally, 200 µl buffer A (80% ACN, 0.5% acetate in ddH₂O) was added and centrifuged. This equilibration step standardizes pH throughout the filter. The sample was then added and centrifuged before the filter was washed in 200 µl buffer A as described above. Samples were eluted in 50 µl buffer B from the filters dropwise by affixing the tip to a 10 ml syringe and applying pressure to the plunger. Samples were then dried using a vacuum centrifuge (65°C, Eppendorf model 07-748-15).

To enrich the samples for phosphorylated peptides, the High-Select^TM Fe-NTA Phosphopeptide Enrichment Kit (Thermo Fisher Scientific) was used according to manufacturer's instructions. Briefly, the 150 µg samples were spun dry in the vacuum centrifuge, then resuspended in the provided wash buffer (corrected to pH<3 with formic acid (FA) to prevent removal of phosphate groups). The samples were then added to equilibrated spin columns containing binding resin and centrifuged at 1000 g for 30 s, washed in wash buffer at 1000 g and eluted in elution buffer. Samples were then dried in a vacuum centrifuge, stored overnight, then resuspended in 0.1% FA to be desalted in C18 columns as described above. Both the total proteome and the phosphoproteome samples were stored in lyophilized form until mass spectrometric analysis.

Liquid chromatography-tandem mass spectrometry (LC-MS/MS) was performed at Bioinformatics Solutions Inc. (Waterloo, Ontario, Canada). Samples were resuspended in 12 µl of 0.1% FA. Then, 6 µl of each sample was injected into the Thermo Orbitrap Exploris 240 (Thermo Fisher Scientific) by nanoflow liquid chromatography using an Ultimate 3000 chromatography system (Thermo Fisher Scientific). Liquid chromatography was performed using a constant flow of 0.25 µl min⁻¹ and a 15 cm reversed-phased column with a 75 µm inner diameter filled with Reprosil C18 (PepSEP, Bruker, Germany). Mobile phase A was 0.1% FA and mobile phase B was 99.9% ACN, 0.1% FA. For total protein, a 120 min solvent gradient was used; the proportion of buffer B was increased from 4% linearly to 32% over the first 110 min and then buffer B was increased to 95% over the following 6 s. This concentration of buffer B was maintained for 5 min to clean the column. Then, buffer B was decreased to 4% for the final 5 min.

Phospho-enriched protein separation was carried out over 60 min as follows; the concentration of buffer B was linearly increased from 4% to 35% over 53 min with an increase to 95% over 5 s where it was maintained for 2 min, 54 s to clean the column. The proportion of buffer B was returned to 4% for the final 4 min.

The mass spectrometry data were acquired on Thermo Orbitrap Exploris 240 in data-dependent mode with a cycle time of 3 s. The first scan data were obtained at 60,000 resolution (at 400 m/z) with a mass range of 400–1600 m/z. The automatic gain control (AGC) was set to standard, with an auto maximum ion injection time. The radio frequency (RF) lens was set to 70%. Isolation for the second scans (MS2) was performed in the quadrupole, with an isolation window of 0.7. MS2 scan data were acquired at a resolution of 15,000 m/z, with a standard AGC target and an auto ion injection time. The scan range of MS2 was also set to auto. Higher energy collisional dissociation (HCD; fixed normalized collision energy: 30%) was used for generating MS2 spectra, with the number of microscans set to 1. The dynamic exclusion was set as 8 and 20 s in the 60 and 120 min gradients, respectively.

The raw mass spectra files were analysed in MaxQuant software (2.2.0.0) (Cox and Mann, 2008) using the integrated Andromeda search engine against the reference Danio rerio (zebrafish, Brachydanio rerio, UP000000437) proteome of 46,696 FASTA sequences (Cox et al., 2011) accessed via Uniprot (16 December 2022). We included trypsin/LysC enzyme specificity with a maximum of two missed cleavage sites, a minimum peptide length of seven amino acids, fixed modification (carbamidomethylation of cysteine) and variable modifications (oxidation of methionine, acetylation of N-terminus). A phospho-specific variable modification (phosphorylation of serine, threonine or tyrosine) was included for the phospho-enriched samples. Protein IDs were adjusted against a false discovery rate (FDR) of 1% with a minimum of one peptide being required between two samples for protein identification (Cox and Mann, 2008). Label-free quantification (LFQ) was used and peptides were matched between runs with a minimum ratio count of 2 (Cox and Mann, 2008).

The difference in expression between the cardiac proteome of cold-exposed and cold-acclimated fish and those of the controls were analysed and visualized in Perseus (v.1.6.15.0) (Tyanova et al., 2016). Using principal component analysis (PCA), we identified one replicate of the control group that differed significantly from the others due to technical variability. This sample was, therefore, excluded from subsequent analysis. All samples were filtered to remove potential reverse peptide matches, peptides identified only by one site and by potential contaminants. Phospho-enriched samples were additionally restricted to those hits with 75% localization probability, which uses the possible peptide fragments to determine the amino acid residues likely to be phosphorylated. Samples were then log₂ transformed and filtered for peptides that occurred in at least 50% of each treatment group to ensure high fidelity in our comparisons. A threshold of 50% was used because it enabled us to maintain a consistent threshold for an odd number of samples. Missing values were interpolated from a normal distribution with a width of 0.3, targeting lowly expressed peptides in the lowest 20% of the distribution. We then compared protein abundance between the control group, the 24 h cold-exposed group and the 1 week cold-acclimated group using a one-way ANOVA (abundance versus treatment). We utilized an FDR of 0.05 to correct for false positives and an S₀ of 1 (which weighs the statistical test against the fold change along with the P-value). However, previous work suggests that FDR correction may be too restrictive for high-throughput biological data (Pascovici et al., 2016); therefore, we considered uncorrected P<0.05 as the threshold for statistical significance (Bárez-López et al., 2023). Previous work has also used a threshold of 1.2 log₂-fold change to delimit biological relevance for proteomic datasets using the more conservative iTRAQ method (Alderman et al., 2021). Such a change represents a 20% increase or decrease in protein abundance. Other work suggests that a 2 log₂-fold change is appropriate for large datasets (Johansen et al., 2023). All proteins identified as significantly different exceeded a 1.2 log₂-fold change. Furthermore, most proteins that were significantly different exceeded a 2 log₂-fold change as well.

All statistical tests were conducted on samples from individual fish. Where multiple technical replicates were utilized (e.g. PCR), these were averaged across all runs from the same individual. Therefore, the sample size in each test corresponds to the number of fish tested within each group.

Proteins identified as being significantly different in abundance between control and cold-acclimated fish were further analysed for GO. Gene names were assigned to each protein from the zebrafish proteome (zebrafish, B. rerio, UP000000437) and run in STRING (NCBI, Danio rerio, https://string-db.org/) to identify proteins with significantly enriched interactions associated with similar physiological pathways. All analyses were conducted independently on the proteins that increased and decreased. All enrichment analyses were run considering significance at P<0.05, adjusted by a false discovery rate of 5% (https://string-db.org, August 2023). The strength of enrichment between pairs of proteins was assessed using k-means so that the network could be grouped into physiological clusters. Where applicable, specific gene function was further explored within the Zebrafish Information Network (ZFIN, https://zfin.org/) and/or the Human Gene Database (Gene Cards, https://www.genecards.org/). Because a change in individual protein abundance can have significant implications for cellular function, we considered it important to report the pathways that significantly affected proteins are a part of even in the absence of GO enrichment. Therefore, while we report only the significantly enriched ‘cellular component’, ‘biological process’, ‘molecular function’ and KEGG terms in our tables, we also used the additional databases that STRING incorporates to evaluate the biological pathways associated with the input protein list. In this way, a parent cluster was characterized first from all proteins that either increased or decreased. Then, k-means were used to identify which proteins contributed to which biological pathways.

The GO analysis in STRING was further used to capture an integrative picture of cold-inducible pathways across the first and second experiments. To do this, we collated the list of proteins whose abundance significantly increased or decreased and re-ran the network analysis. The strength of enrichment was again assessed with k-means and the resultant functional clusters were used to infer important cellular pathways during cold exposure and acclimation in zebrafish. The relative enrichment of GO terms was visualized in R (r-project.org) using the GOPlot package (https://CRAN.R-project.org/package=GOplot).

Mass, length and Fulton's condition factor were compared between the control and cold-exposed groups at 24 h and 1 week using a one-way ANOVA (P<0.05).

Animal husbandry was the same as the first experiment except as follows. Male zebrafish (n=119) were acquired from an in-house breeding colony at the Hagen Aqualab (University of Guelph, Ontario) and were raised to 6 months of age in recirculating tanks (20 l) maintained at 27±1°C. The water parameters were maintained as in the first experiment. At the beginning of the experiment, zebrafish were randomly assigned to one of two groups: controls (n=60) or cold-acclimated (n=59). All fish were transferred to 40 l flow-through aquaria maintained at the same holding conditions, with the exception of temperature. Initially, the tanks of the control and experimental fish were held at 27±1°C. The temperature of the experimental tank was then decreased at a rate of 1°C per day to a final temperature of 20±1°C (Johnson et al., 2014; Klaiman et al., 2011). Both control fish (n=12) and cold-acclimated fish (n=12) were then sampled at 3 days, 1 week, 2 weeks, 4 weeks and 6 weeks. All fish were fed to satiation twice daily. The amount of food provided was weighed before and after feeding to estimate feed intake for each treatment. The food eaten was then standardized to the mass of the fish in the tank (calculated from the average mass of fish sampled at each time point).

During the 6 week experiment, ventilation rates of the zebrafish in both the control and experimental tanks were observed on 2 days mid-week, each week. Fish were observed in the main holding tank in the evening, at least 3 h after their last feeding. This time was chosen as outside the normal feeding schedule, thus reducing stimulation of the fish. Previous work shows that zebrafish fed once have an increase in oxygen consumption for at least 10 h after feeding, temporarily plateauing between 5 and 8 h (Ferreira et al., 2019). Our fish were fed twice daily so they likely were constitutively in a post-prandial state; however, our measurements were ∼8 h after their first feeding, which corresponds with maximum oxygen consumption after feeding (Ferreira et al., 2019). Thus, while we could not account for an increase in oxygen consumption due to specific dynamic action, taking our ventilation measurements at 8 h post first feed and 3 h post latest feed ensures that all fish were in similar (peak) resting physiological state (Fowler et al., 2019). The observer sat stationary for at least 20 min and then zebrafish (n=5–6 per observation) were randomly identified and watched for 5 s. The number of opercular movements in 5 s was recorded and calculated as opercular beats per minute.

Tissue sampling and calculation of body morphometrics as well as protein isolation, quantification and identification were conducted as in the first experiment. In the second experiment, hearts were sampled for RT-qPCR to quantify gene expression of myofilament proteins at 1 and 4 weeks (n=5); as well as for mass spectrometry analysis (n=3–6). Hearts were stored at −80°C until used.

Transcript abundance of serca2a, serca2b, pln1, pln2b, tni2, ctnt and ctnc were measured using RT-qPCR as previously described (Johnston et al., 2013; Johnston and Gillis, 2017). In brief, hearts were solubilized in 500 µl TRIzol on ice using a probe sonicator. Chloroform was then added in 1:2 ratio and RNA was separated using a centrifuge at 4°C. Glycogen was added to the aqueous phase and RNA was precipitated overnight in isopropanol in a 1:2 ratio. The precipitate was pelleted, washed with ethanol, dried and resuspended in 20 µl RNase-free water. A Nano Drop 2000 (Thermo Fisher Scientific) was used to determine the purity and concentration of RNA from each sample. When samples were not pure (A260/280<1.7, A260/230<1.8), samples were reprecipitated using 3 volumes of 95% ethanol and 1:10 volume sodium acetate overnight (Molecular Protocols). cDNA was generated from purified RNA using the High-Capacity cDNA Synthesis Kit (Life Technologies), following standard procedures, except that the samples were incubated after DNase I treatment for 30 min at 37°C. Transcript abundance of the target and the reference genes (ef1a and bac1) were measured in duplicate with RT-qPCR (Bio-Rad CFX2800) using gene-specific primers from the literature or designed from gene sequences (Alderman et al., 2017; Li et al., 2016; Little and Seebacher, 2014). Quantification cycle of each reaction was determined using the default settings and was used to calculate transcript abundance from standard curves generated from serially diluted cDNA. Average gene expression for each individual fish was standardized to the arithmetic mean expression of ef1a and bac1 for that individual (Alderman et al., 2017). All standardized values were then normalized to the mean abundance of the control group at the corresponding time point (Alderman et al., 2017) to evaluate fold change differences at that time. A list of the primers used can be found in Table S1.

Body morphometric and performance data (i.e. ventilation rate, mass, length, Fulton's condition factor, food intake) were evaluated over time using a two-way ANOVA with time as the first factor and treatment as the second. Mass spectrometry data were processed in Perseus the same in the second as in the first experiment. Relative protein abundance and phosphorylation levels in the cold-acclimated fish were compared with the controls using a one-way ANOVA at each time point (2, 4 and 6 weeks) using uncorrected P<0.05 as signifying differential abundance (see experiment 1). Differences in transcript abundance between cold-acclimated fish and controls were compared at 1 and 4 weeks using a two-tailed Student's t-test for each gene. Differences in relative intensity of protein phosphorylation and abundance were compared between cold-acclimated and control groups using a two-tailed t-test (P<0.05). A t-test was used to preserve residuals for the estimation of error owing to low sample size. All tests were set to a significance level of 95%.

There was no difference in the mass or length (P>0.05) of the 24 h cold-exposed, 1 week cold-exposed or control zebrafish in the first experiment (Table 1).

There were 10,776 phosphosites on 3636 proteins characterized in the cardiac proteome of the zebrafish sampled after 24 h. From these we identified 6681 high-quality phosphosites at greater than 75% localization probability. We then refined our search to 1108 phosphosites which were observed in >50% of samples in each of the cold-exposed and control groups (Fig. 1A). At 24 h, ANOVA identified 22 phosphosites on 19 proteins that were significantly more phosphorylated in the cold-exposed group than in the control group (P<0.05; Fig. 1A). We then evaluated the corresponding proteins for enriched GO terms using the zebrafish database in STRING to infer potential pathways affected by cold exposure. The parent cluster matched 18 of the 19 proteins to RNA transport (Table 2). Using the ‘kmeans’ function in STRING to characterize subclusters, two smaller clusters were identified, containing two proteins each. The first cluster was associated with tropomyosin assemblage, the second cluster was the spliceosome and mRNA processing with one significantly enriched GO term.

At 1 week, the 10,776 phosphosites on 3636 proteins were refined to 1174 phosphosites on >50% samples (Fig. 1B). We identified three sites on three proteins at 1 week with significantly increased levels of phosphorylation. These three proteins did not cluster together with enriched terms but the pathways they are involved in are listed in Table 2.

We identified 4960 proteins in the whole proteome of the first experiment. At 24 h, 2052 proteins were found in >50% of samples in each of the cold-exposed and control groups (Fig. 1C). ANOVA revealed 28 differentially abundant proteins (P<0.05); 26 of these differed 1.2 log₂-fold change between the cold-exposed and control group (Fig. 1C). There were eight proteins that were significantly increased in abundance. STRING analysis grouped seven of these proteins into Wnt signalling and protein synthesis pathways (MAPK, RHOAE, TPP2) but no GO terms were enriched (Table 3). There were 21 proteins with significantly decreased abundance at 24 h. Sixteen of these proteins exceeded a 1.2 log₂-fold change, which STRING analysis grouped into apoptotic pathways. Analysis by k-means identified a further two subclusters with two proteins each. The first was associated with but was not enriched for alternative splicing, the second was associated with cell signalling (Table 3).

One week after rapid cooling, there were 2045 proteins identified in the total proteome analysis. Twenty-four of these were differentially abundant (P<0.05); 23 were increased in abundance and also exceeded a 1.2 log₂-fold change (Fig. 1D). STRING analysis identified 23 of the 24 proteins as involved in the synthesis of polyunsaturated fatty acids and lipid metabolism with 6 GO terms enriched. Analysis with k-means identified one subcluster with two proteins enriched for acyl-CoA-dependent β-oxidation. The one protein that decreased in abundance was characterized as an intermediary for cell signal transduction (Table 3).

For full names of the top three proteins at each time point, please refer to Table S2. The three proteins with the greatest decrease in phosphorylated peptides at 24 h were DESMA, TnT3a and KNG1 (Fig. 2A). These are involved with focal adhesion, myofilament Ca²⁺ sensitivity and regulation of cell morphology. The three proteins that increased the most were LDB3A, ATP5MC3A and RHOAE (Fig. 2B), associated with cytoskeletal organization, cell growth and proliferation, respectively. The three proteins that most decreased in abundance were RAB1A, BDH2 and CES2 (Fig. 2B); these proteins are involved with the disruption of mitochondrial function, cytoprotection from reactive oxygen species (ROS) and cell survival/proliferation, respectively.

The proteins with the greatest increase in phosphorylated peptides at 1 week were VTG1, TnT3b and BAIAP2A (Fig. 2C); these proteins represent lipid metabolism and stress response, changes to myofilament Ca²⁺ sensitivity and modulation of cell proliferation. The proteins whose abundance most increased were GFRA3, MuRF-1 and RHAG (Fig. 2D), which are involved in collagen deposition, stretch-induced remodelling and pH regulation in the blood, respectively. The only protein significantly decreased in abundance (YWHABL) is involved in a variety of signalling pathways for mitogenesis and cell growth.

The standard length of the cold-acclimated fish was less than that of the control fish at 2 weeks (–4%, P<0.05) and 4 weeks (–6%, P<0.05). Similarly, the average mass of the cold-acclimated fish was less than that of control fish at day three (–19%, P<0.05) and 4 weeks (–7%, P<0.05; Table 2). Body condition did not differ throughout the experiment as mass and length varied proportionately (Table 1). The cold-acclimated group reached satiation with less food than the control group (Table 1). However, when the amount of food eaten was standardized to the average mass of the fish in the tank (see Materials and Methods), we found that there was no difference in the relative amount of food eaten by either group (Table 1). Ventilation rate was 20% lower in the cold-acclimated group held at 20°C (P<0.05) relative to the control group throughout the experiment (Table 1).

At 2 weeks, we detected 2829 phosphosites on 943 proteins. Of these, 1851 phosphosites were found to have localization probability greater than 75% which were further refined to 365 phosphosites detected in at least 50% of samples in both controls and cold-acclimated fish (Fig. 3A). ANOVA detected differential phosphorylation at 22 phosphosites on 18 proteins (P<0.05). Cold acclimation increased phosphorylation at eight of those sites with five of them exceeding a 1.2 log₂-fold change. These eight phosphosites were enriched for cardiac contraction, phosphatidyl lipid metabolism and cell membrane remodelling (Table 4). There were 14 proteins with significantly decreased phosphorylation, with five exceeding a 1.2 log₂-fold change. Thirteen of these proteins grouped into a single cluster that was associated without enrichment with M-band and sarcomere assemblage as well as changes to cardiac chamber size (Table 4).

At 4 weeks, there were 1473 phosphosites on 491 proteins. Of these phosphosites, 1137 were identified at >75% localization probability and 240 were considered high quality with >50% occurrence in both groups (Fig. 3B). Only one protein had significantly altered phosphorylation, titin (TTN). No enrichment was identified in STRING but TTN is involved in monitoring and responding to mechanical stretch within the myofilament (Table 4).

At 6 weeks, we identified 3351 phosphosites on 1117 proteins. Of these, 2325 phosphosites were detected at greater than 75% localization probability. A total of 416 of these phosphosites were considered high quality (>50% occurrence). ANOVA identified 17 phosphopeptides on 15 proteins that were differentially phosphorylated between the control and the cold-acclimated groups (P<0.05; Fig. 3C). Fifteen of these phosphosites were significantly increased in abundance in the cold-acclimated group and seven exceeded 1.2 log₂-fold change. The 15 proteins were not enriched for any terms but two subclusters comprising two proteins each were identified by k-means. The first enriched for interactions with the actin cytoskeleton while the second was enriched for M-band assemblage (Table 4). The two phosphosites that were found to have a decreased level of phosphorylation did not have any significant enrichment for GO terms but both ARVCFB and TTN are associated with stretch-sensing and myofilament remodelling, respectively (Table 4).

We detected 2059 proteins in the total proteome analysis. 1656 were identified in >50% of both cold-acclimated and control groups at 2 weeks. Of these, we detected 103 high confidence proteins that were differentially abundant (P<0.05; Fig. 3D). Cold acclimation significantly increased the abundance of 51 proteins. Twenty-six of the 51 proteins significantly increased by cold acclimation exceeded 1.2-log₂-fold change. The 51 proteins grouped into a cluster enriched for transcription, translation and alternative splicing as well as protein processing pathways (Table 5). Analysis by k-means identified four subclusters. The first, which had 27 proteins, was enriched for ribosomal assemblage and function and protein biosynthesis (Table 5). The second and third contained three proteins and were enriched for membrane synthesis at the endoplasmic reticulum assemblage and cell signalling via nitric oxide (Table 5). The fourth subcluster contained 2 proteins which were enriched for clathrin-mediated vesicle formation.

We identified 52 proteins that were significantly decreased in abundance after cold acclimation for 2 weeks and 10 of these were decreased by greater than 1.2-log₂-fold. Fifty-one of these were characterized into a cluster for various mitochondrial respiratory processes (Table 5) which relate to the three subclusters identified with k-means. The first contained 36 proteins and was enriched for aerobic metabolism, β-oxidation and the tricarboxylic acid cycle (Table 5). The second and third were did not have any enriched terms but each contained two proteins. The second cluster had proteins associated with peroxisomal function. The third contained proteins involved in endocytosis (Table 5).

After 4 weeks of cold acclimation, we identified 1189 proteins in the total proteome analysis with 635 of these being classified as high confidence (>50% occurrence). The ANOVA identified 23 proteins that were differentially regulated at 4 weeks (Fig. 3E). The seven proteins with increased abundance grouped into a cluster that was enriched for ribosomal processes and regulation of translation (Table 5). The 16 proteins that decreased in abundance were enriched for a variety of mitochondrial respiratory processes (Table 5).

After 6 weeks of cold acclimation, there were 2171 proteins identified in the total proteome analysis. Of these, 1530 proteins were classified as high confidence (>50% representation in both control and cold-acclimated groups). There were 114 proteins that were differentially abundant as detected by ANOVA (P<0.05; Fig. 3F). Of these, 66 proteins significantly increased in abundance (P<0.05) with 43 exceeding a 1.2 log₂-fold change. STRING analysis grouped 61 of these proteins into a cluster that was enriched for proteasomal and phagosomal processes (Table 5). Analysis by k-means identified three subclusters containing 42, 2 and 2 proteins, respectively. The first subcluster was enriched for phagosomal and proteasomal pathways. The second and third subclusters were not enriched but contained proteins that are associated with RNA recognition and nitric oxide signalling, respectively (Table 5).

Of the 48 proteins whose abundance decreased in the cold-acclimated group, 16 exceeded a 1.2-log₂-fold change. STRING analysis grouped 47 of the proteins into a cluster that was enriched for mitochondrial respiration and creatine activity (Table 5). Three subclusters were also identified by k-means with 32, 3 and 2 proteins, respectively. The first was enriched for aerobic and anaerobic pathways. The second and third did not have any enriched GO terms but contained proteins that are associated with micromechanical properties and Z-disc assemblage, respectively (Table 5).

At 2 weeks, the three proteins with the greatest increase in phosphorylated peptides were SORBS2A, PCYT1AA and BAG3 (Fig. 4A). These are associated with glucose metabolism, heat shock and cardiac function and membrane homeoviscous adaptation. The three proteins with the greatest decrease in phosphorylated peptides were SORBS1, CAVIN1A and NDRG4 (Fig. 4A), which are associated with cytoskeletal organization, aerobic and anaerobic metabolism, macromolecular transport and cell signalling. In the total proteome analysis, the top three increased proteins with the greatest log₂-fold increase were GSTAL, LANCL1 and CAP2 (Fig. 4B), which are involved in cytoprotection in antioxidant activity and in cytoskeletal interactions with cAMP The three proteins that decreased the most were CB290, PURAB and TTN1 (Fig. 4B), which are associated with myofilament function, transcription initiation and focal adhesion.

At 4 weeks, the lower levels of phosphorylated TTN2 were concomitant with significantly decreased TTN (P<0.05; Fig. 4C). This phosphosite was found in the stretch-sensitive Titin kinase domain (S31134). The three proteins that increased most in abundance were PTGR1, APOC1 L and RPS9 (Fig. 4D), which are components of the ribosome or serve cytoprotective roles within the plasma membrane. The three proteins found to be decreased in abundance were GSNA, VMHC and ECH (Fig. 4D), which are related to aerobic metabolism, fatty-acid metabolism and anaerobic TCA cycle metabolism, respectively.

At 6 weeks, the three proteins with the greatest increase in phosphorylated peptides were H1L1, DBN2 and MYHC9 L (Fig. 4E), which are associated with localizing intracellular compartments, maintaining DNA integrity and regulating myofilament activity. The two proteins with the greatest decrease in phosphorylated peptides were TTN1 and ARVCFB (Fig. 4E), which are related to modulating cardiac cytoarchitecture and cytoskeletal organization. The three proteins whose abundance increased the most were PLS3, LANCL1 and PTMAB (Fig. 4F), which are involved in cell growth, ROS protection and apoptosis. The proteins whose abundance decreased most were GRPEL, SLC25A3 and TTN1 (Fig. 4F), which are related to ROS cytoprotection, phosphate transport and myofilament tension and protein turn-over.

The changes to the cardiac phosphoproteome and proteome in response to cold exposure and cold acclimation were similar. This included responses that have the potential to alter the utilization of fatty acids, carbohydrate metabolism, mitochondrial function, protein turn-over and myofilament function. This similarity led us to combine the data of the two experiments and repeat the analysis. GO analysis of the results identified enrichment for pathways associated with four main processes: (1) decreased abundance of metabolic enzymes and transporters; (2) increased protein processing and RNA recognition at the endoplasmic reticulum; (3) decreased mitochondrial abundance and increased localization of ATP synthase; and (4) increased utilization of lipids (Fig. S1).

Analysis of the expression of serca2a, serca2b, pln1, pln2b, tnnic, mybpc3 and tnnt2a indicated that cold acclimation increased serca2a expression 3.7-fold at 1 week (P<0.05; Fig. S2). There was no effect of cold acclimation on the transcription of the other genes measured at 1 or 4 weeks. We did not detect changes in the abundance of the proteins associated with any of these genes in the total proteome analysis.

The results of this study indicate that acute cold exposure and cold acclimation lead to significant changes in the phosphoproteome and proteoteme of the zebrafish heart and that these changes are time dependent. It is important to acknowledge that these results reflect the status of individual proteins and through GO analysis, that of multiple components of specific cellular pathways. They do not provide direct functional insight into a cellular capacity or pathway. Keeping this in mind, the results will be put into physiological context using data from studies that have examined the functional response of the fish heart to thermal acclimation and of the cellular pathways that regulate heart function in other vertebrates. By identifying the proteins that were affected by treatment and relating them to the cellular pathways that govern cardiac function in vertebrates, we hope that this work will highlight some areas of research that require greater attention.

The initial changes to the phosphoproteome in response to cold exposure have the potential to stabilize the sarcomere, regulate myofilament function and modulate cytoskeletal structure. This includes changes to the phosphorylation of TTN, TnT and tensin (TNS) at 24 h. Subsequent changes to the proteome would have the potential to increase polyunsaturated fatty acids (PUFAs) content of cellular membranes (ACOX, APOEA and HSD17β4/10) and the remodelling of the cardiomyocyte cytoskeleton (via LDB and ADGRG). One caveat to interpreting these results is that the enzymes that we have found to have altered phosphorylation are also operating at a lower temperature in the cold-acclimated fish. Such a change in temperature would likely decrease enzyme activity. It is not possible to state then, that a change in phosphorylation would increase the activity of the enzyme. This would require functional measurements. However, the reduction in the levels of cytochrome c oxidase (COX), ubiquinol cytochrome c reductase (UQCC) and subunits of ATP synthase (ATP) during cold acclimation do suggest a decrease in metabolic capacity; while the increased levels of ribosomal subunits (e.g. RPL10, RBM24) suggest at least the maintenance of protein synthesis. Whole-animal measurements during the cold acclimation study indicate that changes to the physiological systems of the zebrafish with cold acclimation likely cause a reduction in whole animal metabolic requirements and the maintenance of body condition. The relationship between an apparent decline in metabolic capacity and metabolic requirements requires further functional studies.

The results of the current study suggest that the cellular response to the first 24 h of cold exposure involves an increase in the phosphorylation of proteins involved in stretch-sensing and cellular remodelling. The evidence for this is the elevated phosphorylation of TNS and TTN at its kinase domain, the decreased abundance of YWHABL (or 14-3-3) as well as the increased abundance of stretch-sensitive RHOA and MAPK3. YWHABL is a signalling ligand whose human orthologue inhibits MAPK signalling (Cariolato et al., 2011). Thus, the decrease in YWHABL abundance at 24 h would support an increase in MAPK signalling. One other protein that showed a comparatively large increase in phosphorylation at 24 h was kininogen (KNG1), the protein precursor to the signalling protein, bradykinin (BK) (Howl and Payne, 2003). While targeted studies would be required to confirm the role of KNG1 phosphorylation in BK production, an increase in BK would be consistent with an increase in MAPK signalling (Rex et al., 2022). An increase in BK would also be consistent with a decrease in angiotensin converting enzyme (Xiao et al., 2003), which we saw approaching significance at 24 h and 1 week (P=0.06). A decrease in angiotensin converting enzyme would also likely translate to a decrease in angiotensin II, which would have multiple effects, including a decrease in pro-fibrotic signalling within the heart (Joshi et al., 2021). This is supported in the current study by the STRING analyses that indicate a decrease in angiotensin II synthesis at 24 h (Table S3).

The differential phosphorylation of TTN was observed throughout both experiments. After 1 week of cold exposure, MuRF-1 was elevated alongside a decrease in YWHABL. One role of TTN is to detect stretch and then initiate remodelling of the myofilament via the ubiquitin proteasomal pathway (Puchner et al., 2008; Bogomolovas et al., 2014). In mammals, when the myofilament experiences increased stretch, the kinase domain of TTN becomes exposed and then binds to several components of the autophagosomal pathway in the heart, including Muscle RING-finger protein 1 (MuRF-1) (Bogomolovas et al., 2014). MuRF-1 is a ubiquitin ligase that may be used in zebrafish to recycle damaged proteins after an acute cold exposure (Tamai et al., 2022). The decrease in YWHABL might improve MAPK signalling, which has been shown to enhance MuRF-1 abundance after coronary ligation in rats (Adams et al., 2007). In zebrafish, MuRF-1 colocalizes with the M-line of the sarcomere, where the titin kinase domain is found (Li et al., 2020). The phosphorylation of BAIAP2A at 1 week in the current study would support these proteasomal pathways via G-coupled/Rho-dependent and ERK/MAPK signalling at target tissues (Stephenson et al., 2013). The cellular materials recycled by autophagosomal processes can be used to promote cell growth and proliferation (Mizushima and Levine, 2012; Yu et al., 2018; Yun and Lee, 2018). Functional studies will be required to determine whether these autophagosomal processes are used to fuel cellular remodelling during cold acclimation.

Also at 1 week, protein kinase A (PKA) was significantly increased in abundance. Elevated PKA could enable an increase in the phosphorylation of contractile proteins, such as phospholamban (PLN), TTN, MYBPC and TnT (Sucharov et al., 2011). The increase in the phosphorylation of TnT3b and TTN2 at 24 h, and TnT3b at 1 week, as well as the changes in the phosphorylation of PLN (2 and 6 weeks), TTN2 (4 and 6 weeks) and TTN1 (6 weeks) would suggest that PKA activity is important for the regulation of heart function during cold acclimation. The incorporation of TTN into this pathway for cold acclimation is a novel finding and the functional consequences of this warrants further study.

The cellular pathways most enriched at 2 weeks are associated with homeoviscous adaptation (PCYT1AB), sarcomere alignment (OBSCN, MYOM, TTN) and increased diastolic function (PLN, KCNH). Changes to these cellular processes continued through weeks 4 and 6, with TTN, MYOM, PLN and KCNH. Since the levels of sarcomeric proteins did not significantly change at 2 weeks, it is likely that regulation of myofilament function would have been via the use of the measured PTMs. Obscurin is phosphorylated following periods of maximal exercise and given its sequence homology to TTN, this phosphorylation may have stretch-sensitive functions (Potts et al., 2017). Thus, the dephosphorylated obscurin observed at 2 weeks suggest efforts to compensate for increased blood pressure caused by a decrease in temperature, to decrease cardiac load. However, functional studies will be required to confirm this. The decreased abundance of short TTN at 2 weeks would support the changes in obscurin phosphorylation. Mammals that express a greater proportion of short TTN have been found to have higher myocardial stiffness (Cazorla et al., 2000). This suggests that the comparatively lower levels of the short isoform of TTN in the hearts of the cold-acclimated zebrafish at 2 and 6 weeks may reflect efforts to maintain myocardial stiffness. Another structural component of the myofilament is kinesin, which is thought to be responsible for localizing the mitochondria near the myofilament to facilitate ATP delivery (Paulin and Li, 2004). The decrease in the phosphorylation of this protein at 6 weeks, therefore, suggests changes to the position of the mitochondria, relative to the myofilament.

We quantified elevated levels of RHOA, a Rho-dependent G-protein and MAPK3 at 24 h of cold exposure and an increase in adhesion G protein-coupled receptor 6 at 1 week. These results, coupled with the increase in MAPK3 phosphorylation at 6 weeks of cold acclimation suggest an increased capacity to respond to biomechanical stimulation and increased activation of stretch-sensitive cellular pathways. MAPKs are signalling molecules that mediate cardiac remodelling in mammalian models (Liu and Molkentin, 2016; Ren et al., 2005) and their activity in the heart is sensitive to changes in biomechanical force via membrane receptors, such as G-coupled receptors (Yamazaki et al., 1993) and receptor tyrosine kinase (RTK) (McKay and Morrison, 2007; Rose et al., 2010). RTKs activate Rho-GTPases whose activation enhances the phosphorylation of MAPKs, such as MAPK3 (Rose et al., 2010). Adhesion G protein-coupled receptor 6 also promotes the activity of G-coupled receptors in the heart (Harty et al., 2015). This potential increase in the activity of stretch activated pathways is relevant to cardiac remodelling stimulated by cold acclimation because work by Graham and Farrell (1989) demonstrates that a decrease in temperature causes an increase in the viscosity of blood. Such a change in blood parameters has been proposed to increase the biomechanical forces acting on the heart and as a result activate stretch-sensitive receptors and associated MAPK signalling cascades (Johnston and Gillis, 2022). Importantly, a recent study by Johnston and Gillis (2020) has illustrated that these pathways are stretch activated in cultured trout cardiac fibroblasts, the cells responsible for regulating collagen deposition.

The results of this study also suggest that MAPK may be interacting with the titin kinase stretch-sensing cascade. Previous work indicates that MAPK signalling is important to the response of the zebrafish heart to cold shock and cold acclimation (Long et al., 2013; Yan et al., 2020). In the current study, an increase in titin kinase phosphorylation coincided with an increase in MAPK abundance at 24 h and of MuRF-1 at 1 week of cold exposure. Phosphorylated titin kinase was observed alongside phosphorylated MAPK at 6 weeks of cold acclimation. At 2 and 4 weeks of cold acclimation, there was an increase in the phosphorylation of proteins involved in regulating structural remodelling of the cardiomyocyte, including titin kinase, without changes to the abundance or phosphorylation of MAPK. Therefore, the current study suggests that MAPK signalling is important for the cellular response to an acute change in temperature in adult zebrafish. We also provide the first evidence that titin kinase may contribute to this initial response to cold exposure and plays both an independent and integrated role over time. Further functional work is needed to explore the signalling cascade of TTN and its potential role in initiating remodelling of the heart during cold acclimation in fish.

Cold acclimation was found to decrease the abundance of mitochondrial proteins involved in oxidative phosphorylation including those associated with complex I and complex III, as well as cytochrome c oxidase and ATP synthase (Fig. S3). Simultaneously, apolipoproteins, responsible for lipid transport from the intestinal lumen to target organs and proteins for lipid utilization were increased in abundance in the second experiment (HSD17β4, PTGR, ACOX). These results suggest that cold acclimation causes a reduction in the metabolic capacity of the heart but an increase in fatty acid oxidation in the heart. This suggestion is supported by previous studies demonstrating that cold acclimation causes an elevation of β-oxidation of fatty acids in several fish species and a decrease in glycolytic capacity (Cordiner and Egginton, 1997; Driedzic et al., 1996; Guderley and Gawlicka, 1992; Little et al., 2013; Scott and Johnston, 2012; Seebacher and Simmonds, 2019). Work by Vergauwen et al. (2013) also suggests that cold acclimation of zebrafish causes an increased reliance on lipids to meet energetic demands and a reduction in carbohydrate utilization. Further work is required, however, to determine the relative role of lipids during cold-acclimated metabolism in zebrafish and how this may be affected by diet, as previous work suggests that diet can influence the tolerance of zebrafish to cold-induced ROS production (Lu et al., 2019).

In the current study, we identified increased levels of proteins involved with membrane lipid transport and the biosynthesis of polyunsaturated fatty acids (PUFAs) in the cold-exposed group at 1 week as well as in the cold-acclimated group at 2, 4 and 6 weeks. This result suggests that the cellular membranes are undergoing homeoviscous adaptation, a response stimulated by thermal acclimation that alters the proportion of mono-unsaturated fatty acids (MUFAs) and PUFAs in cellular membranes to regulate membrane fluidity (Grim et al., 2010; Guderley and Gawlicka, 1992; Robertson and Hazel, 1995; Spellner and Hazel, 1982; Gillis and Ballantyne, 1999). With cold acclimation, this results in an increase in the proportion of MUFAs and PUFAs in a membrane and studies by He et al. (2015) using tilapia suggest that these changes can occur within 5 days. An increase in PUFA content in mitochondrial membranes has been found to increase ROS production (Cosgrove et al., 1987; Pamplona et al., 1999) but can also improve oxygen capture at COX (Wodtke, 1981) and this would help combat the production of ROS at low temperatures (Ali et al., 2010).

Previous research indicates that ventilation rate can serve as a predictor of oxygen demand in fish (Frisk et al., 2012; Millidine et al., 2008; Perry et al., 2023). Such measurements in the current study suggest that the oxygen demand of the cold-acclimated fish at 20°C was ∼20% less than that of the control fish at 27°C across all time points. When considering this result, it is important to note that ventilation rate varies with tidal volume controlling water flow across the gills and that a decrease in temperature also increases the oxygen carrying capacity of water. However, the difference in ventilation rate in the zebrafish at 27°C and 20°C is equal to a Q₁₀ of ∼1.35 and this value is similar to that (1.47) reported by Vergauwen et al. (2013) for zebrafish thermally acclimated to 18°C from 26°C. Importantly, in the Vergauwen et al. (2013) study, respiration was measured using an oxygen electrode. This further supports the use of gill opercular movements to estimate respiration. Additionally, the measured decrease in oxygen demand reported here and in the Vergauwen et al. (2013) study is less than would be predicted if the response was due entirely to the impact of temperature change on biochemical rates (Q₁₀=2.2). This result suggests then that modification to relevant cellular pathways was occurring to compensate for the effect of temperature on biochemical rates. Such modifications could include changes to metabolic enzymes, contractile function of the heart and digestive processes. The cold-acclimated zebrafish maintained a consistent body condition despite eating less during the second experiment. This suggests that the resting energetic demand was being met even at low temperature. Therefore, the potential interplay between metabolic suppression, lipid metabolism and homeoviscous adaptation during cold acclimation warrants further study.

Relaxation of the myofilament is dependent upon a decrease in cytosolic Ca²⁺ following myofilament activation. One protein involved in this is SERCA (Smeazzetto et al., 2013). We saw an increase in the transcript abundance of SERCA2A at 1 week of cold exposure which was not observed in the total proteome analysis. However, because of the many steps between mRNA synthesis and protein synthesis, a change in transcription of a gene does not necessarily lead to a change in protein abundance. SERCA activity is also regulated by PLN. PLN is a brake on SERCA activity, but this is released with its phosphorylation (Weber et al., 2021; Said et al., 2003; Zhao et al., 2004). The lower levels of phosphorylated PLN at 2 and 6 weeks would therefore be expected to cause a decrease in the rate at which cytosolic Ca²⁺ returns to basal levels following contraction. Desmin, another phosphorylatable protein, also affects the activity of SERCA and the observed increase in its phosphorylation has the potential to reduce SERCA function (Chen et al., 2020; Kürekçi et al., 2021). Such a reduction in SERCA function would likely prolong and strengthen myocyte contraction as this would increase the time when Ca²⁺ is available to activate the myofilament. Work by Aho and Vornanen (1998) demonstrates that sarcoplasmic reticulum Ca²⁺ uptake is increased in the salmonid heart with cold acclimation, but is decreased in similarly treated crucian carp, Carassius carassius. This suggests that the response of the sarcoplasmic reticulum to cold acclimation is species specific but that the zebrafish might respond in a similar manner as another cyprinid fish.

Cold acclimation has also been shown to cause an increase in the level and rate of ventricular pressure generation in trout as well as the Ca²⁺ sensitivity of the contractile element (Klaiman et al., 2014). These changes in function were proposed to be caused by the observed decrease in TnT phosphorylation (Klaiman et al., 2014). The increase in TnT3b phosphorylation at 24 h and 1 week of cold exposure in the current study, suggests that this protein is also being used to regulate contractile function in the zebrafish heart.

It is important to acknowledge that our methods utilize relatively short reading frames when identifying proteins. It can be difficult to distinguish between highly similar isoforms or splice variants of proteins using only a short reading frame (Shaw et al., 2022). The results from our proteome analysis indicate that there were increases in the abundance of component proteins of the spliceosome which is responsible for generating splice variants. This has been observed previously for cold-acclimated zebrafish (Healy and Schulte, 2019). Therefore, our data may not reflect subtle changes in highly similar isoforms or splice variants that were masked by highly abundant proteins.

The results of this study demonstrate that cold acclimation is an ongoing process that occurs concurrently at multiple levels of biological organization. We saw that protein phosphorylation dominated the response at 24 h and 1 week of cold exposure, while the role of changes to protein abundance increased after 2 weeks of cold acclimation (Fig. 5). At 24 h of cold exposure, the changes in protein abundance and phosphorylation reflected activation of Ca²⁺- and stretch-sensitive pathways. Such changes persisted in the first week and would support homeoviscous adaptation, cardiac remodelling and antioxidant activity. Translational effects were evident by 2 weeks of acclimation when protein pathways involved in protein turn-over, mitochondrial degradation and cellular metabolism were enriched. It is important to consider, based on measurements of gill ventilation, that these changes occur under reduced metabolic rate. The data also suggest that cellular processes, such as the spliceosome, membrane homeostasis and energy production are important both as fast-acting modulators and as ongoing regulators. Other processes, such as the regulation of myofibrillar contractile machinery and sarcomere reorganization were important at specific times, particularly early in acclimation. These variable responses have consequences for the functional capacity of the whole animal as energy must be reallocated to facilitate the acclimation process. Recent work has demonstrated that multiple acclimation/remodelling events in zebrafish translate to reduced body condition (Shaftoe et al., 2023). The impact of physiological remodelling on energy budgets may therefore become more significant as seasonal changes to environmental temperatures become increasingly stochastic. Importantly, acute events leading to energetic limitations are central to the increased mortality characterized in migrating Pacific salmon (Eliason et al., 2017, 2013) and ‘eurythermal’ red-side dace in Ontario (Turko et al., 2020). It will, therefore, be increasingly important to consider how the cost of repeated physiological responses to changes in environmental conditions influence the energy status of natural populations.