Evolutionary Ecology Research, 2019, 20: 213–230

Diversity of gill raker number and diets among
stickleback populations in novel habitats createdby the 2011 Tōhoku earthquake and tsunami
Takuya Hosoki1,2, Seiichi Mori3, Shotaro Nishida3, Manabu Kume4,
5
1,2
Tetsuya Sumi and Jun Kitano
1

Ecological Genetics Laboratory, National Institute of Genetics, Mishima, Shizuoka, Japan,
Department of Genetics, SOKENDAI, Mishima, Shizuoka, Japan, 3Gifu-keizai University,
Ogaki, Gifu, Japan, 4Field Science Education and Research Center, Kyoto University,
Kyoto, Japan and 5Daido University, Nagoya, Aichi, Japan

2

ERRATA
21 January 2020; pp 217, 219

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

Evolutionary Ecology Research, 2019, 20: 213–230

Diversity of gill raker number and diets among
stickleback populations in novel habitats createdby the 2011 Tōhoku earthquake and tsunami
Takuya Hosoki1,2, Seiichi Mori3, Shotaro Nishida3, Manabu Kume4,
5
1,2
Tetsuya Sumi and Jun Kitano
1

Ecological Genetics Laboratory, National Institute of Genetics, Mishima, Shizuoka, Japan,
Department of Genetics, SOKENDAI, Mishima, Shizuoka, Japan, 3Gifu-keizai University,
Ogaki, Gifu, Japan, 4Field Science Education and Research Center, Kyoto University,
Kyoto, Japan and 5Daido University, Nagoya, Aichi, Japan

2

ABSTRACT
Background: Catastrophic disasters often create novel habitats and offer an opportunity toinvestigate how rapidly phenotypic diversification occurs in nature. On 11 March 2011, a large
earthquake and subsequent tsunami struck the downtown area of Otsuchi Town, Iwate Prefecture, Japan. This catastrophic event created new spring water-fed ponds that sticklebacks(genus Gasterosteus) quickly colonized. In 2014, there were 48 newly formed stickleback
habitats.
Questions: Does trophic morphology and ecology vary among stickleback populations in
newly created ponds? Is the morphological variation associated with any ecological or environmental factors in the newly formed habitats?
Methods: We sampled sticklebacks from four native and 17 novel habitats in 2014. We
counted the number of gill rakers (an important trophic trait in fish) and analysed theirstomach contents. We measured salinity, water temperature, pH, dissolved oxygen, and thedensity of benthos in the habitats and conducted principal component analysis on theseenvironmental factors. Then, we tested whether any environmental factors or environmental
principal components are associated with variation in gill raker number.
Results: The number of gill rakers varied significantly among the newly formed populations.
Fish with fewer gill rakers had a higher proportion of benthos in their stomach, suggesting thatgill raker number may contribute to the performance of benthic feeding. Furthermore, we
found that sticklebacks in habitats with higher benthos density had fewer gill rakers.
Keywords: contemporary evolution, diet, empty niche, hybridization, rapid diversification,
trophic ecology.

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

214

Hosoki et al.

INTRODUCTION
The emergence of new habitats provides organisms with opportunities to colonize
empty niches and, potentially, diversify (Simpson, 1955; Schluter, 2000; Losos, 2010). Catastrophic
events, such as volcanic eruptions, earthquakes, and tsunamis, can create new habitats
and therefore play important roles in organismic evolution (Whittaker and Fernández-Palacios, 2007;
Urabe and Nakashizuka, 2016). For example, volcanic activity has created new islands in theoceans, where birds have colonized from the nearby mainland and diversified (Lack, 1947;
Grant and Grant, 2011). Because these catastrophic events are relatively rare and unpredictable,
there are few opportunities to investigate the first stage of colonization of newly formed
habitats.
A large-magnitude earthquake is a catastrophic event that can create new habitats by
inducing landscape and hydrological changes and triggering massive tsunami waves (e.g. Wang
and Manga, 2015; Urabe and Nakashizuka, 2016; Kume et al., 2018). Charles Darwin experienced a large
earthquake in Chile during the voyage of HMS Beagle, which led him to recognize theimportance of crustal movements in organismic evolution (Darwin, 1839). Although recent
studies have highlighted many examples of evolution over contemporary time-scales (Hendry
and Kinnison, 1999; Reznick and Ghalambor, 2001; Bell and Aguirre, 2013), most have occurred in artificiallyintroduced or human-disturbed populations (Hendry and Kinnison, 1999; Reznick and Ghalambor, 2001;
Palumbi, 2002). Therefore, we know little about the effects of earthquakes on contemporary
evolution.
Recently, an excellent example of contemporary evolution in earthquake-created habitats
has been reported (Gelmond et al., 2009; Lescak et al., 2015; Bassham et al., 2018). Uplift of islands by alarge earthquake in Alaska in 1964 created new freshwater ponds in previously marine
areas, where oceanic threespine stickleback (Gasterosteus aculeatus) colonized and diversified (Gelmond et al., 2009; Lescak et al., 2015; Bassham et al., 2018). Genomic analysis of these populations
showed that standing genetic variation in the ancestral oceanic populations enabled rapid
diversification (Bassham et al., 2018).
A more recent large earthquake and subsequent tsunami created new stickleback habitats
in Japan (Kume et al., 2018). On 11 March 2011, a large earthquake with a moment magnitude of
9.0 occurred off the coast of Tōhoku, Japan, where the Pacific Plate is subducting beneath
the Okhotsk Plate (Urabe and Nakashizuka, 2016; Liu and Zhao, 2018). The earthquake triggered amassive tsunami, which struck stickleback habitats in Otsuchi Town, Iwate Prefecture,
Japan (Fig. 1) (Kume et al., 2018; Mori, 2018). Before the tsunami, three rivers in this area contained
sticklebacks: the Gensui River (a tributary of the Otsuchi River) and Teranosawa River(a tributary of the Kozuchi River) provided habitat for freshwater G. aculeatus (Fig. 1),
while the Namaisawa River (another tributary of the Kozuchi River) provided habitat
for both freshwater G. aculeatus and anadromous Japan Sea stickleback (G. nipponicus)(Fig. 1). The huge tsunami and earthquake-induced ground subsidence created new ponds
fed by spring water in a formerly urban area (orange area in Fig. 1). In 2012, we foundsticklebacks in two of these newly created ponds (Fig. 2) (Kume et al., 2018). Although
the source populations in these novel habitats have yet to be identified, previous microsatellite analysis by ourselves showed that these new populations contained hybrids between
G. aculeatus and G. nipponicus (Kume et al., 2018). Subsequently, we identified 48 stickleback
habitats (Fig. 2). However, both variation in environmental factors in these habitats
and stickleback phenotypes among these newly identified populations have yet to beinvestigated.

