microCT imaging of 215 three-spined stickleback

This manuscript (permalink) was generated from habi/sticklebacks-manuscript@5229755 on July 17, 2026.
This version contains 5166 words.

Authors

✉: Correspondence via GitHub Issues or email to Ruslan Hlushchuk.

Abstract

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.

Introduction

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.

Micro-computed tomography

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.

Materials & Methods

Figure 1: Workflow overview

Sample procurement and preparation

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.

micro-CT imaging

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.

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. Further processing of the tomographic dataset was performed with a set of Jupyter (Kluyver Thomas et al., 2016) notebooks (David Haberthür, 2026).

Preview notebook

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.

Figure 2: Maximum intensity projections of one acquired dataset along all three cardinal axes.
Separation notebook

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.

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 verification (see Figure 4).

Figure 4: Mapping lab book notes, photos, and detected regions to fish ID.

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).

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.

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.

Figure 6: Three-dimensional preview of extracted region, automatically thresholded.

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). 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).

Figure 7: Example of Branchial Anatomy with CB1 and CB2 highlighted. Abbreviations: PB = pharyngobranchials, EB = epibranchials, BB = basibranchials, HB = hypobranchials, CB = ceratobranchials.

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.

Machine learning and model training

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).

Landmarking of models

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.

Figure 8: Landmarking shown on dorsal (a) and ventral surfaces (b).

Analysis of shape

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.

Results

micro-CT imaging data

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).

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 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).

Table 1: Estimates of time comparisons between manual and pipeline runs. Note that it would have been completely infeasible to scan each fish separately. On average, we scanned 5.66 fish per scan for a total scan time of 15 days, 23 hours, 50 minutes and 24 seconds. The manual scan time per single fish is thus a calculated 5 times longer than scanning 6 fish at a time in a single scan and separating them after the fact. Also note that the manual splitting does not produce the accompanying log file and figures for double-checking as specified below.
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.

Thresholding

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).

Analysis

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).

Figure 9: PCA of CB1 colored by lake ecotype. Warps are indicated at the extremes of each axis by vectors drawn from a mean shape (red).

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).

Figure 10: PCA of CB2 colored by lake ecotype. Warps are indicated at the extremes of each axis by vectors drawn from a mean shape (red).

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

Figure 11: PCA of CB1 colored by individual lake. Warps are indicated at the extremes of each axis by vectors drawn from a mean shape (red).
Figure 12: PCA of CB2 colored by individual lake. Warps are indicated at the extremes of each axis by vectors drawn from a mean shape (red).

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).

Discussion

Pipeline and efficiency

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.

Findings from Ceratobranchial Analysis

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.

Future improvements and issues

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.

Conclusion

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.

Author Contributions

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

Competing Interests

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

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 the 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.

Analysis scripts

Jupyter; (pre)processing of tomographic data

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.

3D Slicer; extraction

The Python script we used to resample curves is found online at GitHub. The results of the resampling is available online, too.

R; morphometrics

The R script we used to assess the morphometrics is found online at GitHub.

References

