Transition in Population Dynamics of the Intertidal Barnacle Balanus glandula after Invasion: Causes and Consequences of Change in Larval Supply

: To elucidate how the population dynamics of the acorn barnacle Balanus glandula transitioned after its invasion in 2000 along the Paciﬁc coast of Japan, a population census was conducted from 2004 to 2014 at ﬁve shores along 49 km of coastline 144–193 km east outside of the invasion front. Survey areas at each shore consisted of ﬁve paired plots (cleared recruitment plots and control plots). Larval recruitment was ﬁrst detected in 2004 but benthic individuals were not detected until 2 years later. The abundance and occurrence of B. glandula increased until around 2010; abundance then decreased but occurrence remained high (70%) until 2014, suggesting that the metapopulation of this barnacle approached a maximum around 2011. From 2011, the population dynamics of B. glandula changed considerably at two contrasting spatial scales: at a regional scale, the dependency of the number of larvae on stock size decreased, whereas at a local scale, the relative contribution of larval supply as a determinant of local population dynamics decreased. These ﬁndings suggest that the major driving force of population dynamics of the introduced barnacle changed in just a few years after invasion; therefore, population census data from just after an invasion, including larval recruitment monitoring just outside the invasion front, is essential to understanding invasion dynamics by sessile marine organisms.


Introduction
During the early phase of a biological invasion, both the range size and population density of the invading species increase with time [1][2][3], because an invasion event is usually initiated by the arrival of few propagules (i.e., introduced individuals) beyond their native range [4][5][6]. However, such temporal trends in distribution and abundance weaken with time after the invasion and are finally lost after population equilibrium has been reached [7,8]. This is a strong reason to suspect that population dynamics and their underlying processes change during the course of a biological invasion. To deepen our understanding of the nature of such invasions, it is important to reveal the transition that occurs in population dynamics after invasion. However, there have been only a few empirical studies of the early phases of invasion, including propagule supply, in marine habitats (see [9,10]).
Sessile animals are found in a variety of marine habitats and are a major component of marine ecosystems. Most of them have complex life histories that include pelagic and benthic phases linked by larval recruitment [11,12]. The transition between these two contrasting life phases can therefore J. Mar. Sci. Eng. 2020, 8,915 2 of 10 be crucial in determining the population dynamics of marine sessile organisms [13,14]. The number of larvae at a broad spatial scale (i.e., an entire metapopulation scale) will be determined by the regional reproductive output and regional hydrodynamics [15][16][17][18]. On the other hand, at local scales, population dynamics are often strongly affected by recruitment variability, which reflects external larval supply at the locality, especially at low levels of larval supply [19][20][21].
There are at least two fundamental hypotheses to explain the transition of population dynamics of sessile marine species after invasion. First, at a broader spatial scale, the dependency of the number of larvae on stock size, which is a determinant of the regional reproductive output, should decrease with time. This is because dependency of recruitment on stock size tends to be strong when the stock size is exceptionally small [22,23] and the stock size of an introduced species should be inevitably very small just after the invasion [10]. Second, at local scale, the relative contribution of larval supply as a determinant of local population dynamics should decrease with time. This is because the dependency of local population dynamics on larval supply tends to be strong when larval supply and local density are low [24][25][26], and both larval supply and local density should be low just after the invasion and could then increase [10].
Balanus glandula is a rocky intertidal barnacle native to the northwestern coast of North America along the northeastern Pacific Ocean [27,28]. This species sexually matures at 1 year of age [29] and produce around 30,000 larvae per clutch [27]. These planktonic larvae live a pelagic life of 2-4 weeks [30,31]. It is believed that this barnacle arrived on the Pacific coast of Japan in about 1992 [32]. Its distribution subsequently expanded northeastward to eastern Hokkaido [10,32,33].
To elucidate how the population dynamics of B. glandula transitioned after it invaded the eastern Pacific coast of Hokkaido (Figure 1a), northern Japan, population surveys were conducted from 2004 to 2014 at five shores. Surveys at each shore consisted of five paired plots (recruitment and control plots), along 49 km of coastline located 144-193 km east of the 2000 invasion front of this species [33]. While the recruitment of marine benthos is the product of both pre-and post-settlement processes, barnacle recruitment density on cleared suitable substrata is highly dependent on larval supply [34,35]. Therefore, barnacle recruitment density at our recruitment plots can serve as a proxy for larval supply at the locality [20,36]. The specific questions we addressed were (1) how do the distribution and abundance of benthic populations, and larval supply change with time after an invasion; (2) at the broadest spatial scale, did the dependency of the number of larvae on stock size decrease with time; and (3) did the relative importance of larval supply as a determinant of local population growth decrease with time?

