Warming, stochastic diel thermal fluctuations affect physiological performance and gill plasticity in an amphibious mangrove fish Free

Climate change is having dramatic effects on aquatic ecosystems as temperature becomes warmer and more variable. Warming trends are being accompanied by large daily and seasonal variations that compound natural thermal variability (IPCC, 2022; Folguera et al., 2011; Stillman, 2019). For example, marine heatwaves are predicted to become more frequent and have increasing amplitude and duration, thus fuelling the likelihood that shallow-water marine ectotherms will be exposed to temperatures beyond their optima (IPCC, 2022; Oliver et al., 2019; Somero, 2010). Ectothermic animals are particularly vulnerable to increasing temperatures and thermal fluctuations because all aspects of their biology are controlled by their surrounding thermal environment. Moreover, rapid thermal shifts are predicted to pose a greater threat to ectotherms than gradual increases in mean temperature (Vasseur et al., 2014; Verheyen and Stoks, 2019).

Even though most fishes exist in thermally variable environments, much of the research studying the impacts of temperature on performance focuses on the effect of changes in constant environmental conditions over time. A growing number of studies have recognized the need to understand physiological performance under fluctuating temperatures to better represent the diel cycles of natural environments (e.g. Morash et al., 2021; Niehaus et al., 2012; Scheuffele et al., 2021; Schwieterman et al., 2022). Indeed, diel thermal fluctuations and extremes have a considerable effect on an ectotherm's physiology and survival, and must be considered if we are to understand how fish will cope with predicted climate warming scenarios and more intense weather (Burggren, 2019). Incorporating regular diel thermal fluctuations into experimental design is a more accurate representation of what occurs naturally in the wild, specifically in tidal environments with daily fluctuations in water levels, and thus temperature (Drake et al., 2017; Nancollas and Todgham, 2022). However, natural diel temperature variation is generally not regular and predictable (Nancollas and Todgham, 2022; Burggren, 2019; Dillon et al., 2016). Furthermore, it is thought that acclimation/compensation should not occur when daily thermal fluctuations are large or unpredictable (da Silva et al., 2019). Thermally variable environments across all temporal scales are thought to favour thermal generalists that function over a wide range of temperatures and where physiological rates are relatively temperature insensitive (Huey and Hertz, 1984; Somero et al., 1996). Highly variable environments associated with climate warming may thus favour some species and not others, with thermal specialists more at risk (Hayes et al., 2024; Thuiller et al., 2005). Species that cannot buffer reproduction, growth and performance will have difficulty acclimating to irregular diel thermal variation and be negatively affected. Thus, experiments that incorporate realistic, stochastic diel thermal fluctuations will lead to more accurate conservation strategies and predictions on animal thermal limits and their ability to compensate.

Despite this importance, there is a dearth of information on the effects of irregular, stochastic thermal variation on physiological performance in ectotherms. We know that metabolic and cardiac performance significantly varied between predictable and unpredictable daily temperature fluctuations in the California mussel (Mytilus californianus) (Nancollas and Todgham, 2022). Similarly, intertidal limpets (Lottia digitalis) had increased thermal tolerance in daily unpredictable thermal environments compared with predictable thermal cycling (Drake et al., 2017). In terms of molecular responses, predictable and unpredictable thermal environments had independent genomic and proteomic profiles in Drosophila simulans (Sørensen et al., 2020) that suggested distinct selection pressures in these environments. Conversely, there were no differences in growth or thermal tolerance in zebrafish (Danio rerio) acclimated to stochastic thermal variability compared with fish acclimated to predictable thermal cycles (Schaefer and Ryan, 2006). Given the complexity of fluctuating daily thermal environments where variation can exist in mean, amplitude, rate of temperature change, time and predictability, there is a need for more empirical data to provide insight into acclimation capacity in different, ecologically relevant thermally fluctuating conditions.

