High-Throughput Micro-CT and Automated Segmentation for Eco-Evolutionary Studies: A Case Study from Threespine Stickleback

The (anonymized) version v1.0 of this manuscript was submitted to Methods in Ecology and Evolution.

This manuscript (permalink) was generated from habi/sticklebacks-manuscript@9dbbba3 on October 30, 2026 and contains 6468 words. Changes to the submitted version are visualized here on GitHub.

Authors

• David Haberthür
  0000-0003-3388-9187 · habi · @habi@mastodon.social
  microCT research group, Institute of Anatomy, University of Bern, Baltzerstrasse 2, 3012 Bern, Switzerland

• R. Benjamin Sulser
  0000-0002-8750-0942 · sulserrb
  Division of Evolutionary Ecology, Institute of Ecology and Evolution, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland · Funded by Burgergemeinde Bern

• Sheila Christen
  0009-0008-3514-0800 · SheilaChristen
  Division of Evolutionary Ecology, Institute of Ecology and Evolution, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland · Funded by Burgergemeinde Bern

• Catherine L. Peichel
  0000-0002-7731-8944 · cpeichel
  Division of Evolutionary Ecology, Institute of Ecology and Evolution, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland · Funded by Swiss National Science Foundation (TMAG-3_209309/1)

• Ruslan Hlushchuk ^✉
  0000-0002-6722-8996 · RuslanHlushchuk
  microCT research group, Institute of Anatomy, University of Bern, Baltzerstrasse 2, 3012 Bern, Switzerland

Abstract

The threespine stickleback is a well-recognized model system for understanding adaptation to divergent habitats. The populations of benthic stickleback (benthos-feeding) and limnetic stickleback (water-column-feeding) differ in phenotypic traits associated with dietary specialization. Modern investigations of evolutionary change in this organism often require the analysis of hundreds, if not thousands, of individuals. However, analysis of feeding structures—particularly jaws and complex internal branchial anatomy—remains time-consuming and requires considerable expertise, with destructive sampling and fine dissection skills needed for quantitative analysis.

The advent of X-ray microtomography and three-dimensional scanning technologies enables nondestructive imaging and greatly increases the availability of high-resolution morphological data. However, the resulting increase in data complexity and processing demands often makes these approaches impractical for large-scale studies, especially at population scales that require measuring and analyzing hundreds to thousands of samples.

To address these challenges, we developed a rapid and semi-automated segmentation and analysis pipeline based on the Jupyter interactive development environment and the Biomedisa image segmentation platform to investigate three-dimensional morphological adaptation in the threespine stickleback (Gasterosteus aculeatus).

The pipeline includes separation of multi-specimen scans into individual regions of interest, reconstruction of targeted anatomical structures, and morphometric analyses. We applied this pipeline to 38 multi-specimen scans comprising 215 specimens from divergent benthic and limnetic stickleback populations, demonstrating the potential of a high-throughput micro-CT imaging approach for testing ecological and evolutionary hypotheses. These replicable, clade-agnostic methods provide a scalable pipeline that researchers working across diverse systems can adapt for high-throughput micro-CT sampling at scale.

Introduction

The threespine stickleback (Gasterosteus aculeatus) is a well-established model system for understanding the independent evolution of similar traits in similar environments (Bell & Foster, 1994; Reid et al., 2021). This species exhibits marked differences in marine–freshwater, lake–stream, and benthic–limnetic ecotypes (Hendry et al., 2009; Schluter & McPhail, 1992). This study focused on the benthic–limnetic axis, which resolves distinct stickleback ecotypes based primarily on dietary ecology (benthic stickleback forage on macroinvertebrates at the bottom of lakes; limnetic stickleback forage on zooplankton within the water column (WILLACKER et al., 2010)). The samples for this study were collected from the benthic and limnetic source populations used to establish a long-term evolutionary experiment within the Kenai Peninsula of Alaska (USA) (Hendry et al., 2024). This project, the Forward In Time Natural Experimental Study of Selection (FITNESS), aims to study the predictability and repeatability of evolution. Two pools of sticklebacks, one comprising four source populations of limnetic sticklebacks and the other comprising four source populations of benthic sticklebacks, have been placed in eight destination lakes, four of which are small and benthic and four of which are large and limnetic. These new populations have been sampled every year to track the genotypic and phenotypic trajectories of these introduced populations. Understanding the initial variation in the source populations is essential for this project, as it provides the raw material for the evolution in the destination lakes.

