Evolutionary Ecology Research, 2016, 17: 505–523

Trophic niche differentiation and phenotypic
divergence among cryptic species ofJapanese ninespine sticklebacks
Mark Ravinet1,2, Asano Ishikawa1 and Jun Kitano1
1

Division of Ecological Genetics, National Institute of Genetics, Shizuoka, Japan andCentre for Ecological and Evolutionary Synthesis, University of Oslo, Oslo, Norway

2

ABSTRACT
Background: Morphologically similar cryptic species often differ in ecological niches. Three
cryptic species of ninespine stickleback (genus Pungitius), P. tymensis and the freshwater andbrackish types of the P. pungitius–P. sinensis species complex are found in eastern Hokkaido,
Japan. All three co-occur in river catchments, inhabiting the upstream (P. tymensis), midstream(freshwater P. pungitius–P. sinensis type), and downstream (brackish P. pungitius–P. sinensistype).
Questions: Are differences in trophic ecology among these species persistent? Do the three
species differ in geometric body shape? Which ecological and morphological traits best
distinguish these three cryptic species?
Methods: We used stable isotope analysis to estimate trophic position and different source
contributions in order to quantify trophic niche differentiation among three species collectedfrom eastern Hokkaido, Japan. We characterized body shape variation using geometricmorphometrics. We conducted exchangeability analysis to test how well different traits could
distinguish all three species.
Results: Isotopic source and niche overlap estimates clearly indicated a greater contributionof marine prey items to the brackish type than to the other two species and showed little
differences between the freshwater type and P. tymensis. Body shape substantially differed
between the freshwater type and P. tymensis, which had a deeper, elongated shape, whereas thebrackish type was intermediate. Based on morphological and trophic traits, the brackish typeand P. tymensis could be clearly distinguished. In contrast, misclassification of P. tymensis as
the freshwater type based on stable isotope data and trophic morphology was high, indicating
substantial overlap in trophic niche between these species.
Keywords: diet, ecological speciation, exchangeability analysis, resource, sexual dimorphism,
stickleback.

Correspondence: J. Kitano, Division of Ecological Genetics, National Institute of Genetics, Yata 1111, Mishima,
Shizuoka 411-8540, Japan. email: jkitano@lab.nig.ac.jp
Consult the copyright statement on the inside front cover for non-commercial copying policies.
© 2016 Jun Kitano

506

Ravinet et al.

INTRODUCTION
Morphologically similar cryptic species are prevalent across diverse taxa (Bickford et al., 2007).
Because morphological traits are often linked to ecology and foraging behaviour (Schluter,2000), morphologically similar species may exploit similar ecological niches. However, recent
studies have shown that cryptic species often exploit different ecological niches (Rissler andApodaca, 2007; Ashrafi et al., 2011; Fišer et al., 2015; Scriven et al., 2015). In some cases, ecologicaldivergence may be driven by differences in physiological rather than morphological traits(Wellborn et al., 2004; Derycke et al., 2008; DeMeester et al., 2011). Therefore, such ecological information
may help distinguish otherwise phenotypically cryptic species. The extent of ecologicaldivergence between species is generally associated with their divergence time, such that moreclosely related species tend to occupy similar niches (Wiens and Graham, 2005; Losos, 2008; Pyron et al.,2015), although divergence driven by strong competition between closely related sympatricspecies can confound this pattern (Losos et al., 2003). Therefore, investigation of patterns of
ecological differentiation among cryptic species can help us to establish methods for distinguishing them with high confidence and potentially understand better the evolutionary
mechanisms underlying ecological differentiation.
Eastern Asia, in particular the Japanese archipelago, contains a high diversity ofPungitius species, many of which are cryptic (Takahashi et al., 2001; Tsuruta and Goto, 2007; Wang et al.,2015). Three species, P. tymensis and the freshwater type and brackish-water type of theP. pungitius–P. sinensis species complex occur in Hokkaido, Japan (Fig. 1A). Based on
mitochondrial data, divergence between P. tymensis and the other two ninespine sticklebacks is estimated to have occurred 2.0–4.4 million years ago, most probably as a result ofPleistocene or Late Quaternary sea-level change (Takahashi and Goto, 2001; Aldenhoven et al., 2010;
Wang et al., 2015). Because these three species are morphologically similar, it is difficult todistinguish them based only on their appearance (see figure 1a in Ishikawa et al., 2013). Confusion
about the taxonomy of the P. pungitius–P. sinensis species complex is also derived from thediscrepancy between armour plate phenotypes and genotypes at neutral markers. Ninespine
sticklebacks other than P. tymensis in Hokkaido were originally classified into P. pungitius
and P. sinensis based on armour plate morphology (Wootton, 1976). However, recent genetic
studies indicate that lateral plate morphs do not correspond to genetically distinct species(Takata et al., 1987; Ishikawa et al., 2013), so the P. pungitius–P. sinensis species complex is currently
tentatively classified into the freshwater type and the brackish-water type based on theirpreferred habitats (Takata et al., 1987). Despite the difficulty in defining species based onmorphological traits, P. tymensis and the freshwater and brackish-water types mainly inhabit
the upstream, midstream, and downstream of a river respectively, although they sometimes
occur in sympatry (Tsuruta et al., 2008; Ishikawa et al., 2013), suggesting that these three cryptic species
with relatively old divergence times may differ in ecological niches.
Previous work on these three Japanese Pungitius species has shown that they differ inneutral genetic markers, trophic morphology, salinity tolerance, and stomach contents(Ishikawa et al., 2013). Gill raker morphology is strongly divergent among all three species and isalso consistent with species differences in stomach content data. Pungitius tymensis has theshortest and fewest gill rakers and feeds primarily on benthos (i.e. terrestrial insects) andalso zooplankton (e.g. Copepoda); the brackish type has the most numerous and the longest
gill rakers and feeds mainly on planktonic crustaceans such as Mysidacea (Ishikawa et al., 2013).
The freshwater type has intermediate foraging morphology, but is benthic, feeding on items
such as Isopoda. Genetic data showed that only a few hybrids were identified even inEcology of Japanese ninespine sticklebacks

507

sympatry (Takata et al., 1987; Ishikawa et al., 2013), suggesting that there is strong reproductive
isolation among these three species. Habitat isolation, seasonal isolation, and hybrid
incompatibility may contribute to reproductive isolation (Kobayashi, 1959; Takahashi et al., 2005;
Tsuruta et al., 2008; Ishikawa et al., 2013), but their relative contributions to total isolation are
unknown.
The first aim of the present study was to characterize trophic niche differentiation andphenotypic divergence among the three cryptic ninespine stickleback species in eastern
Hokkaido, Japan (Fig. 1). Although we have previously characterized divergence instomach contents (Ishikawa et al., 2013) and stomach content analysis is useful for highresolution identification of prey items, it provides only a snapshot of trophic ecology – thatis, what an individual has ingested in the previous 24 hours (Grey, 2006). This is a particularly
important consideration when ontogenetic shifts occur in foraging behaviour over the life
history of individuals, as these will not be reflected by a dietary snapshot such as stomach
contents (Grey, 2001; Ravinet et al., 2010). Stomach content analysis may also be biased towards
overrepresentation of hard-bodied prey items and does not necessarily reflect energy
intake – that is, gut contents may include items ingested but not assimilated (Hyslop, 1980;
Jackson et al., 1987). Additionally, analyses based on stomach contents alone are often not
sufficient to resolve trophic position with high precision (Vander Zanden and Rasmussen, 1999; Post,
2002). Stable isotopes make it possible to quantify contributions from different dietary
sources and estimate trophic level and an integrated signal of dietary preferences over a