Fig. 1. Location of Otsuchi Town and the epicentre of the earthquake (A) and the southern urban area of Otsuchi Town (orange shaded area in B
and C). New stickleback habitats were formed in the southern urban area (see Fig. 2 for details). The dots in (B) and (C) indicate study sites of thenative habitats (green, Gensui River; light blue, Namaisawa River; dark blue, Teranosawa River). Satellite images used in this study were taken from
Google Earth Pro (v. 7.3.2.5491, Map data: Google Earth).

216

Hosoki et al.

Fig. 2. Landscape of the southern urban area of Otsuchi Town before (2005) and after the tsunami(2012, 2014, and 2017). Blue areas indicate the tsunami-formed pools and channels. Orange dots
denote the sites where stickleback were sampled. Satellite images used in this study were taken from
Google Earth Pro (v. 7.3.2.5491, Map data: Google Earth and Digital Globe).

In the present study, we characterized the patterns of morphological variation among
the newly created habitats, and tested whether the observed morphological variation canbe explained by any environmental factors. We specifically focused on gill raker number.
The number of gill rakers is known to be associated with trophic ecology in fishes: fish
specialized for benthos feeding tend to have fewer gill rakers, while planktivores tend to havea large number of gill rakers (Schluter and McPhail, 1993; Robinson, 2000; Roesch et al., 2013). Therefore, wefirst tested the hypothesis that fish with fewer gill rakers feed more on the benthos than fishwith many gill rakers. Compared with the oceans, freshwater streams and shallow ponds arericher in benthic prey and poorer in planktonic prey, and sticklebacks in these freshwater
environments generally have fewer gill rakers than marine or anadromous sticklebacks(Hagen and Gilbertson, 1972; Gross and Anderson, 1984; Schluter and McPhail, 1993). Therefore, we also tested
the hypothesis that sticklebacks in the habitats with more benthos and/or less access to thesea have fewer gill rakers.

This page corrected 21 January 2020
Stickleback colonization of earthquake-created habitats

217

