Species delimitation can be based on consideration of several different criteria, including differentiation of ecological or functional traits. Two species of Pacific rockfish, the dark rockfish (Sebastes ciliatus) and the dusky rockfish (Sebastes variabilis), appear to represent recently divergent evolutionary lineages. We evaluate evidence for differentiation of these two species in somatic shape using geometric morphometrics in two locations in the northeast Pacific where they occur in sympatry. Somatic shape was significantly different between species, but species’ shape did not vary between the two locations. Sebastes ciliatus had an upturned, relatively smaller head, eye, and jaw, and an elongated mid-body; whereas, S. variabilis had a downturned, larger head, eye, and jaw, and a shorter mid-body. These results suggest that S. ciliatus and S. variabilis are morphometrically differentiated in a similar way in both locations. Somatic shape differentiation between these two sympatric species is similar to genus-wide patterns of somatic shape differentiation.

For this study, we collected S. ciliatus and S. variabilis via hook and line fishing in Frederick Sound near the north end of Kuiu Island, Alaska in late June of 2016, and in Icy Strait, Alaska in June of 2023 (160 km straight line distance between the two sampling locations; Figure 1). We used samples of both species collected in both locations for comparison of divergence in somatic shape. We used a total of 105 specimens; 60 S. ciliatus, 34 from Kuiu (18 males, 16 females) and 26 from Icy Strait (15 males, 11 females); 45 S. variabilis, 17 from Kuiu (8 males, 9 females) and 28 from Icy Strait (12 males, 16 females). Both species are considered sexually monomorphic and we had a relatively even mix of both sexes in each sample. Specimens for both species ranged from 30 – 47 cm in both locations.

Specimens were photographed within a few minutes after capture to ensure that photos reflected fresh, relaxed, and natural shape variation. For each individual in all samples we photographed the right lateral view for shape analysis, and we digitized 13 landmarks (Figure 2) using tpsDig. Landmarks were (1) anterior tip of the lower jaw (i.e., symphyseal knob); (2) anterior extent of the eye; (3) posterior extent of the eye; (4) posterior extent of the operculum; (5) posterior, ventral extent of the maxilla; (6) anterior insertion of the dorsal fin on body outline; (7) posterior extent of scalation at the midline on the caudal peduncle; (8) anterior insertion of anal fin; (9) anterior insertion of pelvic fin; (10) dorsal outline vertical of landmark 3; (11) dorsal outline halfway between landmark 6 and 7; (12) dorsal outline at smallest width of caudal peduncle; (13) ventral outline at smallest width of caudal peduncle. Landmarks 2, 3, 11, 12, and 13 are sliding semilandmarks. All specimens were landmarked by one researcher without respect to predictor variables. The landmarked images were then independently inspected by another researcher to confirm homologous and consistent placement of landmarks. This procedure results in reduced error in digitizing. We used tpsRelW to align specimens via a generalized Procrustes analysis (a superimposition method to remove non-shape variation via rotation, translation, and scaling). We generated shape variables as partial warps and uniform components comprising the weight matrix (i.e., W). The weight matrix is used as input for a principal component analysis and the resulting principal components, termed relative warps [29], are used as shape variables for the subsequent analysis. Like all principal components, relative warps are ordered by the amount of variation they individually explain. Typically, we use all relative warps that account for >1% of variation. In this study, we used the first 10 relative warps (that collectively explain 97.1% of shape variation) to characterize shape.

We used a multivariate linear mixed model to determine if shape differed between S. ciliatus and S. variabilis or between locations. The response variable was shape as represented by the first 10 relative warps. The predictor variables were species (S. ciliatus or S. variabilis), location (Kuiu or Icy Strait), centroid size (a multivariate measure of size), and the index variable (indexing relative warps 1-10, see explanation below). We also included all two-way and three-way interactions with the index variable and other predictor variables. Body size is a common influence on shape among many species of fishes. Although our samples varied little in size, we included centroid size as a covariate to avoid confounding shape caused by size with shape differences between species because of potential sampling differences.

A multivariate linear mixed model assumes a univariate response variable, so we vectorized the shape variables such that each row represented one response variable, but each specimen was represented by multiple rows of data. Thus, the first row represented relative warp 1 for the first specimen, the second row represented relative warp 2 for the first specimen, and so forth until all relative warps were represented in successive rows for the first individual. The same pattern was repeated for all individuals, each with 10 rows. The index variable preserved the order of the relative warps such that comparisons between groups (e.g., S. ciliatus and S. variabilis) were made by matching each relative warp to the same relative warp in each group (e.g., relative warp 1 in one species was compared to relative warp 1 in the other species). Our main goal was to determine how shape differed between species and if shape was different within species in different locations; thus, it is the two-way interaction of species and the index variable and the three-way interaction of species and location and index variable that tested our hypothesis of interest. Main effects by themselves (without the interaction with the index variable) test only for an average effect across all relative warps. Because relative warps are principal components, they have a mean of 0; and more importantly, they have an arbitrary ordination. Thus, a single individual may have a positive score on some relative warps and a negative score on other relative warps so that their mean score across all relative warps may be near 0. It was only by matching relative warps in the same order (by using the index variable as a predictor) that we could accurately test the hypothesis of interest. Specifically, the hypothesis of interest is: does shape differ between species or locations on at least some relative warps (i.e., shape variables)? This vectorization of the multivariate response variables allows parametric testing of multiple and complex interaction effects and has been used successfully to test for shape variation in a variety of systems, among populations, and species. We estimated degrees of freedom using the Kenward and Roger method. We used Proc MIXED in SAS to run this analysis (SAS version 9.4; SAS Institute Inc., Cary, NC, USA).

To visualize the differences between species we plotted mean scores by species for each of the 10 relative warps (i.e., shape variables). In addition, we calculated a divergence vector between S. ciliatus and S. variabilis across all 10 relative warps using methods from Langerhans and Langerhans and Makowicz. To create the divergence vector, we performed a principal components analysis (PCA) of the least squares means for each species derived from the multivariate linear mixed model output. We multiplied values of the first eigenvector (from the PCA) by each of the corresponding relative warp values for each individual and then summed these to create the individual divergence score. We then used the vector of divergence scores as the regressor variable in tpsRegr and the original tps landmark file as the response variable to generate thin-plate spline visualizations of species’ shapes. These visualizations represent overall shape divergence between species across all 10 relative warps simultaneously.