Fig. 1. (A) Map of sampling sites. (B) Landmark positions used in geometric morphometric analysis.
Landmarks are: (1) anterior extent of premaxilla, (2) posterior extent of supraoccipital, (3) anterior
insertion of first dorsal spine, (4) anterior insertion of dorsal fin, (5) dorsal insertion of caudal fin,(6) posterior extent of caudal peduncle, (7) ventral insertion of caudal fin, (8) anterior insertion ofanal fin, (9) insertion point of pelvic spine to pelvic girdle, (10) posteroventral extent of preopercular,(11) anterorventral extent of preopercular, (12) anterior extent of orbit, (13) anterodorsal insertion of
pectoral fin, (14) posterodorsal insertion of pectoral fin, (15) anterior extent of maxilla, (16) posterior
extent of maxilla.

508

Ravinet et al.

longer time period [i.e. 3–4 months with fish muscle tissue (Perga and Gerdeaux, 2005)]. Combining
both stomach content and stable isotope analysis provides a much more robust alternativefor quantifying trophic ecology than relying on a single method (Berner et al., 2008; Bolnick and Paull,
2009; Mathews et al., 2010; Ravinet et al., 2014).
Geometric morphometrics is a powerful method for analysing body shape variation(Zelditch et al., 2012). In Gasterosteus, body shape is substantially associated with habitat (e.g.
Walker, 1997; Aguirre et al., 2008; Sharpe et al., 2008; Hendry et al., 2011). For example, body shape varies
along coastal rivers; downstream fish are typically more slender compared with a deeper
shape seen in fish from upstream (Taylor and McPhail, 1986; Ravinet et al., 2013a). A shallow body depth
may be suitable for planktivorous fish, which perform prolonged swimming to feed in open
water, while deeper body shape can improve manoeuvrability when feeding in more complex
upstream benthic habitats (Taylor and McPhail, 1986; Hendry et al., 2011). Given the fact that cryptic
ninespine stickleback species inhabit different reaches of a river, they may vary in body
shape.
The second aim of the study was to perform exchangeability analysis using these newdata (stable isotope and geometric morphometrics) as well as previously reported data
(genotypes at neutral microsatellite markers and morphological measurements) (Ishikawa et al.,2013) in order to find the traits that can best distinguish these three species. Exchangeability
analysis examines how easily individuals can be exchanged among species based oncomposite traits from multivariate phenotypic values (Hendry et al., 2013). Exchangeability
analysis differs from other methods using standard statistical tests for differences between
trait means or significant variance terms, in that it makes use of the probability distributionof individual classifications to assess exchangeability between populations or species(Hendry et al., 2013). This type of analysis can tell us the extent of phenotypic and niche overlap
among cryptic species as well as the best diagnostic composite trait useful for defining
species.
MATERIALS AND METHODS
Sample collection
Fish captured for this study have been described previously (Ishikawa et al., 2013). Briefly, fish
were collected with minnow traps or hand nets from four different watersheds in May 2011:
Shiomi River, Biwase tidepool, Kiritappu Marsh, and Oboro River (Fig. 1, Table 1). The
number of species present differed between sites with only the brackish type found at
Biwase tidepool and P. tymensis and the freshwater type occurring at Kiritappu Marsh andthe Oboro River (Table 1). All three species were present in the Shiomi River (Table 1) andthis site was sampled at three points: upstream, midstream, and estuary (see the right upper
panel of Fig. 1). Morphological cues such as lateral plate number are not reliable for species
identification because of overlap between species (Takahashi et al., 2001; Ishikawa et al., 2013; Wang et al.,2015). All fish were therefore genotyped with eight microsatellite markers (Pun212, Pun134,
Pun19, Pun68, Pun117, Stn433, Pun171, Pun78) on different linkage groups to confirm
species identification; these results are reported elsewhere (Ishikawa et al., 2013).
Benthic macroinvertebrates were collected from the Shiomi river system and Akkeshi Bay
by hand nets or from local suppliers (for further details, see Ravinet et al., 2014). Samples were taken to
represent putative prey items and also in order to represent the isotopic baseline of benthic
and pelagic food webs for trophic level calculation (Post, 2002) (see ‘Stable isotope analysis’).

Ecology of Japanese ninespine sticklebacks

509

Table 1. Sample sizes used in geometric morphometric and stable isotope analyses
Morphometrics

Stable isotopes

Brackish Freshwater P. tymensis

Brackish Freshwater P. tymensis

Site

Sex

Biwase tidepool

Male
FemaleUnknown
Total

3
12
0
15

0
0
0
0

0
0
0
0

3
12
0
15

0
0
0
0

0
0
0
0

Shiomi River

Male
FemaleUnknown
Total

3
9
0
12

2
8
2
12

5
7
0
12

11
35
34
80

2
8
12
22

7
10
19
36

Kiritappu Marsh Male
FemaleUnknown
Total

0
0
0
0

8
4
2
14

0
0
0
0

0
0
0
0

8
4
0
12

0
0
0
0

Oboro River

0
0
0
0

14
15
1
30

4
10
1
15

0
0
0
0

15
17
1
33

4
10
0
14

Male
FemaleUnknown
Total

Stable isotope analysis
We dissected dorsal muscle from each fish (n = 212 in total: n = 50 for P. tymensis; n = 67 forthe freshwater type; n = 95 for the brackish type; see Table 1 for a full breakdown of sample
sizes for each species at each site) and dried it for 48 hours at 60⬚C. Dried muscle tissue wasthen ground and weighed into micro-tin capsules for analysis. Samples were analysed forδ13C, δ15N, %C and %N on a Carlo Erba Elemental Analyser and a Thermo Finnigan Delta
Plus XL mass spectrometer (Thermo Fisher Scientific, Waltham, MA) at the Duke
Environmental Isotope Laboratory (DEVIL) at Duke University, North Carolina, USA.
Prior to analysis, fish muscle was lipid-normalized using C:N ratios (Kiljunen et al., 2006).
Since we were primarily interested in differences between species, interspecies differencesin δ13C, δ15N, and estimated trophic position were tested using general linear mixed models(GLMMs) with site as a random effect using the nlme package. However, we further investigated intraspecific differences in trophic ecology among sites using standard GLMs. Only asubset of the 212 individuals used for stable isotope analysis were sexed by gonad inspection
and measured for standard length. Therefore, to test for differences in diet between the sexes
and dietary variation with body size, we used this subset (n = 142) in GLMMs with sex andspecies as full effects, standard length as a covariate, and site as a random effect. As with ourmorphometric analyses, we first fit this full model before performing stepwise model selection with AIC using the stepAIC function from the MASS R package (Venables and Ripley, 2002).
Stable isotope analysis can also be used to quantify isotopic niche as a proxy for the
ecological niche that organisms inhabit (Layman et al., 2007; Newsome et al., 2007; Jackson et al., 2011).
However, accurate estimation of isotopic niche should ideally account for sampling andmeasurement error (Jackson et al., 2011). To account for this and to formally test differences in510