An increase in oxygen demand is another implication of rising water temperatures. Fishes may acclimate to a more challenging, warming thermal environment through modifications in their ventilation frequency or gill morphology (Turko et al., 2012; Wood and Eom, 2021). Increases in ventilation rate and/or depth can be used to quickly compensate for decreases in oxygen availability or increases in oxygen demand in warmer water; more water is manually moved over gills, thus picking up more oxygen. However, this mechanism is energetically costly, making it effective only as a short-term coping strategy (Turko et al., 2012). If high temperatures persist, fish may attempt to cope by reversibly modifying their gills through altering the functional gill surface area and/or the blood-to-water diffusion distance. This remodelling alters the rate of oxygen and ion transport across the gill (Wood and Eom, 2021). Gills are sometimes remodelled through changes in perfusion and the growth or regression of histologically unspecialized cells in the interlamellar regions on a gill filament (Gilmour and Perry, 2018; Nilsson et al., 2012). This group of cells is known as the interlamellar cell mass (ILCM) (Nilsson, 2007). The presence of the ILCM poses several benefits including reduced osmoregulatory costs due to a reduction in gill surface area, protection against parasites and pollutants, structural support and decreased evaporation in amphibious fishes (Gilmour and Perry, 2018). The ILCM grows slowly through a balance of mitosis and apoptosis over the course of days to weeks; however, the ILCM can be shed within minutes if oxygen demands quickly increase (Tzaneva and Perry, 2010). This dramatic reversible gill remodelling was first observed in crucian carp (Carassius carassius) in hypoxia (Sollid et al., 2003) and has since been observed in other fish species (e.g. Dhillon et al., 2013; Mitrovic et al., 2009; Tzaneva and Perry, 2010). In the amphibious mangrove rivulus (Kryptolebias marmoratus), air exposure induces an increase in ILCM height, leading to a decrease in lamellar surface area for structural support and to prevent lamellar collapse (Ong et al., 2007; LeBlanc et al., 2010; Turko et al., 2012). In addition to these oxygen-limiting environments, a shift in constant environmental temperature also results in gill remodelling in several teleost species (Sollid et al., 2005; Tzaneva and Perry, 2010; McBryan et al., 2016; Chen et al., 2019). Where measured, there was a regression of the ILCM during exposure to warming, constant temperatures, increasing functional gill surface area. To our knowledge, there are no reports of gill remodelling in response to either regular or irregular, stochastic daily fluctuating temperature.

We capitalized on the amphibious nature and distinct reproductive strategy of the self-fertilizing mangrove rivulus (hereafter ‘rivulus’) to determine the consequences of regular and stochastic thermal variation on growth, fecundity, thermal biology and gill morphology in genetically homozygous lineages of fish. Rivulus originate from mangrove forests along the Central American and northern South American coastlines, where they experience extreme environmental conditions (Avise and Tatarenkov, 2015; Taylor, 2012). This unique fish is tolerant of temperatures ranging from 7°C to 38°C (Taylor and Davis, 1995); however, as a tropical species with a restricted habitat range, it has a narrower thermal tolerance range compared with temperate species (Nati et al., 2021). When water conditions become unfavourable in cases such as habitat drying, low oxygen, predation and/or high temperatures, rivulus may emerge onto land (Turko and Wright, 2015). Rivulus reproduce through selfing where, in theory, all individuals can reproduce independently (Avise and Tatarenkov, 2015). Inbred laboratory populations produce highly homozygous lineages, allowing us to control for inter-individual genetic variation when studying acclimation potential.

Using three distinct isogenic laboratory lineages of rivulus that have been inbred for over 60 generations, we investigated phenotypic changes associated with diel thermal variation. Each lineage is highly homozygous but lineages are genetically distinct (Tatarenkov et al., 2010; Kelley et al., 2016), allowing us to determine whether there is a genetic component in any of our response variables. Our overarching hypothesis was that rivulus exposed to warming, stochastic daily thermal fluctuations would not acclimate as effectively as fish experiencing regular, non-stochastic daily thermal fluctuations. In daily stochastic conditions, fish could spend variable proportions of time above and below their preferred temperatures, possibly leading to diminished performance. We acclimated fish to three different thermal conditions: (1) regular, non-stochastic and (2) irregular, stochastic diel thermal fluctuations with similar heat loads, and (3) a stable, control temperature representing the lower limit of the variation. We measured survival and fecundity, growth rate and measures of gill morphological variation (i.e. ILCM coverage, lamellar surface area). We also used the tendency of this fish to emerge from warming water (Gibson et al., 2015; Currie and Tattersall, 2018) to compare emersion temperature (T_em) and critical thermal tolerance (CT_max). We predicted a lower growth rate, reproductive output and thermal tolerance (CT_max and T_em) in fish exposed to stochastic compared with non-stochastic thermal fluctuations, and that these response variables would be lower in both types of diel fluctuation compared with the stable, lower temperature. We asked whether exposure to fluctuating, warming thermal environments would increase gill surface area through remodelling to ultimately maximize oxygen uptake, and whether the nature of daily thermal variability would influence gill plasticity. We predicted that gill surface area would increase in both warming, daily fluctuations to maximize oxygen uptake, but that we would see differences between stochastic and non-stochastic diel cycling temperatures.

We housed three distinct isogenic lineages of Kryptolebias marmoratus (Poey 1880) in a breeding colony at the Animal Care Facility at Acadia University, Wolfville, NS, Canada: (1) Honduras 11, originally caught in the Bay Islands, Honduras (1996); (2) Twin Cayes, originally caught in Twin Cayes, Belize (1991); and (3) Dangriga, originally caught in Dangriga, Belize (2006). These three lineages are genetically divergent and have been bred in the laboratory for over 60 generations with multiple loci differences (Tatarenkov et al., 2010). In all lineages, the conditions of their original collection site were similar (i.e. Belize and Honduras) and given their generation time in the laboratory, any environmental effects of their origin have likely been removed through successive generations of inbreeding. The main goal for using distinct lineages was to determine whether there may be a genetic component to our response variables and to determine environmental effects independent of genetic variation. We held fish individually in 120 ml cups filled with 80 ml of 15 ppt water made using Instant Ocean (Pets and Ponds, ON, Canada) and reverse osmosis water. The cups were held in a room under controlled conditions (12 h:12 h light:dark photoperiod, 30% humidity, 27±2°C). We fed fish live Artemia 5× per week and changed water every 2 weeks. We used adult hermaphrodites between 6 months and 1 year of age. All experiments were approved by the Acadia Animal Care Committee in accordance with guidelines set by the Canadian Council on Animal Care (protocol no. 01-21).

