A Comparative Analysis of Experimental Selection on the Stickleback Pelvis

The stickleback pelvis offers a compelling case study for understanding evolutionary adaptation. compare.edu.vn delves into a comparative analysis of experimental selection pressures affecting stickleback pelvic morphology, highlighting the role of environmental factors and genetic mechanisms in driving evolutionary change and offering a solution to understanding these variations. Understanding these evolutionary adaptations can provide valuable insights into evolutionary biology, adaptive evolution, and natural selection.

1. Introduction: The Significance of Polymorphism

Polymorphism, the presence of multiple forms of a trait within a population, is a cornerstone of evolutionary biology. Its persistence, especially at frequencies exceeding those expected from simple mutation-selection balance, signifies the influence of balancing selection or heterosis. These mechanisms, where multiple selective forces maintain diversity, offer invaluable opportunities to dissect the genetic and ecological underpinnings that sustain variation. These mechanisms are the key to diversity.

One excellent system for studying evolutionary changes are sticklebacks. Through stickleback studies, we learn about evolutionary processes (Bolnick, 2004; Chan et al., 2010; Colosimo et al., 2005; Hendry et al., 2013; Reimchen, 1994; Schluter, 2000; Schluter et al., 2021; Schluter & McPhail, 1992). The importance of sticklebacks as a model system (Reid et al., 2021) has spurred interest in adaptive polymorphisms within Gasterosteidae, the stickleback family. Spine reduction or loss has evolved separately in four stickleback genera, Gasterosteus, Pungitius, Apeltes, and Culaea, which share a common ancestor approximately 27 million years ago (Blouw & Hagen, 1984a; Chan et al., 2010; Hagen & Blouw, 1983; Nelson, 1969, 1977; Varadharajan et al., 2019). Predation by birds, fish, mammals, and insects has been shown to act as a selective agent on spines in stickleback species (Blouw & Hagen, 1984b; Foster et al., 1988; Hall, 1956; Hoogland et al., 1957; Lescak et al., 2012; Lescak & von Hippel, 2011; Marchinko, 2009; Miller et al., 2017; Reimchen, 1980, 1992, 1994; Reisman & Cade, 1967; Reist, 1980a, 1980b, 1983; Rundle et al., 2003; Vamosi, 2002; Vamosi & Schluter, 2004). Studies demonstrate that gape-limited predators such as trout (Oncorhynchus mykiss) or juvenile pike (Esox lucius) favor increased spine length, whereas predators such as adult pike or larvae of Dytiscus spp. favor reduced spines or spine loss (Reist, 1980a, 1980b). Spine polymorphism is maintained within stickleback populations by selective balance driven by predators (e.g., gape-limited and non-gape-limited) throughout space, seasons, or years (Reimchen, 1994). Spine may be selected against in calcium-deficient environments because of the cost of growing these bony structures (Bell et al., 1993; Giles, 1983; Haines et al., 2023; Spence et al., 2013). Spine polymorphism can be maintained by selection for spines from gape-limited predators and against spines because of physiological costs in low-calcium environments (e.g., in some freshwater lakes).

The ecological correlates of spine polymorphism are characterized in threespine sticklebacks (Bell & Foster, 1994; Haines et al., 2023; Klepaker et al., 2013). These ecological correlates include differences between spine morphs in predation, habitat preference, and diet, which can give rise to additional morphological or behavioral differences between spine morphs. Pelvic spine divergence in some lakes can be associated with morphological and ecological divergence between reproductively isolated limnetic and benthic threespine sticklebacks, where the benthic species often evolves reduced spines (Boughman, 2001; Hatfield, 1997; Hatfield & Schluter, 1999; Peichel et al., 2001; Rundle & Schluter, 2004; Schluter, 1993; Schluter & McPhail, 1992). Even without reproductive isolation, a well-established difference associated with spine reduction in threespine sticklebacks is that spined sticklebacks consume more planktonic food sources and unspined sticklebacks consume more benthic food items, which is associated with a difference in gill raker morphology where fewer, shorter, and thicker gill rakers are found in benthic feeding fish (Matthews et al., 2010; McPhail, 1984; Nagel & Schluter, 1998; Reimchen, 1980; Ridgway & McPhail, 1984; Rundle et al., 2003; Schluter, 1993). It is unclear whether spine reduction that occurred in other stickleback species is associated with ecological and morphological changes, such as changes in diet or gill raker morphology.