MATERIALS AND METHODS
Study sites
Our study sites were 17 locations in newly formed habitats (Fig. 2; Tables 1–3). As some
habitats were difficult to access, environmental and benthos data are not available for allsites. For example, we were unable to conduct water measurements at G131 and W043
(Table 1), or benthos surveys at G019, G131, and W043 (Table 2). Although comparison
with native sites was not the main focus of this paper, we similarly surveyed four native sites(Gensui River, Teranosawa River, and N2 and N4 sites in the Namaisawa River; Fig. 1) andpresent the data here.
Environmental and ecological surveys of stickleback habitats
To characterize the environmental variation among stickleback habitats, we measured watersalinity, dissolved oxygen (mg/L), water temperature (⬚C), and pH (Table 1). At each site,measurements were conducted on the surface water using a multi-parameter water quality
probe (WQC-24, DKK-TOA Co. Ltd., Tokyo, Japan) in August 2014, December 2015, andMay 2016, and a glass electrode pH meter (WM-32EP, DKK-TOA Co. Ltd.) and dissolved
oxygen meter (DO-24P, DKK-TOA Co. Ltd.) in September 2016. As we could not adjust
tidal conditions, we conducted measurements multiple times at each site, except G019, tocover different tidal conditions (Table 1) and used the site-specific averages for statisticalanalysis.
For the analysis of benthos communities, benthic macroinvertebrates were collected by
the Institute of Environmental Ecology IDEA Consultants (Shizuoka, Japan) in August
2014 (Table 2), three times at each site. Each sampling covered a quadrat of 25 cm ×25
cm, to a depth of about 3 cm (1875 cm3). We used a hand-net with a 500 µm mesh forcollecting the benthos; thus, benthos smaller than 500 µm could not be assessed. Benthic
macroinverte-brates were separated from sediments and debris, pooled, and stored in
99.9% ethanol. The total numbers of benthic invertebrates per three replicates (1875 cm3 ×
3 = 5625 cm3) were used for subsequent analysis.
Morphological and stomach content analyses of sticklebacks
For morphological and stomach content analyses, adult fishes were collected with minnow
traps and euthanized using an overdose of MS-222 (tricaine methanesulfonate). We preserved all samples in 99.5% ethanol. Sample sizes and ranges of standard length are shownin Table 3.
Under a stereomicroscope, we counted the number of gill rakers from the first gill arch on
the right side of the fish, except for 10 individuals whose right gill rakers were damaged; forthese 10 fish, the gill rakers on the left side were counted instead. Counting was conductedwithout dissecting the gill or staining it with alizarin red. Gill rakers on the ventral anddorsal arches were summed. Because body size is also known to be an important factor in
determining the diets of sticklebacks (Schluter, 1993), we measured standard length with avernier calliper and used it as a covariate in the statistical analysis (see below). For comparison, we used previously reported data for native populations collected before theearthquake (Kume et al., 2018) and for anadromous sticklebacks from the Bekanbeushi River,
Hokkaido, Japan (Kitano et al., 2007) (Table 3). Because gill raker number was not sexually39.364664
39.352542
39.353724
39.356675
39.356286
39.356939
39.356600
39.356910
39.355697
39.356382
39.355635
39.356672
39.356594
39.356493
39.356124
39.356689
39.355222
39.356674
39.357035

Gensui
Namaisawa N2
Namaisawa N4
Teranosawa
G019
Pool1
Pool3
PoolS4
SAG021
W011
W042
W054
W086
W088
W089
W1034
W1052
W1081
W1088

141.897731
141.897410
141.893194
141.883638
141.902237
141.906903
141.908151
141.905234
141.901815
141.901910
141.902956
141.906668
141.904195
141.903731
141.902999
141.904639
141.903827
141.906280
141.905629

Latitude
2.0
2.5
2.9
4.9
0.7
0.5
0.6
0.5
0.7
0.8
0.8
0.6
0.4
0.6
0.7
0.5
0.8
0.8
0.7

Depth (m)Aug. 2014, Dec. 2015, May 2016
Aug. 2014, Dec. 2015, May 2016
Aug. 2014, Dec. 2015, May 2016
Aug. 2014, Dec. 2015, May 2016
Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May 2016
Aug. 2014, Dec. 2015
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015, May/Sept. 2016
Aug. 2014, Dec. 2015

Date of measurement—
224–297
224–280—226
220–286
260–279
274–313
214–303214–313
210–271
278–282
210–270
210–270
214–277
257–306
212–271
212–306
274–306

Tidal range (m)Note: Tidal range data were taken nearby at Ofunato City (39⬚1⬘N, 141⬚45⬘E), and provided by the Japan Meteorological Agency (https://www.data.jma.go.jp/gmd/
kaiyou/db/tide/suisan).

Longitude

Population

Table 1. Sampling dates and locations

This page corrected 21 January 2020
Stickleback colonization of earthquake-created habitats

219

Table 2. Variation in density of benthos among habitats

Population
Gensui River
Namaisawa N4
Namaisawa N2
Teranosawa River
Pool1
Pool3
PoolS4
SAG021
W011
W042
W054
W086
W088
W089
W1034
W1052
W1081
W1088

Mean depth of sampling (cm)No. individual benthic
organisms/5625 cm3

40
70
40
55
25
30
10
40
30
30
10
25
60
15
5
30
15
25

243
464
360
637
149
828
69
156
292
111
194
167
120
918
266
349
258
246

dimorphic except at W1081 (Wilcoxon rank sum test, W = 73.5, P = 0.0219 for W1081 and
P > 0.05 for other populations), we pooled the two sexes when testing the association
between gill rakers and habitats.
For the stomach content analysis, prey items in the stomach of each fish were sorted into
benthic and planktonic prey (Schluter and McPhail, 1992; Bell and Foster, 1994; Wetzel, 2001). Benthic preyincluded aquatic and terrestrial insects, Isopoda, Tanaidacea, Amphipoda, Ostracoda, andshells. Planktonic prey included Mysida, Cumacea, Copepoda, and Cladocera. For eachfish, we measured the volume of benthic prey, planktonic prey, and others (digested preyand unidentified items) as well as the total volume of the stomach contents.
Statistical analyses
To determine whether the environmental variables, stomach contents, and gill raker number
varied among newly formed stickleback populations, we conducted non-parametric
Kruskal-Wallis tests, because the Shapiro-Wilk test showed violations of the assumption of
normality for several variables. For analysis of gill rakers, including those of native populations, we conducted a Steel-Dwass test for multiple comparisons using both native andnovel populations. As Spearman correlation tests showed that several environmental andecological factors were correlated with each other (Table 4), we condensed these variables
for four native sites and 14 new populations by principal component (PC) analysis using
the prcomp function in R (Table 5). For PC analysis, we used site-specific means of logtransformed environmental variables.