For the diel fluctuating temperatures, we acclimated fish in individual cups in a two-tiered Conviron growth chamber (CMP 5090, model E7/2) programmed with a 4-day thermal profile that repeated throughout the acclimation period (Fig. 1). We generated sinusoidal diel temperature fluctuations by programming the growth chamber to heat the air with consistent ramping to our thermal set points (Fig. 1) and programmed a 12 h:12 h light:dark photoperiod. In one tier of the growth chamber, there was a temperature differential of 2°C between air and water; thus, we programmed the chamber appropriately to allow us to achieve our set thermal profiles.

We monitored both the air and water temperature at least 2× per day with Temp Sticks (tempstick.com; air) and HOBO temperature loggers (HOBO Pendant MX 2201; water) to ensure temperatures were accurate and were within ±1°C. We acclimated our two thermal experimental groups where temperature fluctuated on a 24 h diel cycle for at least 5 weeks: (1) regular diel fluctuation (n=13 for Honduras 11, Twin Cayes and Dangriga), (2) stochastic diel fluctuation (n=13 for Honduras 11, Twin Cayes and Dangriga) and (3) our control group with no fluctuation (n=13 for Honduras 11, Twin Cayes and Dangriga). The first group had temperature fluctuating regularly and predictably around a diel temperature cycle of 27–35°C, with regularly spaced thermal maxima, where the maximum temperature was reached at noon and the minimum temperature at midnight every day. In the second thermal group, temperature fluctuated randomly/stochastically between 27°C and 35°C with irregularly spaced thermal maxima and the maximum and minimum temperature reached at a different time each day (maximum between 06:00 and 18:00 h, minimum between 18:00 and 06:00 h). We set the thermal maxima to vary randomly over the 4-day repeated profile, setting each daily maxima during the light phase of our acclimation (i.e. 06:00, 09:00, 12:00 and 15:00 h). In both diel fluctuating groups, the ΔT was 8°C with a consistent ramp to each programmed set point (Fig. 1). In the regular, non-stochastic fluctuating group, the thermal load, or cumulative degrees, over the 4-day profile was 1064°C (calculated as total °C as in Tunnah et al., 2017). The stochastic diel temperature cycling group had a similar thermal load over the same 4-day time span (1072°C). In the control group, the lower end of the cycle and our colony temperature, the thermal load over 4 days was 972°C. These diel temperature cycles and variations are all environmentally relevant and regularly experienced by these fish in nature (Ellison et al., 2012; Rossi et al., 2019). Water ammonia and nitrate concentration (mg l⁻¹; API test kits) remained in range throughout the experiment. Dissolved oxygen was also measured throughout the course of the experiment (in both the growth chambers and the water bath; see Thermal tolerance) and was never found to be limiting (74.3–96.5%). In their native mangroves, these fish are often found in ephemeral mangrove pools and will experience intermittent hypoxia (Rossi et al., 2019).

During acclimation, fish from both diel fluctuating temperature groups were scattered randomly amongst the acclimation chambers (Conviron growth chambers) to ensure a lack of bias based on location in the chamber (e.g. middle compared to perimeter). The control fish, at a consistent and stable 27°C (rearing temperature), were kept in the fish colony. We changed the water in the cups weekly throughout the duration of the experiment (total of 15 weeks).

We calculated survivorship as the percent of fish living in each lineage in each thermal group per week (number of remaining fish/total number of fish in group×100%).

We checked and recorded the number of embryos laid by each fish daily, removed them to a separate container and monitored the embryos for hatching success. We measured fecundity at the end of the 15-week acclimation period as the number of embryos produced per lineage in each acclimation group. We also documented hatching success; a hatch was considered a success if the juvenile lived beyond 30 days after hatching.

Prior to placement in their respective chamber, we weighed and measured fish to the nearest mg and mm, respectively (Fig. 2). Body length was measured as fork length (tip of the snout to fork of the tail). We weighed and measured fish again at the end of the experiment (15 weeks; Fig. 2). We calculated growth rate as the difference in mass divided by the time span (mg day⁻¹).

