Freshwater influence is associated with differences in bone mineral density and armour configuration in threespine stickleback (Gasterosteus aculeatus)

Adult fish (standard length ≥40 mm) representing four different habitats from the Pender Harbour region, located on the Sunshine Coast of British Columbia, were used in this study. These habitats include the marine Hospital Bay Lagoon (HBL; 49°37′53.4″N, 124°1′48.0″W; n = 24), Bargain Bay Narrows (BBN; 49°36′59.9″N, 124°1′50.3″W; n = 25), Bargain Bay Lagoon (BBL; 49°36′48.6″N, 124°1′46.9″W; n = 25), and the tidally influenced Lily Creek (LYC; 49°37′13.4″N, 124°1′43.2″W; n = 27). A total of 101 specimens were collected in the spring of 2014 using minnow traps, and were euthanized using an overdose of Eugenol (>400 mg/L; Grush et al. 2004; Sigma-Aldrich). All collection and euthanasia procedures were in accordance with the Canadian Council for Animal Care and the British Columbia Ministry of Forests, Lands, and Natural Resources Operations (University of Calgary Animal Care Permit AC13-0040; BC Fish Collection Permit SU14-92793). Specimens were immediately fixed in 10% neutral buffered formalin for 24 h, and then placed in 70% ethanol for long term storage.

All fish were scanned as described below, and subsequently stained with Alizarin Red S (Sigma-Aldrich, Oakville, Ontario) as described by Bell (1979). Lateral plates were then counted on both sides using an Olympus (Olympus Canada, Inc., Richmond Hill, Ontario) dissecting microscope. Fish were assigned to “low” (<10 plates per side), “partial” (10–30 plates per side), or “complete” (>30 plates per side) plate phenotypes, following extensive previous work (e.g., Hagen and Gilbertson 1972; Ziuganov 1983). Plates were counted for both sides twice, and averaged to assign fish to plate phenotypes. Standard length (tip of lower jaw to posterior end of hypural bone) was also measured for each specimen.

Fish were imaged using a Scanco μCT35 instrument (Scanco Medical AG, Brütisellen, Switzerland), using a standardized multiple-scanning protocol at a resolution of 20.00 μm (70 kV, 114 μA, 20.5 mm field of view and a 200 ms integration time). Raw data were thresholded and reconstructed into slice stacks for further analysis using identical parameters for every specimen.

Bone mineral density (BMD) was estimated by estimating hydroxyapatite content (mg HA/cm³) via threshold values obtained from μCT and related to the mineral equivalent using a calibration phantom (calibration was repeated weekly; Bouxsein et al. 2010). BMD was obtained in a grid pattern on the first complete lateral plate (plate five, numbered from the cranial end). This plate was sampled through the anterior margin toward the posterior, taking nine measurements across the plate at seven intervals from dorsal to ventral (Fig. 1). These measurements were repeated three times for each fish, and then a mean value for BMD was calculated for each plate interval for each fish by averaging across slices and repeated measurements. Habitat means for intervals and for whole plates were calculated from these data. BMD data were obtained using IPL software (Scanco Medical AG, Brütisellen, Switzerland).

Analysis of covariance (ANCOVA) was used to test for the presence of significant differences in mean BMD among all four habitats, with fish size (represented by the natural log of standard length) as a covariate. Significant differences in BMD were identified using Tukey post hoc testing, and a Bonferroni adjustment for multiple comparisons was applied where appropriate. Multiple analysis of covariance (MANCOVA) was used to identify differences within each plate interval, again with fish size (represented by the natural log of standard length) as a covariate, and followed by ANCOVA within each interval to identify differences among habitats. Finally, analysis of variance (ANOVA) was used to test for the presence of significant differences in mean BMD among plate intervals. All statistical analyses of bone mineral density were conducted using R version 3.4.2 (R Core Team 2017).

Isosurfaces were created using Amira version 5.4 (FEI Imaging, Thermo Fisher Scientific, Waltham, Massachusetts) from the reconstructed μCT data. Twenty bilateral landmarks were collected in Amira from each individual to represent the shape of the armour apparatus in 3D coordinate space (Fig. 2, Table 1). Landmarks were Procrustes-transformed to remove variance related to translation, rotation, and scale (Dryden and Mardia 1998). Procrustes ANOVA was used to test for the presence of significant effects of size and habitat on armour shape, with effect sizes for these analyses reported as Z-scores (Collyer et al. 2015). All significance testing was conducted using 1000-round permutation tests. All morphometric analyses were conducted using R version 3.4.2 (R Core Team 2017) and the geomorph package version 3.05 (Adams and Otárola-Castillo 2013).

