Site fidelity, size, and morphology may differ by tidal position for an intertidal fish, Bathygobius cocosensis (Perciformes-Gobiidae), in Eastern Australia

Sampling sites

All field work was conducted at Hastings Point, NSW, Australia (−28.36, 153.58) from September 2014 through February 2015. Sampling was conducted in a 150 m zone parallel to the shoreline, which experiences a tidal range of 1.6 m. Random numbers generated by R v3.1.1 (R Core Team, 2014) were attributed to locations determined on Google Earth v7.1. Once on site, discrete rockpools (tide pools) were chosen by walking a straight line from the random point location to the water line. Sampled rockpools had to be at least 10 cm deep with a minimum volume of 40 L as pilot field observations suggested that pools outside of these characteristics were often devoid of fish. We then assigned pools to one of two categories, using “high” and “low” shoreline position. High pools were exposed during every low tide and for most of the tidal cycle. Low pools were only exposed at spring tides and marginally during neap tides when they received intermittent wave splashes. Designated high and pools differ markedly (and consistently) in their biotic characteristics. High pools are largely devoid of macroalgae and above the barnacle zone, whereas low pools maintain extensive macroalgal cover through all seasons and are often contain Pyura stolonifera ascidians (signature species of the low intertidal/shallow subtidal for the region). Similar criteria have been used in other studies (for example, Cox et al., 2011; White, Hose & Brown, 2015) and our categories are consistent with the pool definitions used by White, Hose & Brown (2015) in more southern locations of NSW. (Photographs of tide pool locations and tidepools can be found in the Supplemental information.) As reported for designated high and low shoreline pools in other localities (e.g., Cox et al., 2011; White, Hose & Brown, 2015) fish communities are also quite distinct. There was greater fish species richness in low pools, whereas in high pools fishes consisted mainly of Bathygobius cocosensis and occasionally juvenile Microcanthus strigatus (L Malard, 2015, unpublished data).

Four pools, two high and two low, were selected for this study, based on the aforementioned criteria (shore height and algal cover and sufficient depth and volume). Selected pools were more than 20 m apart to minimize the probability of tagged fish moving between studied tide pools. Previous studies showed that B. cocosensis exhibit strong homing abilities, with 20% of individuals returning to their original pools of capture from translocated distances as great as 20 m (Griffiths, 2003a).

For each sampling event, all fish were captured by completely emptying the pool using a battery-powered bilge pump. While the pool was empty, fish were flushed from crevices using a wash bottle and captured using hand nets to ensure the full sampling of the pool (Griffiths, 2003a; Griffiths, West & Davis, 2003). We made every attempt to capture all fishes from the drained pools and only stopped once a pool was completely drained and five minutes of searching (particularly through macroalgae) without finding any more fish had elapsed. Fishes other than the target species were held in a 20 L container of seawater and were released to their original pool after it was refilled. B. cocosensis individuals were held in a separate container for the mark recapture experiment (see below) and then returned to their home pool.

Site fidelity and survival probability

On the first visit to the site (September 2014), all B. cocosensis captured from each of the four target pools were marked to investigate whether site fidelity and survival differed among high shore and low shore rockpools. Because seasonal recruitment typically ceases by July, all fish encountered in September were juveniles and not fresh recruits; thus, our analyses (tagging and morphometrics) are entirely based on juvenile and adult fish (>18 mm total length), not new recruits. Fish were anaesthetized (0.3 × 10⁻³ mg L^-1 of Aqui-S^®) and tagged on their side with Visible Implant Elastomer (VIE) subdermal tags (Northwest Marine Technologies^®, Inc.). VIE tags are biocompatible fluorescent elastomers injected sub-dermally as a liquid that solidifies after injection (Catalano et al., 2001; FitzGerald, Sheehan & Kocik, 2004). Fish were allowed to recover for at least 10 min in aerated seawater before being returned to the pool from which they were captured. This process was repeated five times, in September, October, November 2014 and January and February 2015. During each sampling event, a different colour of VIE dye was used, in order to identify each cohort and to recognize the previous time of marking for recaptured individuals. As the studied pools were over 20 m apart, migration of marked individuals between sampling pools would be unlikely, and an individual tagging scheme would be difficult to implement given the size of the fish and animal welfare concerns (i.e., an individual marking scheme would necessitate multiple injections per fish per capture event). Newly caught fish were tagged for the first time on their left side. Recaptured individuals were tagged again (on the right side with the date-specific dye) and previous tags were used as indicators of when they were first sampled.