After 5 weeks of acclimation in the chambers, we began testing fish for their emersion/avoidance temperature (T_em) and CT_max. The Honduras 11 lineage was tested between days 33 and 40. The Dangriga lineage was tested between days 75 and 97. We did not measure these variables in the Twin Cayes lineage (see Results). We included acclimation time and body size as co-factors in our statistical analyses and determined that there was no effect of either factor (acclimation time: P=0.33 for T_em and P=0.18 for CT_max; body size: P=0.24 for T_em and P=0.92 for CT_max). We recorded T_em as the temperature at which the fish's gills were completely emerged from the water. To measure T_em, we removed fish from their fluctuating or control conditions and placed fish in an open-top rectangular chamber (7.5×7×4 cm with 2 cm surrounding platform) with ∼27°C water (min. 24.8°C, max. 28.5°C, mean 26.2°C). We built 3 cm ‘walls’ around the side of the chamber to allow the fish to emerge from the chamber but prevent them from emerging completely into the water bath. The chamber was placed in a water bath at 27°C and fish were habituated for 1 h. We then warmed at a constant rate of 1°C min⁻¹ (Gibson et al., 2015; Currie and Tattersall, 2018). The temperature in the chamber was monitored with Logger Pro software (Vernier) for each experiment and the actual rate was 0.83–1.1°C min⁻¹. To minimize any effects of daily rhythms on CT_max measurements (Fangue et al., 2006), all trials took place between 11:00 and 16:00 h and lasted 10–19 min (depending on T_em and CT_max). We recognize that this design introduced a trade-off in that fish were transferred to 27°C water from their fluctuating, experimental conditions at different temperatures. Thus, immediate thermal history, the difference between the experimental temperature and 27°C, varied from 0 to 8°C. Our focus here was on the effect of chronic acclimation to diel cycles; however, we cannot rule out an acute thermal effect from immediate thermal history. If the temperature transfer had a large effect, we would have expected to observe more variability in our data than we did.

We measured CT_max as the temperature at which the fish lost equilibrium for at least 2 s or as the temperature at which the fish fully exposes its ventral surface (Currie and Tattersall, 2018; Lee et al., 2016). We calculated CT_max using the same protocol as T_em, heating the water at a rate of 1°C min⁻¹, as above, but we removed the 3 cm ‘walls’ and placed perforated plastic wrap over the chamber to prevent fish from accessing the water surface and emerging.

To prevent bias towards the order of testing, half of the fish from each temperature group were tested for T_em first, while the other was tested for CT_max first. We did not observe a significant effect of testing order (T_em first versus CT_max first) (t-test; P_Tem=0.714, P_CTmax=0.272); thus, we removed order as a factor from further analyses. The time between T_em and CT_max tests was 48 h for all fish, and they were fed normally during this time and thus not postabsorptive prior to either T_em or CT_max measurements. Fish were tested once for each variable. As with T_em, these trials took place at the same time each day, between 11:00 and 16:00 h. Fish spent ∼2 h 15 min in these trials including habituation to experimental chamber (1 h), T_em/CT_max trial (10−19 min) and recovery (∼1 h) before being returned to their acclimation chamber.

Gill histology was analyzed in the Honduras 11 lineage only, as this lineage consistently showed differences between thermal groups (control: n=11, regular: n=10, stochastic: n=13). Following final mass and length measurements after 15 weeks, we euthanized fish and fixed the tissues by immersing whole fish in 10% neutral buffered formalin (Fisherbrand) for 24 h at 4°C. We decalcified tissues by placing the whole fish in Surgipath Decalcifier II for 1 h at room temperature (20°C) and then transferring them to 70% ethanol at 4°C for 1 week. Before dehydrating, we removed the tail of the fish, cutting at the vent transversely, to decrease fish length and minimize vibrations while sectioning. Samples were dehydrated in 85% ethanol for 3 h, 95% ethanol for 1 h and then 100% ethanol for 30 min and cleared with two changes of toluene (Grade ACS, Bebbington), each for 1 h. Using a vacuum oven (SoTemp, Model 280A), we infiltrated samples three times in paraffin at 56°C (Paraplast Plus), for 20 min each at ∼400 mm Hg. We embedded the samples in paraffin (Thermo electron corporation Histocentre 3) and sectioned them at 6 μm using a rotary microtome (Spencer). Sections were mounted on poly-l-lysine subbed slides, floated using de-gassed reverse osmosis water, and dried overnight (37°C) on a slide warmer. Sections were stained with haematoxylin and eosin (Fisher Scientific).

We were blinded to fish identification for measurements of lamellar length/width and ILCM height. We took measurements (ILCM height, lamellae width and length) from 40 regions of the sections per fish, starting with the gill bars, where the pseudobranchia were first observed (Fig. S1A). To maintain consistency, we only took measurements from the innermost filaments on both the left and right side of the gill chamber (Fig. S1), and only from the first four lamellae or ILCM on a filament (Fig. S1B,C). We excluded all distorted filaments. All measurements were taken to the nearest 1.25 μm. We measured lamellae length starting from the edge adjacent to filament to the base of the filament. ILCM height was measured parallel to the total lamellar length, starting from the edge of the ILCM bordering the filament to the most distal edge of the ILCM from the filament. We measured lamellae width parallel to the filament at the base (Ong et al., 2007). We used a Nikon Y100 compound microscope with an Am-Scope camera to observe and measure ILCM height, lamellae length and width, and calculated the ratio of the ILCM height relative to lamellar height.