All individuals from the marine BBL, BBN, and HBL habitats were completely plated, with the exception of one individual from BBL, whereas the LYC habitat contained all three plate phenotypes (Fig. 3). The “partial” individual from the BBL locality was found to have 32 plates on the right side, but only 27 on the left, giving an average of 29.5 plates. Other fish with left–right variation varied by only one or two plates. Only three individuals in the LYC habitat were found to be completely plated.

Fish from two of the three marine habitats were found to have significantly greater mean plate BMD than fish from the LYC habitat (BBL-LYC and BBN-LYC p < 0.01), and fish from the three marine habitats did not vary significantly in mean plate BMD among themselves (p > 0.05 for all comparisons). Although size was a significant source of variance, there was no interaction between size and BMD across the plate and therefore an interaction term was not included in the final model (Tables 2 and 3).

BMD varied across the plate, with the densest region located in the middle of the plate, represented by interval 4, which was significantly denser than both the superior and inferior regions of the plate (Fig. 4). BMD was relatively reduced following a similar pattern on either side of interval 4 in all four habitats, such that intervals 1 and 7, 2 and 6, and 3 and 5 were similar to each other and density tapered toward the superior and inferior aspects of the plate.

MANCOVA indicated significant differences in all seven intervals between the grouped marine and freshwater habitats (p < 0.05; Table 4). Detailed hypothesis testing using ANCOVA for each plate interval revealed significant differences among LYC fish and those from at least one marine habitat for all but the first and fourth intervals, and significantly lower BMD than all fish from all three marine habitats in the third interval (p < 0.05). Once again, although size was a significant source of variance, there was no interaction between size and BMD within plate intervals, and therefore interaction terms were not included in these models.

Procrustes ANOVA revealed that size (as represented by standard length) had a significant effect on armour shape, and that in the case of armour shape there was an interaction between size and habitat. We, therefore, included an interaction term in the final model for armour shape, and found a significant effect of habitat on armour shape after controlling for size (Table 5). Principal component analysis of an “allometry-free” data set consisting of regression residuals added back to the mean landmark configuration for each group revealed that shape differences between armour configurations differentiate fish from the LYC habitat from fish captured in the marine habitats (Fig. 5). Fish from the LYC habitat clustered toward the negative end of both PC1 and PC2 axes. This region of morphospace represents a broader armour configuration and reduction of the distance between the ascending branch of the pelvic girdle and the point of articulation of the plate with the vertebral column, and also includes a shortening of the distance between the anterior portion of the pectoral girdle and the tips of the pelvic spines. LYC fish are also more variable than those from the marine habitats, occupying a larger region of morphospace than the marine specimens.

We used traditional staining methods to confirm that in the system described here, as in most other marine–freshwater stickleback systems, stickleback inhabiting a tidally influenced habitat display a notable reduction in armour plate number (Fig. 3). It is interesting to note that in the tidally influenced habitat considered here, we found both a group of low-plated individuals (n = 13) and a group of partially-plated individuals (n = 11) with at least 10, but fewer than 30, plates. Although we do not have genetic information for these habitats, and it is certainly possible that gene flow is influencing our results, these data indicate that the “low-plated” allele is not fixed in LYC. We also note that we observed a reduction in plate size dorsoventrally among fish from the LYC habitat, which is consistent with previous work (Leinonen et al. 2012; Wiig et al. 2016) and lends support to the paedomorphosis hypothesis proposed as a mechanism for plate reduction (Bell 1981), whereby slowing of the developmental program for plates produces smaller plates in adult fish. It is possible that the production of this phenotype is mediated by cis-regulation of GDF6 (Indjeian et al. 2016), and further work is required to elucidate the gene regulatory networks that may underlie this observation. Alternatively, it is also possible that the results we observe are due to a plastic response to salinity in these habitats. Thus, further study using a common garden framework will help to clarify the relative roles of plasticity and selection in this system.