Among sticklebacks and other bony fish, differences in jaw structures are directly related to functional and kinematic differences between ecotypes (Haines et al., 2020). Benthic sticklebacks have modified jaws for enhanced suction force and hypertrophied epaxial muscles to aid in foraging on benthic invertebrates, while limnetic sticklebacks have modifications for larger jaw protrusions and quick strikes during ram feeding (McGee et al., 2013). The internal hyoid arch–branchial arch complex is an important structure involved in diet and feeding ecology (BERNER et al., 2008; Schluter & McPhail, 1992). Although the shape and arrangement of the paired ceratobranchial and pharyngobranchial bones within this complex aid in food processing and water vortex generation during feeding (Brooks et al., 2018), the shape of these bones has received relatively little attention in relation to other aspects of dietary anatomy. This is likely due to the flattening and destructive sampling used in traditional raker counting methods, which dissect and deform these structures to make them visible for manual measurements (Ellis & Miller, 2016). Moreover, these structures are difficult to study without full cranial dissection and the corresponding distortion of the branchial anatomy. X-ray microtomography-based three-dimensional (3D) imaging preserves these features at a high resolution and allows for nondestructive analysis. However, imaging large numbers of stickleback specimens has until now been prohibitive due to the limitations of acquisition and processing time.

X-Ray Microtomography

High-resolution X-ray tomography (micro-CT) imaging is an indispensable tool to gain nondestructive insights into the inner structure of highly diverse samples, specifically for specimens in biomedical research (Rawson et al., 2020). Micro-CT imaging is ideally suited for nondestructive anatomical assessments in fish (Ford et al., 2023), capturing internal structures and delicate skeletal elements that are otherwise difficult to quantify without destructive preparation.

While in theory, each specimen could be scanned individually, performing hundreds of scans of individual specimens would represent a large drain of time and resources. Instead many investigators and tomography institutes choose to scan several specimens at once in a “multiscan” format. Although efficient in terms of beam time and number of individuals per scan, multiscans result in considerably larger datasets that must be partitioned correctly with respect to each individual specimen prior to downstream analysis. This creates a bottleneck for large-scale research where sizable multiscans bloat processing time, in addition to potential risks regarding errors with specimen splitting and specimen identification.

Furthermore, once partitioned, small structures can be manually rendered by a skilled investigator using segmentation tools such as 3D Slicer (Kikinis et al., 2013; Rolfe et al., 2021), yet the time and resources required per sample make this approach, too, infeasible for population-level and evolutionary studies. Accessing these structures instead typically requires destruction of the mandibular and cranial anatomy, preventing subsequent analyses of these features.

This project aims to address these gaps, demonstrating a novel pipeline for automatic splitting and rendering of multi-specimen scans for mass sampling, creating a dataset with consistent parameters that can be used as input for downstream machine learning approaches (Lösel et al., 2020) to aid in the segmentation of individual bony structures in each scan. Once a Biomedisa model is trained, the entire pipeline runs from multi-specimen input to reconstructed anatomical structures for each specimen in a fraction of the time and resources required for traditional analysis.

Materials & Methods

Sample Procurement and Preparation

