Outstanding performance of an invasive alien tree Bischofia javanica relative to native tree species and implications for management of insular primary forests

Study site

The oceanic Ogasawara Islands are located in a subtropical region of the Pacific Ocean (between 24°14′N and 27°44′N, and 140°52′E and 142°16′E). The resident biota contains a high percentage of endemic species (Shimizu, 2003). Haha-jima Island is one of the two inhabited islands in the archipelago. It covers 20 km² and has a maximum elevation of 463 m above sea level. The island’s central mountains are covered by mesic forests that consist of relatively tall trees (about 15 m in height) compared with that of other forests in the Ogasawara Islands. The Sekimon mesic forests cover uplifted limestone in the northeastern corner of Haha-jima. The uplifted limestone has a doline-like central depression. Relatively thick sedimentary soil (Okamoto et al., 1995) and protection from wind by the walls of the doline have favored the growth of dense, tall forest on the base of the doline. This environment provides habitat for many plant species that the distributions are restricted to the Sekimon (Abe, Tanaka & Shimizu, 2018). B. javanica was introduced to the Ogasawara Islands for the silvicultural purpose before 1905 (Toyoshima, 1938; Shimizu, 2003). Although there is no record of planting B. javanica in the Sekimon in the forest management ledger, a participator attested that B. javanica had been planted before 1935 (Toyoda, 2003). In 1997, the seaward edge of the doline collapsed (Fig. 1A) and, subsequently, many trees have been exposed to salt-bearing onshore wind, causing desiccation and salt damage to the trees.

This area was struck by a strong typhoon in late 2006. Typhoon 0614 YAGI was spawned on 19 September in the northwestern Pacific (20.3°N, 159.2°E), about 1,800 km southeast of the Ogasawara Islands. The typhoon was closest to Haha-jima Island on 22 and 23 September when it passed about 100 km west of Haha-jima. At that time, the atmospheric pressure decreased to 930 hPa, the maximum wind velocity attained 45 m s⁻¹, and the 170 km radius of storm area experienced a wind velocity ≥ 25 m s⁻¹ estimated by the Dvorak method (Japan Meteorological Agency, 2018).

Field survey

We selected a survey area in the central portion of the primary forests in the Sekimon area and established two 2-ha census plots (100 m × 200 m) because there is a steep limestone ridge difficult to traverse between the two plots. We surveyed all trees with diameter at breast height (DBH) ≥ 10 cm in 2006 and described the status of each individual’s crown in terms of whether it formed part of the forest canopy or understory. We defined canopy trees as individuals in which more than half of the crown surface was exposed to direct sunlight (i.e., not shaded by neighboring trees); for individuals classified as an understory tree, we recorded the tree species that covered the largest proportion of its crown. This judgement was conducted by eyesight, aided by observation using binoculars when necessary. In 2008, we conducted a second census following the same method of the first census. The abbreviations shown in Table 1 were used for the species names used in the figures and tables in this paper.

Shimizu (1994) surveyed a portion of our study site in 1987 using a 100 m × 50 m plot (Fig. S1). The southern portion of this plot disappeared in a landslide in 1997 (Fig. 1A). The present study plot included the remaining portion (60 m × 50 m) of the Shimizu plot in the southeastern part of the western plot. Our reconstruction of the Shimizu plot was based on a tree-by-tree map drawn in 1987 (Shimizu, 1994). We checked the position of characteristic large trees (e.g., Melia azedarach) and old stumps of Morus boninensis that had been cut about 130 years previously but had not decomposed because of the strong, decay-resistant wood (Yoshida & Oka, 2000). The 1987 data enabled us to analyze changes in species composition in terms of the number of stems and basal area. However, we could not analyze individual mortality and growth since 1987 because Shimizu (1994) did not label individual trees.

To detect the impacts of typhoon 0614 YAGI, we surveyed the damage soon after the first tree census (November and December 2006). We recorded the types of damage for individual trees with DBH ≥ 10 cm in the northern half of the western plot (1 ha, N = 2675). The damage to each tree was classified as defoliated, snapped, uprooted, or trapped (under one or more uprooted trees). Among these damaged trees, the stems that died at the 2008 survey were judged to have died due to typhoon damage, and the mortality rate was defined as the number of the dead stems in 2008 divided by the number of stems in 2006 damage survey.

Field survey was approved for the Ogasawara National Park by the Ministry of Environment (No. 0606328007, No.080507006) and for the Ogasawara National Forest by the Forest Agency (No.18-2-50 and No.20-1-32).

Statistical analyses