Ravinet et al.

isotopic niche overlap, we used a Bayesian method to estimate standard ellipse area (SEAB)
implemented in siber (Stable Isotope Bayesian Ellipses in R) (Jackson et al., 2011).
We additionally estimated the relative dietary source contributions of putative marine
and freshwater prey items to each of the ninespine stickleback species. To do so, we used a
Bayesian isotopic mixing model implemented in siar (Stable Isotope Analysis in R) (Parnell etal., 2010). Assuming isotopic values represent a composite signal of diet, the mixing model
approach makes use of known fractionation rates and isotopic concentrations in order toestimate the proportion of diet made by each source. We ran siar using standard trophic
fractionation rates for fish muscle (mean ‰ ± SD: δ13C, 1.63 ± 0.63; δ15N, 3.54 ± 0.74)(Minagawa and Wada, 1984; France, 1995), both of which work well for stable isotope analysis on
sticklebacks (Ravinet et al., 2013b, 2014). Similar to siber, the siar model produces posterior probability distributions of estimates that can account for error. Furthermore, Bayesian analysis
of source contributions and niche size make it possible to perform exact tests for differences
between groups by directly comparing posterior probability distributions.
Geometric morphometrics
Pictures were taken from the left sides of euthanized fishes. A total of 16 landmarks werechosen to capture overall body shape (see Fig. 1B). Digitized coordinates were placed on
photographs using tpsDig2 (Rohlf, 2006). Initial shape analysis was performed in MorphoJ(Klingenberg, 2011). We quantified shape by performing a Procrustes fit using all individuals
before correcting any incorrectly placed or swapped landmarks from any individuals exhibiting extreme deviations from the mean shape across the entire data set using the MorphoJ‘Find Outliers’ feature (http:www.flywings.org.uk/MorphoJ_guide/). Outliers are identified
by Procrustes distance deviations in the observed data set from an expected multivariate
normal distribution (http:www.flywings.org.uk/MorphoJ_guide/). Following manual
correction for erroneously placed landmarks, we performed principal component analysis(PCA) to summarize shape variation among all individuals. For all morphometric analyses,
centroid size measured from the digital images was used as a proxy for body size. Shapevariation between the sexes and species for all principal components was tested using a
multivariate analysis of covariance (MANCOVA) with centroid size as a covariate.
Following this, univariate general linear mixed models (GLMMs) were used to test PCs
individually explaining >10% of the total variance. For each GLMM, we first fit a fullmodel with sex and species as fixed factors, centroid size as a covariate, and sampling site asa random effect using the package nlme in R (Pinheiro et al., 2016). We then used the function
stepAIC from the R package MASS to perform stepwise model selection based on AIC
scores (Venables and Ripley, 2002); in short, terms were removed from the model sequentially tofind the nested model with the lowest AIC score.
Exchangeability analysis
To determine what measurements best distinguish Pungitius species, we used exchangeability analysis – that is, how easily individuals can be exchanged among species based ontrait values (Hendry et al., 2013). In their paper proposing this method, Hendry et al. (2013)
focused on conspecific exchangeability, but here we used the method to assess which trait
measurements best distinguish each of the three species studied: P. tymensis, the freshwatertype and brackish type.

Ecology of Japanese ninespine sticklebacks

511

We examined the extent of exchangeability between each of the three species based on
five categories of measured values: diet (using δ13C and δ15N values), defence traits (firstdorsal spine length, second dorsal spine length, the most posterior dorsal spine length, anal
spine length, pelvic spine length, and pelvic girdle length), trophic morphology (gill raker
length, number of gill rakers on upper gill arch, and number of gill rakers on lower gillarch), body shape (Procrustes coordinates), and neutral genetic markers (Pun212, Pun134,
Pun19, Pun68, Pun117, Stn433, Pun171, and Pun78). Morphological measurements andneutral microsatellite genotypes have been reported previously (Ishikawa et al., 2013), and othertraits were taken from the current study (see ‘Results’).
For each trait category, we performed discriminant analysis on principal components(DAPC) (Jombart et al. 2010). This straightforward method first requires a principal component
analysis and then a discriminant function analysis on principal components in order toclassify individuals back to the species level. For diet, defence, trophic and body shape, PC
and discriminant function analyses were conducted using the prcomp and lda functions in
R v.3.2.2 (R Development Core Team 2015), while neutral microsatellite markers were analysed using
the dapc function in adegenet (Jombart et al. 2010). Following Hendry et al. (2013), we calculated
assignment probability and cross-classification probability. Assignment probability is theprobability of classification to the correct species out of all three included in the analysis.
For the cross-classification analysis, the self-classification category was omitted and theprobability of classification to the other two species was calculated. Then, both assignmentand cross-classification probabilities were divided by the random expectations (33.3% forassignment probability and 50% for cross-classification probability). Our approach madeone small departure from the original DAPC and exchangeability formulation (Jombart et al.,2010, Hendry et al., 2013). Namely, we did not substitute missing measurements with mean trait
values as some measurements such as body shape were only performed on a subset ofindividuals; this also meant that the number of samples used in each category differed. To
account for this, we calculated proportions assigned to each individual trait category rather
than focusing on the absolute number of individual assignments.
RESULTS
Trophic niche differentiation estimated by stable isotope analysis
Mean δ C and δ15N values of macroinvertebrates did not differ between the Akkeshi Bay
and Biwase tidal sites (GLMMs with species as random factor: δ13C, P = 0.66; δ15N,
P = 0.44); thus values for invertebrates sampled at these locations were pooled and treatedas a marine signal. Freshwater prey showed a strong depletion of δ13C (GLM: R 2 = 0.96,
F3,34 = 6.66, P = 0.001; Fig. 2A), indicating a transition from a marine to freshwater benthic
carbon source. Similarly, δ15N values were lower in freshwater prey relative to marine
food items, suggesting a more depleted δ15N baseline (R 2 = 0.96, F3,34 = 11.17, P < 0.0001;
Fig. 2A).
Mean δ13C differed among the three Pungitius species (R 2 = 0.81, F2,206 = 363.28,
P < 0.0001; GLMMs with site as random effect and size as a covariate). The brackish typehad a less depleted mean δ13C value (15.25 ± 3.19‰; mean ± SD) than either the freshwater type or P. tymensis (−26.90 ± 2.02 and 25.52 ± 2.43‰, respectively; Tukey HSD
13
P < 0.001; Table 2, Fig. 2A). The similar, more depleted δ C signal in both P. tymensis andthe freshwater type suggested predominately freshwater foraging in these species.
13