The samples used for this study were collected from six of the eight source lakes used for the FITNESS project in the region of Cook Inlet, Alaska. Fish were collected using unbaited minnow traps in two separate field seasons, the first from May 26–June 10, 2023, and the second from May 25–June 11, 2024. Specimens were collected from a random sample of up to 30 fish from each lake, under Alaska Department of Fish and Game (ADFG) permits SF2023-030 and P-24-015 for 2023 and 2024, respectively. The fish were euthanized with MS-222, photographed, labeled, and preserved in 10% formalin in individual bags, all under McGill University Animal Use Protocol (AUP) MCGL-8265. At the end of each field season, samples were shipped from Anchorage (AK, USA) to Bern (BE, CH), where they were stored until scanning. The total number of fish for each lake was: 52 from Watson Lake (benthic), 48 from Finger Lake (benthic), 46 from Spirit Lake (limnetic), 35 from South Rolly Lake (limnetic), 20 from Walby Lake (benthic), and 14 from Tern Lake (benthic) for a total of 215 specimens. As there is evidence that sexual dimorphism may be present in the stickleback feeding structures (McGee & Wainwright, 2013), caudal fin clips were additionally collected from fish sampled during the 2024 field season for molecular sex identification through genotyping. DNA was extracted from fin tissue following Meeker et al. (Meeker et al., 2007), and subsequent PCR genotyping with LRRc61 primers was performed, following Archambeault et al. (Archambeault et al., 2020).

Due to their inherent contrast difference with the surrounding tissue, the structures of interest in this study (teeth and bones, i.e., jaws and skull) are well visualized in unstained samples, so no further preparation of the fish was necessary.

Micro-CT Imaging

In a pilot study, we identified optimal scanning parameters to meet the constraints on total scanning time, resolution, and sample handling.

To optimize for these constraints, we scanned all sticklebacks in batches of six fish in a custom-made 3D-printed sample holder per scan. This holder was designed in OpenSCAD (OpenSCAD Developers; version 2021.01) and is available online, either directly as an STL file for printing or as a (parameterized) OpenSCAD file for adaptation to other classes of samples. Both files are part of a library of 3D-printable sample holders for tomographic imaging (Haberthür, 2019).

Fish were taken out of single, labeled Falcon tubes, wrapped in melamine resin foam (virtually radiotransparent, but stable) and put head-first into the labeled compartments of the sample holder. The foam was moistened with formalin and the top of the sample holder was sealed with Parafilm to avoid drying of the fish during the tomographic imaging.

Tomographic imaging was performed on a Bruker SkyScan 2214 (Bruker microCT, Kontich, Belgium) at the Institute of Anatomy, University of Bern, Switzerland. In total, we performed 38 scans, each scan usually containing up to six fish in the sample holder.

The relevant details of each scan are summarized in a table in the Supplementary Materials; a short overview of the scanning parameters is given below. The X-ray source was set to a voltage of 60 kV and a current of around 110 µA for all scans except one, where we used a source voltage of 49 kV and 159 µA due to operator error. Setting different voltage and current for one scan had no discernible influence on the downstream pipeline, since the segmentation of the features of interest was performed semi-automatically, as explained later. For each scan, we recorded a set of 3601 projections per stack (each projection with a size of 3000 by 2000 pixels) at 0.1° intervals over a 360° sample rotation. Every single projection was exposed for about a second. Because of the length of the fish, we had to acquire so-called stacked scans; on average, we scanned three fields of view along the rotation axis of the sample holder. This resulted in an average scan time of 10 hours for each scan. The projection images were subsequently reconstructed into stacks of 8-bit png images with NRecon (Bruker microCT, Kontich, Belgium; version 2.2.0.6), without applying any ring artifact or beam hardening correction. The isotropic voxel sizes in the resulting datasets vary from 15 to 17.5 µm.

The result of one such multi-specimen scan is visualized in Figure 2, which shows the six fish contained in this scan. The 3D-printed labels on the holder are highlighted in yellow in Figures 2 and 4–6 (image processing in notebooks), and visible (upside down) in Figure 7.

Data Analysis

Preparation and Handling of Tomographic Datasets

After acquisition, a simple script was used to copy the relevant data to both archival storage and storage accessible by all coauthors. Further processing of the tomographic datasets was performed with a set of Jupyter (Kluyver Thomas et al., 2016) notebooks (David Haberthür, 2026).

Preview Notebook