220

Hosoki et al.

Table 3. Sample sizes for morphological (Nmorph) and stomach content (Nstomach) analyses, and themean and the range of standard lengths (SL) of stickleback fishes
Population

Sampling dateNmorph

Nstomach

Mean (range) of SL (mm)

29
5
20
20
43
20
20
19
4
20
20
20
23
5
5
20
3
20
27
10
25
21
34
24
27
26

—
—
20
20
—
20
20
19
4
20
20
20
23
5
5
20
3
20
27
10
25
21
34
24
—
—

50.88 (44.19–68.34)
61.63 (39.83–77.81)
56.77 (50.02–63.44)
52.56 (45.67–61.00)
55.89 (35.38–78.69)
53.21 (45.40–58.93)
55.12 (40.35–64.62)
47.61 (42.36–53.09)
56.87 (51.48–62.69)
46.14 (40.27–55.52)
47.76 (41.54–59.79)
45.67 (39.93–52.75)
52.02 (37.37–61.48)
51.03 (39.99–55.77)
49.91 (44.22–54.27)
50.62 (33.35–64.88)
57.71 (53.96–61.64)
36.18 (30.06–46.00)
47.28 (32.24–56.83)
46.57 (32.54–53.83)
40.80 (33.26–60.54)
46.46 (30.24–64.14)
45.50 (29.81–57.83)
43.49 (32.22–66.08)
79.17 (70.58–84.47)
64.40 (60.29–75.29)

Gensui River
December 1998
Gensui River
June 2010
Gensui River
February 2014
Namaisawa N2
February 2014
Namaisawa N3
June 2010
Namaisawa N4
February 2014
Teranosawa River
February 2014
G019
July 2014
G131
August 2014
Pool1
February 2014
Pool3
February 2014
PoolS4
February 2014
SAG021
July 2014
W011
July 2014
W042
March 2014
W043
March 2014
W054
August 2014
W086
March 2014
W088
March/July 2014
W089
March 2014
W1034
February 2014
W1052
March 2014
W1081
February/March 2014
W1088
March 2014
Akkeshi G. aculeatus
May 2003
Akkeshi G. nipponicus
May 2003

Table 4. Correlations among environmental and ecological variables

Salinity
Dissolved oxygen
Water temperature
pH
Benthos density

Salinity

Dissolved
oxygenWatertemperature

pH

Benthosdensity

—
0.838
0.079
0.682
0.066

−0.052
—
0.046
0.004
0.064

0.425
0.476
—
0.058
<0.001−0.104
0.643
0.455
—
0.276

−0.442−0.445−0.779−0.271—

Note: The matrix above the diagonal shows the Spearman’s correlation coefficients, while the lower matrix shows
P-values. Site-specific mean values were used for the analysis.

Stickleback colonization of earthquake-created habitats

221

To test the effects of gill raker number on stomach contents, we used a generalized linear
mixed model (GLMM) with stomach contents (either the volume of benthos or plankton in
the stomach) as a response variable, gill raker number and standard length as explanatory
variables, the total volume of stomach contents as an offset, and sampling site, sex, andsampling date as random factors. An offset term specifies that the volume of prey items inthe stomach is proportional to the total volume of stomach contents (Zuur et al., 2009; Kubo, 2012).
For the GLMM, we used the lme4 statistical package in R (Bates et al., 2015). For statisticalanalysis of stomach contents, fish with empty stomachs were excluded.
To determine whether any environmental factors influenced gill raker number, we used a
generalized linear model (GLM) with an environmental factor (either one of the PCs orsite-specific mean of each environmental factor) as an explanatory variable and site-specificmean number of gill rakers as a response variable: each environmental factor was tested
independently. A GLM with a gamma distribution and log-link function was used. For
testing the effects of habitat benthos density on stomach contents, we used habitat benthos
density as an explanatory variable and the site-specific mean of the proportion of either
benthos or plankton in the stomach contents as a response variable in the GLM with agamma distribution and log-link function.
RESULTS
Environmental variation among novel habitats
Novel habitats differed in salinity (Fig. 3; Kruskal-Wallis test, χ = 32.87, P = 0.003) anddissolved oxygen (χ2 = 26.89, P = 0.020) but not in water temperature or pH (P > 0.05). The
density of benthos also showed variations among novel habitats (Table 2). These environmental variations can be summarized into two major principal components: PC1 and PC2
explain 48.6% and 33.0% of the variation, respectively (Table 5; Fig. 4). Larger PC1 values
reflect higher values of water temperature, dissolved oxygen, pH and salinity, and lower
density of benthos, while higher PC2 values reflect higher values of salinity and lower valuesof pH and dissolved oxygen (Table 5). PC1 was significantly higher in the new habitats thanin the native habitats (Fig. 4) (Kruskal-Wallis test, χ2 = 5.46, P = 0.020), while PC2 did not
differ between habitats (Kruskal-Wallis test, χ2 = 0.18, P = 0.671).
2

