Temperature, species identity and morphological traits predict carbonate excretion and mineralogy in tropical reef fishes

Animal collection and holding for this project was conducted under Marine Research Permit RE-19–28 issued by the Ministry of Natural Resources, Environment, and Tourism of the Republic of Palau (10.03.2019), Marine Research/Collection Permit and Agreement 62 issued by the Koror State Government (08.10.2019), Queensland Government GBRMPA Marine Parks Permit G14/36689.1, Queensland Government DNPRSR Marine Parks Permits QS2014/MAN247 and QS2014/MAN247a, Queensland Government General Fisheries Permit 168991, Queensland Government DAFF Animal Ethics approval CA2013/11/733, approval by The Bahamas Department of Marine Resources, approval by the Animal Care Officer of both the University of Bremen and the Leibniz Centre for Tropical Marine Research (ZMT), and in accordance with UK and Germany animal care guidelines.

Sample collection

We collected fish carbonate samples at four study locations across three tropical and subtropical regions: Eleuthera (24°50’N, 76°20’W), The Bahamas, between 2009 and 2011^27,37; Heron Reef (23°27’S, 151°55’E) and Moreton Bay (27°29’S, 153°24’E) in Queensland, Australia, in 2014 and 2015²⁸; and Koror (7°20’N, 134°28’E), Palau, during November and December 2019. These are located within four distinct marine biogeographic provinces and three realms (Tropical Atlantic, Central Indo-Pacific, and Temperate Australasia)⁴³. At each location fish were collected using barrier nets, dip nets, clove oil or hook and line, and immediately transferred to aquaria facilities at the Cape Eleuthera Institute, Heron Island and Moreton Bay Research Stations, and the Palau International Coral Reef Center. Fish were held in a range of tanks (60, 400, or 1400 L in the Bahamas, 10, 60, 100, 120, or 400 L in Heron Island and Moreton Bay, and 8, 80, 280, or 400 L in Palau) of suitable dimensions for different fish sizes (<1 g to 11 kg), either individually or, for particularly social species, in small groups of similar sized individuals of the same species. All tanks were supplied with flow-through locally-drawn filtered (1–5 µm) natural seawater, except in Moreton Bay where we used locally-drawn filtered natural seawater in a recirculation system, and maintained at ambient conditions. Food was withheld throughout the sampling period (typically three days but sometimes shorter or longer) and for at least 48 h prior to ensure that sample material comprised only carbonate precipitated within the intestine from imbibed seawater calcium, rather than from dietary sources. Additionally, each tank was fitted with a false mesh bottom to prevent further disturbance of excreted carbonate pellets or potential ingestion by fish. Carbonate pellets were collected from the bottom of the tanks using a siphon or disposable Pasteur pipette at 24 h intervals in The Bahamas and Australia (except a few non-scarine Labridae that were sampled at 4 h intervals—see ref. ²⁸), and at 8 h intervals in Palau. Samples were rinsed three times with deionised water (centrifuging each time for 3 min at 2655 x g) to remove saltwater and excess salts and soaked in sodium hypochlorite (commercial bleach; < 4% available chlorine) for 6–12 h to disaggregate organic material⁸³. All traces of bleach were removed with further rinses with deionised water before drying the samples for 24 h at 50 ˚C. Full details of carbonate collection in The Bahamas and Australia are described in refs. ^27,28,37.

Sample analysis

Carbonate composition

Samples were characterised for their morphological and mineralogical composition using scanning electron microscopy (SEM), energy-dispersive X-ray spectroscopy (EDX), and Fourier-Transform Infrared spectroscopy (FTIR). Detailed procedures for samples collected in The Bahamas and Australia are described in refs. ^28,37. The same methodology was applied to samples collected in Palau. Morphological and chemical composition were analysed using a Tescan Vega 3 XMU SEM with integrated Oxford Instruments X-MAX EDX detector. Dry samples were mounted on adhesive carbon tape and covered with a 20 nm conductive coating (Au). Morphological observations were made on at least five pellets per sample and electron microscope images were acquired at accelerating voltages between 5 and 15 kEV and working distances of 7–12 mm using either secondary electron or backscatter detectors. A minimum of 30 EDX scans were performed on each sample, incorporating all present particle morphotypes, using an accelerating voltage of 20 kEV, a working distance of 15 mm, and acquisition time of > 40 s. Scans were only performed on particles surrounded by others of similar morphology, or on particles of sufficient size to ensure that data were representative of a particular morphology.