The preview notebook was used to identify scanning and reconstruction issues. For this, we read all the relevant scan and reconstruction parameters from the log files of each scan. Afterwards, we efficiently loaded the reconstruction png images from disk with the dask_image.imread.imread function (Dask Development Team, 2016). This approach allowed us to map all the reconstructions to memory and quickly generate maximum intensity projections (MIP) of each scan (see Figure 3 for an example) for both quality control and further processing.

Separation Notebook

The separation notebook processes all the acquired scans to extract each individual fish from each scan that encompasses six fish in total. As in the preview notebook, we efficiently load all the pngs from disk with dask (Dask Development Team, 2016).

A connected region labeling (skimage.measure.label) of the previously extracted MIP images was used to detect the single fish in their compartments (see right panel of Figure 4). The 3D-printed compartment labels of the sample holder were detected in the same way and are highlighted in yellow in the left panel of Figure 4. Since skimage numbers the detected regions in raster order (top to bottom, left to right), this numbering depends on the exact position of each fish within its compartment. Due to a slight tilt of the sample holder in the scanner and inconsistent size of the single fish, the detected numbering does not consistently follow the 3D-printed compartment numbering. We thus assigned each detected region to a compartment based on the angle of its centroid (skimage.measure.regionprops) around the rotation axis of the scan, which coincides with the center of the sample holder. Each region is assigned to the compartment with the closest nominal direction (spaced by 60°, counted anticlockwise from compartment 1); empty compartments are skipped.

The bounding boxes of the individual fish in the scan (see right panel of Figure 4) were used to extract small cropped regions from the large multi-specimen dataset. This extraction process is fully reproducible and is tailored to the arrangement of the samples in the custom sample holder.

We then mapped the compartment numbers to the ID numbers of the scanned fish, as noted in the lab book during sample mounting, and presented these images together with photos of the lab book and sample tubes for visual verification (see Figure 5).

We extracted each region of each fish separately from the large reconstructions (with a configurable border buffer; see Figure 6) and wrote these extracted regions to disk in discrete folders for efficient further analysis. As a first step, we wrote the regions of the single fish to disk in zarr (Miles et al., 2020) format, which is well suited to storing n-dimensional arrays on disk. In addition to this, we also wrote a log file for each extracted region, containing all relevant information to redo the cropping step manually in any image processing software. All these log files are also available in the processing repository; an example of such a log file (from the region shown in Figure 7) is shown online here.

Writing the regions as zarr files made it possible to efficiently work with the image data of each extracted fish and to convert those data to any desired format for further analysis. For this further analysis, we also wrote stacks of png images and nrrd files for each fish region in both cropped and cropped-and-binarized forms. These binarized regions were segmented into bone and background based on a simple multi-Otsu thresholding algorithm (廖炳松 et al., 2001) (with four classes, using only the middle threshold). Providing the regions as nrrd files helped to efficiently work with the datasets as specified in the following sections.

Using either K3D-jupyter (Overview, n.d.) or napari (Sofroniew et al., 2026), we implemented a quick way to view any of the extracted regions directly in the Jupyter notebook. K3D-jupyter allows saving an interactive version of the preview. One such interactive preview is available online. napari is well suited to load the intermediate *.zarr files; Figure 7 shows such a preview.

Extraction of Features of Interest

After separation, the cropped image files were checked and rendered using 3D Slicer (Kikinis et al., 2013) and the SlicerMorph extension (Rolfe et al., 2021). A subset of 51 specimens was chosen to create an initial training dataset. For each of these samples, individual elements of the branchial apparatus were rendered and cleaned by hand using a combination of thresholding and “Split Islands” tools to separate the pharyngobranchials, epibranchials, basibranchials, hypobranchials, and ceratobranchials as distinct labels (see Figure 8).

Once segmented, the labels corresponding to each set of bones were exported as a colored labelmap along with the nrrd file from which they were segmented to pass to Biomedisa.

Machine Learning and Model Training