Table 5. Loading values and percentage of variance explained by principal components (PCs) ofenvironmental and ecological variables
Environmental variable

PC1

PC2

PC3

PC4

PC5

Salinity
Dissolved oxygen
Water temperature
pH
Benthos density% Variance explained

0.284
0.389
0.580
0.345
−0.55948.6

0.613−0.5180.173−0.535−0.19833.0

0.448−0.376−0.056
0.726
0.357
8.3

0.557
0.662
−0.190−0.217
0.411
6.0

0.182
0.040
−0.7700.141−0.5934.0

222

Hosoki et al.

Fig. 3. Boxplots showing variation in water parameters among study sites. Environmental data werecollected in 2014–2016. For comparison, we show the data from the native habitats.

Diversity in gill raker number and diets
Gill raker number varied significantly among novel habitats (Kruskal-Wallis test,
χ2 = 86.70, P < 0.001) (Fig. 5). The proportions of benthos and plankton in the stomach
contents also showed significant variation among novel habitats (Kruskal-Wallis test,
χ2 = 87.04, P < 0.001 and χ2 = 107.41, P < 0.001 for benthos and plankton, respectively)(Fig. 6).
Fish with fewer gill rakers consumed higher proportions of benthic prey (GLMM:
P = 0.017, coefficient estimate ± S.E. = −3.72 ± 1.55) (Fig. 7). We also found that larger

Stickleback colonization of earthquake-created habitats

223

Fig. 4. Scatter plot of principal components (PC1 and PC2) of environmental and ecologicalvariables. The percentage of variance explained by each principal component is shown in parentheses.
Orange dots indicate novel habits, while other colours indicate native habitats (green, Gensui River;light blue, Namaisawa River; dark blue, Teranosawa River). The orange and grey circles indicate
95% confidence eclipses of novel and native sites, respectively.

fish foraged more benthos (GLMM: P < 0.001, coefficient estimate ± S.E. = 2.15 ± 0.35)(Fig. 7). These data are consistent with the idea that fish with fewer gill rakers and larger
body sizes have higher performance concerning benthos feeding (Schluter and McPhail, 1992; Schluter,1993). The proportion of planktonic prey was not significantly associated with gill raker
number (GLMM: P = 0.486, coefficient estimate ± S.E. = −0.37 ± 0.54) or standard length(GLMM: P = 0.854, coefficient estimate ± S.E. = −0.02 ± 0.12).
Next, we tested whether gill raker number is associated with any environmental factors.
We found that fish in habitats with higher PC1 values had higher numbers of gill rakers(GLM: P = 0.0175, coefficient estimate ± S.E. = 0.019 ± 0.007) (Fig. 8). PC2–PC5 had no
significant effects on gill raker number. When we analysed each of the environmental
factors, we found significant negative effects of benthos density (GLM: P = 0.0347,
coefficient estimate ± S.E. = −0.0001 ± 4 × 10−5) and marginally significant positive effectsof water temperature (GLM: P = 0.0511, coefficient estimate ± S.E. = 0.009 ± 0.004) on gillraker number.
Finally, we tested whether stomach contents were associated with any environmental
factors. Except for the effect of environmental PC5 on the proportion of benthos, wefound no significant effects of any environmental variable on the proportions of benthos or224

Hosoki et al.

Fig. 5. Box plots showing variation in gill raker number. For comparison, we show the data for sticklebacks in the native habitats, if available, and anadromous Gasterosteus aculeatus and G. nipponicus
collected from Hokkaido, Japan. Study sites with different letters indicate a statistically significant
difference in gill raker number between them (Steel-Dwass test, P < 0.05).

Fig. 6. Box plots of the proportions of prey items in stomach contents. Data are shown for benthic(upper panel) and planktonic prey (lower panel). All data were taken from fish collected in 2014. For
comparison, we show the data from the native habitats.

Stickleback colonization of earthquake-created habitats

225

Fig. 7. Relationship between the proportion of benthic prey in the stomach contents and morphology(gill raker number in the upper panel and standard length in the lower panel, N = 237) in fish fromnewly formed habitats. The lines show linear regression relationships with 95% confidence intervals.
Each dot represents a single fish. Note that the regression line is for visualization purposes only, andGLMM was used for statistical analysis (see Materials and Methods).