Mineralogical composition was assessed using Attenuated Total Reflection FTIR (ATR-FTIR). Spectra were obtained by the co-addition of 32 repeated scans performed at a resolution of 2 cm⁻¹ using a Nicolet 380 FTIR spectrometer coupled with a Thermo Scientific SMART iTR ATR sampler equipped with a diamond reflecting cell. To ensure spectra were representative, analyses were performed on at least three sub-samples (each comprising 2–3 pellets) per species. Spectra were then compared against an extensive spectral database for the identification of carbonate phases (see ref. ³⁷).

Finally, each particle morphotype produced by a species was assigned to a carbonate polymorph based on compositional and mineralogical data. The relative abundances of different particle morphotypes were then estimated visually for every sample observed using SEM. These were converted into relative abundances of carbonate polymorphs, which were then averaged for each species.

Carbonate excretion rates

The amount of carbonate excreted by fish in The Bahamas and Australia was quantified using a double-titration approach^18,27,28. Samples were homogenised in 20 mL of distilled water, titrated with HCl to below pH 4.0 and titrated back to the starting pH with NaOH, while continuously aerating with CO₂-free air to remove all HCO₃⁻ and CO₃²⁻ as gaseous CO₂. Concentrations of 0.001–0.1 N were used for both HCl and NaOH as appropriate for sample size. Titrations were performed using a Metrohm Titrando autotitrator and Metrohm Aquatrode pH electrode (Australian samples), or manually with combination pH electrodes (Radiometer PHC 2401) and handheld pH meters (Hanna HI 8314 and Russell RL 200), with acid and base delivery via 2 mL micrometer syringes (Gilmont Instruments, Barrington, USA) with a precision of ± 1 ({{{{{rm{mu }}}}}})L (Bahamian samples). The amount of HCl used minus the amount of NaOH required to return to the starting pH corresponds to the amount of bicarbonate equivalents (i.e., HCO₃⁻ + 2CO₃²⁻) in the sample. Therefore, the molar amount of (Ca, Mg)CO₃ in the sample was calculated as:

assuming that each carbonate molecule yields two bicarbonate equivalents.

Due to laboratory constraints, a slightly different approach was applied to the samples collected in Palau. Specifically, carbonate alkalinity was determined by single end point titration using the mixed indicator Bromocresol Green-Methyl Red⁸⁴. Samples were suspended in 5 mL of distilled water and sonicated, then 50 ({{{{{rm{mu }}}}}})l of mixed indicator was added to the solution which turned blue (pH > 5). Each sample was titrated with 0.01–0.5 N HCl (with continuous aeration with CO₂-free air) until the end point (grey-lavender; pH~4.80) was reached and stable for at least 10 min. If the sample was over-titrated (pink), 0.01–0.1 N NaOH was added to titrate back to the end point and the amount of base used was subtracted from the amount of acid. Acid and base were added using an electronic multi-dispenser pipette (Eppendorf Repeater ®E3X, Eppendorf, Hamburg, Germany) with a precision of  ± 1 ({{{{{rm{mu }}}}}})L. Additionally, the pH of several samples was monitored using a pH microelectrode (Mettler Toledo InLab Micro) to ascertain the correctness of the colorimetric end point. The amount of carbonate in the sample was then calculated using Eq. (1). The method was validated using certified reference material (Alkalinity Standard Solution, 25,000 mg/L as CaCO₃, HACH) and the accuracy in the determination of solid samples was verified using certified CaCO₃ powder (Suprapur, ≥ 99.95% purity, Merck) samples (60–500 ({{{{{rm{mu }}}}}})g) and resulted in 96.53 ± 1.94% accuracy (mean ± SE; n = 8).