The hand-labeled dataset of 51 samples (including nrrd and label files) was input into Biomedisa (Lösel et al., 2020) to train a segmentation model. We allowed a 180° rotation to account for possible specimen variability, with an 80/20 split between training and validation data. The model was trained with a batch size of 24 for 50 epochs, using a network architecture of 32-64-128-256-512. The final model performs well, with a Dice score of 0.9159 on the validation dataset. We only performed manual touch-ups on the resulting segmentation output where the adjacent bones were in close proximity (causing their appearance to be “stuck” in the final render; this is also an issue with manual segmentation and not exclusive to the approach described here).

Landmarking of Models

To demonstrate the effectiveness of this tool and the importance of 3D morphometry to answer evolutionary questions, we quantified the shape differences of the ceratobranchial bones. Once trained, we applied the Biomedisa segmentation model to the remaining 160 specimens, resulting in segmentation labels for all structures in the branchial apparatus.

As a test and for subsequent analysis, only the first and second right ceratobranchials were chosen for comparison between all samples. The segmentation labels corresponding to these bones were imported together for landmarking using Stratovan Checkpoint (Stratovan Corporation, n.d.) and 3D Slicer (Kikinis et al., 2013).

Type II landmarks were placed at the ends of each bone, with semi-landmarks between them to cover the curvature axes along the bone (see Figure 9).

In total, seven landmarks and four semi-landmark curves (two containing 20 semi-landmarks, two containing 15) were placed on the first ceratobranchial (CB1), and five landmarks and three semi-landmark curves (one containing 20 semi-landmarks, two containing 15) were placed on the second ceratobranchial (CB2). Fixed landmarks were transferred using the “Import/Apply Template” function in Stratovan Checkpoint, and manually inspected and adjusted prior to export. Once curves and fixed points were established, equal distances of semi-landmark curves were ensured using the resample_curves function in 3D Slicer. While this particular workflow was chosen for time-efficiency, accuracy, and user familiarity, any preferred 3D landmarking program should be able to replicate these results.

Analysis of Shape

All subsequent analyses were run using R (version 4.4.1, (R Core Team, 2021)) and the geomorph package (Adams et al., 2012; Baken et al., 2021; M. L. Collyer & Adams, 2018; M. Collyer & Adams, 2018). Landmarks corresponding to each bone were split and analyzed separately after generalized Procrustes analysis (GPA) using the gpagen() function, with principal component analysis (PCA) and linear models run with gm.prcomp() and procD.lm(), respectively. Linear fits were further investigated using the pairwise() function to analyze differences in pairwise statistics.

Results

Micro-CT Imaging Data

The acquisition and reconstruction of fish datasets were successful and enabled high-throughput processing. A total of 215 unique specimens were scanned in 38 different scans with a total scan duration of nearly 16 days (spread over 209 calendar days). We acquired 403312 projections, reconstructed into a total of 154645 slices, resulting in approximately 4000 reconstruction png files per scan (N = 38), covering (on average) a length parallel to the fish of 65 mm.

Fish Separation

Our method reproducibly extracts each of the six fish scanned simultaneously in one scan. The custom-made sample holder aligns each fish along the vertical axis around the rotation axis of the tomographic scan. The extraction based on the MIP image along the rotation axis is deterministic and reproducible. The assignment of detected regions to compartments is computed from their position around the rotation axis, and the only manual input is the mapping of compartment number to fish ID from the lab book.

The resulting image stack, containing six fish in one multi-specimen scan, creates a large output dataset that can be unwieldy to manage even for experienced users. Depending on the available hardware, it may not even be possible to load the resulting image stack into image processing software to manually perform the cropping. The reconstructions of one exemplary dataset (Sticklebucket_H) are stored as losslessly compressed 8-bit PNG files and take up 8.5 GB on disk, but require 42.8 GB of memory when fully loaded in Fiji (Schindelin et al., 2012) (3064 × 3064 pixels × 4891 slices). Large stacks can be loaded as virtual stacks and cropped in XY with low memory requirements, either manually or with a macro implementing steps similar to our pipeline. Done manually, this would need to be repeated six times per scan, followed by assigning each cropped stack to its fish ID in an error-prone manual way.