512

Ravinet et al.

Fig. 2. Stable isotope analysis. (A) Isotopic biplot showing δ13C and δ15N variation among ninespine
stickleback species. Black diamonds indicate mean values for prey items, horizontal and vertical error
bars represent standard deviation around the mean. Boxplots showing species distributions ofestimated trophic position (B) and isotopic niche size measured as Bayesian standard ellipse area
(SEAB) (C).

Table 2. Mean stable isotope values, trophic position, and estimates of standard ellipse area (SEAB)
for three ninespine stickleback species found in the Shiomi River, Hokkaido
δ13C

δ15N

Trophic position

SEAB

Species

Mean

SD

Mean

SD

Mean

SD

Median L95 HPD U95 HPD

Brackish type
Freshwater typeP. tymensis−15.26−25.52−26.91

3.20
2.44
2.02

16.23
11.91
12.08

1.54
1.56
1.68

3.70
2.20
2.22

0.51
0.49
0.51

11.09
9.99
9.56

9.12
7.93
7.33

13.67
12.80
12.73

Note: L95 HPD = lower 95% highest posterior density, U95 HPD = upper 95% highest posterior density (HPD).

Variation in δ15N values showed a similar pattern of enrichment in the brackish type but
lower values in the freshwater type and P. tymensis (R 2 = 0.67, F2,206 = 171.98, P < 0.0001;
see Table 2, Fig. 2A), while the freshwater type and P. tymensis had very similar values. A
more enriched mean δ15N value in the brackish type (16.22 ± 1.54‰) suggests foraging at
a higher trophic level. Estimates of trophic level following baseline correction showed that
this was indeed the case (R 2 = 0.71, F2,206 = 205.23, P < 0.0001); the brackish type had amean trophic level ( ± SD) of 3.70 ± 0.50, a level above that of P. tymensis and the freshwater type (2.20 ± 0.49 and 2.22 ± 0.51, respectively; P < 0.0001; Table 2, Fig. 2B). Analysis

Ecology of Japanese ninespine sticklebacks

513

on a subset of individuals that were sexed and measured for standard length (n = 142)
13
15
indicated that δ C, δ N, and trophic position did not differ significantly between males andfemales (GLMMs with site as a random effect, size as a covariate, and sex and species as
fixed effects; P = 0.09 for the effects of sex on δ13C; sex term removed from the best-fit
model for δ15N and trophic position). Size was a significant overall positive term for all three
dietary models (δ13C: R 2 = 0.89, F1,133 = 99.31, P < 0.0001; δ15N: R 2 = 0.78, F1,135 = 67.29,
P < 0.0001; trophic position: R 2 = 0.81, F1,135 = 87.15, P < 0.0001) but was only significant
within species for δ15N and trophic position (P < 0.0001 in both cases); there were no
species × size interactions. In short, δ15N and trophic position increased with standard
length.
Although our main focus was to investigate isotopic differences between species, someintraspecific variation between sites was observed. For the freshwater type, mean δ13C
was slightly less depleted for individuals from Kiritappu (−23.97 ± 2.92‰) comparedwith the Oboro River and Shiomi River (−26.21 ± 2.94 and −25.63 ± 1.54‰, respectively;
R 2 = 0.07, F2,64 = 3.61, P = 0.033), suggesting more freshwater foraging or a more depleted
15
freshwater baseline in the latter two populations. The mean δ N value for P. tymensis was
slightly higher in the Oboro River (13.54 ± 1.01‰) compared with the Shiomi River(11.51 ± 1.54‰; R 2 = 0.29, F1,48 = 20.64, P = 3.74 × 10−5) and was reflected by a slightincrease in estimated trophic level in the Oboro River (2.66 ± 0.30 vs. 2.05 ± 0.48 in Shiomi
River).
Isotopic niche size did not differ significantly between the species based on posterior
comparisons of SEAB estimates (P > 0.05 based on posterior comparisons; Fig. 2C), withthe distribution of difference between pairwise estimates always encompassing zero. There
was no niche overlap at all between the brackish type and either of the other two species.
However, niche overlap was high between the freshwater type and P. tymensis, with a mean
overlap of 52% for the two species (51% and 53% of total isotopic niche area of eachspecies, respectively).
Estimated dietary source contributions suggested a predominately putative marine dietfor the brackish type (Fig. 3): putative marine dietary source contribution was greater forthe brackish type than the other two species (P < 0.001 in both cases, based on posterior
comparisons). Estimated dietary source contributions of P. tymensis were almost entirely
freshwater (Fig. 3). Despite no obvious difference in diet between the freshwater type andP. tymensis using standard stable isotope analysis (see previous section), dietary source
contributions indicated a greater proportion of putative marine dietary source in the freshwater type than in P. tymensis (P < 0.002).
Body shape variation revealed by geometric morphometrics
Analysis of centroid size revealed that the freshwater and brackish types showed larger bodysize than P. tymensis (GLMM with site as a random effect: R 2 = 0.39, F2,96 = 7.62,
P = 0.0008 for species term; Fig. 4); centroid size also differed consistently between the sexes
in all three species, with females being larger than males (same model, F1,96 = 8.47,
P = 0.0045; Fig. 4). Multivariate analysis on all shape principal components revealed
significant and large effects of species, body size, and sex on shape variation (see Table 3).
Three major axes accounted for >10% of shape variation each, cumulatively explaining
57.2% of the total variance explained (PC1 = 22.9%, PC2 = 19.2%, PC3 = 14.5%; Fig. 5)
and we therefore focus on these primarily here. Variation along PC1 contained some

514

Ravinet et al.

Fig. 3. Posterior probability distributions of mean dietary source contributions to each ninespine
stickleback species estimated using siar.
Table 3. MANCOVA on all body shape principal components with centroid size as a covariate
Factord.f.

Wilks’ λFPη2Body size
Sex
SpeciesSize × Sex
Size × Species
Sex × Species
Size × Species × Sex

1,28
2,56
2,56
2,56
2,56
3,84
2,56

0.241
0.164
0.008
0.387
0.417
0.139
0.448

7.62
3.56
24.19
1.48
1.33
2.27
1.20

<0.0001<0.0001<0.0001<0.05
0.09<0.0001
0.2

0.76
0.99
0.58
0.55
0.84
0.61
0.86

Note: η2 represents the partial variance explained by each factor and is calculated as 1 – Wilks’ λ. Note that for
multivariate models, partial variance explained can sum to greater than 1 and is intended as a relative means ofcomparing factors.

specimen-bending artefacts but also represented an increase in depth in the caudal
peduncle, elongation of the snout, and elongation of the premaxilla, as well as a slight
posterior shift in the position of the pectoral fin (lower right panel in Fig. 5). PC2 clearly
captured an overall deepening of the body and a slight elongation of both the snout
and caudal peduncle (lower right panel in Fig. 5); PC3 also showed evidence of posterior
deepening of the body and the caudal peduncle (lower right panel in Fig. 5).
GLMMs with site as a random effect indicated that centroid size contributed to shapevariation on all three axes and significant shape differences occurred between species forPC1 and PC3 but not PC2 (Table 4). Nonetheless, there was a significant size × species
interaction for PC2, indicating that shape variation with size for this axis differed between
species (Table 4). Significant sex effects were present only for PC1 and PC2, although asignificant sex × species interaction was present for PC3; the sex × species interaction was