We measured ventilation frequency by habituating Honduras 11 fish to smaller 50 ml cups with 30 ml of 15 ppt seawater at 27°C under a dissecting scope for 5 min until they were swimming normally. We then counted the number of times the operculum opened and closed per 15 s period. We did not control the temperature during these measurements but treated each group the same and reasoned that the water would not have changed appreciably in this time frame.

We completed all statistical analyses for physiological measurements using RStudio version 1.1.463 with an alpha value of 0.05. To determine whether the errors were normally distributed, we plotted residuals and fitted values; all data were homoscedastic, thus, the statistical test was accepted. To determine whether there were differences in body size among lineages or treatments, we used a two-way ANOVA followed by a Tukey HSD test to determine specific differences. To compare differences in growth rate and thermal biology (T_em and CT_max) we used a two-way ANOVA with two fixed factors (temperature and lineage) and acclimation time or body size as a co-factor. There was no statistical interaction between thermal condition and lineage in any of the traits (P>0.05); however, we observed significant differences between lineages and experimental groups. When significant differences were detected among groups and/or lineages (P<0.05), we used a Tukey HSD test to determine specifically which groups or lineages were statistically different from one another.

We analysed gill histology with a generalized linear mixed model (GLMM; SPSS 28) to compare the lamellar surface area, ILCM height, lamellar length and ILCM:lamellae height ratio among treatments (α=0.05). Subjects were individual fish, and the 40 measurements per variable taken from each fish were analysed as repeated measures. Fixed effects were treatment group (i.e. control, predictable, unpredictable). Significant differences (α=0.05) amongst groups were estimated with a post hoc pairwise comparison. The data for gill morphology were not normally distributed, as visible in the violin plots (Fig. 6). We analysed ventilation frequency with a one-way ANOVA followed by a Tukey HSD test to determine which acclimation groups were statistically different from one another. We created all graphs using Prism 9 (version 9.4.1).

Survivorship was influenced by experimental temperature. One hundred percent of fish survived acclimation to control and regular, non-stochastic diel fluctuating temperatures, and from the Honduras 11 and Dangriga lineages at all temperatures. The only group experiencing mortality was the Twin Cayes lineage from the stochastic diel fluctuating condition. After 4 weeks, survivorship began to decline and after 15 weeks, survivorship dropped precipitously to 62% and we terminated the experiment for the Twin Cayes lineage.

At the beginning of the experiment, Honduras 11 lineage was significantly larger than the Dangriga and Twin Cayes lineages (two-way ANOVA, P<0.001), and these two lineages were not different from each other. We did not detect body size differences among thermal treatments. Fecundity, measured as the number of embryos laid per individual, and hatching success, the percent of embryos laid that hatch and mature to 30 days, varied with thermal acclimation (Fig. 3). We did not find any embryos from fish acclimated to stochastic temperature fluctuations in any lineage or at any of the acclimation temperatures in the Twin Cayes lineage throughout the 15-week acclimation. Like Twin Cayes, the Dangriga lineage did not produce embryos under control conditions. We only observed hatching success in the Honduras 11 lineage under control and regular diel fluctuating conditions (Fig. 3).

Thermal variability had a significant effect on growth rate in all three lineages (Fig. 4). Interestingly, growth rate was affected by both temperature (two-way ANOVA, P<0.001) and lineage (two-way ANOVA, P<0.001), with no detectable interaction (P>0.5) (Fig. 4). In the Honduras 11 lineage, growth rate at the stable control temperature was significantly greater than during stochastic, diel temperature variation (Tukey HSD, P<0.001; Fig. 4A). Fish acclimated to regular diel thermal cycles were not significantly different from fish acclimated to stochastic diel cycles or control temperature (Tukey HSD P>0.05). In the Twin Cayes lineage, there was also a significant difference in growth rate between control and stochastic temperature fluctuation (Tukey HSD, P<0.01; Fig. 4B), with fish experiencing stochastic temperature cycles having the lowest growth rate. Fish acclimated to regular temperature cycles were not significantly different from fish acclimated to either control or stochastic temperatures (Tukey HSD, P>0.05). Finally, in the Dangriga lineage, growth rate was significantly different between fish acclimated to control temperatures and stochastic diel temperature cycles (Tukey HSD, P<0.01; Fig. 3C), but as in the other two lineages, we did not observe differences between fish acclimated to regular temperature cycles compared with either control temperatures or stochastic diel cycles (Tukey HSD, P>0.05).

We investigated thermal tolerance in only the Honduras 11 and Dangriga lineages given the unexpected mortality we observed in the Twin Cayes lineage. Thermal acclimation had a modest but insignificant effect on T_em (two-way ANOVA, P=0.067; Table 1), with stochastic temperature cycles higher than control temperature (P=0.06). There was no effect of lineage on T_em (two-way ANOVA, P>0.1). We observed significant differences in CT_max among temperature acclimation groups (two-way ANOVA, P<0.001) and among lineages (two-way ANOVA, P<0.001). In the Honduras 11 lineage, fish acclimated to regular temperature cycles had a significantly higher CT_max (40.9°C) than fish at control temperatures (39.8°C; Tukey HSD, P<0.01; Table 1). In fish acclimated to stable control temperatures, CT_max in the Honduras 11 lineage was significantly lower than that of the Dangriga lineage (Tukey HSD, P<0.01; Table 1). We did not observe differences between fish acclimated to stochastic and non-stochastic fluctuating temperatures in either T_em or CT_max.