In addition, such a manual process is not reproducible in the sense that it cannot be consistently replicated by others using the same data, since the manual cropping is operator-dependent. Algorithmically/automatically cropping the large datasets based on the axial MIP image yields reproducible cropped regions and uses operator time efficiently by eliminating manual cropping steps (see Table 1).

Note that the automated extraction writes intermediate files, which facilitate the handling of the data (this step takes about 3 minutes per scan). These files are technically not necessary for the process, but we still took into account the time spent to write them. In addition, human-readable log files documenting the cropping position in the encompassing dataset and the crop extent, as well as images for double-checking the process, are written to disk, which the manual process does not provide reproducibly. This enables reproducible double-checking and confirmation of the process after the fact (see this direct link for one such log file and one such image).

The automated extraction produced 215 nrrd files with a total size of 385 GB (median 1.65 GB). The exported regions varied in physical size, resulting in datasets with a median anteroposterior length of 4250 px and median x- and y-dimensions of 650 px (due to the rotation of the fish, the x- and y-dimensions do not map to anatomical directions). These cropped regions include the region of the originally acquired data where the 3D-printed label is visible, aiding in end-to-end reproducibility of the extraction process.

We have uploaded further cropped datasets of each analyzed fish to MorphoSource (see Supplementary Materials). These further cropped datasets include only the region of the fish in anteroposterior direction, drop the 3D-printed label region and have a median length of approximately 3000 pixels; the x- and y-directions are the same as for the previously mentioned datasets.

Thresholding

The separated fish were segmented based on a simple multi-Otsu thresholding algorithm. This simple segmentation was sufficient to extract all of the features we analyzed further, and we did not have to employ more advanced thresholding methods in our separation pipeline. Selection and individual rendering of branchial structures takes between 10 and 15 minutes; the segmentation of an average Biomedisa render takes 15 seconds once trained (see Table 1).

Analysis

The speed and quality of these data allow us to study the internal branchial anatomy on a large scale and in situ, without the need for fine dissection.

After GPA alignment, we quantified the shape differences among all scanned fish for this project. Changes due to allometry (using centroid size or standard length as the size metric) were significant but slight, explaining only a small fraction of shape variation (CB1: p = 0.001, \(R^2\) = 0.03246; CB2: p = 0.001, \(R^2\) = 0.06220) in both bones. Both linear models and PCA results suggest that the lakes themselves—and not overarching categories of ecotype or sex—drive most of the shape variation in these bones (CB1: p = 0.001, \(R^2\) = 0.09463; CB2: p = 0.001, \(R^2\) = 0.11039). However, the variation in the first ceratobranchial was significantly associated with different ecotypes, although the effect size was small (p = 0.009, \(R^2\) = 0.02056), with substantial overlap in the resulting shape space (see Figure 10).

On the other hand, the second ceratobranchial bone showed equally small but significant changes associated with the ecotype (p = 0.001, \(R^2\) = 0.0377). The differences between benthic and limnetic specimens are, for this bone, clearly divergent in shape space (see Figure 11).

These ecological patterns were further examined at the level of individual lakes (see Figures 12 and 13).

The differences in the second ceratobranchial appear to be driven by divergence in the South Rolly population, supported by significant pairwise differences between this lake and all other lakes in CB2 but not in CB1 (see Supplementary Materials).

Discussion

Data Acquisition

The main driver for developing our pipeline was efficient use of both machine time and human time in this project.

At the voxel size needed to unambiguously resolve the branchial skeleton (less than 20 µm, as determined in pilot experiments), the time needed for acquisition of a tomographic scan is determined by the number of projections and exposure time, not by the number of fish in the field of view. Efficiently “filling” the field of view of the scanner with as many fish as possible in one scan enabled us to scan a total of 215 unique specimens in 38 different scans with a total scan duration of 16 days. If we would have scanned each specimen individually, we would have required about 90 days of scan time; while not impossible, this would have increased the instrument time and its costs severalfold.