We evaluated the annual diameter growth rate in 2-year period as ((DBH in 2008)-(DBH in 2006))/(survey interval months)*12/(DBH in 2006)*100 for each tree species. The morality rate of each tree species was defined as the number of dead stems in the 2008 survey divided by the number of stems in 2006 survey. The population growth rate was defined as the period growth rate of the number of stems: (N in 2008)-(N in 2006)/(N in 2006)*100, where N is the number of stems. Generally, trees have a trade-off relationship between growth and survival (Grubb, 1977; Hubbell & Foster, 1992; Wright et al., 2003), but B. javanica on Hahajima Island seemed to have good performance for both. To confirm this, the Pearson’s product-moment correlation coefficient between the annual diameter growth rate and the population growth rate was examined when all tree species were used and when only B. javanica was removed.

Differences between B. javanica and native trees for typhoon damages and stem dynamics (mortality and recruitment) were examined by a Tukey’s HSD multiple comparison after generalized linear model (GLM) analyses using the multcomp package in R ver. 3.3.2 (R Core Team, 2016). The GLMs of typhoon damage were conducted independently for each type of damage and mortality assuming a binomial error distribution with the number of damaged stems as a responsible variable and the tree species as an explanatory variable. The GLMs of population growth were conducted assuming a binomial error distribution with the number of recruited stems or the number of dead stems as a responsible variable and the tree species as an explanatory variable, respectively. We examined the effects of crown position on diameter growth of understory tree stems using two types of analysis: the effect of the canopy tree species on a given understory species and the growth differences among the understory tree species under a given canopy species. Both analyses used a general linear model (GLM) with a Gaussian link function and a multiple-comparison test using R. The responsible variable was the annual diameter growth rate of understory tree stems in both GLM analyses. The explanatory variable was understory tree species in the comparison among understory species under a given canopy species and was canopy tree species in the comparison among canopy species over a given understory species.

In the tree invasion process, it is effective to cover the understory trees with a wide crown in addition to the fast growth. Even if individual understory stems are likely to die sooner or later, there are always many stems under the wide canopy in the process of development of canopy trees, and conversely there would be only fewer stems with more than 10 cm DBH under the narrow canopy. Since we did not directly measure individual crown widths, we used, simply assuming that there are many stems under the wide crown, the following formula to index the crown area (CW) of each tree species: ${\displaystyle CW=NS∕NC}$ where NS is the number of stems covered by the crown of the canopy species and NC is the number of canopy stems of the species.

Prediction of increase in B. javanica occupancy

It is preferable to use highly accurate models, such as a population matrix, to predict the population dynamics of an invasive tree species (e.g., Buckley, Briese & Rees, 2003). However, we could not use such a model in the present analysis because we surveyed the young trees less than 10 cm in DBH including seedlings only once (Abe, Tanaka & Shimizu, 2018). Instead, we used a simple logistic curve (Radosevich, Stubbs & Ghersa, 2003; Webster & Wangen, 2009) to predict future population growth of B. javanica in terms of the number of stems and basal area. Given that it can be assumed that the spread of an invasive tree species is random and continuous within the forest, a simple model prediction is considered to be sufficiently practicable (Frappier et al., 2003). The model represented the proportion of B. javanica (D_BJ) with an upper limit of 1.0 for the proportion, as follows: ${\displaystyle D_{\mathrm{BJ}}=1∕\left\{1+a\times exp\left(-b\times t\right)\right\}}$ where t represents the number of years since 2006. The coefficients a and b were determined based on the data from the 1987 measurements in the Shimizu (1994) plot and the 2006 measurements in Abe, Tanaka & Shimizu (2018) (Table S1). Although the two plots were separated for convenience because of the cliff between them, the vegetation of both plots is considered to be homogeneous. Accordingly, we applied these parameters to the prediction of B. javanica dynamics in both plots.

We predicted the time required for B. javanica to attain 30% and 50% of the number of stems and basal area for the western plot and eastern plot, using logistic regression models. The lower percentage (30%) was based on the guideline of the National Forest that restricts the proportion of tree removal less than 30% of the total volume to prevent soil erosion. The higher percentage (50%) was based on data from the forests on Mt Kuwanoki (Haha-jima Island), where the former forest type had been identical to that at the Sekimon but now resembles a B. javanica forest stand with more than 40% occupancy of the total basal area (Shimizu, 1988). In addition, as a property of the logistic model, the estimated year tends to include a smaller error in the central portion of the logistic curve (e.g., between 30% and 70% occupancy) than that at each extreme (i.e., the first year of invasion and the end of the simulation period). Therefore, forecast years reaching 30% and 50% occupancy are expected to be most accurate and robust.

Survival, growth, and typhoon damage