We used the pseudobranch (Fig. S2) to mark where in the gill chamber we started measuring lamellar length/width and ILCM height. There was an overall decrease in respiration surface area for fish in both thermally fluctuating conditions when compared with control (Fig. 5). This surface area was reduced by both a decrease in lamellar surface area (Fig. 6A; see below) and an increase in ILCM (Fig. 6B). The regular thermal cycling condition had the smallest lamellar surface area, and the stochastic diel fluctuations had the highest ILCM height; therefore, there was no overall difference in overall gill surface area between stochastic and non-stochastic diel fluctuations (Figs 5, 6). Lamellar surface area was significantly different among temperature conditions (P<0.001; Fig. 6A; Table S2) and lower in both regular (P<0.001) and stochastic (P<0.001) cycling conditions compared with the control. Comparing the two thermally fluctuating groups, the lamellar surface area was significantly higher in the stochastic compared with the non-stochastic temperature cycling condition (P=0.007; Fig. 6A).

ILCM height was also significantly affected by temperature acclimation (P<0.001; Fig. 6B; Table S3). When compared with the control condition, ILCM height was significantly greater in both the warm, regular (P<0.001) and stochastic (P<0.001) fluctuating temperatures. ILCM height was significantly greater in the stochastic cycling condition when compared with non-stochastic, regular thermal cycling (P<0.01).

Lamellar length was significantly different among the three temperature conditions (P<0.001; Fig. 5C; Table S4;). Lamellae were significantly shorter in the regular, stochastic diel fluctuations when compared with control (P<0.001) and shorter in the non-stochastic compared with the stochastic cycling condition (P<0.001).

We expressed ILCM height relative to lamellar height as a ratio to determine ILCM coverage and found that it significantly differed among groups (P<0.001; Fig. 6D; Table S5). ILCM:lamellae ratio was significantly greater in both the stochastic and non-stochastic diel thermal cycling conditions when compared with the control (P<0.001; Fig. 6D). The two thermally fluctuating conditions were not significantly different from each other.

Ventilation frequency in the Honduras 11 lineage was significantly higher in both stochastic and non-stochastic thermal fluctuations when compared with the stable control temperature (one-way ANOVA, P<0.001; Fig. 7). As above, there was no difference between the two thermally fluctuating conditions.

The goal of our study was to investigate the effects of diel temperature variability on measures of fitness, thermal tolerance and gill morphology independently of changes in mean temperature. We reasoned that rivulus acclimated to irregular, stochastic diel thermal cycles would not compensate as well as rivulus in regular, non-stochastic conditions, given the likely varying proportions of time spent outside their preferred temperature and the greater vulnerability of tropical species to variable and extreme warming (Nati et al., 2021). Thus, we predicted lower growth rate, reproductive output and thermal tolerance in fish in the stochastic thermal cycle compared with the regular fluctuations. Our hypothesis was supported in that fish exposed to stochastic diel thermal cycles had lower reproductive output than fish in regular, diel thermal fluctuations. Stochastic fluctuating diel temperatures significantly reduced growth rate compared with the control temperature, and thermal tolerance modestly increased only in regular, warming diel fluctuations. We were also interested in how thermal fluctuations would affect gill morphology given that acclimation to warm, constant temperatures significantly increases gill surface area (McBryan et al., 2016; Potts et al., 2021; Sollid et al., 2005; Tzaneva and Perry, 2010). We predicted that both warming fluctuating conditions would lead to an increase in gill surface area, with differences between stochastic and non-stochastic cycling temperatures. In contrast to our prediction, acclimation to both warming conditions significantly reduced gill surface area, possibly for structural support in these amphibious fish (see below). Stochastic diel thermal fluctuation significantly increased ILCM height compared with regular, non-stochastic temperature fluctuation, suggesting that the nature of the thermal cycle affects gill remodelling.

Our experimental design used ecologically relevant warming, diel fluctuations; thus, we did not expect to see any effects on survival in any of our lineages. However, the Twin Cayes lineage experienced mortality in the stochastic thermal fluctuations. This lineage-specific response suggests that there is a genetic component in coping with these warming conditions. We know that there is variability in stress tolerance in rivulus. For example, emergence following exposure to hydrogen sulphide (H₂S) was influenced by the fish's genetic lineage, with the Honduras 11 lineage having a greater sensitivity to H₂S when compared with Twin Cayes (Martin and Currie, 2020). Additionally, fish originating from the Honduras 11 lineage tolerate air exposure (emersion) significantly longer than fish originating from Belize (such as the Twin Cayes lineage), suggesting that rivulus from different lineages have genetically based differences in acclimatory capabilities and/or strategies (Turko et al., 2019).