Measurement of hydroxyapatite content across the entire plate revealed both a significant effect of size on bone mineral density as well as a conserved plate bone density profile independent of fish size (Fig. 4). This profile describes a structure that is densest near the dorsoventral midpoint of the plate (interval 4) and tapers to its lowest densities at its dorsal and ventral extremes (intervals 1 and 7). This profile was recovered in fish from all four habitats, and although the relative density of the plate in LYC individuals was significantly lower at every measurement interval except the first and fourth, the shape of the profile remained unchanged, with the densest point identified at interval 4 and tapering to the least dense at intervals 1 and 7. This result suggests that although these animals’ skeletons show changes in response to novel selective pressures in their new habitat, the functional complex to which the armour plates belong may be maintained.

Plate interval 4 corresponds to the point at which the plate articulates with the vertebral column via the epiplural rib as well as with the preceding and following plate, and the plate measured in this study (plate 5) also articulates with the anterior edge of the pelvic girdle (Fig. 1). Armour plates both alone and in association with spines clearly have a protective function (Reimchen 1992; Reimchen 1994; Reimchen 2000). Although fish occupying fresh water reduce both the number and density of their plates, those that retain some plates, such as the fish described here, always retain the anterior plate complement, as has been documented elsewhere (Reimchen 1983). This retention is due, at least in part, to the role of the anterior plates in supporting the spines (Reimchen 1983), but other possible biomechanical functions for plates are less well understood. Swimming performance and plate number have been shown to be negatively correlated (Bergstrom 2002), but the biomechanical interactions between plates and the axial skeleton remain undocumented. Our results indicate the conservation of a mineralization profile that we hypothesize may be responsible for maintenance of structural integrity in the face of reduction of total plate number. Once again, further work is required to elucidate the functional role(s) of lateral plates and the ways in which these might be altered in changing environments. It is important to note that plate ossification begins from the neuromast at the center of the plate and moves toward the dorsal and ventral edges (e.g., Igarashi 1964; Igarashi 1970; Pistore 2018); an alternative hypothesis would be that the pattern observed merely reflects that plates are developmentally “older” in their central regions and have had longer to accumulate bone. Rejection of the latter hypothesis would require comparison of bone mineral density in a sample of adult stickleback of known age to determine if the pattern observed here is lost in older specimens. A balance between bone deposition and removal maintains skeletal form throughout life, and the fact that we obtain a similar plate density pattern across three representative marine habitats in a sample containing adult fish is an argument against this pattern being a feature restricted to younger specimens. Once again, further work is required to test this hypothesis explicitly. In addition, our study included only three fully plated LYC individuals, which precluded direct comparison between these and their fully plated marine counterparts. Important future work will involve comparisons of bone mineral density among different LYC plate morphs as well as between fully plated morphs of LYC and marine individuals.

We found that allometric size is a major determinant of shape differences among armour phenotypes. There remains a significant effect of habitat on shape after controlling for allometry (Fig. 5), which is consistent with previous work documenting substantial differences in cranial shape among occupants of these habitats (Jamniczky et al. 2015a, 2015b). Fish from the LYC habitat occupy a larger portion of morphospace indicating greater variability in phenotype, but generally display a more compact armour configuration that is both dorsoventrally and anteroposteriorly reduced but does not include any apparently major, possibly functionally important, configuration changes. This is again consistent with both our bone mineral density findings and with previous work indicating an overall reduction in the elaboration of the skeleton in fresh water habitats (Leinonen et al. 2012; Wiig et al. 2016), likely by paedomorphosis (Bell 1981). Future studies should focus on the functional role(s) of lateral plates as they relate to trophic morphology and swimming performance, as well as previously documented defensive functions, to better elucidate the nature of the possible functional constraints hinted at in our results. Further, it is likely that sexual dimorphism is influencing our results to some extent, given that dimorphism has been found in other studies of plate phenotype (e.g., Moodie and Reimchen 1976; Reimchen et al. 2016), as well as in many other stickleback traits (e.g., Kitano et al. 2007; Aguirre et al. 2008; Spoljaric and Reimchen 2008; Aguirre and Akinpelu 2010; Leinonen et al. 2011; McGee and Wainwright 2013; Reimchen et al. 2016), but sex data are unfortunately not available for this data set. Further work should certainly test specific hypotheses about the relationship of sex to plate phenotype in G. aculeatus.