To compare values obtained with the two titration methods we further analysed 12 samples collected at Lizard Island, Australia, in February 2016. Samples were collected at 24 h intervals from one individual of Lethrinus atkinsoni (f. Lethrinidae, body mass: 245 g), a group of five Lutjanus fulvus (f. Lutjanidae, mean body mass: 21 g), and an individual of Cephalopholis cyanostigma (f. Serranidae, body mass: 295 g), following the procedures described above. During sample collection water temperature ranged from 29.1 °C during the night to 32.6 °C during the day, with an average of ~31 °C, mean salinity was 35.4, and pH_NBS ranged from 8.13 to 8.21. To compare the amount of carbonate measured by the two methods we added carbonate samples to 20 ml ultrapure water and disaggregated crystals via sonication. We then used a Metrohm Titrando autotitrator and Metrohm Aquatrode pH electrode to measure initial pH of the suspension of carbonates, then titrated each sample of carbonate in two stages. Firstly, they were titrated down to pH 4.80 using 0.1 M HCl, adding 20 µl increments of acid until this was sufficient to keep pH below 4.80 for 10 min whilst bubbling with CO₂-free air. This first stage was comparable to the single end point titration used for samples collected in Palau. Secondly, whilst continuing to bubble with CO₂-free air, further acid was added to the sample until it reached pH 3.89 and was stable for 1 min. Then 0.1 M NaOH was added to the samples to return them to the initial pH. For all samples the first end point titration (to pH 4.80) yielded slightly higher values for carbonate content than the second double titration. The ratio between the two methods (single end point/double titration) was 1.08 ± 0.01 (mean ± SE; range: 1.04–1.14; Supplementary Table 2). As we found a small but consistent difference between the two methods, all following analyses were initially performed on the actual data obtained with the double titration for samples from Australia and The Bahamas, and the single end point titration for samples from Palau. Then, to assess the robustness of the results, we repeated the analyses after applying a correction factor of 1.08 to the excretion rates of Palauan fishes (that used the single end point titration method). All results were consistent and robust to the measured difference between the titration methods (Supplementary Figs. 8, 9).

Finally, measurements of multiple samples from each individual collected over periods of 18–169 h (median: 64 h) were combined to produce an average individual excretion rate in ({{{{{rm{mu }}}}}})mol h⁻¹. For fish held in groups, carbonate excretion rates per individual (of average biomass) were obtained by averaging the total excretion rate of the group across the sampling period and dividing it by the number of individuals in the tank. Excretion rates obtained from fish groups thus evened the intraspecific variability within tanks, and are therefore more robust than those directly obtained from fish held individually. This aspect was considered in our models by fitting weighted regressions (see the “Statistical modelling” section). In total, we measured the carbonate excretion rates of 382 individual fishes arranged in 192 groups (i.e., independent observations), representing 85 species from 35 families across three tropical regions (180 individuals from 29 species in Australia, 90 individuals from 10 species in the Bahamas, and 112 individuals from 46 species in Palau; Supplementary Table 1).

We assume that during the sampling of carbonates fishes were close to their resting metabolic rate and that their carbonate excretion rates are representative of fish at rest. Although the ratio of tank volume to fish volume in our study (median ~660; inter-quartile range ~180–1700) typically greatly exceeds the guideline ideal range for measuring resting metabolic rate (20–50)⁸⁵, fishes were fasted prior to and throughout sampling, and in most instances their movement was somewhat constrained by tank volume. Fasting reduces metabolic rate in all animals, including fish, as they do not undergo energy-intensive digestive processes and use energy reserves to support vital processes, triggering metabolic changes in many tissues and reducing activity levels^86,87. Additionally, other than the carbonate syphoning (< 2 min), no stressors were present. Fishes were not engaged in foraging activity, they experienced no predator-prey interactions, many were held individually so did not engage in social interactions, and social species were held in groups to minimise stress. Therefore, although spontaneous activity likely occurred, fishes were placid throughout the sampling period and the constrained space, minimal disturbance, and fasting suggest that they had very low activity levels.

Selection of families

To assess the main determinants of variability in fish carbonate excretion and composition we considered only families with at least three independent observations (i.e., three individuals or groups of fish). This removed 14 families characterised by only one or two observations, thus reducing our dataset to 175 independent observations from 352 individuals representing 21 families. In our analyses, we considered parrotfishes (f. Labridae: tribe Scarini) separately from other labrids following ref. ²⁸, as preliminary data showed that they excrete distinct carbonate products, possibly due to their distinct trophic ecology and unique gut chemistry⁸⁸.

Explanatory variables

Fish carbonate production is the result of both extrinsic environmental conditions and intrinsic species traits. To analyse the major factors determining carbonate excretion rate and composition we considered a suite of potential variables, while accounting for taxonomic relationships. Some are known to influence fish carbonate excretion rates, such as salinity, body mass, and temperature^8,25,27,63. Others, such as CO₂ and AR, are likely to indirectly affect carbonate excretion by influencing acid–base regulation³¹ and activity level⁵⁵, respectively. We did not have data for seawater pCO₂ during the sampling period, however, due to good aeration of tanks, we assumed that it was close to atmospheric equilibrium (~400 ({{{{{rm{mu }}}}}})atm) at all locations, and thus would not have been a relevant factor in our analysis.