In brook sticklebacks (Culaea inconstans Kirtland, 1840), within–population pelvic spine polymorphism does not involve reproductively isolated morphs (Lowey et al., 2020) and the ecological and morphological changes associated with pelvic spine reduction are not consistent among populations (Willerth et al., 2022). To date, no study has addressed whether pelvic spine reduction is associated with diet differences in brook sticklebacks similarly to the differences in diet in pelvic-reduced relative to spined threespine sticklebacks. The extent of ecological or phenotypic parallelism among stickleback species associated with pelvic spines is unknown. Answering this question on parallelism among stickleback species will address questions about phenotypic parallelism across taxa (Bolnick et al., 2018; De Lisle & Bolnick, 2020).

This study seeks to understand whether pelvic spine reduction correlates with ecological and morphological changes in brook sticklebacks. We tested the hypothesis that pelvic spine polymorphism is associated with diet differences in brook sticklebacks, and we predicted that pelvic-reduced fish eat a more benthic diet, similar to threespine sticklebacks. To address our hypothesis, we used metabarcoding of gut contents to compare diet among spined and pelvic-reduced brook sticklebacks in two polymorphic populations. Because gill raker morphology differs between ecomorphs of threespine sticklebacks and other fish species (see, for example, Mee et al., 2015), we hypothesized that gill raker morphology would differ between spined and pelvic-reduced brook sticklebacks.

2. Methods: Studying Brook Sticklebacks

2.1. Sample Collection and Study Site

Brook sticklebacks were collected in Alberta, Canada, from Muir Lake (53.627659, -113.957524) and Shunda Lake (52.453899, -116.146192) in the summer of 2017 and 2019 using unbaited minnow traps. To sample fish from the littoral zone, traps were set near the shore (maximum 4–5 m from shore) at a depth of 0.5–2 m. To sample fish from the limnetic zone, traps were set at least 50 m from shore and suspended from floats at a depth of 1–2 m. All traps were retrieved within 1–12 h after being set. In each year, all brook sticklebacks were retained until we had collected 40 (in 2017) or 30 (in 2019) spined individuals. Once we captured 30 or 40 spined fish, we retained only pelvic-reduced individuals until we had achieved a balanced sample of spined and pelvic-reduced morphotypes (Table 1). Pelvic morphology was assessed by close inspection of the pelvic region of each fish with forceps, and sex was identified based on male nuptial colouration or the presence of eggs. Each fish was euthanized by immersion in an overdose mixture of lake water and eugenol. In 2019, the entire gut of each fish was removed and preserved in 95% ethanol. The body of each fish was preserved in 70% ethanol. As described in Willerth et al. (2022), we stained each fish with Alizarin red to observe calcium-containing osteocytes more clearly and confirmed pelvic morphology in the lab using magnified digital images.

Table 1: A summary of samples analyzed in this study. Numbers in parentheses indicate the samples used in the gut metabarcoding analysis.