(n.d.). ANSI/NISO Z39.104-2022, CRediT, Contributor Roles Taxonomy. NISO. https://doi.org/10.3789/ansi.niso.z39.104-2022
Adams, D. C., & Otárola‐Castillo, E. (2013). geomorph: an<scp>r</scp>package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution, 4(4), 393–399. https://doi.org/10.1111/2041-210x.12035
Archambeault, S. L., Bärtschi, L. R., Merminod, A. D., & Peichel, C. L. (2020). Adaptation via pleiotropy and linkage: Association mapping reveals a complex genetic architecture within the stickleback<i>Eda</i>locus. Evolution Letters, 4(4), 282–301. https://doi.org/10.1002/evl3.175
Bell, M. A., & Foster, S. a. (Eds.). (1994). The Evolutionary Biology of the Threespine Stickleback. Oxford University PressOxford. https://doi.org/10.1093/oso/9780198577287.001.0001
BERNER, D., ADAMS, D. C., GRANDCHAMP, A. ‐C., & HENDRY, A. P. (2008). Natural selection drives patterns of lake–stream divergence in stickleback foraging morphology. Journal of Evolutionary Biology, 21(6), 1653–1665. https://doi.org/10.1111/j.1420-9101.2008.01583.x
Brooks, H., Haines, G. E., Lin, M. C., & Sanderson, S. L. (2018). Physical modeling of vortical cross-step flow in the American paddlefish, Polyodon spathula. PLOS ONE, 13(3), e0193874. https://doi.org/10.1371/journal.pone.0193874
Dask Development Team. (2016). Dask: Library for dynamic task scheduling. https://dask.org
David Haberthür. (2026). habi/sticklebacks: Maintenance release (Version 0.99) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.18257528
Eckert, L., Bolnick, D. I., Derry, A. M., Haines, G. E., Heckley, A. M., Lind, Å. J., Peichel, C. L., Roth, A. M., Steinel, N. C., Vlahiotis, K., Weber, J. N., Hendry, A. P., & Barrett, R. D. H. (2026). Intrinsic fitness differences outweigh environmental matching in shaping introduction outcomes in nature. openRxiv. https://doi.org/10.64898/2026.02.04.699496
Ellis, N. A., & Miller, C. T. (2016). Dissection and Flat-mounting of the Threespine Stickleback Branchial Skeleton. Journal of Visualized Experiments, (111). https://doi.org/10.3791/54056
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
Haberthür, D. (2019). TomoGraphics/Hol3Drs: A release (Version v0.618) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.2587555
Haines, G. E., Stuart, Y. E., Hanson, D., Tasneem, T., Bolnick, D. I., Larsson, H. C. E., & Hendry, A. P. (2020). Adding the third dimension to studies of parallel evolution of morphology and function: An exploration based on parapatric lake‐stream stickleback. Ecology and Evolution, 10(23), 13297–13311. https://doi.org/10.1002/ece3.6929
Hendry, A. P., Barrett, R. D. H., Bell, A. M., Bell, M. A., Bolnick, D. I., Gotanda, K. M., Haines, G. E., Lind, Å. J., Packer, M., Peichel, C. L., Peterson, C. R., Poore, H. A., Massengill, R. L., Milligan‐McClellan, K., Steinel, N. C., Sanderson, S., Walsh, M. R., Weber, J. N., & Derry, A. M. (2024). Designing eco‐evolutionary experiments for restoration projects: Opportunities and constraints revealed during stickleback introductions. Ecology and Evolution, 14(6). https://doi.org/10.1002/ece3.11503
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
Kikinis, R., Pieper, S. D., & Vosburgh, K. G. (2013). 3D Slicer: A Platform for Subject-Specific Image Analysis, Visualization, and Clinical Support. In Intraoperative Imaging and Image-Guided Therapy (pp. 277–289). Springer New York. https://doi.org/10.1007/978-1-4614-7657-3_19
Kluyver Thomas, Ragan-Kelley Benjamin, P&eacute;rez Fernando, Granger Brian, Bussonnier Matthias, Frederic Jonathan, Kelley Kyle, Hamrick Jessica, Grout Jason, Corlay Sylvain, Ivanov Paul, Avila Dami&aacute;n, Abdalla Safia, Willing Carol, & Jupyter Development Team. (2016). Jupyter Notebooks &amp;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
McGee, M. D., Schluter, D., & Wainwright, P. C. (2013). Functional basis of ecological divergence in sympatric stickleback. BMC Evolutionary Biology, 13(1), 277. https://doi.org/10.1186/1471-2148-13-277
Meeker, N. D., Hutchinson, S. A., Ho, L., & Trede, N. S. (2007). Method for Isolation of PCR-Ready Genomic DNA from Zebrafish Tissues. BioTechniques, 43(5), 610–614. https://doi.org/10.2144/000112619
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) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.3773450
Overview — K3D-jupyter documentation. (n.d.). Retrieved July 17, 2026, from https://k3d-jupyter.org/
R Core Team. (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.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
Reid, K., Bell, M. A., & Veeramah, K. R. (2021). Threespine Stickleback: A Model System For Evolutionary Genomics. Annual Review of Genomics and Human Genetics, 22(1), 357–383. https://doi.org/10.1146/annurev-genom-111720-081402
Rolfe, S., Pieper, S., Porto, A., Diamond, K., Winchester, J., Shan, S., Kirveslahti, H., Boyer, D., Summers, A., & Maga, A. M. (2021). SlicerMorph: An open and extensible platform to retrieve, visualize and analyse 3D morphology. Methods in Ecology and Evolution, 12(10), 1816–1825. https://doi.org/10.1111/2041-210x.13669
Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D. J., Hartenstein, V., Eliceiri, K., Tomancak, P., & Cardona, A. (2012). Fiji: an open-source platform for biological-image analysis. Nature Methods, 9(7), 676–682. https://doi.org/10.1038/nmeth.2019
Schluter, D., & McPhail, J. D. (1992). Ecological Character Displacement and Speciation in Sticklebacks. The American Naturalist, 140(1), 85–108. https://doi.org/10.1086/285404
Stratovan Corporation. (n.d.). Stratovan checkpoint (Version 2025.06.02.09:11 WIN x64). Retrieved https://www.stratovan.com/products/checkpoint
Takata, Y., & Sasaki, K. (2005). Branchial structures in the Gasterosteiformes, with special reference to myology and phylogenetic implications. Ichthyological Research, 52(1), 33–49. https://doi.org/10.1007/s10228-004-0251-5
Vandewalle, P., Parmentier, E., & Chardon, M. (2000). The branchial basket in Teleost feeding. Cybium, 24(4). https://orbi.uliege.be/handle/2268/14698
WILLACKER, J. J., VON HIPPEL, F. A., WILTON, P. R., & WALTON, K. M. (2010). Classification of threespine stickleback along the benthic-limnetic axis. Biological Journal of the Linnean Society, 101(3), 595–608. https://doi.org/10.1111/j.1095-8312.2010.01531.x
廖炳松, 陳澤生, & 詹寶珠. (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