In mark-recapture experiments such as this one, an individual can disappear because they lose their tag, because they die, or because they leave the sampling region. VIE tags do not detectably impact growth and mortality rates (FitzGerald, Sheehan & Kocik, 2004) and typically have high retention rates over 6 months (Catalano et al., 2001; Olsen & Vøllestad, 2001; Josephson et al., 2008). Of particular relevance here, a previous study of tag retention and mortality in B. cocosensis demonstrated high tag retention and low mortality of VIE tagged fish (for White & Brown, 2013: 100% retention and 0% mortality for 42 days, and for Griffiths (2002): 77 ± 19% tag retention and 7% mortality compared to 11% mortality in a non tagged, control group after 90 days). Tag loss introduces measurement error into mark–recapture analyses, but here we have no reason to expect tag loss to occur differentially between high and low shoreline pools such that this source of error would bias data interpretation. We did not directly observe any mortality events in this study, either as fish deaths during capture, VIE tagging, or during photography. Thus, although we expect some fish to have died between monthly sampling events, we assume that these mortalities are independent of handling effects. Untagged fish were detected during monthly sampling events; they were presumed to be immigrants and were subsequently tagged.

Given our repeat sampling design, the Cormack–Jolly–Seber (CJS) (Cormack, 1964; Jolly, 1965; Seber, 1965) model is appropriate for estimating survival and site fidelity (White & Burnham, 1999; Cooch & White, 2006). CJS models are based on the main assumption that marked and unmarked individuals have the same capture probability. This method also uses instances where tagged individuals reappear after being absent for one or more sampling events to separately estimate survival and emigration. The CJS method was implemented in MARK v8 (White & Burnham, 1999) and was used to estimate the survival (ϕ) and encounter probabilities (p) (Cooch & White, 2006). An encounter history, which recorded encounters of individuals with regards to sampling occasion and pool ID, was input into MARK and the resultant models estimated the survival (ϕ) probability between each recapture event and encounter probabilities (p) at each event. Both parameters can be either constant (⋅), time-dependent (t), have a group effect (g, in this case high vs. low rockpool) or have a group and time-dependent effect (g^∗t). (Notation as per MARK instruction and literature, Cooch & White (2006)). All sixteen possible models were tested, from the simplest model ϕ(.)p(.) to the most general with the most parameters ϕ(g^∗t)p(g^∗t). From these sixteen possible models, the most likely model was selected using the corrected Akaike Information Criterion (AIC_c) (Burnham & Anderson, 2002). The model with the lowest AIC_c was considered the best model and models with a ΔAIC_c >2 were not considered further. The goodness-of-fit (GOF) was tested using the median $\hat{c}$ function (Cooch & White, 2006). The median $\hat{c}$ was estimated using 30 replicates at each design point. This process was repeated three times for the selected models and the results were respectively averaged. If $1<\hat{c}<3$ then the model fits the data and the closer $\hat{c}$ is to 1, the better the fit. If $\hat{c}$ is outside of this range, the model does not fit and should be revised. Finally, estimates of the survival probability (ϕ) between each recapture event and the encounter probability (p) at each event were obtained by a weighted average of the best models and the overall ϕ and p were calculated by averaging pool estimates by habitat.

Size and morphology

Each captured B. cocosensis individual (once anaesthetised, see previous) was photographed using a NIKON Coolpix AW120 camera against a laminated graph paper background to provide uniform size scale; fish were photographed immediately after tagging, and any tagged fish that was re-captured was not photographed again, ensuring that our dataset does not contain repeated measures. Morphology was characterised by readily visible landmarks (Lele & Richtsmeier, 2001) (Fig. 1), localized on each photograph in IMAGEJ 1.48, after setting the scale using the background 1 mm grid (ImageJ, National Institutes of Health, Bethesda, Maryland, USA). Fish size was determined by total length (TL, distance between landmarks 1 and 6 in Fig. 1). In geometric morphometrics, centroid size, the square root of the sum of squared distances of each landmark to their centroid, is typically used to describe specimen size. Here, centroid size was highly correlated with TL (Pearson’s product moment correlation coefficient = 0.999, d.f. = 131, P < 0.0001), and we only analysed TL as it is readily comparable across studies.