Ecology of Japanese ninespine sticklebacks

515

Table 4. General linear mixed models for differences in bodymorphology summariszd by three major principal components (PC)PC axisFactor

d.f.

F

P

PC1

Size
Sex
SpeciesSize × Species

1,93
1,93
1,93
1,93

5.64
28.10
13.81
2.30

0.020
<0.0001<0.0001
0.106

PC2

Size
Sex
SpeciesSize × Sex
Size × Species
Sex × Species

1,90
1,90
1,90
1,90
1,90
1,90

11.94
13.75
2.06
1.32
4.27
2.93

<0.001<0.001
0.133
0.253
0.017
0.059

PC3

Size
Sex
SpeciesSize × Species
Sex × Species

1,91
1,91
1,91
1,91
1,91

7.80
1.24
79.98
1.76
3.24

0.006
0.269
<0.0001
0.178
0.044

Note: Terms differ because in all cases a full model was fitted (i.e. size as
covariate, sex and species as fixed factors, all interactions permitted) and thenthe final model was chosen using stepwise analysis of AIC.

close to significant for PC2. Taken together, these data suggest that morphological differences between the sexes are not consistent among species (Fig. 5, Table 4). PC2 appears
to reflect primarily sex differences with an increased body depth in males versus females(Fig. 5). Focusing on PC1 and PC3 as the main axes describing shape variation between thespecies (Fig. 5), there is considerable morphological divergence between P. tymensis and the
freshwater type: P. tymensis individuals are characterized by an elongated (larger PC1) andtypically deeper body shape (larger PC3), whereas the freshwater type has a shorter (smaller
PC1) and thinner body shape (smaller PC3) (Fig. 5). In contrast, the brackish type wasintermediate in body shape. In all species, males tend to have larger PC1 (more compressedbody and larger head) and larger PC3 (deeper body shape) than females.
Species classification based on exchangeability analysis
Classification to the correct species (i.e. self-classification) was typically high (i.e. 2–3 times
the random expectation) for all three species for all five trait categories (Fig. 6, upper panel),
except that P. tymensis had the lowest rates of correct classification for both diet and trophic
morphology (1.4 and 1.7 times the random expectation, respectively). Pungitius tymensis
showed a high rate of misclassification to the freshwater type for these categories (1.25–1.5
times above the random expectation; see Fig. 6, upper panel), suggesting that foraging
ecology offers a poor means of distinguishing these species.
For cross-classification analyses, the freshwater type showed a high rate of crossclassification to the brackish type for all categories (2–3 times the expected amount), except
for diet and trophic morphology, for which cross-classification was highest to P. tymensis

Fig. 4.

Fig. 5.

Ecology of Japanese ninespine sticklebacks

517(2.1–2.7 times the expectation; see Fig. 6, lower panel). Both brackish type and P. tymensis
have high rates of cross-classification to the freshwater type for all traits, except thatP. tymensis was cross-classified to not only the freshwater type but also the brackish typebased on body shape.
DISCUSSION
Trophic niche differentiation amongst the three species
Our results demonstrate that the three Pungitius cryptic species show trophic niche
divergence, but that there is a considerable niche overlap between the freshwater type andP. tymensis. Overall, our stable isotope data are consistent with our previous stomach
content data that these three cryptic ninespine stickleback species exploit divergent preyresources along a marine–freshwater gradient (Ishikawa et al., 2013). Source contributions
suggested a greater proportion of putative marine prey in the brackish type than the othertwo species. Standard stable isotope analysis and Bayesian estimates of niche overlap suggested a similar benthic, freshwater diet for both the freshwater type and P. tymensis.
However, the results of our mixing model analysis suggest subtle source contributiondifferences between the freshwater type and P. tymensis – specifically, a greater freshwater
dietary contribution for the latter. Previous stomach content data also suggest some subtle
niche differentiation between these two species; P. tymensis fed on benthic larvae, terrestrial
insects, and some zooplankton, whereas the freshwater type fed primarily on benthic
isopods in the midstream (Ishikawa et al., 2013).
Stable isotope analysis is complementary to stomach content analysis. The long-term
integration of diet is a potentially useful property of stable isotope analysis: for δ13C andδ15N, fish muscle tissue typically reflects dietary assimilation over the previous 3–4 months(Hesslein et al., 1993; Grey, 2000). This may be extended even further if sampling is conducted after
overwintering, as growth typically ceases and nutrients are used to sustain basal metabolic
processes (Perga and Gerdeaux, 2005). As the fish included in this study were collected in May,
stable isotope signals likely mainly reflect growth from the previous summer. While
it is possible that our isotopic data do not fully reflect trophic ecology after overwintering,the isotopic signal observed here might represent dietary preference over a considerable
proportion of the life history of these species – an important consideration given that manyfish species exhibit ontogenetic dietary shifts (Grey, 2001; Ravinet et al., 2010), including sticklebacks(Sillett and Foster, 2000; Ravinet et al., 2014). Indeed, a positive relationship between standard lengthand trophic position suggests some dietary shifts occur with growth in all three species.
Furthermore, stable isotope analysis can give us information about trophic position. Higher
δ15N in the brackish type indicates foraging at a higher trophic level than either of the otherFig. 4. Body size (measured as centroid size) differences among ninespine stickleback species.
Fig. 5. Principal component (PC) analysis of size-corrected Procrustes residuals produced by geometric morphometrics results in three major axes of divergence (PC1, PC2, and PC3). Boxplots show
variation among species and sexes – full-colour boxes to left are female, opaque boxes to right are
male. The lower right panel indicates outlines showing how body shape varies along each axis. Blackoutlines represent shapes with a lower value (−0.5); red outlines represent shapes with a higher value(+0.5); grey dots indicate the positions of consensus body shape.

518

Ravinet et al.

Fig. 6. Exchangeability analysis. Ratio of proportion assigned to each category to the expected
proportion assigned by chance is shown. Upper panels show assignment to correct species (i.e. selfclassification) or misclassification to one of the two alternative species; here, random expectation is0.33. Lower panels show cross-classification where individuals are assigned to the next highestprobability following correct classification; here, random expectation is 0.5. Individual panels
represent species assigned to; bar colours indicate ‘true’ species.

two species. The difference in trophic levels likely reflects the presence of cryptic trophic
levels within the marine planktonic food-web (Harrod and Grey, 2006; Santer et al., 2006).
It should be noted that our Bayesian mixing model analysis suggests some contributionof putatively marine prey resources to P. tymensis (median 21%, 14–27% highest posterior
density). We argue that this is mostly due to the fact that the putatively marine prey resource
we used in this analysis does not represent a complete marine signal – this would also
explain why the majority of brackish individuals are more δ13C enriched than their putative
marine prey (Fig. 2A).