We also observed lineage-specific differences in fecundity. Neither the Dangriga nor Twin Cayes lineage deposited embryos at the control temperatures, suggesting that these fish were not fecund at that time. Thus, we cannot accurately conclude how thermal fluctuations affected reproductive output in these two lineages. Previous studies have shown that, in general, K. marmoratus follow a temporal reproductive pattern that is not strictly seasonal (Lomax et al., 2017). The Honduras 11 lineage was significantly larger than the other two lineages, which may partly explain the higher fecundity. In this lineage, regular diel temperature cycling led to significantly fewer embryos compared with control temperatures, and in stochastic diel thermal cycling reproduction was inhibited. When California mussels (M. californianus) were exposed to unpredictably fluctuating thermal environments, they had significantly more glycogen stores than their predictable counterparts, suggesting the possibility of different strategies for energy allocation (Nancollas and Todgham, 2022). Given the lower fecundity in warm, diel thermal cycling and the inhibition of reproduction in stochastic diel cycles, it is possible that during times of environmental or thermal uncertainty, energy stores are mainly used for maintaining critical life functions.

In all three lineages, only fish acclimated to stochastic diel cycles had significantly reduced growth rates compared with the stable control temperature of 27°C. At high temperatures, the metabolic cost of surviving increases; therefore, energy available for growth decreases (Baudron et al., 2014; Ikpewe et al., 2021; Neuheimer et al., 2011). Thus, our reduced growth rate with warm, diel cycling may not be surprising. Both diel cycles were warmer than the control with a mean of 31°C, yet only stochastic diel cycling affected growth rate. Reduced growth rates have been reported in fish acclimated to thermal fluctuations compared with constant temperature conditions (e.g. Chadwick and McCormick, 2017; Morash et al., 2018), but these growth effects can be species and population dependent and will be influenced by the magnitude of thermal variation and the thermal optimum. Our experimental design was focused on comparing the effects of stochastic and non-stochastic diel thermal variability, rather than temperature constancy. In contrast to our results, growth rates in stochastic diel thermal cycles were similar to those in regular, variable diel thermal fluctuations in zebrafish; however, the average daily cycles in these groups were different (28±3°C and 28±4°C, respectively; Schaefer and Ryan, 2006). Like our study, Pisano et al. (2019) compared growth rate in brook trout (Salvelinus fontinalis) acclimated to regular, non-stochastic and stochastic thermal fluctuations (both high and low stochasticity). They determined that fish were longer in the low-stochasticity conditions compared with non-stochastic cyclical temperatures. These authors concluded that the influence of diel thermal fluctuations (stochastic and non-stochastic) on fitness are not readily predictable, at least partly because thermal performance curves are plastic (Pisano et al., 2019). Given the reduced fitness metrics we observed with stochastic diel temperature variation (i.e. reproduction and growth rate), it is possible that fish spent reduced time at their T_opt (the optimum temperature where performance is maximized) in the stochastic thermal profile compared with the regular fluctuating profile.

We used the temperature at which individuals emerged from warm water (T_em) as a proxy for pejus temperature when performance begins to decline (Pörtner et al., 2006; Melanson et al., 2023). Our fluctuating diel cycles had no significant effect on pejus temperature and modestly increased CT_max in the regular diel fluctuating environment. Increases in critical thermal tolerance with warming (to a limit) is not surprising and has been shown for other fishes (e.g. killifish, Fangue et al., 2006; salmon, Anttila et al., 2015; zebrafish, Morgan et al., 2020). Notably, both fluctuating cycles had a similar thermal load, or cumulative degrees, but we only observed a significant increase in CT_max in the regular diel cycle. In the fingered limpet (Lottia digitalis), thermal tolerance did not vary between predictably and unpredictably fluctuating temperatures with the same thermal load (Drake et al., 2017). Stochastic diel temperature fluctuations enhanced upper thermal tolerance of M. californianus as determined by cardiac performance (Nancollas and Todgham, 2022). In fluctuating temperature conditions, the recent thermal history of the organism may play a key role in determining thermal tolerance. For example, limpets in unpredictable cycles that ended with high temperatures had different upper thermal limits than limpets exposed to cycles ending in lower temperatures (Drake et al., 2017). Our experimental design was focused on chronic acclimation to diel thermal cycles and did not allow us to determine the effects of immediate thermal history on our variables. This would be an interesting area for future study.

