Understanding decay in marine calcifiers: Micro‐CT analysis of skeletal structures provides insight into the impacts of a changing climate in marine ecosystems

Calcifying organisms and their exoskeletons support some of the most diverse and economically important ecosystems in our oceans. Under a changing climate, we are beginning to see alterations to the structure and properties of these exoskeletons due to ocean acidification, warming and accelerated rates of bioerosion. Our understanding has grown as a result of using micro‐computed tomography (µCT) but its applications in marine biology have not taken full advantage of the technological development in this methodology. We present a significant advancement in the use of this method to studying decalcification in a marine calcifier. We present a detailed workflow on best practice for µCT image processing and analysis of marine calcifiers, designed using coral skeletons subjected to acute, short‐term microbial bioerosion. This includes estimating subresolution microporosity and describing pore space morphological characteristics of macroporosity, in perforate and imperforate exoskeletons. These metrics are compared between control and bieroded samples, and are correlated with skeletal hardness as measured by nanoindentation. Our results suggest that using subresolution microporosity analysis improves the spatiotemporal resolution of µCT data and can detect changes not seen in macroporosity, in both perforate and imperforate skeletons. In imperforate samples, the mean size and relative number of pores in the macroporous portion of the images changed significantly where total macroporosity did not. The increased number of pores and higher microporosity are both directly related to a physical weakening of the calcareous exoskeletons of imperforate corals only. In perforate corals, increased macroporosity was accompanied by an overall widening of pore spaces though this did not correlate with sample hardness. These novel techniques complement traditional approaches and in combination demonstrate the potential for using µCT scanning to sensitively track the process of decalcification from a structural and morphological perspective. Importantly, these approaches do not necessarily rely on ultra‐high resolution (i.e. single micron) scans and so maintain the accessibility of this technology. The continued optimization of these tools for a variety of marine calcifiers will advance our understanding of the effect of climate change on marine biogenic calcified structures.