Study Area
The study area is situated along 49 km of coastline 144-193 km east of Hiroo, which was known to be the eastern front of area invaded zone by B. glandula in 2000 [33] (Figure 1a). The study area is considered a sub-arctic zone [37,38] influenced by the cold Oyashio (Kurile Current), which flows westward along the shoreline, and where ice scour occurs once every several years. Low spring tides occur during the day during the warm season (from late March to early September) and during the night during the cold season (from late September to early March). In the study area, the mid-tidal zone is dominated by the native barnacle Chthamalus dalli [39]; the barnacle Semibalanus cariosus and perennial seaweeds, such as the brown alga Analipus japonicus, the red algae Chondrus yendoi and Gloiopeltis furcate, and the crustose coralline alga Corallina pilulifera are also present [40,41]. Invertebrate predators include the whelk Nucella lima and starfish Leptasterias ochotensis; the former preys upon barnacles [36], whereas the latter species is rare and restricted to the lower tidal zone [41,42].

Census Design
We applied a hierarchical nested sampling design [43] for the layout of the study sites; five sites were nested within each of the five shores along the Pacific coast of eastern Hokkaido (Figure 1b). Within each shore, we chose five sites located on steep rock walls at semi-exposed locations; distances between neighboring pairs of plots ranged from 10 to 15 m. At each site, permanent anchors drilled into roughly vertical rock marked adjacent recruitment (30 cm wide and 30 cm high) and control (50 cm wide and 100 cm high) plots, which were separated from each other by several tens of centimeters. For both recruitment and control plots, their vertical midpoints corresponded to the mean tidal level.
In recruitment plots, all surface organisms were cleared each May and August from 2004 to 2014 during low tide events by burning and scraping the rock surface with a wire brush. Censuses of recruitment plots were performed twice per year approximately 90 days after clearing, during the major recruitment period of B. glandula [10], by photographing 5 × 5 cm quadrats using a Canon IXX Digital 320 camera. Twelve replicate quadrat photographs were taken from each recruitment plot. The annual recruitment density of B. glandula was obtained as a summation of recruits in August and October at each site.
Control plots were observed each August from 2004 to 2014 to assess the occurrence (i.e., presence or absence in a plot) of B. glandula. In addition, estimates of B. glandula coverage were obtained each May from 2004 to 2014 by using a grid overlain on control plots with 200 grid points at evenly spaced intervals (5 cm vertically and 5 cm horizontally). The proportion of grid points occupied by B. glandula was then used to estimate the total coverage of B. glandula for each control plot at each census.

Dependence of Regional Larval Supply on Stock Size
As a surrogate for the amount of the regional larval supply for each year, we utilized the annual means of annual recruitment density obtained from recruitment plots. As a surrogate for the regional stock size for each year, we utilized the annual mean coverage of B. glandula obtained from control plots. To evaluate the temporal changes in dependency of regional larval supply on stock size, we used a moving window regression. For this analysis, we employed a simple linear regression in which the log of annual mean recruitment density each year and the log of mean coverage each year were treated as the dependent and independent variables, respectively, in a moving window (years of census), with equal weights assigned to all observations within the moving window. In a movingwindow regression, a narrow window can contain excessive noise, whereas a wide window is liable