Scanning up to six fish at once makes efficient use of the scanner, but does shift the effort to the processing of the data, since each scan has to be separated into single-fish datasets. Since our pipeline automates the separation (see Fish Separation) reproducibly, this shift of effort is cost-effective.

Pipeline and Efficiency

Once all components of the pipeline are combined, running a simple script enables automatic separation, thresholding, and segmentation of stickleback specimens. Requiring minimal active user input, all automated pipeline steps are substantially faster than equivalent manual processing by an expert operator. This reproducible pipeline allows high-throughput sampling and population-scale analysis of stickleback specimens. Adding more (multi-specimen) scans to the pipeline does not increase the time users have to spend on manually separating these scans into single fish regions.

The imaging and separation step is easy to generalize to either other micro-CT setups or other species and specimens. As long as the specimens show up as distinct blobs in the MIP (see Figure 4), our pipeline can reproducibly separate these blobs into a cropped dataset for any further processing, enabling a general, reproducible, and high-throughput micro-CT workflow in biology and other fields. We have already applied the same pipeline to multi-specimen scans of other fish species, namely Spinachia spinachia and several Lake Victoria and Neotropical cichlid species (manuscripts in preparation).

Findings from Ceratobranchial Analysis

The first and second ceratobranchials differ remarkably in morphology, size, and breadth, suggesting that these structures may respond differently to changes in diet, even within the same feeding apparatus.

CB1 showed statistically significant, albeit small, differences between benthic and limnetic ecotypes based on linear and pairwise analyses (PC4, accounting for 4.33% of the variation, appears to capture these differences; see supplementary figure 14). However, we caution that the greater number of benthic lakes and specimens available for this study may contribute to these patterns, as the two groups (benthic and limnetic) exhibit statistically different variances. This pipeline will allow for more extensive analysis from future sampling years to confirm these findings.

Within CB2, limnetic fish (and particularly those from South Rolly Lake) appear to have narrower, less keeled bones than benthic fish. The muscles that attach to the ceratobranchials (m. adductor branchialis, m. abductor filament, and m. obliquus ventralis) are attached along the lateral surfaces of these bones. The increased surface area observed in benthic fish may provide a greater area for muscle attachment, which could enhance the ability of fish to abduct these structures during water filtration (Takata & Sasaki, 2005; Vandewalle et al., 2000). Although the dietary analyses of these fish are still ongoing, these findings suggest that the South Rolly population may have unique dietary specializations and would be expected to feed differently than the fish from the other lakes in this study. In terms of reintroduction, populations from this lake might be expected to fare better than others in terms of limnetic specialization. Indeed, fish with South Rolly heritage outperform Spirit Lake fish in every transplant experiment in which both limnetic populations are included as part of the FITNESS project (Eckert et al., 2026).

The different responses of the first and second ceratobranchials raise the possibility of a modular response to dietary changes within the branchial basket. Studies treating the branchial basket as a single structure, focusing only on the first ceratobranchial, or investigating only external morphology could potentially miss significant changes in shape and size of these structures, and this pipeline provides a wealth of data with which to conduct a follow-up study.

Future Improvements and Issues

As with other multi-specimen scanning projects, the scanning parameters can be optimized individually for each scan but not for each individual specimen. Additionally, atypically large or dense specimens cause an issue for the holder and the reproducibility between scans. As with most machine learning frameworks, the training dataset must encompass the full range of anatomical variation to prevent segmentation errors. And while the separation step outlined in this manuscript is independent of the scanned specimens, the segmentation model has to be retrained for other species or anatomical structures.

Conclusions

The provided pipeline enables repeatable, high-throughput analysis of 3D shape in stickleback specimens. Although applied here to stickleback specimens, the described methods can readily be applied to large-scale sampling efforts of multiple taxonomic groups, with data acquired with diverse imaging setups (not only micro-CT) and different sample holders, due to the use of simple region detection and reproducible logging. Extracting individual specimens from multi-specimen scans is efficient; a custom 3D-printed holder and automated splitting pipeline minimize operator time and maximize throughput.