These lakes are well-suited for brook stickleback pelvic spine polymorphism studies because they contain relatively high abundances of spined and pelvic-reduced morphs (Lowey et al., 2020; Nelson, 1977). The pelvic-reduced morphs include unspined individuals that completely lack spines or any vestige of a pelvic girdle, but also individuals with vestigial (a.k.a. intermediate) pelvic structures that lack some pelvic element (e.g. one spine missing, half the pelvic girdle missing, etc.; see Klepaker et al. (2013) for a detailed description of variation in pelvic morphology across stickleback species). In Muir Lake, approximately 40% of the brook sticklebacks are unspined, with very few (<2%) intermediate morphs (Lowey et al., 2020, and personal observation). In Shunda Lake, approximately 20% of the brook sticklebacks have some sort of pelvic reduction (unspined or intermediate morphs), and approximately 25% of these pelvic-reduced morphs are completely unspined (Lowey et al., 2020, and personal observation). Recent surveys and stocking records (Alberta Environment and Parks, 2021) suggest that salmonids, such as brown trout (Salmo trutta Linnaeus 1758) and rainbow trout [Ocorhychus mykiss (Walbaum 1792)], are present in both lakes, although we also observed other predators, such as loons [Gavia immer (Brunnich 1764)], dragonfly nymphs (Gomphidae and other unidentified families), giant water bugs [Lethocerus americanus (Leidy 1847)], and backswimmers (Notonectidae) in both lakes (Willerth et al., 2022). Both lakes contain abundant Daphniid zooplankton, and littoral macroinvertebrate sampling revealed Trichoptera larvae, amphipods, leeches, and Chironomidae larvae (Willerth et al., 2022). Sampling permits were issued by the Government of Alberta, and fish handling protocols were approved by the Animal Care Committee at Mount Royal University (Animal Care Protocol ID 101029 and 101795).

2.2. Diet Analysis