Body shape divergence among the three species
Our geometric morphometric analysis showed considerable body shape variation amongninespine stickleback species. Compared with the patterns of body shape variation in threespine sticklebacks, we found a somewhat mixed pattern of convergence in body shapeEcology of Japanese ninespine sticklebacks

519

between Pungitius and Gasterosteus. Benthic foraging P. tymensis has an increased body
depth and a compressed caudal peduncle compared with the planktivorous brackish type,whereas the similarly benthic freshwater type has a much more shallow body depth and amore slender, elongated caudal peduncle than the brackish type. In threespine sticklebacks,
a shallower, more slender body shape is associated with open-water, planktivorous foraging,whereas a deeper bodied shape is associated with benthic feeding (Schluter, 1995; Berner et al., 2008;
Willacker et al., 2010). Increased body depth suggests more manoeuvrability in complex environments, which is likely linked to greater foraging efficiency on benthic prey (Walker, 1997;
Hendry et al., 2011). Our data therefore indicate that there are some differences in patterns ofmorphological variation between threespine and ninespine sticklebacks.
In addition to morphological differences among species, we also found significant sex
differences in body shape within species. In all cases, males tended to have larger heads anddeeper body shapes, while females tended to have longer bodies. This pattern of sexual
dimorphism is similar to that observed in threespine sticklebacks (Kitano et al., 2007, 2012; Spoljaric
and Reimchen, 2008; Aguirre and Akinpelu, 2010). However, body shape sexual dimorphism appears to
be more pronounced in threespine sticklebacks than ninespine sticklebacks: there is very
little overlap between male and female body shape in the former (Kitano et al., 2007, 2012; Aguirre
et al., 2008; Cooper et al., 2011), whereas substantial overlap between the sexes is observed in thelatter. Although adaptation to divergent resources between the sexes is one of the possible
driving forces of the evolution of sexual dimorphism in threespine sticklebacks (Bolnick andLau, 2008; Reimchen et al., 2008; Cooper et al., 2011), we did not find any significant differences instable isotopes between the sexes in the ninespine stickleback species analysed. These results
suggest that the patterns and driving forces of sexual dimorphism may differ between
threespine and ninespine stickleback lineages despite the fact that directions of shapedifferences are often shared.
Evolutionary inferences from exchangeability analysis
Our exchangeability analysis indicated that multiple trait measurements enable us to distinguish these three cryptic species with a high degree of confidence. However, both diet andtrophic morphology proved to be relatively poor distinguishing traits between the freshwater type and P. tymensis. Despite the niche similarity between the freshwater type and P.
tymensis, no hybrids were found where these two species are sympatric (Ishikawa et al., 2013),
suggesting that trophic divergence is unlikely to play an important role in maintaining
reproductive isolation. However, there is habitat isolation, with P. tymensis inhabiting more
upstream habitats than the freshwater type (Tsuruta et al., 2008; Ishikawa et al., 2013). Divergence insome other physiological traits, including salinity tolerance (Ishikawa et al., 2013), may contributeto habitat isolation.
Reflecting the differences in preferred habitats (the most upstream for P. tymensis andthe most downstream for the brackish type), misclassification of the brackish type andP. tymensis into each other was generally low except for body shape. The brackish type
and the freshwater type could also be confidently classified based on the majority of traits,
although some brackish individuals were misclassified into the freshwater type based on
trophic morphology and defence morphology. Because the frequencies of hybrids were very
low between the brackish type and the freshwater type (Ishikawa et al., 2013), a combination ofisolating barriers, including ecological divergence (Ishikawa et al., 2013 and this study) and hybrid
male sterility (Takahashi et al., 2005), likely contribute to total reproductive isolation. It should be520

Ravinet et al.

noted that we analysed only wild-caught fish in this study, so we cannot exclude the possibility that phenotypic divergence may be caused by phenotypic plasticity (Wund et al.,2008; Svanbäck and Schluter, 2012). Further work is necessary to characterize physiologicaldivergence and isolating barriers among these three species and also understand their
genetic basis.
Finally, the high level of accuracy exhibited by our exchangeability analysis suggests thismethod, which incorporates variation at the individual level (Hendry et al., 2013), could be used
to identify cryptic Pungitius species. Although identification of species was conducted using
neutral genetic markers in this study, our data set could be feasibly used as a training set for
further analyses, removing the need for genotyping new samples. By distinguishing cryptic
species non-lethally at low cost (i.e. using body shape and defence trait measurements only),
exchangeability analysis could be used in conservation efforts to preserve rare ninespine
stickleback populations.

ACKNOWLEDGEMENTS
We thank M. Kume, K. Yoshida, and other field assistants for help with sampling and two anonymous
reviewers and Andrew Hendry for their comments on an earlier version of the manuscript. This
research was supported by Grant-in-Aid for Scientific Research on Innovative Areas from the
Ministry of Education, Science, Sports, and Culture to J.K. (23113007 and 23113001). M.R. was
supported by the Japanese Society for the Promotion of Science on the Short-term and the Standard
Postdoctoral Fellowship for Foreign Researchers and the NIG Cooperative Research Program (2013A33). A.I. was supported by the Japanese Society for the Promotion of Science Research Fellowshipfor Young Scientists (PD).

REFERENCES
Aguirre, W.E. and Akinpelu, O. 2010. Sexual dimorphism of head morphology in three-spined
stickleback Gasterosteus aculeatus. J. Fish Biol., 77: 802–821.
Aguirre, W.E., Ellis, K.E., Kusenda, M. and Bell, M.A. 2008. Phenotypic variation and sexual
dimorphism in anadromous threespine stickleback: implications for postglacial adaptive
radiation. Biol. J. Linn. Soc., 95: 465-478.
Aldenhoven, J.T., Miller, M.A., Corneli, P.S. and Shapiro, M.D. 2010. Phylogeography of ninespine
sticklebacks (Pungitius pungitius) in North America: glacial refugia and the origins of adaptive
traits. Mol. Ecol., 19: 4061–4076.
Ashrafi, S., Beck, A., Rutishauser, M., Arlettaz, R. and Bontadina, F. 2011. Trophic niche partitioning of cryptic species of long-eared bats in Switzerland: implications for conservation. Eur.
J. Wildl. Res., 57: 843–849.
Berner, D., Adams, D.C., Grandchamp, A.C. and Hendry, A.P. 2008. Natural selection
drives patterns of lake–stream divergence in stickleback foraging morphology. J. Evol. Biol., 21:
1653–1665.
Bickford, D., Lohman, D.J., Sohdi, N.S., Ng, P.K.L., Meier, R., Winker, K. et al. 2007. Cryptic
species as a window on diversity and conservation. Trends Ecol. Evol., 22: 148–155.
Bolnick, D.I. and Lau, O.L. 2008. Predictable patterns of disruptive selection in stickleback inpostglacial lakes. Am. Nat., 172: 1–11.
Bolnick, D.I. and Paull, J.S. 2009. Morphological and dietary differences between individuals are
weakly but positively correlated within a population of threespine stickleback. Evol. Ecol. Res.,
11: 1217–1233.