A Generalized Procrustes Alignment (GPA), implemented in MorphoJ (Klingenberg, 2011), was used to align landmarks. The GPA removes scale (size), position and orientation (rotational) variation from the landmark coordinates (Bookstein, 1996; Adams, Rohlf & Slice, 2004; Zelditch, Swiderski & Sheets, 2012). This alignment generates collinearity among Procrustes aligned coordinates, reducing the dimensionality of the data from 2p (where p is the number of landmarks, which in our case was 7) to 2p-4 (Rohlf & Slice, 1990; Zelditch, Swiderski & Sheets, 2012). To remove this collinearity and to reduce the dimensionality of the data, a Principal Component Analysis (PCA) was conducted on the covariance matrix of the aligned landmark coordinates in MorphoJ (Klingenberg, 2011). The first seven (of 2p-4 = 10) Principal Components (PCs) explained 96.2% of the variance in body shape, and were retained for further analyses. All PCs were standardised (mean = 0 and standard deviation = 1) prior to analysis.

We also used Mahalanobis distance (De Maesschalck, Jouan-Rimbaud & Massart, 2000) to identify individuals that were multivariate outliers and outlier values were removed from the data, resulting in a total of seven missing records. Sample sizes were uneven among sites: there were a total of 41 low shore pool samples (19 and 22 from the two replicate pools) and 82 (35 and 47) high shore pool samples. If uneven sample sizes combine with unequal dispersion of samples, significance testing will be biased. The standard test for differences in variance between habitats is a Box’s M test, where a significance level of P ≤ 0.001 is generally taken as sufficient variance to bias sampling (Tabachnick & Fidell, 2007); in this instance the value indicated uneven variance (Box’s M = 48.02, F_28,23603.08 = 1.59, P = 0.025), but insufficient to substantially bias analyses. Moreover, inspection of individual traits indicated that only PC1 differed in variance between habitats. When PC1 was excluded from the analysis (PCs 2–7 included), Box’s M no longer supported heterogeneous dispersion (Box’s M = 31.24, F_21,24735.10 = 1.39, P = 0.107). Therefore to be conservative in our statistical tests, PC1 was excluded from the multivariate analysis.

Our experimental design corresponds to a mixed model nested ANOVA: (1)${\displaystyle y_{ijk}=\mu +{\mathrm{Habitat}}_i+{\mathrm{Pool}}_{j\left(i\right)}+{\epsilon }_{ijk}}$ where an individual’s trait value (y) was modelled as a function of the population mean (μ) and tide zone habitat (Habitat), which were fit as a fixed effects, and the replicate pool within each habitat (Pool) and the residual variance (ε), which were fit as a random effects. For body shape (PCs 2–6) the multivariate version of model Eq. (1) was fit, where y was a response vector of the six traits. Because there were only two replicate pools per habitat, we were unable to estimate an unconstrained covariance matrix for this random effect in the MANOVA. We therefore fit a diagonal covariance structure for this effect, estimating among pool variance for each trait, but not allowing covariance among traits. Using AIC we compared model fit of this covariance structure to more complex factor analytic covariance structures (Littell et al., 2006), which estimate reduced dimensions of variance. The diagonal covariance structure was the best-fit (smallest AIC), and therefore used in subsequent analyses.

Analyses were implemented using PROC MIXED in SAS (v. 9.4, SAS Institute Inc., 2012) under a maximum likelihood framework. To test for an effect of habitat on body shape (PCs 2–7) or size (TL), model 1 (above) was fit including or excluding the main effect of habitat, and a log-likelihood ratio test was used to determine if excluding habitat resulted in a significantly worse fit, with the degrees of freedom corresponding to the difference in the number of parameters fit (Wright, 1998). For the analysis of body shape, to determine which traits (PC2–PC6) contributed most to body shape differences between high and low tide pools, we fit habitat as a random effect, specifying a covariance structure with a single dimension as with only two levels of habitat (high and low tide pools) there can be only a single dimension of divergence between them. Specifically, we fit a first order factor analytic covariance structure (FA0(1): Littell et al., 2006) using restricted maximum likelihood. This analysis is analogous to a discriminant function analysis and identifies the trait combination that explained the variance between habitats when variation within habitats (among and within replicate pools) was appropriately taken into account.