Census Design
We applied a hierarchical nested sampling design [43] for the layout of the study sites; five sites were nested within each of the five shores along the Pacific coast of eastern Hokkaido ( Figure 1b). Within each shore, we chose five sites located on steep rock walls at semi-exposed locations; distances between neighboring pairs of plots ranged from 10 to 15 m. At each site, permanent anchors drilled into roughly vertical rock marked adjacent recruitment (30 cm wide and 30 cm high) and control (50 cm wide and 100 cm high) plots, which were separated from each other by several tens of centimeters. For both recruitment and control plots, their vertical midpoints corresponded to the mean tidal level.
In recruitment plots, all surface organisms were cleared each May and August from 2004 to 2014 during low tide events by burning and scraping the rock surface with a wire brush. Censuses of recruitment plots were performed twice per year approximately 90 days after clearing, during the major recruitment period of B. glandula [10], by photographing 5 × 5 cm quadrats using a Canon IXX Digital 320 camera. Twelve replicate quadrat photographs were taken from each recruitment plot. The annual recruitment density of B. glandula was obtained as a summation of recruits in August and October at each site.
Control plots were observed each August from 2004 to 2014 to assess the occurrence (i.e., presence or absence in a plot) of B. glandula. In addition, estimates of B. glandula coverage were obtained each May from 2004 to 2014 by using a grid overlain on control plots with 200 grid points at evenly spaced intervals (5 cm vertically and 5 cm horizontally). The proportion of grid points occupied by B. glandula was then used to estimate the total coverage of B. glandula for each control plot at each census.

Dependence of Regional Larval Supply on Stock Size
As a surrogate for the amount of the regional larval supply for each year, we utilized the annual means of annual recruitment density obtained from recruitment plots. As a surrogate for the regional stock size for each year, we utilized the annual mean coverage of B. glandula obtained from control plots. To evaluate the temporal changes in dependency of regional larval supply on stock size, we used a moving window regression. For this analysis, we employed a simple linear regression in which the log of annual mean recruitment density each year and the log of mean coverage each year were treated as the dependent and independent variables, respectively, in a moving window (years of census), with equal weights assigned to all observations within the moving window. In a moving-window regression, a narrow window can contain excessive noise, whereas a wide window is liable to smooth over significant changes of interest. We therefore chose window widths of 5, 6, 7, and 8 years after performing preliminary analyses at a range of widths.

Importance of Larval Supply as a Determinant of Local Population Growth
To examine temporal changes in the relative importance of larval supply as a determinant of spatial variation in local population growth, we calculated Spearman's rank correlation coefficients between the annual recruitment density at year t − 1 at each site and the amount of change in coverage between year t − 1 and year t at each site. These calculations were performed for four year-groups: were pooled because the sample sizes were too small (n < 9) to perform separate calculations for each year. These calculations were not performed before 2007 because recruitment density at year t − 1 was zero at each site during the period. to smooth over significant changes of interest. We therefore chose window widths of 5, 6, 7, and 8 years after performing preliminary analyses at a range of widths.

Importance of Larval Supply as a Determinant of Local Population Growth
To examine temporal changes in the relative importance of larval supply as a determinant of spatial variation in local population growth, we calculated Spearman's rank correlation coefficients between the annual recruitment density at were pooled because the sample sizes were too small (n < 9) to perform separate calculations for each year. These calculations were not performed before 2007 because recruitment density at year t − 1 was zero at each site during the period.

Change in Dependence of Regional Larval Supply on Stock Size
The temporal relationship between annual recruitment density (a surrogate for regional larval supply) and mean coverage (a surrogate for stock size) yields an asymptotic up pattern when plotted, showing saturation of annual recruitment density at log mean coverages greater than around 0.1 (Figure 3). The dependence of annual recruitment density on mean coverage tended to decrease with time after invasion ( Table 1). Regardless of the window size, the regression coefficient and R 2 , the