We know that alterations in ventilation frequency can help mitigate the short-term effects of warming temperatures as oxygen demand increases. Indeed, we observed a significant increase in ventilation frequency in the Honduras 11 lineage with both stochastic and non-stochastic warming cycles. Given these increases with thermal fluctuations in this lineage, we predicted a regression in the ILCM between lamellae with warm temperature would underpin these increases, to increase gill surface area and ultimately oxygen uptake capacity, as has been shown in several fish species (e.g. killifish, McBryan et al., 2016; carp, Sollid et al., 2005; goldfish, Tzaneva and Perry, 2010). Notably, the experimental design in these studies acclimated fish to stable temperatures before measuring gill morphology. In contrast to our prediction, we observed an increase in ILCM and a decrease in lamellar surface area after exposure to warming fluctuations when compared with the stable, cooler control temperature. It is important to note that oxygen was not limiting in our study, despite an increased oxygen demand. At this point, it is not clear why we see this surprising decrease in gill surface area with warm diel fluctuations, but we suggest a few possibilities.

Increases in temperature increase metabolic oxygen demand. The fish gill responds through increases in effective gill permeability to enhance oxygen uptake (Randall et al., 1972). According to the osmo-respiratory compromise, these changes to gill oxygen permeability will also increase permeability to water and ions (Gilmour and Perry, 2018; Wood and Eom, 2021), and there is significant empirical evidence of these temperature-induced changes on gill permeability (e.g. Evans, 1969; Giacomin et al., 2017; Onukwufor and Wood, 2022). Studies on the role of thermal fluctuations on the osmo-respiratory compromise are scarce. It is possible that a thermally fluctuating, high temperature environment exacerbates gill water and ion permeability beyond changes observed under stable warm environments, such that covering the gills with ILCM and reducing surface area is energetically necessary and a trade-off for oxygen uptake. Thus, it may be beneficial for some fish to decrease their gill surface area, as we observed here, to mitigate temperature-induced ion imbalances and water loss. Mangrove rivulus can also rely on ionocytes in their skin for osmoregulation (LeBlanc et al., 2010), so covering the gills (and possibly gill ionocytes) is not as critical in this amphibious fish as in aquatic fishes. Given their ability to use cutaneous respiration, an alternate possibility is that the osmorespiratory compromise is negated, and will not occur, in amphibious fish such as mangrove rivulus, that osmoregulate across the skin.

It may be that an increase in ILCM with warming is a consequence of the amphibious nature of and is a preparatory response for emersion. Prior to and during emersion, rivulus experience reversible gill remodelling when they emerge onto land and switch to cutaneous respiration (Ong et al., 2007; Turko et al., 2012). These structural changes occur even when emersion is of short duration. Challenges associated with emergence include structural support of the gills (i.e. no longer in neutral buoyancy) and evaporative water loss. Indeed, rivulus decrease their gill surface area during emergence through a filling in of the ILCM between lamellae (Ong et al., 2007; Turko et al., 2012). The ILCM acts as a protective structure to support lamellae, aiding in aerial respiration (i.e. prevents gill collapse) and preventing desiccation (Ong et al., 2007). Other strategies for preventing lamellar collapse in amphibious fish include short, widely spaced lamellae in the mudskipper Boleophthalmus chinensis (Tamura et al., 1976), thick, long and enlarged gills in the rockskipper Mnierpes macrocephalus (Graham, 1970), and permanent fusions between lamellae in the mudskipper Periophthalmodon schlosseri (Kok et al., 1998; Wilson et al., 1999). To our knowledge, these structural modifications are not observed in rivulus, thus these fish would rely on the growth of ILCM for support and possibly stouter lamellae, as we observed. The reduced surface area would also reduce evaporative water loss on land (Ong et al., 2007). Here, rivulus acclimated to stochastic and non-stochastic warming fluctuations were in environments approaching their emersion threshold; thus, it is possible that they remodelled their gills to better suit impending aerial respiration.

The goal of our study was to compare the physiological effects of warming, stochastic and non-stochastic diel temperature cycles in a mangrove fish, which would experience such fluctuations in nature. When considering how performance is affected by temperature variation, if the upper and lower limits of the fluctuations reach the boundaries of the their thermal scope, then organisms could be become temporarily stressed and performance can be negatively affected (Deutsch et al., 2008; Folguera et al., 2011). Alternatively, heating and cooling phases of fluctuations could cancel each other out (Salachan and Sørensen, 2017), with no deleterious effects on performance. Our data suggest that amphibious fish acclimated to stochastic, warm thermal cycling do not compensate as well as fish in regular diel cycles. There is a decrease in reproductive output and growth rate compared with regular, thermal fluctuation. Not surprisingly, thermal tolerance increases with warming, but only in regular diel cycles. To our knowledge, this is the first study to investigate gill remodelling with diel thermal cycling, and our observed increase in ILCM and consequent decrease in gill surface area with warm, temperature fluctuation in this amphibious fish suggest unique gill plasticity to warming and/or thermal fluctuations. At this point, it is unclear whether such gill remodelling is adaptive in these warm, temperature cycles. We have also determined that there is likely a genetic component to how fish cope with thermal stochasticity given the differences we observed among lineages in some traits. Understanding the fitness consequences of warming, diel fluctuations will allow us to more accurately predict the capacity of fish to cope with thermal heterogeneity and increasing temperature challenges and extremes.

We are grateful to Dawn Miner and several undergraduate students at Acadia University for assistance with animal care and to Dr T. Avery for statistical assistance.