Ecology of Japanese ninespine sticklebacks

521

Cooper, I.A., Gilman, R.T. and Boughman, J.W. 2011. Sexual dimorphism and speciation on two
ecological coins: patterns from nature and theoretical predictions. Evolution, 65: 2553–2571.
DeMeester, N., Derycke, S., Bonte, D. and Moens, T. 2011. Salinity effects on the coexistence of
cryptic species: a case study on marine nematods. Mar. Biol., 158: 2717–2726.
Derycke, S., Remerie, T., Backeljau, T., Vierstraete, A., Vanfleteren, J., Vincx, M. et al. 2008. Phylogeography of the Rhabditis (Pellioditis) marina species complex: evidence for long-distancedispersal, and for range expansions and restricted gene flow in the northeast Atlantic. Mol. Ecol.,17: 3306–3322.
Fišer, Ž., Altermatt, F., Zakšek, V., Knapič, T. and Fišer, C. 2015. Morphologically cryptic amphipod species are ‘ecological clones’ at regional but not at local scale: a case study of four Niphargus
species. PLoS One, 10: e0134384.
France, R.L. 1995. Carbon-13 enrichment in benthic compared to planktonic algae: foodweb
implications. Mar. Ecol. Progr. Ser., 124: 307–312.
Grey, J. 2000. Trophic fractionation and the effects of diet switch on the carbon stable isotope‘signatures’ of pelagic consumers. Verh. Int. Ver. Theor. Angew. Limnol., 27: 3187–3191.
Grey, J. 2001. Ontogeny and dietary specialization in brown trout (Salmo trutta L.) from Loch
Ness, Scotland, examined using stable isotopes of carbon and nitrogen. Ecol. Freshw. Fish, 10:
168–176.
Grey, J. 2006. The use of stable isotope analyses in freshwater ecology: current awareness. Polish
J. Ecol., 54: 563–584.
Harrod, C. and Grey, J. 2006. Isotopic variation complicates analysis of trophic relations within
the fish community of Plußsee: a small, deep, stratifying lake. Arch. für Hydrobiol., 167: 281–299.
Hendry, A.P., Hudson, K., Walker, J.A., Räsänen, K. and Chapman, L. 2011. Genetic divergencein morphology–performance mapping between Misty Lake and inlet stickleback. J. Evol. Biol.,
24: 23–35.
Hendry, A.P., Kaeuffer, R., Crispo, E., Peichel, C.L. and Bolnick, D.I. 2013. Evolutionary inferencesfrom the analysis of exchangeability. Evolution, 67: 3429–3441.
Hesslein, R.H., Hallard, K.A. and Ramlal, P. 1993. Replacement of sulfur, carbon and nitrogenin tissue of growing broad whitefish (Coregonus nasus) in response to change in diet traced by
δ34S, δ13C and δ15N. Can. J. Fish. Aquat. Sci., 50: 2071–2075.
Hyslop, E.J. 1980. Stomach contents analysis – a review of methods and their application. J. Fish
Biol., 17: 411–429.
Ishikawa, A., Takeuchi, N., Kusakabe, M., Kume, M., Mori, S., Takahashi, H. et al. 2013. Speciationin ninespine stickleback: reproductive isolation and phenotypic divergence among cryptic speciesof Japanese ninespine stickleback. J. Evol. Biol., 26: 1417–1430.
Jackson, A.L., Inger, R., Parnell, A. and Bearhop, S. 2011. Comparing isotopic niche widths among
and within communities: SIBER – Stable Isotope Bayesian Ellipses in R. J. Anim. Ecol., 80:
595–602.
Jackson, S., Duffy, D.C. & Jenkins, G. 1987. Gastric digestion in marine vertebrate predators: in vitro
standards. Funct. Ecol., 1: 287–291.
Jombart, T., Devillard, S. and Balloux, F. 2010. Discriminant analysis of principal components: a
new method for the analysis of genetically structured populations. BMC Genet., 11: 94.
Kiljunen, M., Grey, J., Sinisalo, T., Harrod, C., Immonen, H. and Jones, R.I. 2006. A revised model
for lipid-normalizing δ13C values from aquatic organisms, with implications for isotope mixing
models. J. Appl. Ecol., 43: 1213–1222.
Kitano, J., Mori, S. and Peichel, C.L. 2007. Sexual dimorphism in the external morphology of thethreespine stickleback (Gasterosteus aculeatus). Copeia, 2007: 336–349.
Kitano, J., Mori, S. and Peichel, C.L. 2012. Reduction of sexual dimorphism in stream-resident
forms of three-spined stickleback Gasterosteus aculeatus. J. Fish Biol., 80: 131–146.
Klingenberg, C.P. 2011. MORPHOJ: an integrated software package for geometric morphometrics.
Mol. Ecol. Resour., 11: 353–357.

522

Ravinet et al.