Change in Dependence of Regional Larval Supply on Stock Size
The temporal relationship between annual recruitment density (a surrogate for regional larval supply) and mean coverage (a surrogate for stock size) yields an asymptotic up pattern when plotted, showing saturation of annual recruitment density at log mean coverages greater than around 0.1 (Figure 3). The dependence of annual recruitment density on mean coverage tended to decrease with time after invasion ( Table 1). Regardless of the window size, the regression coefficient and R 2 , the coefficient of determination, were highest in the first or second regression intervals starting in 2004.
After that, they decreased with time. Both parameters were higher before the peak of coverage (2011) than after. coefficient of determination, were highest in the first or second regression intervals starting in 2004. After that, they decreased with time. Both parameters were higher before the peak of coverage (2011) than after.

Change in Importance of Larval Supply as a Determinant of Local Population Growth
There was no significant spatial association between local population growth and annual recruitment density throughout the entire study period (Spearman's rank correlation coefficient = 0.26, p = 0.083) (Figure 4). The spatial association between local population growth and annual recruitment density changed with time after invasion ( Table 2). The Spearman's rank correlation coefficients were positive before the time interval of 2010-2011, with a significant value for 2010-2011. Thereafter, the coefficients were negative but statistically insignificant.

Change in Importance of Larval Supply as a Determinant of Local Population Growth
There was no significant spatial association between local population growth and annual recruitment density throughout the entire study period (Spearman's rank correlation coefficient = 0.26, p = 0.083) (Figure 4). The spatial association between local population growth and annual recruitment density changed with time after invasion ( Table 2). The Spearman's rank correlation coefficients were positive before the time interval of 2010-2011, with a significant value for 2010-2011. Thereafter, the coefficients were negative but statistically insignificant.

Temporal Changes in Coverage, Occurrence, and Recruitment
Balanus glandula, which invaded the study area [10], increased its coverage and occurrence in this area until around 2010, and then considerably decreased in coverage and occurrence until 2013. Similar "boom-and-bust" dynamics have been documented from the early phase of invasion of many introduced species [44,45], but the mechanism underlying these dynamics has not been understood in most of these cases [45] (but see [46,47], including the case of B. glandula introduced to eastern Hokkaido). It is likely that the decline of B. glandula during 2011-2013, which occurred across the entire 49 km of coastline in the study area, was caused by some degradation of the abiotic environment occurring at a large spatial scale. Ice scour often heavily damages intertidal barnacle populations in the sub-arctic zone [48], but it did not occur in this region after 2008 [49]. Also, this region did not experience any extremely hot summers from 2011 to 2013; such summers are another major cause of drastic declines in the abundance of sub-arctic intertidal barnacles [50,51]. None of the three indices for the hotness of summer-the average daily maximum temperature per month, the maximum temperature per month, or the amount of solar radiation-showed notably high values from 2010 to 2013 [49]. Finally, the 2011 earthquake off the Pacific coast of Tohoku (Mw 9.0), which caused a tsunami in this region, did not coincide with the decline of B. glandula during 2011-2013.

Temporal Changes in Coverage, Occurrence, and Recruitment
Balanus glandula, which invaded the study area [10], increased its coverage and occurrence in this area until around 2010, and then considerably decreased in coverage and occurrence until 2013. Similar "boom-and-bust" dynamics have been documented from the early phase of invasion of many introduced species [44,45], but the mechanism underlying these dynamics has not been understood in most of these cases [45] (but see [46,47], including the case of B. glandula introduced to eastern Hokkaido). It is likely that the decline of B. glandula during 2011-2013, which occurred across the entire 49 km of coastline in the study area, was caused by some degradation of the abiotic environment occurring at a large spatial scale. Ice scour often heavily damages intertidal barnacle populations in the sub-arctic zone [48], but it did not occur in this region after 2008 [49]. Also, this region did not experience any extremely hot summers from 2011 to 2013; such summers are another major cause of drastic declines in the abundance of sub-arctic intertidal barnacles [50,51]. None of the three indices for the hotness of summer-the average daily maximum temperature per month, the maximum temperature per month, or the amount of solar radiation-showed notably high values from 2010 to 2013 [49]. Finally, the 2011 earthquake off the Pacific coast of Tohoku (Mw 9.0), which caused a tsunami in this region, did not coincide with the decline of B. glandula during 2011-2013.
Although the occurrence of B. glandula decreased from 2011 to 2013, it was still at a relatively high level (e.g., 70%) during 2013 and 2014. This suggests that the metapopulation of B. glandula in the region had been approaching its maximum rather than continuing in an exponentially growing phase after 2011.

Change in Dependence of Regional Larval Supply on Stock Size
For B. glandula introduced to eastern Hokkaido, the dependence of annual recruitment density (a surrogate for regional larval supply) on mean coverage (a surrogate for stock size) decreased with time after invasion, exhibiting especially weak dependency after the value of mean coverage reached its peak in 2011 (Table 1). This result suggests that, at a broader spatial scale encompassing the entire metapopulation, the dependency of the number of larvae on stock size decreased with time after invasion, because the number of larvae returning to the near-shore zone reached a maximum set by limitation through factors other than the stock size. It is unlikely that the upper limit was related to the availability of open space, which is often a major factor limiting larval recruitment [52][53][54]. This is because the proportion of free space on recruitment plots (i.e., the space unoccupied by sessile organisms) consistently exceeded 70% during 2010-2014. Alternatively, food availability, which often limits the growth and survival of barnacle larvae [55,56], might explain the peak in larval recruitment on the recruitment plots. This is because many laboratory experiments have demonstrated that the planktotrophic larvae of several marine benthic species exhibit density-dependent mortality, with food being suggested as a factor limiting larval survival during their pelagic lives [57][58][59].

Change in Importance of Larval Supply as a Determinant of Local Population Growth
For B. glandula introduced to eastern Hokkaido, the spatial association between local population growth and annual recruitment density was stronger before 2011 than after. This suggests that, at a local scale, the relative contribution of larval supply as a determinant of local population dynamics decreased with time after invasion. This may be explained by increases in local density but not in larval supply. This is because the spatial association between local population growth and annual recruitment density was strongest in 2011, when larval supply was highest during the study period. In addition, multiple regression analysis, in which annual population growth on plots was treated as the dependent variable, and the initial cover of B. glandula and annual recruitment rates were treated as independent variables (Table S1), showed an inverse relationship between initial coverage and population growth.

Conclusions
The invading barnacle B. glandula increased in coverage and occurrence from 2004 until around 2010, and then declined in coverage but with a high level of occurrence (e.g., 70%) until 2014, suggesting that the metapopulation of B. glandula in the region approached its peak around 2011. After 2011, the population dynamics of B. glandula exhibited considerable change at two contrasting spatial scales: at a regional scale, the dependency of number of larvae on stock size decreased, whereas at a local scale, the relative contribution of larval supply as a determinant of local population dynamics decreased.
Our results indicate that the major driving force of the population dynamics of an introduced barnacle changed in just a few years after invasion. This finding highlights the importance of understanding the early phase of invasion dynamics, during which the population has not reached asymptotic equilibrium [10,60,61]. Such research, however, may be easier said than done, because there are only limited opportunities for this kind of study, especially for marine benthos with life histories similar to that of B. glandula. It is therefore crucial to conduct long-term population censuses, including monitoring of larval recruitment, in the areas surrounding the invasion fronts of these organisms.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-1312/8/11/915/s1, Table S1: Multiple regression analysis, in which annual population growth on plots was treated as the dependent variable and initial cover of B. glandula and annual recruitment rates were treated as independent variables. Funding: This research was funded by grants-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology Japan to T.N. (nos. 20570012, 24570012, 15K07208, and 18H02503).