Typhoon 0614 YAGI was situated closest to Haha-jima Island on 22 and 23 September 2006. The typhoon defoliated all standing stems (Fig. 1B), and snapped, uprooted, and trapped trees accounted for 6.9%, 2.6%, and 0.2% of the total, respectively (Table 2). There was no significant difference in the proportion of stems of these types of typhoon damage between native species and B. javanica. Pioneer trees (sun-lit trees growing rapidly in the early stage of succession or in the gaps) exhibited relatively high mortality (Zanthoxylum ailanthoides var. inerme at 16.7%, Trema orientalis at 33.3%, and Cyathea mertensiana at 21.4%), as did some later-successional species (Ochrosia nakaiana at 50.0% and Psychotria homalosperma at 21.4%). B. javanica showed low mortality (1.9%) in response to the typhoon disturbance.

The number of stems decreased between 2006 and 2008 among the most frequent tree species (more than 30 stems in the plots) except for B. javanica (7.4% increase) (Fig. 2). The increment in B. javanica was the result of recruitment of 44 individuals to the DBH ≥ 10 cm size class and the death of 10 individuals. Species that showed the greatest decrease in number of stems were an endemic pioneer, Z. ailanthoides var. inerme (−43.3%), and an endemic tree fern, Cyathea mertensiana (−34.8%). The proportion of the number of recruitments into the stem size class DBH ≥ 10 cm was largest for the alien species B. javanica (8.8%) followed by Callicarpa subpubescens (6.9%) and Ficus boninsimae (6.6%). Some native species had a significantly higher proportion of the number of dead stems and significantly less proportion of the number of recruitments than B. javanica (Fig. 2). Annual diameter growth rate (Fig. 3) was largest in B. javanica (3.1 ± 0.1%, mean ± SE) followed by three pioneers, C. mertensiana (2.1 ± 0.4%), Z. ailanthoides var. inerme (2.1 ± 0.3%), and C. subpubescens (2.0 ± 0.3%). The diameter growth rates of dominant native species were less than half that of B. javanica (e.g., Ardisia sieboldii at 0.8 ± 0.0%, Elaeocarpus photiniifolius at 1.0 ± 0.1%, and Pisonia umbellifera at 1.3 ± 0.1%). Annual diameter growth rate was negatively correlated with population growth rate when the data for B. javanica were omitted from those for the most frequent tree species (Pearson’s product-moment correlation, r = −0.635, t = −3.182, df = 15, p = 0.006), but no significant relationship was observed when the data for B. javanica were included (r = −0.225, t = −0.922, df = 16, p = 0.370).

Effects of crown shading

The number of trees in which more than half of the crown was shaded by the crown of a neighboring tree in 2008 was 2761 (39.9% of all stems, Fig. 4); the number was largest for A. sieboldii (1956), P. umbellifera (301), and B. javanica (105). The most frequent canopy species were E. photiniifolius (793), B. javanica (685), and Celtis boninensis (219).

The mean annual diameter growth of understory trees was significantly less than that of canopy trees (GLM with a Gaussian link function; estimate = 0.059, t = 8.32, P < 0.001). The canopy of B. javanica significantly decreased the diameter growth of several understory tree species: diameter growth was significantly decreased for A. sieboldii than under E. photiniifolius and under Z. ailanthoides var. inerme, and for P. umbellifera than under A. sieboldii (Fig. 5). On the other hand, understory individuals of B. javanica exhibited superior growth compared with that of native understory tree species, regardless of the canopy tree species (Fig. 6). Although the CW index was much larger in M. azedarach (CW = 5.3) and C. boninensis (4.9) compared with that of all other species (Fig. 7), the largest values of CW among dominant species (i.e., those with ≥ 100 canopy individuals) were for E. photiniifolius (2.2), followed by B. javanica (1.9) and Planchonella obovata var. obovata (1.1). The most frequent dominant species, A. sieboldii, showed a small CW index (<0.1).

Prediction of invasion by B. javanica

In the Shimizu plot, B. javanica increased substantially in both the number of stems (176.4%) and basal area (177.8%) for the 19-year period (Table S1). We applied these changes for B. javanica to estimate the coefficients of logistic curves (Fig. 8). The coefficients of the logistic model were a = 36.214 and b = 0.038 based on the number of stems, and a = 36.155 and b = 0.051 based on the basal area. The model predicted that in the eastern plot, B. javanica will account for 30% of the number of stems in 2033 and 30% of the basal area in 2017. In the eastern plot, B. javanica will account for 30% of the number of stems in 2087 and 30% of the basal area in 2057. In the eastern plot, B. javanica will account for 50% of the number of stems in 2056 and 50% of the basal area in 2034. In the western plot, B. javanica will account for 50% of the number of stems in 2109 and 50% of the basal area in 2074.