plankton in the stomach, regardless of whether any PC or single environmental factor wasused for the analysis (P > 0.05). Fish from habitats with higher PC5, which reflects lower
values of water temperature and benthos density (Table 5), tended to have more gill rakers(GLM: P = 0.0184, coefficient estimate ± S.E. = 0.8152 ± 0.299), although PC5 explained
only 4.0% of the environmental variations (Table 5).
DISCUSSION
We found diversity in gill raker number among new stickleback populations created by the2011 earthquake and subsequent tsunami in Japan. Gill raker number is highly heritable(Hagen, 1973; Hermida et al., 2002; Aguirre et al., 2004; Glazer et al., 2014) and rarely influenced by phenotypic
plasticity in sticklebacks (Day et al., 1994; Oke et al., 2016). Therefore, the observed diversity in gillraker number is likely to have a genetic basis. What drove the diversification of gill raker
number in the newly established populations?

226

Hosoki et al.

Fig. 8. Relationship between gill raker number and PC1 of the environmental factors (upper panel)and benthos density (lower panel). The lines show linear regression relationships with 95% confidence
intervals. Each dot represents a single novel site.

One possibility is that this variation may reflect adaptation to diverse environments. Fish
with fewer gill rakers were found in habitats with smaller PC1 (Fig. 8), which reflects ahigher density of benthos (Table 5). Furthermore, we found that fish with fewer gill rakers
had a higher proportion of benthos in their stomach contents (Fig. 7). These data areconsistent with previous findings that benthos-feeding sticklebacks generally have fewer gillrakers than plankton-feeding sticklebacks (Gross and Anderson, 1984; Schluter, 2000; Kitano et al., 2007;
Berner et al., 2009; Ravinet et al., 2014). However, we did not find a significant association between thedensity of benthos in a habitat and stomach contents. This is contrary to the idea that thedensity of benthos in a habitat influences the diets of sticklebacks residing there. Thus, itremains unclear how much the density of benthos in a habitat contributes to the evolutionof gill raker number. The results of the present study are based on few sampling time points.
Further detailed long-term surveys across multiple seasons will be necessary to reach a clear
conclusion on the relationships among gill raker number, diet, and habitat.
A second possible explanation is that different immigration rates of anadromous sticklebacks may cause variations in gill raker number. Anadromous G. nipponicus, which inhabits
the sea around Otsuchi Town, has a larger number of gill rakers than native freshwater

Stickleback colonization of earthquake-created habitats

227

G. aculeatus (see Fig. 5) (Kitano et al., 2007; Kume et al., 2018). Therefore, migration of the anadromous G. nipponicus into new habitats and hybridization with the fish there will result in
an increase in gill raker number in habitats accessible to G. nipponicus. Our preliminary
microsatellite analysis of two novel habitats in 2012–2013 showed that newly establishedpopulations contained hybrids between G. aculeatus and G. nipponicus (Kume et al., 2018),
suggesting that hybridization occurred in some of the new habitats. Genomic analysis ofadditional individuals in multiple newly formed populations will enable us to determine
when and how much gene flow has occurred from anadromous G. nipponicus and whether
gill raker number is associated with the rate of introgression.
The third possible explanation is that variation in gill raker number results from selection
on other phenotypic traits. In North American populations, one of the genetic loci controlling gill raker number was found on Chromosome 4 (Miller et al., 2014); the same chromosome
contains genes important for other morphological, physiological, and behavioural traits(Colosimo et al., 2005; Barrett et al., 2009; Greenwood et al., 2013; Peichel and Marques, 2017). Therefore, selection onother traits might change the allele frequency of genes regulating gill raker number throughhitchhiking effects (Maynard Smith and Haigh, 1974; Barton, 2000). Using quantitative trait loci (QTL)
mapping (Peichel et al., 2001; Peichel and Marques, 2017), it will be possible to identify genetic lociimportant for variation in gill raker number and other adaptive traits in the Japanesestickleback populations.
Finally, gill raker divergence may occur by chance. The genetic composition of founderpopulations may differ among newly colonized habitats because of the random sampling ofa limited number of alleles from the ancestral source populations (Whittaker and FernándezPalacios, 2007). Even after colonization, genetic drift will continue to fix different alleles indifferent populations for every generation (Kimura, 1983). If we can identify the genetic lociunderlying variation in gill raker number, we will be able to investigate the change in allele
frequency at QTL and get better insight into the roles of founder effects and genetic drift ingill raker number diversification.
In conclusion, we have highlighted a clear case of a catastrophic earthquake creating new
habitats for sticklebacks, which then diversified in gill raker number. Variation in gill raker
number among habitats was associated with environmental factors. This is a good model
system for investigating the ecological and genetic mechanisms by which phenotypic
diversification occurs in the earliest stage of colonization of newly formed habitats overcontemporary time-scales.
ACKNOWLEDGEMENTS
We thank Ken Sasaki, Masao Himura, and many Otsuchi citizens for help with sampling and fieldwork; participants at Stickleback 2018 for discussion; Yo Yamasaki, Ryo Kakioka, Asano Ishikawa,Andrew Hendry, and two anonymous reviewers for comments on the manuscript; Google Earth Pro
for the satellite images; the Japan Meteorological Agency for the tidal data; and the City of Ono and
ROIS for support of Stickleback 2018. This study was supported by JSPS Kakenhi 15H02418 and
Future Investment Project 2018 of ROIS to J.K., and the Environment Research and Technology
Development Fund (ZD-1203) of the Ministry of Environment Japan to S.M., J.K., and T.S.