Kobayashi, H. 1959. Cross-experiments with three species of stickleback, Pungitius pungitius (L.),
Pungitius tymensis (Nikolsk), and Pungitius sinensis (Guichenot), with special reference to theirsystematic relationship. J. Hokkaido Gakugei Univ. B, 10: 363–384.
Layman, C.A., Arrington, D.A., Montaña, C.G. and Post, D.M. 2007. Can stable isotope ratios
provide for community-wide measures of trophic structure? Ecology, 88: 42–48.
Losos, J.B. 2008. Phylogenetic niche conservatism, phylogenetic signal and the relationship between
phylogenetic relatedness and ecological similarity among species. Ecol. Lett., 11: 995–1007.
Losos, J.B., Leal, M., Glor, R.E., De Queiroz, K., Hertz, P.E., Rodríguez Schettino, L. et al. 2003.
Niche lability in the evolution of a Caribbean lizard community. Nature, 424: 542–545.
Mathews, B., Marchinko, K.B., Bolnick, D.I. and Mazumder, A. 2010. Specialization of trophic
position and habitat use by sticklebacks in an adaptive radiation. Ecology, 91: 1025–1034.
Minagawa, M. and Wada, E. 1984. Stepwise enrichment of 15N along food chains: further evidence
and the relation between δ15N and animal age. Geochim. Cosmochim. Acta, 48: 1135–1140.
Newsome, S., Martínez Del Rio, C., Bearhop, S. and Phillips, D.L. 2007. A niche for isotopic
ecology. Front. Ecol. Environ., 5: 429–436.
Parnell, A.C., Inger, R., Bearhop, S. and Jackson, A.L. 2010. Source partitioning using stable
isotopes: coping with too much variation. PLoS One, 5: e9672.
Perga, M.E. and Gerdeaux, D. 2005. ‘Are fish what they eat’ all year round? Oecologia, 144:
598–606.
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. and R Development Core Team. 2016. nlme: Linear
and Nonlinear Mixed Effects Models. R package v.3.1-124 [http:CRAN.R-project.org/
package = nlme].
Post, D.M. 2002. Using stable isotopes to estimate trophic position: models, methods, andassumptions. Ecology, 83: 703–718.
Pyron, R.A., Costa, G.C., Patten, M.A. and Burbrink, F.T. 2015. Phylogenetic niche conservatism
and the evolutionary basis of ecological speciation. Biol. Rev., 90: 1248–1262.
Ravinet, M., Syväranta, J., Jones, R.I. and Grey, J. 2010. A trophic pathway from biogenic methane
supports fish biomass in a temperate lake ecosystem. Oikos, 119: 409–416.
Ravinet, M., Prodöhl, P.A. and Harrod, C. 2013a. On Irish stickleback: morphological diversification in a secondary contact zone. Evol. Ecol. Res., 15: 271–294.
Ravinet, M., Prodöhl, P.A. and Harrod, C. 2013b. Parallel and nonparallel ecological, morphological and genetic divergence in lake–stream stickleback from a single catchment. J. Evol. Biol.,
26: 186–204.
Ravinet, M., Takeuchi, N., Kume, M., Mori, S. and Kitano, J. 2014. Comparative analysis ofJapanese three-spined stickleback clades reveals the Pacific Ocean lineage has adapted to
freshwater environments while the Japan Sea has not. PLoS One, 9: e112404.
R Development Core Team. 2015. R: A Language and Environment for Statistical Computing.
Vienna, Austria: R Foundation for Statistical Computing.
Reimchen, T.E., Ingram, T. and Hansen, S.C. 2008. Assessing niche differences of sex, armour andasymmetry phenotypes using stable isotope analyses in Haida Gwaii sticklebacks. Behaviour,
145: 561–577.
Rissler, L.J. and Apodaca, J.J. 2007. Adding more ecology into species delimitation: ecological nichemodels and phylogeography help define cryptic species in the black salamander (Aneides
flavipunctatus). Syst. Biol., 56: 924–942.
Rohlf, F.J. 2006. tpsDig2 2.05. Stony Brook, NY: Department of Ecology and Evolution, State
University of New York.
Santer, B., Sommerwerk, N. and Grey, J. 2006. Food niches of cylopoid copepods in eutrophic
Plußsee determined by stable isotope analysis. Arch. für Hydrobiol., 167: 301–316.
Schluter, D. 1995. Adaptive radiation in sticklebacks: trade-offs in feeding performance and growth.
Ecology, 76: 82–90.
Schluter D. 2000. The Ecology of Adaptive Radiation. Oxford: Oxford University Press.

Ecology of Japanese ninespine sticklebacks

523

Scriven, J.J., Woodall, L.C., Tinsley, M.C., Knight, M.E., Williams, P.H., Carolan, J.C. et al. 2015.
Revealing the hidden niches of cryptic bumblebees in Great Britain: implications for conservation. Biol. Conserv., 182: 126–133.
Sharpe, D.M.T., Räsänen, K., Berner, D. and Hendry, A.P. 2008. Genetic and environmental
contributions to the morphology of lake and stream stickleback: implications for gene flow andreproductive isolation. Evol. Ecol. Res., 10: 849–866.
Sillett, K.B. and Foster, S.A. 2000. Ontogenetic niche shifts in two populations of juvenile threespine
stickleback, Gasterosteus aculeatus, that differ in pelvic spine. Oikos, 91: 468–476.
Spoljaric, M.A. and Reimchen, T.E. 2008. Habitat-dependent reduction of sexual dimorphism ingeometric body shape of Haida Gwaii threespine stickleback. Biol. J. Linn. Soc., 95: 505–516.
Svanbäck, R. and Schluter, D. 2012. Niche specialization influences adaptive phenotypic plasticity inthe threespine stickleback. Am. Nat., 180: 50–59.
Takahashi, H. and Goto, A. 2001. Evolution of East Asian ninespine sticklebacks as shown by
mitochondrial DNA control region sequences. Mol. Phylogenet. Evol., 21: 135–155.
Takahashi, H., Takata, K. and Goto, A. 2001. Phylogeography of lateral plate dimorphism in the
freshwater type of ninespine sticklebacks, genus Pungitius. Ichthyol. Res., 48: 143–154.
Takahashi, H., Nagai, T. and Goto, A. 2005. Hybrid male sterility between the fresh- and brackishwater types of ninespine stickleback Pungitius pungitius (Pisces, Gasterosteidae). Zool. Sci., 22:
35–40.
Takata, K., Goto, A. and Yamazaki, F. 1987. Biochemical identification of a brackish water typeof Pungitius pungitius, and its morphological and ecological features in Hokkaido, Japan. Jpn.
J. Ichthyol., 34: 176–183.
Taylor, E.B. and McPhail, J.D. 1986. Prolonged and burst swimming in anadromous and freshwater
threespine stickleback, Gasterosteus aculeatus. Can. J. Zool., 64: 416–420.
Tsuruta,T. and Goto, A. 2007. Resource partitioning and asymmetric competition between
sympatric freshwater and Omono types of ninespine stickleback, Pungitius pungitius. Can. J.
Zool., 85: 159–168.
Tsuruta, T., Machida, Y. and Goto, A. 2008. Nesting habitat use and partitioning of three sympatricninespine sticklebacks (genus Pungitius): implications for reproductive isolation. Environ. Biol.
Fish., 82: 143–150.
Vander Zanden, M.J. and Rasmussen, J.B. 1999. Primary consumer δ13C and δ15N and the trophic
position of aquatic consumers. Ecology, 80: 1395–1404.
Venables, W.N. and Ripley, B.D. 2002. Modern Applied Statistics with S. New York: Springer.
Walker, J.A. 1997. Ecological morphology of lacustrine threespine stickleback (Gasterosteus
aculeatus) L. (Gasterosteidae) body shape. Biol. J. Linn. Soc., 61: 3–50.
Wang, C., Shikano, T., Persat, H. and Merilä, J. 2015. Mitochondrial phylogeography and cryptic
divergence in the stickleback genus Pungitius. J. Biogeogr., 42: 2334–2348.
Wellborn, G.A. and Cothran, R.D. 2004. Phenotypic similarity and differentiation among sympatriccryptic species in a freshwater amphipod species complex. Freshw. Biol., 49: 1–13.
Wiens, J.J. and Graham, C.H. 2005. Niche conservatism: integrating evolution, ecology, andconservation biology. Annu. Rev. Ecol. Evol. Syst., 36: 519–539.
Willacker, J.J., von Hippel, F.A., Wilton, P.R. and Walton, K.M. 2010. Classification of threespine
stickleback along the benthic–limnetic axis. Biol. J. Linn. Soc., 101: 595–608.
Wootton, R.J. 1976. The Biology of Sticklebacks. London: Academic Press.
Wund, M.A., Baker, J.A., Clancy, B., Golub, J.L. and Foster, S.A. 2008. A test of the ‘flexible stem’
model of evolution: ancestral plasticity, genetic accommodation, and morphological divergencein the threespine stickleback radiation. Am. Nat., 172: 449–462
Zelditch, M.L., Swiderski, D.L. and Sheets, H.D. 2012. Geometric Morphometrics for Biologists:
A Primer. London: Academic Press.