Furthermore, diet is expected to strongly influence carbonate excretion rate and composition, as fish obtain large amounts of calcium and magnesium from their food, which are then likely largely precipitated and excreted as carbonates²³. However, our data were collected on fasting fishes, thus we cannot account for the direct effect of diet on carbonate excretion and composition in our analyses. Nevertheless, we accounted for the indirect effect of diet (i.e., the adaptation of the gut morphology to the typical diets each species has evolved to consume) by including RIL as an additional potential variable in our analyses. Indeed, RIL is strictly related to diet in reef fishes⁶⁰. As we did not have direct measurements of intestinal length of fishes used to collect carbonate samples, we predicted species-level RIL using a Bayesian phylogenetic model trained with the largest available dataset of intestinal length of reef fishes⁸⁹ (see Supplementary Methods).

Statistical modelling

Predictors of total carbonate excretion

Before modelling carbonate excretion rates, we used bivariate correlations to identify potential multicollinearity among all explanatory variables, including two covariates related to our methodologies (total sampling period and titration method). The titration method was strongly correlated with all environmental variables as one protocol was applied to all samples from Palau and the other to all samples from Australia and The Bahamas. Titration method was therefore initially excluded from our models. Conversely, salinity was relatively strongly correlated with RIL (r = −0.60) in our dataset. Therefore, these variables were alternatively included in our models (i.e., the same models were fitted twice including either RIL or salinity as a covariate).

We fitted a series of Bayesian regression models to predict carbonate excretion rates based on the selected traits and environmental variables. Let y_ijk be the carbonate excretion rate of the i^th individual of the j^th species, belonging to the k^th family. We assumed that each observation of the response variable (y_ijk) was t-distributed:

with degrees of freedom (nu), scale (sigma), and observation-specific locations µ_ijk defined as

where (gamma)₀ is the average model intercept, u_k is the random variation in (gamma)₀ based on taxonomic family, and (beta)_x is a vector of regression coefficients of the fixed effects x.

We fitted a series of 36 linear and multilevel models starting from an intercept-only model and increasing in complexity. All models were fitted by weighting the response variable based on whether fish were kept individually or in groups. Although some fish were kept in relatively large groups (up to 13 individuals), most were kept individually (61% of tanks). Therefore, to avoid overweighting observations from groups larger than two individuals, we gave a weight of two to all observations derived from fish kept in groups of two or more individuals. We built linear models by first including body mass which was the known major predictor. We then added either RIL or salinity, which in our exploratory data analysis showed the strongest correlation with the response after accounting for body mass. Lastly we included temperature, AR, and total sampling period, either alone or in combination. This procedure resulted in 18 linear models which were then refitted including taxonomic family as a group-level effect. Model selection was performed through leave-one-out cross-validation (LOO-CV) (Supplementary Table 3). All multilevel models had a better fit than the corresponding linear models highlighting the importance of including the fish family as a group-level effect. Similarly, all models including RIL performed better than the same models where RIL was replaced by salinity. The selected model included the following set of covariates:

where M is the body mass (in kg), RIL and AR are the species-level relative intestinal length and caudal fin aspect ratio (the latter obtained from FishBase⁷⁸), respectively, and T is the average water temperature (in °C) during the sampling period.

To investigate whether there was some unexplained variance in the response that could be attributed to the excluded explanatory variables, we tested for correlations between model residuals and average salinity during the sampling period, titration method, and total sampling period. No residual correlation was observed. However, residual variance was related to the titration method used to quantify carbonate excretion rate (Levene’s test, F_1,173 = 22.82, p < 0.001), with a larger residual variance in samples analysed through single end point titration. Therefore, to account for this, we refitted the selected model as a distributional regression where we allowed the scale parameter (sigma) of the t-distribution to vary with respect to the titration method used:

where ({gamma }_{sigma }) is the intercept and ({beta }_{sigma }) is the regression coefficient for the titration method. Thus, (sigma=exp left({gamma }_{sigma }right)) for observations obtained through the reference method (i.e., double titration), while (sigma=exp left({gamma }_{sigma }+{beta }_{sigma }right)) for observations obtained through single end point titration.