REFERENCES
Aguirre, W.E., Doherty, P.K. and Bell, M.A. 2004. Genetics of lateral plate and gillraker phenotypes
in a rapidly evolving population of threespine stickleback. Behaviour, 141: 1465–1483.

228

Hosoki et al.

Barrett, R.D.H., Vines, T.H., Bystriansky, J.S. and Schulte, P.M. 2009. Should I stay or should I go?
The Ectodysplasin locus is associated with behavioural differences in threespine stickleback. Biol.
Lett., 5: 788–791.
Barton, N. 2000. Genetic hitchhiking. Phil. Trans. R. Soc. Lond. B: Biol. Sci., 355: 1553–1562.
Bassham, S., Catchen, J., Lescak, E., von Hippel, F.A. and Cresko, W.A. 2018. Repeated selection ofalternatively adapted haplotypes creates sweeping genomic remodeling in stickleback. Genetics,
209: 921–939.
Bates, D., Maechler, M., Bolker, B. and Walker, S. 2015. Fitting linear mixed-effects models using
lme4. J. Stat. Softw., 67: 1–48.
Bell, M.A. and Aguirre, W.E. 2013. Contemporary evolution, allelic recycling, and adaptive
radiation of the threespine stickleback. Evol. Ecol. Res., 15: 377–411.
Bell, M.A. and Foster, S.A., eds. 1994. The Evolutionary Biology of the Threespine Stickleback.
Oxford: Oxford University Press.
Berner, D., Grandchamp, A.C. and Hendry, A.P. 2009. Variable progress toward ecological
speciation in parapatry: stickleback across eight lake–stream transitions. Evolution, 63:
1740–1753.
Colosimo, P.F., Hosemann, K.E., Balabhadra, S., Villarreal, G. Jr., Dickson, M., Grimwood, J. et al.
2005. Widespread parallel evolution in sticklebacks by repeated fixation of Ectodysplasin alleles.
Science, 307: 1928–1933.
Darwin, C. 1839. The Voyage of the Beagle: Journal of Researches into the Natural History andGeology of the Countries Visited During the Voyage of H.M.S. Beagle Round the World. London:
Henry Colburn.
Day, T., Pritchard, J. and Schluter, D. 1994. A comparison of two sticklebacks. Evolution, 48:
1723–1734.
Gelmond, O., Von Hippel, F.A. and Christy, M.S. 2009. Rapid ecological speciation in three-spined
stickleback Gasterosteus aculeatus from Middleton Island, Alaska: the roles of selection andgeographic isolation. J. Fish Biol., 75: 2037–2051.
Glazer, A.M., Cleves, P.A., Erickson, P.A., Lam, A.Y. and Miller, C.T. 2014. Parallel developmental
genetic features underlie stickleback gill raker evolution. EvoDevo, 5: 19. Available at: https://doi.org/10.1186/2041-9139-5-19.
Grant, P.R. and Grant, B.R. 2011. How and Why Species Multiply: The Radiation of Darwin’s
Finches. Princeton, NJ: Princeton University Press.
Greenwood, A.K., Wark, A.R., Yoshida, K. and Peichel, C.L. 2013. Genetic and neural modularity
underlie the evolution of schooling behavior in threespine sticklebacks. Curr. Biol., 23: 1884–
1888.
Gross, H.P. and Anderson, J.M. 1984. Geographic variation in the gillrakers and diet of European
threespine sticklebacks, Gasterosteus aculeatus. Copeia, 1984: 87–97.
Hagen, D.W. 1973. Inheritance of numbers of lateral plates and gill rakers in Gasterosteus aculeatus.
Heredity, 30: 303–312.
Hagen, D.W. and Gilbertson, L.G. 1972. Geographic variation and environmental selection in
Gasterosteus aculeatus L. in the Pacific Northwest, America. Evolution, 26: 32–51.
Hendry, A.P. and Kinnison, M.T. 1999. Perspective: The pace of modern life: measuring rates ofcontemporary microevolution. Evolution, 53: 1637–1653.
Hermida, M., Fernández, C., Amaro, R. and Miguel, E.S. 2002. Heritability and ‘evolvability’ ofmeristic characters in a natural population of Gasterosteus aculeatus. Can. J. Zool., 80: 532–541.
Kimura, M. 1983. The Neutral Theory of Molecular Evolution. Cambridge: Cambridge University
Press.
Kitano, J., Mori, S. and Peichel, C.L. 2007. Phenotypic divergence and reproductive isolation
between sympatric forms of Japanese threespine sticklebacks. Biol. J. Linn. Soc., 91: 671–685.
Kubo, T. 2012. Introduction to Statistical Modeling for Data Analayis: Generalized Linear Model,
Hierarchical Bayesian Model, and MCMC. Tokyo: Iwanami Press.

Stickleback colonization of earthquake-created habitats

229

