This manuscript
(permalink)
was generated
from habi/sticklebacks-manuscript@5229755
on July 17, 2026.
This version contains 5166 words.
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
Ben Sulser
0000-0002-8750-0942
·
sulserrb
Evolutionary Ecology Group, Institute of Ecology and Evolution, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland
· Funded by Burgergemeinde Bern
Sheila Christen
0009-0008-3514-0800
·
SheilaChristen
Evolutionary Ecology Group, Institute of Ecology and Evolution, University of Bern, Baltzerstrasse 6, 3012 Bern,
· 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
✉: Correspondence via GitHub Issues or email to Ruslan Hlushchuk.
Can we use high-throughput imaging to understand evolution?
The three-spined stickleback (Gasterosteus aculeatus) is a well-recognized model system for understanding adaptation to divergent habitats. Benthic (benthos-feeding) and limnetic (water-column-feeding) stickleback populations 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, analyses of feeding structures—particularly the jaws and complex internal branchial anatomy—remain time-consuming and require considerable expertise, with destructive sampling and fine dissection skills needed for quantitative analysis.
The advent of X-ray microtomography and 3D-scanning technologies enables non-destructive 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 studies at an eco-evolutionary scale.
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 three-spined stickleback.
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 high-throughput imaging approaches for testing ecological and evolutionary hypotheses.
The three-spined stickleback (Gasterosteus aculeatus) is an oft-studied 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 (WILLACKER et al., 2010). This study will focus on the benthic-limnetic axis, using samples from a long-term evolutionary experiment currently investigating divergent populations of limnetic and benthic stickleback 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 made from four source populations of limnetic sticklebacks and four source populations of benthic sticklebacks—have been placed into 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 in order to track the genotypic and phenotypic trajectories of these introduced populations. Understanding the initial variation in the source populations is essential to this project, as this initial variation would be expected to reflect which phenotypes are associated with each ecotype under study.
Among sticklebacks and other bony fishes, differences in jaw structures are directly related to functional and kinematic differences between different 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, whereas 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 implicated in diet and feeding ecology (BERNER et al., 2008; Schluter & McPhail, 1992). While the shape and arrangement of paired ceratobranchial and pharyngobranchial bones within this complex aid food processing and water vortex generation during feeding (Brooks et al., 2018), the shape of these bones has received comparatively little attention relative 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 render them visible for manual measurements (Ellis & Miller, 2016). These structures are, however, difficult to study without full cranial dissection and corresponding distortion of the branchial anatomy. 3D imaging preserves these features at a high resolution and allows for non-destructive analysis. This work is embedded within the Genomics axis of the Alaska Stickleback Restoration Project, with which Katie Peichel, Ben Sulser, and Sheila Christen are affiliated.
X-ray micro-computed tomography (micro-CT) imaging is an indispensable tool to gain non-destructive insights into the inner structure of highly diverse samples, specifically for specimens studied in the biomedical sciences (Rawson et al., 2020). Microtomographic imaging is ideally suited to non-destructively assess the morphology of different fish species (Ford et al., 2023), including the internal anatomy and small structures difficult to quantify without additional preparation.
While these structures can be rendered by hand by a skilled investigator with segmentation tools such as 3DSlicer (Kikinis et al., 2013)-(Rolfe et al., 2021), the time and cost required per specimen is inefficient for the scale required via eco-evolutionary study and requires destruction of the mandibular and cranial anatomy of the specimen.
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.
The specimens used for this study were collected from source lakes as part of the FITNESS project in the region of Cook Inlet, Alaska. Fish were collected using unbaited minnow traps in two separate field seasons, the first taking place from May 26–June 10, 2023 and the second taking place 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. Fish were euthanized with MS-222, photographed, labeled and preserved in 10% formalin in individual bags, all under 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 are: Watson lake (52), Finger lake (48), Spirit lake (46), South Rolly lake (35), Walby lake (20), and Tern lake (14) for a total of 215 specimens. As an additional step, fish collected in the 2024 field season had DNA fin clips taken via a small piece of the caudal fin in order to sex each individual. Samples underwent hotshot DNA extraction protocol (Meeker et al., 2007) and subsequent PCR followed Stickleback LRRc61 Sexing with primers for this purpose (see Archambeault et al., 2020 (Archambeault et al., 2020)).
Due to their inherent contrast difference to the surrounding tissue, the structures of interest in this study (teeth and bones, i.e., jaws and skull) are well visualized in unstained samples, hence no further preparation of the fish was necessary.
In a small pilot study, we determined the optimal scanning parameters to meet the constraints on total scanning time, resolution, and sample handling. To optimize for these constraints, we scanned all the sticklebacks in batches of six fish in a custom-made 3D printed sample holder in a single scan. This holder was designed in OpenSCAD (OpenSCAD Developers. Version 2021.01) and is available online, either directly as STL file for printing or as (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).
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 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 but one scan where we used a source voltage of 49 kV and 159 µA due to operator error. For each sample, we recorded a set of 3601 projections of approximately 3000 x 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 3 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 then subsequently reconstructed into stacks of 8bit PNG images with NRecon (Bruker microCT, Kontich, Belgium. Version: 2.2.0.6), without applying any ring artefact or beam hardening correction. The isometric voxel sizes in the resulting datasets vary from 15 to 17.5 µm.
After acquisition, a simple script was used to copy the relevant data to both archival storage and storage accessible by all co-authors. Further processing of the tomographic dataset was performed with a set of Jupyter (Kluyver Thomas et al., 2016) notebooks (David Haberthür, 2026).
The preview notebook is used to identify issues with the scanning.
For this, we read all relevant scanning 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 2 for an example) for both quality control and further processing.
The separation notebook processes all the acquired scans to extract each individual fish from each scan encompassing six fish in total.
As in the preview notebook, we efficiently load all the PNGs from disk with dask (Dask Development Team, 2016).
Based on the previously extracted MIP images and a simple labeling of these images (skimage.measure.label), we extract both the labels in the custom-made sample holder and the positions of individual fish in the scan (skimage.measure.regionprops) (see Figure 3).
This extraction is completely reproducible and well-adapted to the custom-made sample holder.
Based on a simple mapping of the detected region to the ID numbers of the scanned fish, we labeled the resulting images and presented these images together with photos of the lab book and sample tubes for verification (see Figure 4).
The skimage.measure.regionprops function used for labeling not only returns the positions of each detected fish, but also the extent of the bounding box of each detected region.
We extracted each region of each fish separately out of the large reconstructions (with a configurable border buffer, see Figure 5) and wrote these extracted regions to disk in discrete folders for efficient further analysis.
In a first step, we wrote the regions of the single fish to disk in zarr (Miles et al., 2020) format, which is a preferred format to store 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 completely manually (an example of such a log file is shown as part of the processing repository).
Writing the regions as zarr files made it possible to efficiently work with the image data of each extracted fish and to convert that data to any desired format for further analysis.
For this further analysis, we also wrote stacks of PNG images and, additionally, 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-level Otsu thresholding method (廖炳松 et al., 2001).
Providing the regions as NRRD files helped to efficiently work with the datasets as specified in the following sections.
Using K3D-jupyter (Overview — K3D-Jupyter Documentation, n.d.) we implemented a quick way to view any of the extracted regions directly in the Jupyter notebook (see Figure 6).
An interactive version of this figure is available online.
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). The individual elements of the branchial apparatus were rendered using a combination of thresholding and ‘Split Islands’ tools to separate the pharyngobranchials, epibranchials, basibranchials, hypobranchials and ceratobranchials (see Figure 7).
Once rendered, these bones were exported as a colored labelmap alongside the NRRD file from which they were segmented to pass to the Biomedisa program.
As a group, a dataset of 51 specimens (including NRRD and .label files) was passed to Biomedisa (Lösel et al., 2020) to train a segmentation model.
We allowed for rotation of 180° to account for possible specimen variability, and an 80/20 split between training and validation data.
The model was trained with a batch size of 24 and 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.
Manual touchups were only needed and performed where bones were extremely close together (causing their appearance to be “stuck” in the final render; this is an issue with manual segmentation as well).
To demonstrate the effectiveness of this tool and the importance of 3D morphometrics for answering eco-evolutionary questions, we have run a demonstration quantifying the shape differences of the ceratobranchial bones.
Once trained, we applied the Biomedisa segmentation model to the remaining 160 specimen volumes and landmarked the final results using Stratovan Checkpoint (Stratovan Corporation, n.d.).
As a test and for subsequent analysis, the first and second right ceratobranchials were chosen for comparison across all specimens.
Type II landmarks were set on the ends of each bone, with semilandmarks in-between each to cover axes of curvature along the bone (see Figure 8).
In total, 7 landmarks and 4 semilandmark curves (two containing 20 semilandmarks, two containing 15) were placed on the first ceratobranchial (CB1), and 5 landmarks and 3 semilandmark curves (one containing 20 semilandmarks, two containing 15) were placed on the second ceratobranchial (CB2).
Equal distances were ensured using the resample_curves function in 3D Slicer.
All subsequent analyses were run using R (version 4.4.1, (R Core Team, 2021)) and the geomorph package (Adams & Otárola‐Castillo, 2013).
Both bones 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 via the pairwise() function to analyze differences in pairwise statistics.
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 scanning duration of nearly 16 days. We acquired 136838 projections, reconstructed into a total of 154622 reconstructions, resulting in approximately 4000 reconstruction PNG files per scan (N=38).
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 completely automated and very robust, since the detected fish ‘regions’ do not overlap in the resulting image.
Depending on the available hardware, it may not even be possible to load the full stack of each scan into software to manually perform the cropping, such as Fiji (Schindelin et al., 2012).
Large stacks of images (in other words larger than the available RAM of the available machine) can be loaded as virtual stacks, but to manually crop the region of each fish from the large scan with the Crop (3D) function, one needs to load the full dataset.
Since one (exemplary) dataset (Sticklebucket_10) is 7 GB on disk and reported as 35.4 GB when loaded in Fiji, using the 3D cropping function on an uncropped single dataset is not possible without a powerful workstation.
Extracting individual fish from the encompassing dataset would thus be a two-step manual process, e.g. cropping the full dataset loaded as virtual stack and then cropping it down further before writing out the cropped stack. For each encompassing scan this would need to be repeated 6 times (for each of the 6 fish in each of the encompassing scans). 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 leads to both reproducible cropped regions and efficiently uses the operator time by eliminating manual cropping steps (see Table 1).
| Task | Est. Manual Time | Pipeline Time | Speed-up |
|---|---|---|---|
| Scanning single fish | ~10 hours | 1 hour, 45 minutes | ~5.5 x |
| Splitting scans into single fish | 17 minutes | 9 minutes | ~2 x |
| Rendering volumes | 5 minutes | 15 seconds | 20 x |
| Segmentation | 10-15 minutes | 15 seconds | ~60 x |
Note that the automated extraction process writes intermediate files during the extraction process which facilitate the handling of the data (this process takes about 3 minutes per scan). These files are technically not necessary for the process, but we still accounted for 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 extraction and sampling process led to a total of ~64 GB of NRRD files, which were assessed as specified before.
The separated fish were segmented based on a simple multi-level Otsu thresholding method. This relatively simple segmentation was sufficient to extract all 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 the branchial structures takes between 10-15 minutes; the average Biomedisa render takes 2.5 minutes once trained (see Table 1).
The speed and quality of these data allow us to study the internal branchial anatomy at scale and in situ, without the need for fine dissection.
Numerous studies have shown the relationships between gill rakers (bony protrusions arising from the branchial complex) and diet (BERNER et al., 2008; Schluter & McPhail, 1992).
While the shape and arrangement of the ceratobranchials and the corresponding bony gill rakers are hypothesized to work in tandem for food processing and water vortex generation during suspension feeding (Brooks et al., 2018), the shape of these bones has received comparatively little attention.
This is likely due to the flattening and destructive sampling used in traditional raker counting methods, which dissect and deform these structures to render them visible for manual measurement. 3D imaging preserves these features at high resolution.
After GPA alignment, we quantified the shape differences among all scanned fish for this project. Changes due to allometry (using the metric of centroid size or standard length of the fish) were significant but slight, explaining only a small fraction of shape variation 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.03246, CB2: p = 0.001, R^2 = 0.06220). Variation in the first ceratobranchial was significantly associated with lake of origin, although the effect size was small (p = 0.009, R^2 = 0.02056), and with a substantial overlap in the resulting shape space (see Figure 9).
The second ceratobranchial bone, on the other hand, shows equally small yet significant shifts associated with the ecotype (p = 0.001, R^2 = 0.0377). The differences between benthic and limnetic gill rakers are, for this bone, clearly divergent in shape space (see Figure 10).
These ecological patterns were further examined at the level of individual lakes (see Figures 11 and 12).
The differences in the 2nd ceratobranchial appear to be driven by divergence in the South Rolly population, supported by significant pairwise differences observed between this lake and all other lakes observed in CB2 and not in CB1 (see supplementary information).
Once all components of the pipeline are combined, running a simple script enables automatic reconstruction, thresholding, and segmentation of stickleback specimens. All automated pipeline steps are substantially faster than equivalent manual processing by an expert operator, with minimal active user input. This reproducible pipeline allows for high-throughput sampling and population-scale analysis of stickleback specimens.
The imaging and separation step is easy to generalize to either other micro-CT setups or other species/specimen. As long as the specimen show up as distint blobs in the MIP (see Figure 3) our pipeline can reproducibly separate these blobs into cropped dataset for any further processing, enabling a general, reproducible and high-throughput micro-CT workflow in biology and many other fields.
The first and second gill rakers differ remarkably in morphology, size, and breadth, suggesting that these structures may respond differently to shifts in diet—even within the same feeding apparatus. CB1 showed statistically significant differences based on linear and pairwise analyses (and PC4, accounting for 4.33% of the variation, appears to capture these differences; see supplementary figures). However, we caution that the greater number of lakes and specimens available for this study may contribute to these patterns, as the two groups exhibited 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 greater area for muscle attachment, which may enhance the ability of fish to abduct these structures during water filtration (Takata & Sasaki, 2005; Vandewalle et al., 2000). While 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 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 populations are included as a part of the FITNESS project (Eckert et al., 2026).
The different response of the first and second ceratobranchial raises the possibility of a modular response to dietary shifts within the branchial basket. Studies treating the branchial basket as a single structure, focusing on the first ceratobranchial, or investigating external morphology might potentially miss significant changes in shape and size of these structures—and this pipeline provides a wealth of data with which to perform a follow-up study.
As with many multi-scan projects, the scanning parameters can be optimized individually for each scan but not for each individual specimen. In addition, atypically large or dense specimens cause an issue for the holder and the reproducibility across scans. As with most machine learning approaches, it is also important to ensure that the full range of variations across the entire dataset are represented in training to avoid erroneous segmentations.
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 on many different micro-CT machines and different sample holders, due to the use of simple region detection and reproducible logging. The extraction of individual specimens from multi-specimen scans is efficient, using a custom 3D-printed sample holder and automated splitting procedure to reduce the need for active operator time and maximize automated processing.
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 performs on average 60 times faster than manual segmentation by a skilled operator, and without the inter-operator variability inherent to manually segmenting large numbers of specimens. This brings virtual, non-destructive dissection of internal stickleback anatomy up to parity with hand-dissected methods.
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 possible through destructive sampling approaches.
Contributor Roles Taxonomy (CRediT), as defined in (n.d.):
| Author | Competing Interests | Last Reviewed |
|---|---|---|
| David Haberthür | Nothing to declare | 2026-01-14 |
| Ben Sulser | Nothing to declare | |
| Sheila Christen | Nothing to declare | 2026-07-13 |
| Catherine L. Peichel | Nothing to declare | 2026-01-19 |
| Ruslan Hlushchuk | Nothing to declare | 2026-01-19 |
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.
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.
All Jupyter scripts to process the acquired tomographic data as described in the text are found online at GitHub, and can be easily previewed online, too.
The Python script we used to resample curves is found online at GitHub. The results of the resampling is available online, too.
The R script we used to assess the morphometrics is found online at GitHub.