Model comparison through LOO-CV showed that modelling(,sigma) as a function of the titration method improved model fit (Supplementary Table 4). Moreover, the distributional model showed that the parameter (sigma) was actually different between methods (mean and 95% CI: 0.39 [0.31, 0.47] and 0.80 [0.63, 0.98], for double and single end point titration, respectively). This model was therefore selected to draw conclusions on the relationship between carbonate excretion rate and the explanatory variables.

Predictors of carbonate composition

Five major carbonate polymorphs are produced by fish: LMC, aragonite, HMC, MHC, and ACMC, in order of increasing expected solubility. To investigate the factors determining the excretion of the different carbonate minerals by fish we modelled the excretion rates of individual polymorphs. This approach has two major strengths: (1) it facilitates the investigation of both what drives the probability of a polymorph being excreted, and the predictors of the polymorph-specific excretion rates, and (2) it allows direct predictions of polymorph-specific excretion rates.

To obtain the excretion rate of individual polymorphs, the carbonate excretion rate of each fish was multiplied by the species-level relative carbonate composition. Then, we modelled these excretion rates using a Bayesian multivariate hurdle-lognormal model. We used a multivariate model (i.e., a model with multiple response variables) because it accounts for the correlation among polymorphs within taxonomic group, while allowing the use of different sets of predictors for each response. As all response variables contained zeros (14 to 83% of the observations), we opted for a hurdle-lognormal model, which is a two-part model that combines a logistic regression for the probability that the outcome is zero or not, with a lognormal model for the non-zero responses:

where (theta) is the probability of zero outcome (i.e., no excretion), (1 – (theta)) is the probability of positive outcome (i.e., excretion), and logN_CDF is the cumulative distribution function for the lognormal distribution of the non-hurdle part.

The hurdle probability of each carbonate polymorph (({theta }_{{ijk}}^{m})), i.e., the probability that the i^th individual of the j^th species, belonging to the k^th family, did not excrete the m^th polymorph, was estimated using a multilevel logistic regression as:

Where ({gamma }_{0}^{m}) is the average intercept for the m^th polymorph, ({v}_{k}^{m}) is the random variation in ({gamma }_{0}^{m}) based on taxonomic family, and ({beta }_{1:2}^{m}) are the regression coefficients of RIL and T, respectively. While, for each response, we modelled the mean of the lognormal distribution (({mu }_{{ijk}}^{m})) according to Eqs. (3) and (4) and the standard deviation (({sigma }_{i{jk}}^{m})) according to Eq. (5). Body mass and AR were only included as predictors of the excretion rates, but not as predictors of the probability of excretion of different carbonate polymorphs, because no mechanistic link is described or expected for these variables. Conversely, we used fish family, temperature, and RIL to predict the probability of excretion of the polymorphs because fish carbonate mineralogy is generally consistent within families^28,38, a potential thermal effect has been suggested³⁸, and fish with long intestines have long gut residence times^64,65 likely affecting the precipitation of different polymorphs.

To account for both the between-family variance (({tau }^{2})) and covariance ((rho)) we modelled the group-level effects as correlated. Therefore, we assumed the family-specific intercepts of both the hurdle (({v}_{k}^{m})) and non-hurdle part (({u}_{k}^{m})) of the model to follow a multivariate normal distribution with zero means and covariance matrix Σ with 2 m(2 m + 1)/2 components:

with

Finally, we fitted a second model specifying a different formula for the mean of the lognormal distribution of each polymorph (({mu }_{{ijk}}^{m})). This was achieved by removing the fixed effects with relatively large errors (i.e., those with an estimated error greater than the mean estimate). Specifically, we removed the effects of RIL, AR, and temperature on LMC and the effect of temperature on MHC. Model comparison through LOO-CV showed no difference in model fit between the two models, therefore, we selected the more parsimonious model to create figures and interpret results.

All analyses were performed with the software program R (version 4.1.3⁹⁰) and all models were fitted with the R package brms (version 2.15.0⁹¹). Linear and multilevel models were run for 4 chains, each with 4000 iterations and a warm-up of 1000 iterations, whereas hurdle models were run for 3 chains, each with 4000 iterations and a warm-up of 2000 iterations. All models were examined for evidence of convergence using trace plots and Gelman–Rubin statistics and we used posterior predictive distributions to check for models’ fit.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com