Kume, M., Mori, S., Kitano, J., Sumi, T. and Nishida, S. 2018. Impact of the huge 2011 Tohoku-oki
tsunami on the phenotypes and genotypes of Japanese coastal threespine stickleback populations. Sci. Rep., 8: 1684. Available at: https://doi.org/10.1038/s41598-018-20075-z.
Lack, D. 1947. Darwin’s Finches. Cambridge: Cambridge University Press.
Lescak, E.A., Bassham, S.L., Catchen, J., Gelmond, O., Sherbick, M.L., von Hippel, F.A. et al. 2015.
Evolution of stickleback in 50 years on earthquake-uplifted islands. Proc. Natl. Acad. Sci. USA,112: E7204–E7212.
Liu, X. and Zhao, D. 2018. Upper and lower plate controls on the great 2011 Tohoku-oki earthquake. Sci. Adv., 4: eaat4396. Available at: https://doi.og/10.1126/sciadv.aat4396.
Losos, J.B. 2010. Adaptive radiation, ecological opportunity, and evolutionary determinism. Am.
Nat., 175: 623–639.
Maynard Smith, J. and Haigh, J. 1974. The hitch-hiking effect of a favorable gene. Genet. Res.,
23: 23–35.
Miller, C.T., Glazer, A.M., Summers, B.R., Blackman, B.K., Norman, A.R., Shapiro, M.D. et al.
2014. Modular skeletal evolution in sticklebacks is controlled by additive and clustered quantitative trait loci. Genetics, 197: 405–420.
Mori, S. 2018. Utilization of environmental water resources in the reconstruction of Otsuchi
Town after the 2011 tsunami. In The Water–Energy–Food Nexus (A. Endo and T. Oh, eds.),
pp. 175–193. Singapore: Springer.
Oke, K.B., Bukhari, M., Kaeuffer, R., Rolshausen, G., Räsänen, K., Bolnick, D.I. et al. 2016. Does
plasticity enhance or dampen phenotypic parallelism? A test with three lake–stream stickleback
pairs. J. Evol. Biol., 29: 126–143.
Palumbi, S.R. 2002. The Evolution Explosion: How Humans Cause Rapid Evolutionary Change. New
York: Norton.
Peichel, C.L. and Marques, D.A. 2017. The genetic and molecular architecture of phenotypic
diversity in sticklebacks. Phil. Trans. R. Soc. Lond. B: Biol. Sci., 372: 20150486. Available at:https://doi.org/10.1098/rstb.2015.0486.
Peichel, C.L., Nereng, K.S., Ohgi, K.A., Cole, B.L.E., Colosimo, P.F., Buerkle, C.A. et al. 2001. The
genetic architecture of divergence between threespine stickleback species. Nature, 414: 901–905.
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. Available at:https://doi.org/10.1371/journal.pone.0112404.
Reznick, D.N. and Ghalambor, C.K. 2001. The population ecology of contemporary adaptations:what empirical studies reveal about the conditions that promote adaptive evolution. Genetica,
112/113: 183–198.
Robinson, B.W. 2000. Trade offs in habitat-specific foraging efficiency and the nascent adaptive
divergence of sticklebacks in lakes. Behaviour, 137: 865–888.
Roesch, C., Lundsgaard-Hansen, B., Vonlanthen, P., Taverna, A. and Seehausen, O. 2013. Experimental evidence for trait utility of gill raker number in adaptive radiation of a north temperate
fish. J. Evol. Biol., 26: 1578–1587.
Schluter, D. 1993. Adaptive radiation in sticklebacks: size, shape, and habitat use efficiency. Ecology,
74: 699–709.
Schluter, D. 2000. The Ecology of Adaptive Radiation. Oxford: Oxford University Press.
Schluter, D. and McPhail, J.D. 1992. Ecological character displacement and speciation in sticklebacks. Am. Nat., 140: 85–108.
Schluter, D. and McPhail, J.D. 1993. Character displacement and replicate adaptive radiation.
Trends Ecol. Evol., 8: 197–200.
Simpson, G.G. 1955. Major Features of Evolution. New York: Columbia University Press.
Urabe, J. and Nakashizuka, T. 2016. Ecological Impacts of Tsunamis on Coastal Ecosystems: Lessons
from the Great East Japan Earthquake. Tokyo: Springer Japan.

230

Hosoki et al.

Wang, C.Y. and Manga, M. 2015. New streams and springs after the 2014 Mw6.0 South Napa
earthquake. Nat. Commun., 6: 7597. Available at: https://doi.org/10.1038/ncomms8597.
Wetzel, R.G. 2001. Limnology: Lake and River Ecosystems. Oxford: Gulf Professional Publishing.
Whittaker, R.J. and Fernández-Palacios, J.M. 2007. Island Biogeography: Ecology, Evolution, andConservation. Oxford: Oxford University Press.
Wilson, G.A. and Rannala, B. 2003. Bayesian inference of recent migration rates using multilocus
genotypes. Genetics, 163: 1177–1191.
Zuur, A., Ieno, E.N., Walker, N., Saveliev, A.A. and Smith, G.M. 2009. Mixed Effects Models andExtensions in Ecology with R. New York: Springer.