Fish guts collected from 95 fish in 2019 (Table 1) were shipped to the Canadian Centre for DNA Barcoding at the University of Guelph for membrane–based DNA extraction and amplicon sequencing (Ivanova et al., 2006). This metabarcoding approach used five PCR primer sets designed to amplify a barcode region of the cytochrome c oxidase subunit I (COI) gene in arthropods, mollusks, annelids, amphipods, and microalgae (details of COI amplification can be found at http://ccdb.ca/site/wp‐content/uploads/2016/09/CCDB_Amplification.pdf). A second round of PCR amplification added Ion Torrent sequencing adapters and a unique multiplex identifier (MID) sequence to the 5′ end of amplicons for each sample. Amplified COI fragments were then pooled and single-end sequenced on an Ion Torrent PGM sequencing platform (Thermo Fisher). The resulting sequence data were automatically de-multiplexed by the PGM Torrent Browser (Thermo Fisher). We removed primer and adapter sequences from the de-multiplexed reads, and discarded reads that consisted of only primer or adapter sequences using CUTADAPT v2.3 (Martin, 2011). We used DADA2 (Callahan et al., 2016) in R v4.2.2 (R Core Team, 2022) to trim and filter reads, estimate error rates for the filtered and trimmed amplicon dataset, combine identical reads into unique sequences, and infer biologically meaningful (as opposed to spurious) amplified sequence variants (ASVs). In the filtering and trimming step, we removed reads containing any unknown basepairs, truncated reads to 200 bp, discarded reads <200 bp, truncated reads at the first instance of a quality score <2, discarded reads that match the phiX genome, and discarded reads with greater than 2 “expected errors” after truncation (where the number of “expected errors” is calculated from the per-basepair quality scores for each read). In the inference of ASVs, because Ion Torrent sequencing tends to make homopolymer errors and produces higher numbers of indels (relative to Illumina sequencing), we applied a homopolymer gap penalty of −1 and applied a more restrictive threshold for the number of insertions in one sequence relative to another (BAND_SIZE = 32). We then created an ASV abundance table (with the abundance of each ASV for each fish) and removed chimeric sequences using the “consensus” method. We exported the list of ASVs across all samples to a fasta file and used the JAMP pipeline in BOLDigger (Buchner & Leese, 2020) to identify the best-fitting taxonomic identification for each ASV.

We used PHYLOSEQ (McMurdie & Holmes, 2013) to create a data object in R that allowed us to associate the abundance of each taxon identified in each fish’s gut with other sample data (i.e. lake, sex, size, and pelvic phenotype). We removed any taxa identified as brook sticklebacks or bacteria (the latter being more indicative of gut microbiome than diet). Because many taxa were not identified to genus or species, we agglomerated taxa at the rank of family, while also removing taxa that were not identified to the level of family, using the tax_glom function in phyloseq. We analyzed individual-level variation in gut contents by performing an ordination using non-metric multidimensional scaling (NMDS) and Bray-Curtis distances based on family-level taxonomic assignment.

To test the hypothesis that brook stickleback pelvic spine reduction is associated with a difference in diet, we used two analyses. First, to investigate an overall difference in diet, we conducted a PERMANOVA using the adonis2 function in the vegan package (Oksanen et al., 2022) with 10,000 permutations to evaluate whether centroids of the Bray-Curtis distances (calculated using phyloseq, as described above) differed among groups. Second, to investigate differences between groups for specific diet items, we used DESeq2 (Love et al., 2014) to conduct Wald tests of the Log2 fold-differences in the normalized abundances of each eukaryote family in the diets of brook sticklebacks (normalized by individual fish) between lakes and between pelvic phenotypes (with p-values corrected for multiple testing with the Benjamini-Hochberg method; Benjamini & Hochberg, 1995). Because we had previously found a significant difference in length between fish collected in near-shore “littoral” traps and fish collected in “limnetic” traps suspended >50 m from shore (Willerth et al., 2022), which would lead to collinearity in statistical models investigating the effects of fish length and trap location on diet, we looked for correlations among predictor variables using generalized linear models with a Gaussian error distribution (for length—a continuous numerical variable) or a binomial error distribution (for factor variables—pelvic phenotype, sex, and trap location). Given the results of these tests (details below in the Section 3) showing a significant correlation between length and habitat (but between no other predictor variables), we used only one of either habitat or length in our analyses. In the comparison between lakes (for the PERMANOVA and Wald tests), we specified a model that included lake, pelvic morphology, sex, and scaled total length (i.e. scaled to mean = 0 and variance = 1), but not habitat because no fish were captured in “limnetic” traps in Shunda Lake. Because we expect diet might be affected by allometry (Aguirre et al., 2008; Reimchen et al., 2016), we also included interactions with length in all our diet models. Because of concerns with over-fitting given our sample size, we did not include other interaction terms in our diet models. For the Wald tests, we subsequently conducted separate analyses for each lake, wherein we specified models (one per lake) that included pelvic phenotype (spined or pelvic-reduced), sex, scaled total length (or habitat in the case of Muir Lake), and two-way interactions with length.

2.3. Gill Rakers

Fish that had been bleached and stained with alizarin red for morphological analyses (Willerth et al., 2022) were used to investigate variation in gill raker morphology associated with different pelvic phenotypes. We removed the first branchial arches from both sides of each fish with forceps and scissors, cleared away all the gill filaments, and photographed the gill rakers in a plastic weighing boat using a Nikon SMZ-75T Greenough-type stereo microscope with a mounted Imaging Source camera and NIS-Elements D Software (©2021 Nikon Corporation). Following Kaeuffer et al. (2012), we counted the number of gill rakers on the left and right gill arches, and we measured the length of the second through fourth gill raker from the epibranchial-ceratobranchial joint on the ceratobranchial, which were consistently the three longest gill rakers on each arch. We used generalized linear models in R with a Gaussian error distribution to model the effects of fish length, sex, year of capture, lake, and pelvic phenotype on variation in gill raker length or number. As described above, we used length as a proxy for habitat to avoid collinearity between the length and habitat variables. We set contrasts among factors using the contr.sum function, and we used the CAR package (Fox & Weisberg, 2019) to generate an ANOVA table with type III sums of squares. Given the larger sample size for our gill raker dataset and the reasonable possibility that the effect of pelvic reduction on gill raker morphology differs between sexes, between lakes, or between years, we included all two-way interactions in our initial models, and then used a reverse-stepwise approach for model selection, removing any interaction terms that were not significant (p > .05).

3. Results: Diet and Gill Raker Analysis

A list of all taxa identified via metabarcoding in Muir Lake and Shunda Lake is available in the supplementary material (Tables S1 and S2). At the population level, brook sticklebacks in Muir Lake and Shunda Lake had taxonomically diverse diets. We identified 13 Eukaryote phyla among the Amoebozoa, Archaeplastida, Chromalveolata, Fungi, and Metazoa. The Metazoa was the most diverse taxon in both lakes where we identified six different phyla, the most abundant and diverse of which were the arthropods. We identified six arthropod classes: Arachnida, Branchiopoda, Copepoda, Insecta, Malacostraca, and Ostracoda. Insecta was the most abundant and diverse class (among all phyla) in the diets of brook sticklebacks in both lakes, with six orders and 16 families. We focused on family-level classifications because many taxa were not resolved to genus or species. Prior to statistically evaluating differences in diet between lakes and between pelvic phenotypes, we evaluated correlations among predictor variables. Using a type II sums of squares ANOVA (i.e. evaluating a main effect conditional on all other main effects without significant interaction effects), we found a significant effect of habitat (limnetic or littoral) on fish length in Muir Lake (df = 1, Wald’s Χ 2 = 4.7997, p = .028), and a significant effect of fish length on habitat in Muir Lake (df = 1, Wald’s Χ 2 = 4.8788, p = .027). No other predictor variables were correlated.

Individual-level differences in diet among brook sticklebacks obscured group-level differences in diet between lakes or between pelvic morphologies in an ordination plot (Figure 1). A PERMANOVA test of overall differences in diet did reveal that we can reject the null hypothesis that pelvic spine morphs have the same centroid, indicating that diet does differ significantly between pelvic spine morphs (Table 2). Wald tests of the Log2 fold-differences in the abundances of each eukaryote family in the diets of brook sticklebacks between lakes suggested that there were significant differences in the diets of brook sticklebacks in Muir Lake compared to Shunda Lake attributable to several taxa (Figure 2). Several families found in Muir Lake were not found in Shunda Lake and vice versa: Unionicolidae, Culicidae, Muscidae, Chaoboridae, Hydrozetidae, Ascarididae, Succineidae, Planorbidae, Crambidae, Candonidae, Naididae, Gelechiidae, Chlorellaceae, Scirtidae, and Cyclopidae were not found in Shunda Lake, whereas Simuliidae, Baetidae, Tephritidae, Sellaphoraceae, Saprolegniaceae, Leptolegniaceae, Spionidae, Cyprididae, and Bucculatricidae were not found in Muir Lake. The diet differences between lakes were not driven entirely by the absence of taxa in one lake—several taxa that were present in both lakes (e.g. Caenidae) were consumed in greater abundance in one lake relative to the other (Figure 2).

Figure 1: Ordination plot of Non-metric Multidimensional Scaling (NMDS) axes based on sample-wise Bray-Curtis distances in brook stickleback diet (with diet items classified to the level of family). Ellipses show the 95% probability interval for each group, assuming a normal distribution of points. The larger points indicate the center of each ellipse.

Table 2: Results of a PERMANOVA test of diet differences.

Figure 2: Log2 fold change between Shunda Lake and Muir Lake in the abundance of eukaryotic families in the diets of brook sticklebacks.

To further address the hypothesis that brook stickleback pelvic spine reduction (or loss) is associated with a difference in diet, we evaluated the differential abundance of each taxon identified in the stomachs of spined and pelvic-reduced brook sticklebacks in Muir Lake and Shunda Lake. When fish length was included as a factor in analyses of diets in Muir Lake, Daphniidae (daphnia, or water fleas) were significantly more abundant in the diets of spined fish, whereas Caenidae (squaregill mayflies), Gammaridae (amphipods commonly known as “scuds”), and Hyalelidae (amphipods) were significantly more abundant in the diets of pelvic-reduced fish (Figure 3). When habitat (instead of length) was included as a factor in analyses of diets in Muir Lake, only Caenidae were significantly more abundant in the diets of pelvic-reduced fish (Figure 4). In Shunda Lake, Lumbricidae (earthworms) were significantly more abundant in the diets of spined fish (Figure 3). Diets also differed significantly based on fish length and between sexes (Figures 3 and 4).

Figure 3: Log2 fold change in the abundance of eukaryotic families in the diets of brook sticklebacks in Muir Lake and Shunda Lake with length as a co-factor.

Figure 4: Log2 fold change in the abundance of eukaryotic families in the diets of brook sticklebacks in Muir Lake with habitat as a co-factor.

To investigate a potential change in gill raker morphology between spined and unspined brook sticklebacks, we compared gill raker morphology between individuals with spined and reduced pelvic phenotypes in Muir Lake and Shunda Lake. We conducted separate analyses where we looked for an effect on maximum gill raker length (i.e. the length of the longest gill raker in each individual), mean gill raker length, the length of a specific gill raker (e.g. the second gill raker from the epibranchial-ceratobranchial joint on the ceratobranchial), the number of gill rakers on the left or right gill arch, or the mean number of gill rakers on both arches. We also conducted analyses involving the effect of pelvic morphology with three categories (i.e. spined, vestigial, and unspined) instead of two categories (i.e. spined and pelvic-reduced). Although the effects of some factors varied in magnitude among these analyses, no analysis revealed a significant effect of pelvic morphology on gill raker size or number. For clarity, we report only the details of our analyses involving the effect of pelvic reduction (with two categories of pelvic morphology) on maximum gill raker length (Table 3, Figure 5) and the mean number of gill rakers (Table 4, Figure 6). Gill raker length and number did not differ significantly between spined and pelvic-reduced brook sticklebacks. Gill raker length and number were significantly affected by fish size, with larger fish having more and longer gill rakers. Brook sticklebacks in Muir Lake had more gill rakers and longer gill rakers than those in Shunda Lake. Male brook sticklebacks had significantly longer gill rakers than female brook sticklebacks, especially at larger sizes. The effect of fish size on gill raker length was weaker for samples collected in 2019 than for samples collected in 2017 (see Tables 3 and 4 for details of statistical results).

Table 3: ANOVA table for the generalized linear model of the effects of pelvic reduction (spined or pelvic-reduced), lake (Muir Lake or Shunda Lake), sex, sampling year (2017 or 2018), and body length (total length in mm) on maximum gill raker length in brook sticklebacks.

Figure 5: Comparison of maximum gill raker lengths between spined and pelvic-reduced brook sticklebacks in Muir Lake and Shunda Lake.

Table 4: ANOVA table for the generalized linear model of the effects of pelvic reduction (spined or pelvic-reduced), lake (Muir Lake or Shunda Lake), sex, sampling year (2017 or 2018), and body length (total length in mm) on mean number of gill rakers in brook sticklebacks.

Figure 6: Comparison of mean number of gill rakers between spined and pelvic-reduced brook sticklebacks in Muir Lake and Shunda Lake.

4. Discussion: Ecological Consequences of Spine Reduction

The factors maintaining the pelvic spine polymorphism within brook stickleback populations are unclear, it is possible, if not likely, that spined and unspined fish differ ecologically. A previous study of stable isotope signatures revealed that spined and pelvic reduced brook sticklebacks are associated with different habitat types (Willerth et al., 2022), suggesting that the pelvic spine polymorphism may have ecological correlates. In this study, we tested the hypothesis that spined and pelvic-reduced brook sticklebacks have different diets. Consistent with our hypothesis, we found significant differences in the abundance of some taxa in the diets of spined and pelvic-reduced brook sticklebacks.

Pelvic spine reduction in brook stickleback and threespine stickleback populations has evolved in parallel, likely driven by similar selective forces. The extent to which we should expect parallelism to extend beyond a particular phenotype of interest (in this case, pelvic spine reduction), resulting in parallel ecological, behavioral, and morphological consequences, is unknown. We predicted that, like threespine sticklebacks, pelvic spine reduction in brook sticklebacks should be associated with a more benthic diet. Our observations did not align with this prediction, which implies non-parallelism in the ecological correlates of spine reduction. Diet differences between spined and pelvic-reduced brook sticklebacks did not indicate a clear or consistent trend towards a benthic or limnetic diet in