Site fidelity and survival probability

Site fidelity differed substantially with shore height for the four pools considered in our study (Table 1). A total of 222 individuals were tagged at least once, corresponding to 158 individuals collected in high shore pools and 64 individuals collected in low shore pools over the five-month study period. Fifty-three individuals were encountered in their original pool again (tagged twice minimally), corresponding to 40 individuals from high shore pools and 13 individuals from low shore pools (Table 2). The average density of B. cocosensis was somewhat greater in low pools (low pools, A: 0.088 ± 0.059; B: 0.010 ± 0.003; high pools, C: 0.032 ± 0.007; D: 0.033 ± 0.025 B. cocosensis per liter). Of the sixteen models for encounter probability (i.e., site fidelity) and survival, the two best models were indistinguishable based on AIC_c scores (299.74 and 300.39) with the remaining models having substantially poorer fits to the data (ΔAIC_c > 2; Table S1) and were not considered further. The two best models included: (i) encounter probability (p) that differed by group (i.e., high vs. low shore pools) and survival probability (ϕ) that did not differ by groups or time (best model: ϕ(.)p(g)) and, (ii) survival probability that differed by group and encounter probability that differed by group (second best model: ϕ(g)p(g)). Both models fitted the data well, with the average $\hat{c}=1.25$ for model (i) and $\hat{c}=1.19$ for model (ii). Estimates for the survival probability ϕ and the encounter probability p were obtained by averaging these two models together to obtain the most realistic estimates (Table S2). Thus, encounter probability in the high shore habitat was more than twice that in the low shore habitat (Table 1). This habitat-level difference in site fidelity was strongly supported by the inclusion of the p(g) term in both of the best-fit models. Survival was higher for both low shore pools than for the high shore pools, with average survival probability of low shore individuals 1.2 times that of high shore individuals (Table 1). Although this result is suggestive, the lack of difference in model fit between models (i) and (ii) indicates that the present data do not support a significant difference in survival between the habitats.

Size and morphology

Fish size and morphology differed significantly between high and low shore pools. Fish from the two low shore pools (mean TL = 48.25 mm ± 1.46 SE) were significantly bigger than fish from the two high shore pools (mean TL = 40.00 mm ± 1.00) (χ² = 6.53, d.f. = 1, P = 0.011; Fig. 2). However, the size range within both habitats overlapped considerably, with fish captured from low shore pools spanning a wider range of sizes (18.6–68.3 mm) than fish from the high shore pools (25.4–59.9 mm) (Fig. 3). There appeared to be greater variance in size both between and within the two replicate low shore pools than for the replicate high shore pools (Figs. 2 and 3) but this pattern was not statistically supported (log-likelihood ratio test, LRT, comparing model with environment-specific between pool variance to a model with the same variance in both environments: χ² = 1.26, d.f. = 1, P = 0.261).

There was a significant difference in body shape between high and low shore pools (LRT of model with Habitat fit as a fixed effect vs. Habitat not fit in the model: χ² = 16.07, d.f. = 6, P = 0.013). PC6 contributed most strongly to the body shape divergence between habitats (normalised vector loading 0.78). This axis of body shape variation was associated with variation in snout length and eye diameter (Fig. 1). On average, fish from the high tide zone pools had higher PC6 scores than did low shore pool individuals (Fig. 4), corresponding to longer snouts and larger eyes for their size. It is worth noting that PC6 explained only 3.38% of the total phenotypic variation across the dataset, suggesting that differences in habitat generated relatively subtle differences in morphology. To verify that PC6 was not affected by body size, we regressed PC6 on TL (not significant: F_1,121 = 3.16, R² = 1.7%, P = 0.078) and used thin plate spline exploration to check for a non-linear relationship, of which there was no evidence. For each trait (PCs 1–10), we used the mixed model to partition variation among habitats, among pools within habitats and to the residual variation among individuals within a pool (including measurement error). Over 96% of the total variance in body shape was attributed to variation within pools.