A decrease in exoskeletal durability and strength has broad consequences for ecosystem structure and function. For pelagic and benthic calcifiers, alterations to their shell structure will affect their ability to capture light, and protect themselves from predation and UV damage (Enríquez, Méndez Eugenio, Hoegh-Guldberg, & Iglesias-Prieto, 2017;Meng et al., 2018;Monteiro et al., 2016;Orr et al., 2005). Reef-building organisms will suffer a loss of structural stability, compromising key ecosystem services such as commercial fishing and coastal protection (Fantazzini et al., 2015;Graham & Nash, 2013;Meng et al., 2018). Impairment to these critical functions has the potential to affect primary productivity over large areas (e.g. the Caribbean sea; Allgeier, Valdivia, Cox, & Layman, 2016). Therefore, the nature and progression of climate change impacts in these marine systems are directly related to changes in the physical and morphological properties of biogenic CaCO 3 structures.
Despite this recognized link, we still lack a comprehensive knowledge base to predict changing structural properties of calcifiers in the Anthropocene. This is, in part, due to limitations in the tools we use to study decalcification (i.e. the degradation or dissolution of CaCO 3 structures; see Table 1).
However, the application of this method in marine biology stands to benefit from improvements in µCT technology in other disciplines. Recently, we applied advanced µCT methodology to estimate the subresolution porosity (termed microporosity) of coral skeletons subjected to accelerated microbial bioerosion by resident endoliths during a marine heatwave (Leggat et al., 2019). Endolithic microbes are near ubiquitous in biogenic carbonates (Tribollet, 2008;Tribollet & Golubic, 2011) and are highly effective bioeroders (Carreiro-Silva, McClanahan, & Kiene, 2005;Chazottes, Campion-Alsumard, & Peyrot-Clausade, 1995;Perry et al., 2014). Shortterm microbial bioerosion is an ideal process off which we can design sensitive methods for measuring decalcification because the small-scale changes are hard to detect and because reef-wide microbial bioerosion is expected to significantly increase as a result of climate change (Leggat et al., 2019;Reyes-Nivia et al., 2013). Here we present a major advance in the application of tomographic (μCT) analyses, adapting three geological techniques that can be applied to investigate decay and dissolution in calcified structures. We also provide a case study of quantifying imperforate and perforate pore space morphology in two species of scleractinian coral at a high spatiotemporal resolution. Included is a step-by-step protocol ( Figure 1) that provides a guide for tomographic image processing and analyses. of decalcification from a structural and morphological perspective. Importantly, these approaches do not necessarily rely on ultra-high resolution (i.e. single micron) scans and so maintain the accessibility of this technology. The continued optimization of these tools for a variety of marine calcifiers will advance our understanding of the effect of climate change on marine biogenic calcified structures.  Fantazzini et al. (2013Fantazzini et al. ( , 2015 and Anovitz and Cole (2015) X-ray radiography A 2D X-ray radiograph is taken of a cut surface or slice of the sample. Growth bands can be used to measure calcification and in sclerochronology, erosion scars can be used to measure dissolution/ degradation Can be used to measure inter-annual growth patterns and identify historical stress events. Cost-effective for high biological replication Destructive of whole sample. Measurement is dependent upon where the 2D image is taken. Images artefacts can arise from 3D variation in density Barnes, Lough, and Tobin (1989), Le Tissier, Clayton, Brown, and Davis (1994)  (b) washed with deionized water; (c) and dried at a maximum temperature of 50°C for 48 hr to prevent crystalline phase shifts (Caroselli et al., 2011). Skeletons were packed in 50 ml centrifuge tubes with cotton wool to prevent damage during storage and transport. For scanning, two samples were mounted on top of each other inside a 50 ml falcon tube, separated by a disk of clear plastic.

| Scanning (Step 2, Figure 1)
Samples were scanned using a HeliScan MicroCT system with an optimized space-filling trajectory (Kingston et al., 2018) to yield sharp images (Latham, Varslot, & Sheppard, 2008; with an isotropic voxel size of 26.2 µm (  (Barrett & Keat, 2004). To correct for this, the polychromatic X-ray source was passed through a 0.25-mm stainless steel filter. This removes this artefact without sacrificing the improved resolution associated with polychromatic versus monochromatic X-ray sources (Barrett & Keat, 2004).

| Reconstruction (Step 3, Figure 1)
Image reconstruction refers to the process of recombining twodimensional information contained in sinograms into a three-dimensional tomogram ( Table 2). The optimization of reconstruction algorithms requires analytical expertise to ensure a theoretically exact reconstruction; in this case, a Katsevich 1PI inversion formula was used (Kingston, Sakellariou, Varslot, Myers, & Sheppard, 2011;Varslot, Kingston, Latham, et al., 2011;. This was performed in Mango image analysis software.

| Image masking (
Step 4, Figure 1) All tomographic image analyses were performed using WebMango, an online implementation of Mango image analysis software; however, ImageJ could also be used. The software NCViewer was used to examine 2D slices from tomograms in netcdf format; any netcdf viewer can be used in its place.
One challenge in quantifying skeletal morphology in tomograms of irregular objects is defining the region of interest by excluding external space outside of the object (Roche et al., 2010). To achieve this, we used a three-step process. First, we manually 'shrank' the external space to be subsequently masked algorithmically (first panel, Step 4, Figure 1). In practice, this involves choosing a 2D slice of the tomogram on the z-axis and defining the smallest possible circle about the sample that excludes the maximum empty space without cropping the sample.
We then cycle through the rest of the z-slices, growing the circle as necessary to ensure none of the tomogram is cropped. We then applied a flood fill algorithm (FFM) in order to mask the empty external space within our theoretical cylinder (second panel, Step 4, Figure 1).
The algorithm grows from user-defined voxels with a particular intensity value (Table 2), replacing these seeded voxels with a pre-defined replacement intensity. Here, we replaced zero-intensity external voxels with the max value 65,335, the maximum value on the 16-bit scale (Table 2), as no part of our sample had this intensity, and we wanted to partition internal air space (i.e. porosity) from external air space. The process is the same as using the 'fill' function in digital image analysis/ design programs to replace one colour with another across a whole image. The FFM grows along connections between neighbouring voxels with the same intensity; therefore, how you define a connection between two voxels affects the outcome. There are two variations based on whether we consider voxels touching at the corners as connected: eight-way (faces and corners) and four-way (faces only). Here we used eight-way connections and constrained the 'invasion radius' of the FFM to six voxels (157.2 µm) wide. The FFM tests if it can grow at least six voxels in each direction before making any replacements; this prevented it from incorrectly masking pore space that was connected externally. Finally, we shrank the resulting mask by two voxels (third panel, Step 4, Figure 1); this was a trade-off between excluding external space and removing some small features of the skeleton.

| Image segmentation (Step 5, Figures 1 and 2)
Image segmentation, also known as thresholding ( (Silbiger et al., 2016) or using Otsu's thresholding (Otsu, 1979) which can automatically segment images based on maximizing variance between the classes of voxels (e.g. pore and skeleton). Binary objective segmentation is valuable because it overcomes a phenomenon known as 'partial voxel filling', in which voxels in a tomogram are neither wholly solid nor air (Barrett & Keat, 2004;DeCarlo, 2017;Silbiger et al., 2016). This occurs at the boundaries between the two phases or due to subresolution porosity, where pore spaces are smaller than one voxel. It also allows us to convert greyscale values

Sobel gradient
The gradient of a slope representing an increase or decrease in greyscale intensity as you move from one voxel to its neighbours. The bigger the difference between neighbours, the greater the gradient Sobel and Feldman (1968) Isotropic voxel size The length of a cubic 3D pixel. Akin to the image resolution of our tomogram (though resolution and voxel size can be different; see source) du Plessis, Broeckhoven, Guelpa, and le Roux (2017) Greyscale intensity The intensity of X-ray attenuation by a voxel or group of voxels. Related to sample density but using a variable scale from 8-bit (0 to 255) to 16-bit as presented here (0 to 65,335)

Ketcham and
Carlson (2001) Macroporosity The porosity defined as air that is 'visible' at the chosen scan resolution (i.e. resolvable porosity) Leggat et al. (2019) Microporosity Subresolution porosity due to pore spaces that are smaller than the isotropic voxel size but affect X-ray attenuation  (Table 1; Fantazzini et al., 2015). Otsu's thresholding is a straightforward and objective method but its performance depends upon how well separated the peaks in the intensity histogram are and how deep the 'valley' is between them (Kittler & Illingworth, 1985).
Binary thresholding addresses the problem of partial voxel filling at phase boundaries; but partial voxel filling does not always represent an artefact of the method but might also represent variation in subresolution porosity (Figure 2b-d). Therefore, the accuracy of objective thresholding may be context-dependent and it limits inferences about subresolution characteristics of the sample.
A challenge with biological samples is pore size variability relative to voxel size, creating a need to correct for boundary effects while retaining information about subresolution porosity.
Our proposed solution was to apply three-phase segmentation, a technique designed for the geosciences that can be used to simultaneously define and 'grow' phases from multiple materials with variable density in polymineralic samples (Prodanović, Mehmani, & Sheppard, 2015;Sheppard, Arns, et al., 2006;Sheppard, Sok, & Averdunk, 2004). We defined three phases in our monomineralic sample: air, solid and an air-solid mixed F I G U R E 2 Three-phase segmentation of a single Acropora aspera skeleton. (a) Voxels in an image are seeded or undecided based on their greyscale intensity and sobel gradient ( This phase is akin to that used to estimate microbial bioerosion over long time frames (>1 year) by  and Silbiger et al. (2016), which was referred to as variable skeletal density. Three-phase segmentation uses a combination of greyscale intensity and sobel gradients ( Figure 2a; Table 2). Using sobel gradients includes additional information about the average change in density in all directions of 3D space, which indicates a possible phase boundary (i.e. high gradient), to distinguish between true microporosity and partially filled voxels at the boundary between air and solid. First, we set the seeded voxels and the 'undecided' voxels ( Figure 2a,b). The seeded voxels were then used to grow each of the three phases into the 'undecided' fractions, using a converging active contours (CAC) algorithm (Sheppard et al., 2004). The speed at which this algorithm grows each phase is defined by the sobel gradient value.
Near a true phase boundary, the algorithm will slow down due to a higher intensity difference between adjacent voxels. This generates a sharp, well-defined image.
Segmentation steps: 1. The operator defined two intensity thresholds, T1 and T4. All voxels below T1 are assigned as air; all voxels above T4 are assigned as solid (Figure 2a 5. Once these thresholds are set, then the CAC algorithm grows all three phases simultaneously into the undecided fractions of the tomograms (Figure 2a,b).
The operator needs to examine the segmented image and compare it in detail to the original greyscale tomogram. In order to optimize the segmentation process, we needed to go through multiple iterations of steps 1-5 during which thresholds are adjusted one at a time. As a result, three-phase segmentation is a time-consuming process. It is essential that the same bias (underor overestimation) be applied to all the samples. It is useful for final segmentations to be checked against the original tomogram by a third party, with a focus on phase boundaries, in order to identify potentially consistent operator bias. At the end of this process, we could quantify the macroporosity (i.e. resolvable porosity) of our samples as the volume of voxels in the air phase ( Figure 1). While this method of segmentation uses thresholds that are initially defined by the user which can introduce subjectivity, it also allows us to make inferences about subresolution characteristics of the sample.

| Microporosity analysis (Step 6, Figures 1 and 3)
We used microporosity analysis to estimate the total volume of air contained within all of subresolution pore spaces in the intermediate phase (Sheppard, Arns, et al., 2006;Varslot, Kingston, Myers, et al., 2011). This can then be summed with macroporosity to obtain a more precise value for skeletal porosity. Microporosity  (Sheppard, Arns, et al., 2006).
Here that is a reasonable assumption because our samples are, for practical purposes, monomineralic. It is well known that calcite is present in small amounts in coral skeletons but its absolute density (2.71 g/cm 3 ) is close to the dominant CaCO 3 isomorph, aragonite (2.93 g/cm 3 ), and not distinguishable through µCT scanning at this resolution.
Microporosity steps: 1. The operator defined two thresholds, T m 1 and T m 2, for microporosity based upon the thresholds used in, and results from, the preceding three-phase segmentation (T1 < T m 1 < T2 and T3 < T m 2 < T4; Figure 3a).
2. To visualize the range of greyscale intensities captured by T m 1 and T m 2, we set the voxels below T m 1 and above T m 2 to a neutral colour such as black ( Figure 3b) and compared this to the original and segmented tomograms. We identified originally 'uncertain' areas, noted how the CAC algorithm assigned these voxels and compared whether these agreed with our microporosity thresholds.
3. We then categorized the greyscale values captured by the microporosity thresholds into 100 bins. Based upon the linear relationship between greyscale intensity and microporosity, these 2.1.7 | Image transformations (Steps 7 and 8, Figure 1) The analyses in steps 9 and 10 are approaches to describing the morphology of the resolvable pore space. Prior to these, we necessarily merged the solid and intermediate phases to produce a twophase tomogram (Figures 4b and 5c). At this stage you can also calculate the internal surface area of the resolved void space using the methods described by Ohser and Mücklich (2000). Following this, we applied a Euclidean distance transformation (Step 8, Figure 1) to voxels in the air phase. Each voxel is revalued based upon the straight-line distance to the nearest phase boundary in 3D space which represents its position within the pore space. The resultant 'distance maps' appear as in Figures 4c and 5d. 2.1.8 | Particle size distribution analysis for imperforate corals (Step 9, Figures 1 and 4) TD-NMR has been previously used to describe the size-frequency distribution curves of pores within coral skeletons in three dimensions (Table 1; Fantazzini et al., 2013Fantazzini et al., , 2015. But it only describes the shape of, and shifts, in distribution (Lawrence & Jiang, 2017), as opposed to quantifying the volumetric size of individual pores. Here we applied a common technique from the geosciences, particle size distribution (PSD) analysis (Lawrence & Jiang, 2017 Working upon the transformed tomogram, our steps for applying PSD to biogenic carbonates were as follows: 1.  (2000), we also calculated pore-wise surface area to volume ratio as a proxy for pore shape. Finally, by normalizing the number of pores to the total volume of the sample, we calculated pore density (i.e. no. of pores/mm 3 ).

| Maximal inscribed sphere analysis (Step 10, Figures 1 and 5)
We described the pore geometry of complex perforate corals through the maximum inscribed sphere technique (Fanwua, Chunguang, Haiming, & Juan, 2011;Sheppard, Arns, et al., 2006;. This approach assigns a value to each pore voxel based on the radius of the largest possible theoretical sphere that can be formed around it while still fully contained within the pore space (Figure 5a). In isolation, the greatest sphere size for each voxel has a radius equal to its Euclidean distance value (Figure 5a). Each voxel within this 'inscribed' sphere is labelled with this value. When applied to all the voxels in our sample simultaneously, this produces a series of overlapping spheres.
The problem then is that voxels can take multiple values; either that of the sphere formed about them, or of the spheres formed by neighbouring voxels that encompass them. This is resolved by a size-dependent hierarchy, wherein voxels that exist in an overlap take the value of the larger sphere (Figure 5b), hence the term 'maximum inscribed sphere' (MIS). It allows us to describe both maximum and minimum geometric features of pore space. The output is in the form of a relative proportion of all voxels that are in each sphere size class, expressed as binned radii measured in voxels ( Figure 8).

| Statistical analyses
All statistical analyses were conducted using r 3.6.0 (R Core Team, 2019). The statistical tests and packages used were based upon the formats and distributions of the data output by the above analyses (Table 3; Table S1). Almost certainly the data generated will vary by taxa, species, resolution and possibly by the

| Beta regression
This type of regression is designed to examine data in the form of relative proportions (or percentages), bounded between 0 and 1 (0% and 100%; Table S5; Cribari-Neto & Zeileis, 2009;Ferrari & Cribari-Neto, 2004). It was the most common type of regression used here (Table 3). To perform this, we used the betareg package (Cribari-Neto & Zeileis, 2009). This model cannot take zero values which were present in low numbers in our microporosity phase and MIS data. Therefore we applied the following transformation recommended by Cribari-Neto and Zeileis (2009), where n is the number of samples in your analysis:

| Regression with gamma distribution
The PSD and pore 'shape' data were non-negative and extremely skewed, suggesting that a gamma distribution would be the optimum choice for modelling these data (Table 3; Table S5; Figures S1 and S2). We used the fitdistrplus (Delignette-Muller & Dutang, 2015) and DescTools (Signorell, 2016) to compare the goodness-of-fit for gamma, Weibull and log-normal models ( Figure S2). In particular, the Weibull distribution was a good candidate given its historical application to PSD data (Vesilind, 1980). However, here we judged Gamma regression to be a similar, if not superior, fit to the data and easier to interpret.
We use the shape (α) and rate (β) parameters to describe distributions. α relates to the skewness of the data (skewness = 2/√α) while β describes the 'rate of decay', essentially the gradient of the slope between the peak and the tail. Low values of α correspond to high skewness; high values of β indicate a narrow peak. A total of 174,679 pores were analysed for size and shape.
To avoid pseudoreplication and account for the unequal samplewise contribution to model variance, we included sample ID as a random effect.

| Correlating with hardness
To elucidate the relationship between these metrics and mechanical properties, we used the data on sample hardness (n = 5 per species per treatment) from TP 1 , as measured by nanoindentation in Leggat et al. (2019). Total sample hardness is represented on a scale between 0 and 100 and is an average of three locations in the skeleton: outer surface, inner skeleton and internal pore surface. We used multiple imputation by chained equations (Buuren & Groothuis-Oudshoorn, 2010) to account for the unequal sample sizes, generating 20 imputed datasets for each component of hardness. We then ran beta regression models on each dataset and pooled results. For P. damicornis, we included macroporosity, microporosity, median pore size and pore density as covariates. Pore shape was excluded due to its high variance inflation factor (>10) and strong correlation with median pore size (r 2 = 0.92). For A. aspera, we included the sphere size class mode and median as well as macroporosity and microporosity.

| RE SULTS
Here we provide the results of our case study, in which we applied this µCT workflow to two scleractinian corals, P. damicornis and A. aspera. We compare the analysis results of heat-stress impacted, eroded corals to healthy, control corals. We used imperforate and perforate coral growth forms to demonstrate the variations in applying and interpreting μCT analysis depending on exoskeletal structure and also exposure to microbial bioeroders.
Outliers within each analysis were identified by their Cook's distance relative to a cut-off of 4/n (0.222; Table S5; Cook, 1977) and the model outcomes with and without these outliers are presented.

| Microporous phase distribution
At time point 0 (TP 0 ), only the 83.5% microporous bin had significantly higher relative proportion for heat-exposed corals compared to controls (Figure 6a; z = −2.055, p = 0.040). In contrast, control and treatment groups had markedly different distributions at time point 1 (TP 1 ; Figure 6b). Voxels between 71.5% and 83.5% microporous in the treatment group were proportionally greater compared to the control (p < 0.05; Table S2). Concurrently, there was a significant (p < 0.05) decrease in the relative proportion of voxels in the range 32.5%-18.5% microporous for the heat-exposed corals compared to control corals.

| Total surface area to volume ratio
Graphically, it appeared that the ordinary least squares (OLS) model violated the assumption of homoscedasticity with respect to time point. This is likely due to different sample sizes for TP 0 and TP 1 (n = 6 vs. n = 12). We applied a Breusch-Pagan test (Breusch & Pagan, 1979) for heteroscedasticity which was nonsignificant (χ 2 = 1.183; p = 0.277). At neither TP 0 nor TP 1 were the two groups significantly different (TP 0 : t = 2.097, p = 0.055; TP 1 : t: −0.401, p = 0.695). However excluding an outlier in the control group at TP 1 (Cook's distance = 0.375; surface area to y = (n − 1) + 0.5 n .

| Pore density
A potential deviation from normality arose due to two outliers, both from TP 1 : one in the treatment group (Cook's distance = 0.303; density = 10.300 pores/mm 3 ) and the other in the control group (Cook's distance = 0.285; density = 8.595 pores/ mm 3 ). When included in the OLS model, there were no significant differences between control and heat exposed at either time point

| Pore size-frequency distribution
The number of pores per sample ranged from 3,027 to 20,890. The size-frequency distribution curve of these was extremely right skewed. In all cases, the shape parameter α was <0.3 (Figure 7c,d) which means the curve was effectively asymptotic to the y-axis.

| Pore-wise surface area to volume ratio
Graphically, it appeared that high volume pores in TP 1 treatment corals had a lower ratio compared to controls and better approximated the ratio of a sphere (Figure 7e,f). However, the mixed effects gamma regression revealed this apparent change to be non-significant.

F I G U R E 6
The microporosity phase distributions for Pocillopora damicornis (a and b) and Acropora aspera (c and d), at TP 0 (a and c) and TP 1 (b and d). Green curves represent treatment/eroded corals; black curves are controls. In c and b, *denotes a significant (p < 0.05) divergence in distribution. Bin-wise values and contrasts are in Tables S2 and S3

| Hardness versus porosity
Pore density was significantly negatively related to inner skeletal hardness in all of our samples (z = −3.048, p < 0.05). Outer skeletal hardness was not significantly related to any of the porosity metrics used here. Finally, the hardness of the inner surface of pores was negatively correlated with sample microporosity (z = −3.516, p < 0.05) with an estimated 3.8% decrease in hardness per % increase in microporosity.

| Microporous phase distribution
At TP 0 , there were considerable differences between the microporosity phase distributions of control and heat-exposed corals ( Figure 6c). Forty-four of the 100 microporosity bins were significantly different between the two groups ( Table S3). Voxels that were 93.5%-89.5% or 37.5%-21.5% were proportionally more common in control compared to heat-exposed corals (p < 0.05). At TP 1 , these differences were no longer evident (Figure 6d). Voxels between 94.5% and 79.5% were of greater relative proportion in control compared to heat-exposed corals (p < 0.05).

| Total surface area to volume ratio
The assumptions of the OLS model were confirmed graphically and results indicated that there was no significant difference in surface area to volume ratio at TP 0 (t = 0.312, p = 0.7595). But at TP 1 , the surface area to volume ratio in treatment corals was almost double than that in control corals (12.8 vs. 6.62; t = −2.284, p = 0.039).

| MIS distribution
The distribution of sphere radii was, as with PSD, right skewed though to a lesser extent (Figure 8). At TP 0 , there were no significant differences in the relative proportions of each bin between heat exposed and control (Table S4; p > 0.05). In contrast, seven of the bins were significantly different at TP 1 , following microbial bioerosion (p < 0.05; Leggat et al., 2019). The relative proportions of spheres with radii 1.25 and 4.75-6.25 voxels decreased in treatment groups while the relative proportions in bins 2.75 and 3.25 increased. This latter bin in particular accounted for 17.13% of all air voxels compared to 12.84% to control corals (Table S4).

| Hardness versus porosity
We found no significant correlations between measures of skeletal hardness, porosity and pore space morphology, likely reflecting the lack of significant changes in hardness for this species (Leggat et al., 2019).

F I G U R E 7
Results from the particle size distribution analysis (black represents the control group in all cases; red/orange represents the treatment group). In all cases, *denotes a significant difference (p < 0.05) between treatment groups. (a and b) Changes in mean pore density (pores/mm 3 ), including (a) or excluding (b) two outlier (one in each treatment group). (c and d) Mean volumes of individual pores in control and treatment groups at TP 0 (c) and TP 1 (d). α and β refer to the shape and rate parameters, respectively, of the gamma distributions for each group. Low values of α correspond to highly skewed data; high values of β correspond to a steep slope between the distribution peak and tail. (e and f) Porewise surface area to volume relationships for control and treatment groups at time points zero (e) and one (f). The blue dotted line is the surface area to volume ratio of a sphere

| Interpretation, P. damicornis
Using a single density threshold to distinguish between air and solid overcomes artefacts such as partial voxel filling but invariably ignores information contained in the tomogram by partitioning and assigning voxels with variable density to one phase or the other (Roche et al., 2010;Scherf & Tilgner, 2009). By using three-phase segmentation, we make an attempt to quantify this information by considering this range of densities as a distinct phase we call microporosity. In this case study, this is particularly useful as we are studying biological behaviour occurring below scan resolution. The results of P. damicornis highlight this. The macroporosity of heattreated corals showed no significant increase in an 8-week period while the microporosity of our samples doubled. Interestingly, this overall increase was accompanied by a shift in the relative distribution of voxels towards a higher voxel-wise percentage of microporosity ( Figure 6). We find evidence of this crossover when we use PSD analysis. The decrease in mean pore size (Figure 7c,d) in eroded samples is counterintuitive to that we might expect but our working hypothesis is that it can be explained by an increased number of small pores that were previously too small to appear at our resolution.
There are no significant changes in pore shape despite apparent shift towards a relationship characteristic of a sphere. As this change seemed most evident in larger pores, we re-ran the analysis on pores ≥0.5 mm 3 but this did not affect the results. It is possible that pores undergo alternating roughening and smoothing in relation to our scan resolution. Subtle features first eroded, leading to an apparent smoothing of pore surfaces; this may also explain the non-significant difference in whole sample surface area to volume ratio at TP 1 compared to TP 0 , in that model which excluded outliers. The 'smoother' surface might then be subsequently roughened; this is supported by the increase in surface area to volume ratio in A. aspera whose macroporosity increased more than that of P. damicornis. Different measures of pore shape should be explored in the future as this is a known determinant of the mechanical properties, including hardness, of manufactured materials (Speirs, Humbeeck, Schrooten, Luyten, & Kruth, 2013;Torres-Sanchez, McLaughlin, & Bonallo, 2018).
We found that the hardness of internal pore surfaces was significantly correlated with microporosity but not macroporosity, highlighting the potential for this measure to predict changes in exoskeletal mechanical properties where macroporosity cannot. This is consistent with the findings of both geological (Torres-Sanchez et al., 2018) and medical (Wu, Adeeb, & Doschak, 2015;Zhang, Fan, Dunne, & Li, 2018) studies. In nanoporous materials, Esqué-de Los Ojos, Pellicer, and Sort (2016) show that for a fixed degree of porosity, material deformation by nanoindentation was greater when the material was composed of many small pores than those with fewer big pores. The inner pore surface is effectively the reaction substrate for CaCO 3 dissolution given that microbial bioerosion is generally accepted to occur due to decreases in pore water pH resulting from microbial respiration (Garcia-Pichel, 2006;Tribollet, 2008). Further we show that increased pore density, concurrent with decreased mean pore size, correlated with decreased hardness inside the coral exoskeletons.
These results suggest that pore density and microporosity could be useful avenues for better predicting a loss in skeletal hardness compared to macroporosity.

| Interpretation, A. aspera
In addition to the increase in macroporosity reported by Leggat et al. (2019), we found a significant change in sample microporosity in our extended analysis. Unexpectedly, bioeroded corals had lower F I G U R E 8 Results from maximum inscribed sphere (MIS) analysis. Shown for TP 1 only as there were no significant differences at TP 0 . It is the distribution of voxels across different sphere radii, normalised to a relative proportion for each size class. *Denotes a statistically significant (p < 0.05) difference between groups in a particular size class microporosity. One possible explanation is that the increase in microporosity, which presumably precedes increased macroporosity in P. damicornis, happened earlier in the experimental period in A. aspera.
We are possibly seeing a later stage of the same process, hence the flattening of the distribution for A. aspera at TP 1 (Figure 6d). This has interesting implications for combining micro-and macroporosity to measure fine-scale rates of carbonate dissolution in different species/organisms. The overall decrease in microporosity also indicates that voxels in this phase were not 'replaced' by the erosion of solid carbonates (i.e. the emergence of a new peak at the solid-intermediate threshold). We suggest this is due to a nonlinear relationship between the rates of microbioerosion of solid and microporous carbonates. Microbioerosion is a chemical process driven by the acidification of the microbes' external environment (Tribollet & Golubic, 2011). Comparing solid and microporous carbonates, the latter has a greater surface area for carbonate dissolution, suggesting that the rate of dissolution increases exponentially with microporosity. But this is potentially limited to very small pore spaces due to higher water volumes in larger pores, given that we did not see an increase in macroporosity in P. damicornis despite higher microporosity.

| Limitations
A potential confounding factor is the presence of a scaffolding proteinaceous organic matrix (OM) in calcified exoskeletons (Weiner, Traub, & Lowenstam, 1983). Soaking in bleach is intended to remove OM from the skeleton, but the success of this depends on how well the bleach saturates the skeleton. For the same reason that buoyant weight and TD-NMR underestimate porosity (Table 3), it is likely that some intra-crystalline OM remains within the skeleton. This will attenuate low energy X-rays and so may appear as very low-density areas of skeleton (i.e. high microporosity). Fortunately, OM tends to constitute <1% of the total skeletal weight (Caroselli et al., 2011). Nonetheless, it remains a possible source of variation between organisms and/or species.
We find that the microporous phase distribution does not precisely match the shape of the original intensity histogram for the same sample (Figure 3). Specifically, in the original intensity histogram we see the curve rising towards the air/solid peaks while the microporosity distribution falls to zero. This is due to the 'uncertain' areas assigned by the CAC algorithm. Microporosity analysis ignores uncertain voxels that were assigned as solid or air by the CAC algorithm and only displays those assigned as 'intermediate'. As the greyscale intensity approaches either peak, the algorithm is more likely to assign an uncertain voxel as air or solid. Thus the curve approaches zero at the tail ends of the microporosity distribution.
In some cases, the location of the T2 and T3 thresholds from the segmentation is visible in the microporosity distribution as sudden jumps in the curve ( Figure 3b); however, in Figure 6 they are less clear due to averaging across samples. Given that these areas were previously uncertain, we consider the bin-wise distribution here to retain some of that uncertainty.

| Further applications and future directions
The biomechanical strength of carbonate structures built by marine calcifiers underpins the stability of the physical habitats these organisms build (Byrne & Fitzer, 2019), and directly affects population scale dynamics such as survivorship and productivity (Davidson et al., 2018;Orr et al., 2005). It is therefore vital that we are able to precisely quantify how environmental change in the Anthropocene is affecting the strength and decay rate of calcified structures. µCT scanning is a powerful tool that has so far been used to advance our understanding of the micromorphology of calcified exoskeletons, but its full potential in ecological research remains untapped (Gutiérrez et al., 2018). It has been employed for very high-resolution (i.e. single microns) studies of a range of calcifying taxa including pteropods, gastropods, rhodophytes, bivalves and annelids, quantifying density, shell thickness and volume in the context of biomechanical change (Chatzinikolaou, Grigoriou, Keklikoglou, Faulwetter, & Papageorgiou, 2016;Li et al., 2014;Meng et al., 2018;Oakes & Sessa, 2019). But a lack of advanced analytical techniques to complement this means that certain features of carbonate structures, such as pore space morphology, are still predominantly measured using 2D microscopy such as SEM (Byrne & Fitzer, 2019;Li et al., 2014;Meng et al., 2018).
The results of published studies which have used µCT scanning strongly suggest that the techniques presented in this study are generally applicable to studying the internal structure of other marine calcifiers. For example, a 2D cross-section of the shortspined sea urchin Heliocidaris erythrogramma from a 3D tomogram displays a highly connected and complex pore space (Johnson et al., 2020), similar to that of Porites cylindrica (Video S1), which would make it amenable to MIS analysis. The SEM images of a study. This suggests that microporosity analysis might be a viable alternative to ultra-high resolution SEM, at least in identifying losses in biomechanical strength. Further developing these methods signifies that µCT could replace SEM as the foremost method for studying biomechanical change, which would reduce the cost of analysis compared to using both methods to study 2D and 3D changes respectively.
The methods described here could also be used to study the nature of biotic interactions between ubiquitous endolithic microbes and their calcified hosts. These microbes play an important role in primary production and biogeochemical cycling (Mason et al., 2009;Tribollet, Langdon, Golubic, & Atkinson, 2006), including the exchange of nitrogen and carbon with their host (Barille et al., 2017;Pernice et al., 2019). Different skeletomorphs will impose restrictions on the biomass and composition of this community, through limitations in physical space and/or hydrodynamic diffusion of pore water; this variation can be modelled using these µCT analyses. For example, the permeability of abiogenic carbonates has been successfully quantified in carbon sequestration research and involves the application of microporosity and pore size analysis (Shah, Yang, Crawshaw, Gharbi, & Boek, 2013;Sheppard, Arns, et al., 2006). Understanding the relationship between skeletal structure and the endolithic microbiome could improve our knowledge of the ecological consequences of evolutionary divergence in calcification patterns.

| CON CLUS IONS
Our ever-improving ability to effectively measure change in marine ecosystems is an important driver of our growing understanding of the natural world. µCT is a tool that is growing in its power and popularity in the medical and geological sciences and has the potential to provide new insights through its application to new areas of research and new ecosystems. Here we have described in detail three analytical techniques that are common to the geosciences and adapted them for the first time to studying the impact of climate change upon marine calcifiers. The successful transfer of these methods across disciplines suggests that these tools, as well as others, can be applied to the study of marine calcifiers in general due to commonalities in the features of all carbonate structures (Prodanović et al., 2015;Shah et al., 2013). By relating these measures to the mechanical properties of our samples, we highlight their potential for predicting change in the functional ecology of marine calcifiers.
These tools can advance our understanding of how the physical structures built by marine calcifiers are changing in the Anthropocene as oceans become warmer and more acidic. Combined with further investigations that link the morphological and mechanical properties of these structures, we can begin to model shifts in their ecological function and how this will affect the ecosystems they reside.
Interdisciplinary approaches such as these allow ecological researchers to benefit from the technological advancement in a seemingly unrelated field.

ACK N OWLED G EM ENTS
We thank two anonymous reviewers and Dr Thomas M. DeCarlo for their comments that greatly improved the quality, clarity and context of the manuscript. This collaborative project was funded by an Australian Research Council Discovery Project Grant (DP180103199) and International Society for Reef Studies Graduate Fellowship. All sample collection was approved by the Great Barrier Reef Marine Park Authority (Permit # G17/39039.1). The authors declare no conflict of interest.

AUTH O R S ' CO NTR I B UTI O N S
A.J.F., T.D.A. and W.L. conceived the study; L.B. and M.T. optimized and performed the µCT scanning; A.J.F. and L.K. designed and implemented the image analysis workflow; A.J.F. performed the statistical analyses; A.J.F., L.K., T.D.A. and W.L. wrote and finalized the manuscript.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/2041-210X.13439.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and associated r code (Fordyce, 2020)