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 Bern Burgergemeinde
Sheila Christen
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-0001-2345-6789
· RuslanHlushchuk
microCT research group, Institute of Anatomy, University of Bern, Baltzerstrasse 2, 3012 Bern, Switzerland
✉ — Correspondence possible via GitHub Issues
or email to
Ruslan Hlushchuk <ruslan.hlushchuk@unibe.ch>.
Abstract
Can we predict evolution?
The three-spined stickleback (Gasterosteus aculeatus) is a well-recognized system for understanding adaptation to divergent habitats.
Populations of benthic and limnetic stickleback differ in a number of phenotypic traits that are associated with shifts in dietary specialization.
However, analyses of the structures required for feeding – especially the jaws and complex internal branchial anatomy – requires considerable time and expertise, with destructive sampling and fine dissection skills needed for quantitative analysis.
The advent of µCT and 3D-scanning technology affords non-destructive sampling and an increase the resolution of data available for study, but at the substantial cost of increasing complexity and processing time for each specimen.
To address these concerns, we developed a rapid and semi-automated segmentation and analysis pipeline based on both the Jupyter interactive development environment and the Biomedisa image segmentation platform for investigating three-dimensional morphological adaptation within the three-spined stickleback.
The pipeline includes splitting 44 multi-specimen scans into regions of interest for each specimen, reconstruction of targeted anatomy, and morphometric analyses.
We then applied this pipeline to a sampling effort encompassing hundreds of samples of divergent benthic and limnetic stickleback populations (N=216), showcasing the possibility of using high-throughput scanning data to provide tests of ecological and evolutionary hypotheses.
Introduction
The threespine stickleback (Gasterosteus aculeatus) is an oft-studied organism for understanding the independent evolution of similar traits in similar environments.
This species exhibits marked differences in marine-freshwater, lake-stream, and benthic-limnetic ecotypes4.
This study will focus on the benthic-limnetic axis, using samples from a long-term evolutionary experiment currently in-process studying divergent populations of limnetic and benthic stickleback within the Kenai peninsula of Alaska (USA).
This project, the Forward In Time Natural Experimental Study of Selection (FITNESS), aims to study the predictability and repeated of evolution.
Two pools of sticklebacks — one made from four source populations of limnetic 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 endeavor, as this initial variation would be expected to reflect which phenotypes are associated with each ecotype under study.
Among these and other bony fish, the internal hyoid arch-branchial arch complex is an important structure implicated in diet and feeding ecology.
The pharyngobranchial bones at the posterior margin of this complex (and their corresponding ceratobranchials) aid in dietary processing.
The “gill rakers” on the ceratobranchials and the “internal jaws” of the are strong indicators of different foraging styles, and have been shown to be related to genotypic and phenotypic changes in different species, and even among different populations.
These structures are, however, difficult to study without full cranial dissection and distortion of the branchial anatomy.
Why are we studying what we are studying in these fish?
Micro-computed tomography
X-ray microtomography is an indispensable tool to gain non-destructive insights into the inner structure of highly diverse samples, namely 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, large and small (Ford et al., 2023), including the internal anatomy.
While these structures can be rendered by-hand by a skilled investigator, the time and cost required per-specimen is inefficient for the scale required via eco-evolutionary study.
This project aims to address these gaps, demonstrating a novel pipeline for rendering and auto-splitting of a multispecimen scan for mass sampling, creating a dataset with consistent parameters that can be fed to donwstream machine-learning techniques [Biomedisa] 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 rendered structures for each specimen in a fraction of the time and resources used in traditional analysis.
Materials & Methods
Figure 1: Workflow overview
Sample procurement and preparation
The specimens used for this study were collected from source lakes as a 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.
Specimen collections were taken from a random sample of 30 fish from each lake.
Fish were euthanized with MS-222, photographed, and preserved in formalin in a bag with a specific label.
At the end of each field season, samples were shipped from Anchorage (AK, USA) to Bern (BE, CH) where they were stored until scanning time.
Due to their inherent contrast difference to the surrounding tissue, the structures of interest we touch upon in this manuscript (teeth and bones, i.e. jaws and skull) are well visualized in unstained samples, hence no further preparation of the fish was necessary.
microtomographic imaging
In a small pilot study we determined the optimal scanning parameters to fit 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 generated with OpenSCAD 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 at the Institute of Anatomy, University of Bern, Switzerland.
In total we performed 44 scans, each of the scan usually containing 6 fish in the sample holder.
The relevant details of each scan are collated 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 every 0.1° over a 360° sample rotation.
Every single projection was exposed for about a second (depending on the sample).
Due to 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 scan times between 3 to 5 hours.
The projection images were then subsequently reconstructed into a 3D stack of images with NRecon (Bruker microCT, Kontich Belgium, Version: 2.1.0.1 or 2.2.0.6).
The isometric voxel sized in the resulting datasets vary from 15 to 19 µm.
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 co-authors at the same time.
Further processing of the tomographic dataset was performed with a set of Jupyter (Kluyver Thomas et al., 2016)notebooks(David Haberthür, 2026).
The scripts are freely available online under the MIT License and may be freely used, modified, and redistributed for research, teaching, and other non-military purposes.
Preview notebook
The preview notebook is used for surfacing issues with the scanning.
For this, we read all relevant scanning and reconstruction parameters from the log files of each scan.
Afterwards, we efficiently loading the reconstruction PNG images from disk with the dask_image.imread.imread function (Dask Development Team, 2016).
Like so, we can map all the generated 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.
Figure 2: Maximum intensity projections of one acquired dataset along all three cardinal axes.
Separation notebook
The separation notebook processes all the performed scans to extract the single fish out from each scan, where 6 fish have been scanned.
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 the image (skimage.measure.label), we extract both the labels in the custom-made sample holder and the positions of single fish in the scan (skimage.measure.regionprops) (see Figure 3).
Figure 3: Automatically detected regions based on maximum intensity projection along the rotation axis of the tomographic scan.
The regions are numbered consecutively from the top left to the bottom right.
These numbers are mapped to the correct fish ID in the next step.
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 double-checking (see Figure 4).
Figure 4: Mapping lab book notes, photos and detected regions to fish ID.
The skimage.measure.regionprops function we used for labeling returns not only the positions of the detected fish, but also the extent of the bounding box of the region of the fish shown in the original image.
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).
Figure 5: Double-checking crop extent and fish ID.
The top row shows the extracted regions from the previously calculated MIP of the full scan, the bottom row shows the MIP images of the extracted regions.
Both rows must show exactly the same region.
Saving out 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 wrote out stacks of PNG images and additionally, as nrrd files for each fish region as a simple crop out of the original dataset and as binarized regions, which are segmented into bone and background based on a simple multi-level Otsu thresholding method (廖炳松 et al., 2001).
We actually scanned N=216 fishes (unique Fish IDs)
The total scanning duration was 18 days, 12 hours and 6 minutes
In total, we acquired 158444 projections
These were reconstructed into a total of 177749 reconstructions
Which is about 4040 files per scan (N=44)
Analysis
3D shape variation on internal pharyngobranchial bone.
Only possible to get this information in 3D.
Discussion
Both repeatable and reproducible research
Automated cropping out of single fish from combined scan very efficiently uses machine time.
Several fish can be scanned together, splitting is performed after the fact, reproducibly and without manual input.
Combination of methods cuts down on time a lot.
Biomedisa makes more “extraction” possible.
Other biological questions can be answered, too.
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 the fishes
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.
A copy of each log file containing all scanning parameters is available in a folder in the data processing repository.
Ford, K. L., Albert, J. S., Summers, A. P., Hedrick, B. P., Schachner, E. R., Jones, A. S., Evans, K., & Chakrabarty, P. (2023). A New Era of Morphological Investigations: Reviewing Methods for Comparative Anatomical Studies. Integrative Organismal Biology, 5(1). https://doi.org/10.1093/iob/obad008
Himmelstein, D. S., Rubinetti, V., Slochower, D. R., Hu, D., Malladi, V. S., Greene, C. S., & Gitter, A. (2019). Open collaborative writing with Manubot. PLOS Computational Biology, 15(6), e1007128. https://doi.org/10.1371/journal.pcbi.1007128
Kluyver Thomas, Ragan-Kelley Benjamin, Pérez Fernando, Granger Brian, Bussonnier Matthias, Frederic Jonathan, Kelley Kyle, Hamrick Jessica, Grout Jason, Corlay Sylvain, Ivanov Paul, Avila Damián, Abdalla Safia, Willing Carol, & Jupyter Development Team. (2016). Jupyter Notebooks &ndash; a publishing format for reproducible computational workflows. In Positioning and Power in Academic Publishing: Players, Agents and Agendas. IOS Press. https://doi.org/10.3233/978-1-61499-649-1-87
Lösel, P. D., van de Kamp, T., Jayme, A., Ershov, A., Faragó, T., Pichler, O., Tan Jerome, N., Aadepu, N., Bremer, S., Chilingaryan, S. A., Heethoff, M., Kopmann, A., Odar, J., Schmelzle, S., Zuber, M., Wittbrodt, J., Baumbach, T., & Heuveline, V. (2020). Introducing Biomedisa as an open-source online platform for biomedical image segmentation. Nature Communications, 11(1). https://doi.org/10.1038/s41467-020-19303-w
Miles, A., Kirkham, J., Durant, M., Bourbeau, J., Onalan, T., Hamman, J., Zain Patel, Shikharsg, Rocklin, M., Dussin, R., Schut, V., Andrade, E. S. D., Abernathey, R., Noyes, C., Sbalmer, Pyup.Io Bot, Tran, T., Saalfeld, S., Swaney, J., … Banihirwe, A. (2020). zarr-developers/zarr-python: v2.4.0 (Version v2.4.0). Zenodo. https://doi.org/10.5281/zenodo.3773450
Overview — K3D-jupyter documentation. (n.d.). Retrieved March 13, 2026, from https://k3d-jupyter.org/
Rawson, S. D., Maksimcuka, J., Withers, P. J., & Cartmell, S. H. (2020). X-ray computed tomography in life sciences. BMC Biology, 18(1). https://doi.org/10.1186/s12915-020-0753-2
廖炳松, 陳澤生, & 詹寶珠. (2001). A Fast Algorithm for Multilevel Thresholding. Journal of Information Science and Engineering, 17(5). https://doi.org/10.6688/jise.2001.17.5.1
0000-0003-3388-9187
· 
habi
· 
@habi@mastodon.social