The reproducible scans and their consistent quality rapidly provide a large amount of similar data, ideal for training machine learning models. In its current implementation, Biomedisa is on average 60 times faster than manual segmentation by a skilled operator, and without the interoperator variability inherent to manually segmenting large numbers of specimens. This puts virtual, nondestructive dissection of internal stickleback anatomy on par with, if not beyond, hand dissection.

Finally, the 3D analysis step of the pipeline allows for insights from 3D data that cannot be obtained from traditional dissection, including complex shapes and arrangements not accessible through destructive sampling approaches.

Authors’ Contributions

Contributor Roles Taxonomy (CRediT), as defined in (n.d.):

• Conceptualization: R. Benjamin Sulser, Sheila Christen, Ruslan Hlushchuk
• Data curation: David Haberthür, R. Benjamin Sulser, Sheila Christen
• Formal analysis: David Haberthür, R. Benjamin Sulser, Sheila Christen
• Funding acquisition: R. Benjamin Sulser, Catherine L. Peichel, Ruslan Hlushchuk
• Investigation: David Haberthür, R. Benjamin Sulser, Sheila Christen
• Methodology: David Haberthür, R. Benjamin Sulser, Sheila Christen, Ruslan Hlushchuk
• Project administration: David Haberthür, R. Benjamin Sulser, Catherine L. Peichel, Ruslan Hlushchuk
• Resources: R. Benjamin Sulser, Ruslan Hlushchuk
• Software: David Haberthür, R. Benjamin Sulser
• Supervision: R. Benjamin Sulser, Catherine L. Peichel, Ruslan Hlushchuk
• Validation: David Haberthür, R. Benjamin Sulser, Sheila Christen
• Visualization: David Haberthür, R. Benjamin Sulser, Sheila Christen
• Writing – original draft: David Haberthür, R. Benjamin Sulser
• Writing – review & editing: David Haberthür, R. Benjamin Sulser, Sheila Christen, Catherine L. Peichel, Ruslan Hlushchuk

Use of AI tools

No generative AI was used to generate the scientific content of the manuscript.

Large language models were used to: - identify spelling and grammatical errors - to suggest wording for individual sentences during the revision - to assist with writing and debugging the analysis code, - and to obtain feedback from the Stanford Agentic Reviewer.

The authors reviewed all suggestions, selectively incorporated them, and take full responsibility for the content of the manuscript and the code.

Competing Interests

Acknowledgments

We are grateful to the Microscopy Imaging Center of the University of Bern for their infrastructural support. We also thank the manubot project (Himmelstein et al., 2019) for facilitating collaborative writing of this manuscript.

Supplementary Materials

Parameters of Tomographic Scans of All Fish

The CSV file ScanningDetails.csv gives a tabular overview of all the (relevant) parameters of all the scans we performed. This file was generated with the data processing notebook and collates the relevant data read from all the log files of all the scans we performed for this study. A copy of each log file containing all scanning parameters is available in a folder in the data processing repository of this project.

Cropped Tomographic Data of Each Extracted Fish

The cropped datasets for analysis include the region with the 3D-printed label of the multi-specimen holder to keep the process reproducible in an end-to-end fashion. We further cropped the datasets along the anteroposterior axis with a helper notebook. These datasets do not include the region of the 3D-printed label, but are minimized to only the fish region. A png stack of each such cropped fish region is available in a collection at MorphoSource and can be downloaded on request.

Analysis Scripts

Jupyter; (Pre)processing of Tomographic Data

All Jupyter notebooks used to process the acquired tomographic data as described in the text are available on GitHub and can be previewed with nbviewer.

3D Slicer; Extraction

The Python script we used to resample curves is found online on GitHub. The result of the resampling is also available online.

R; Morphometrics

The R script we used to assess the morphometrics is available online on GitHub.

Supplementary Figure

References