APP下载

Investigating factors driving phytoplankton growth and grazing loss rates in waters around Peninsular Malaysia

2021-02-22JoonHaiLIMChoonWengLEEChuiWeiBONG

Journal of Oceanology and Limnology 2021年1期

Joon Hai LIM , Choon Weng LEE , Chui Wei BONG

1 Laboratory of Microbial Ecology, Institute of Biological Sciences, University of Malaya, Kuala Lumpur 50603, Malaysia

2 Institute of Ocean and Earth Sciences, University of Malaya, Kuala Lumpur 50603, Malaysia

3 Institute for Advanced Studies, University of Malaya, Kuala Lumpur 50603, Malaysia

Abstract In tropical waters where temperatures are relatively stable, we investigated whether the relationship between phytoplankton growth and grazing loss rate across different habitats around Peninsular Malaysia can be attributed to its eutrophication states. We measured phytoplankton growth ( μ) and grazing loss ( g) rates in waters off Bachok Marine Research Station (BMRS), located northeast of Peninsular Malaysia. Chlorophyll- a (chl-a) concentration ranged from 2.90 to 15.78 μg/L and was dominated by nano- and micro-phytoplankton (>2 μm in size). Using the Landry and Hassett dilution method, μ at BMRS ranged from 1.02 to 1.58/d whereas g varied from 0.07 to 0.88/d. Grazing accounted for 35% of the primary production at BMRS. A systematic review of available data in waters around Peninsular Malaysia, revealed how μ fluctuated over a wide range (0.01-1.80/d) and correlated with chl a distribution ( R 2=0.181, P <0.001). However, the relationship was only significant at <9 μg/L chl a for mesotrophic waters and <16 μg/L chl a for eutrophic waters. In contrast, g ranged from 0.00 to 1.01/d, and correlated with μ at all locations. The g/ μ slope ranged from 19% to 84%, and was generally similar for waters around Peninsular Malaysia. However, all the g/ μ slopes had a positive y-intercept except for BMRS, and this seemed to suggest the availability of alternative prey supporting grazing at the other stations.

Keyword: phytoplankton growth; grazing loss; grazing pressure; Bachok Marine Research Station (BMRS)

1 INTRODUCTION

In the oceans, primary production is mainly carried out by phytoplankton (Duarte and Cebrián, 1996; Field et al., 1998), and phytoplankton biomass that is produced is transferred to higher trophic level via grazing (Vargas et al., 2007; Selph et al., 2016). Across different regions, microzooplankton grazing accounts for an estimated 44% to 79% of daily primary production (Schmoker et al., 2013). Phytoplankton production not grazed is either lost through horizontal export or sinking (Turner, 2002; Stukel et al., 2011). As primary production and grazing loss rates are affected by different environmental drivers (Boyce et al., 2010; López-Urrutia and Morán, 2015; Zhang et al., 2018), their values are distributed over a wide range and are sitespecific. It is important to measure these biological processes as the balance between both primary production and grazing rates affects the distribution of primary producers (Zheng et al., 2015; Cermeño et al., 2016).

In a cross-compilation study of 110 datasets by Schmoker et al. (2013), only 24% were from tropical waters. Data from Sunda Shelf were clearly lacking even though two additional studies were published from the Sunda Shelf waters i.e. Strait of Malacca (Lim et al., 2015, 2018) and Strait of Singapore (Schmoker et al., 2016). Sunda Shelf waters are highly productive tropical waters (Ooi et al., 2013), and clear differences in terms of eutrophication exist in waters around Peninsular Malaysia (Lim et al., 2018).

Fig.1 The study area and sampling locations

Although temperature is an important environmental driver of biological processes (Boyce et al., 2010), in tropical waters where temperatures are stable, most biological processes are already functioning at near optimum (Pomeroy and Wiebe, 2001). Therefore, the balance between primary production and grazing loss rates in these waters are primarily determined by the productivity or trophic state of the waters as grazing loss is reported to be higher in productive waters (Schmoker et al., 2013). In this study, we hope to determine whether the statement holds true by adopting the Landry and Hassett (1982) dilution method to measure concurrent phytoplankton growth and grazing loss rate. Despite the possible limitation of this approach, i.e. nutrient limitation, grazing saturation and non-linear response (Ayukai, 1996; Paterson et al., 2008; Schmoker et al., 2016), the rates obtained via this method have generally been comparable to the more common light-dark dissolved oxygen and14C uptake methods (Landry et al., 2011; Lim et al., 2015; Selph et al., 2016).

In this study, we measured phytoplankton growth and grazing loss rates in waters northeast of Peninsular Malaysia (BMRS), and then analysed these rates from waters around Peninsular Malaysia: west of Peninsular Malaysia i.e. Port Klang (PK) and Port Dickson (PD) (Lim et al., 2015); northwest (NW) and southeast (SE) of Peninsular Malaysia (submitted); and south of Peninsular Malaysia i.e. Strait of Singapore (SS) (Schmoker et al., 2016). In general, NW and SE stations were shelf waters located >5 km away from the shoreline whereas BMRS and PD were stations located next to sandy beaches. PK was at Klang River estuary and SS was an enclosed bay south of Singapore. For nearshore habitats such as BMRS, PD, PK and SS, they are commonly influenced by semidiurnal tide where two episodes of equally high and low waters are recorded each day.

At present, we aim to address the following questions: What are the growth and grazing loss rate of primary producers in the northeast coast of Peninsular Malaysia and can the differences in phytoplankton growth and grazing loss rate across different habitats around Peninsular Malaysia be attributed to its trophic states? Other than filling in an obvious data gap, the short-term sampling at BMRS would also allow us to better understand the fluctuation of phytoplankton growth and grazing loss rates throughout the day. These short-term data are unavailable from published reports in this region (Lim et al., 2015; Schmoker et al., 2016). As primary production and grazing loss rates can respond rapidly to environmental change (Goes et al., 2016), the results obtained here could also direct us on possible long-term research strategies.

2 MATERIAL AND METHOD

2.1 Sampling and physiochemical parameters

Sampling was carried out at Bachok Marine Research Station (BMRS) (6°00′N, 102°25′E), located northeast coast of Peninsular Malaysia (Fig.1). Twelve samplings were carried out over a three-day period where about 10-L sample was collected every time. Temperature, salinity (YSI Pro-30, USA), and pH (Metrohm 827, Switzerland) were measured in-situ whereas dissolved oxygen (DO) concentration was determined via Winkler titration method (Grasshoff et al., 1999).

Seawater sample was filtered through glass fiber filters (GF/F) that were pre-combusted at 500 °C for 3 h. The filtrate was then kept at -20 °C before dissolved inorganic nutrient analyses (ammonium: NH4; phosphate: PO4; nitrate+nitrite: NO3+NO2; and silicate: SiO4) whereas the filters were kept for chlorophyll a (chl a) and total suspended solid (TSS). Dissolved inorganic nutrients were measured by the colorimetric method with a spectrophotometer (Hitachi U-1800, Japan) (Parsons et al., 1984). NH4was measured by the blue indophenol method whereas NO3was reduced to NO2with a cadmium reduction column before NO3+NO2was detected as pink diazo compound. For PO4, it was determined through the molybdenum blue method whereas for SiO4, it was first reduced to molybdosilicic acid before its concentration was determined by the molybdenum blue method. TSS was measured by gravimetry method as the net weight increase of the filter after drying at 50 °C for seven days. The same filter was later combusted at 500 °C for 3 h and the weight loss was calculated as particulate organic matter (POM).

C hlorophyll a retained on the precombusted GF/F was extracted in 90% (v/v) ice-cold acetone at -20 °C for approximately 16 h. The fluorescence intensity was then detected at the excitation/emission wavelength of 440/680 nm with a spectrofluorometer (Perkin Elmer LS55 spectrofluorometer, USA). Standardization was carried out with a commercially available chl a standard (Sigma-Aldrich, USA) (Lohrenz et al., 2003). Seawater sample was also passed through the 2-μm and 20-μm (pore size) polycarbonate membrane filters, and chl a in the filtrate as collected on the GF/F filters represented the <2-μm and <20-μm fraction, respectively. The >20-μm fraction was calculated as the difference between total chl- a and the <20-μm fraction whereas the 2-20-μm fraction was calculated as the difference between the <20-μm and the <2-μm (Arin et al., 2002).

2.2 Phytoplankton production and grazing loss rate

In the Landry and Hassett (1982) dilution method, seawater sample was filtered through a 0.2-μm filter to prepare particle-free seawater. The particle-free seawater was then used to dilute the sample to 0.2, 0.4, 0.6, 0.8, and 1.0-strength (or undiluted) microcosms. These microcosms were then incubated under 12 h light (340 μmol/(m2·s)): 12 h dark photoperiod at in-situ temperature. The increase in chl- a concentration for each fraction ( μi) was calculated as μi=ln chl at-ln chl a0, where chl a0was the initial chl- a concentration and chl atwas the chl- a concentration at the end of the incubation. The μifor each microcosm was then plotted against its dilution factor, and a linear regression analysis was carried out. From the linear regression slope, the y-axis intercept was measured as μ whereas the gradient of the slope was measured as g. In this study, g was total grazing as we did not distinguish between micro- and meso-zooplankton grazing. However, it is generally accepted that microzooplankton grazing is the major factor for the loss of primary producers (Calbet and Landry, 2004; Schmoker et al., 2013).

2.3 Systematic review of waters around Peninsular Malaysia

In order to understand the drivers for phytoplankton growth and loss rates in waters around Peninsular Malaysia, we carried out a systematic review of published and available data from the NW ( n=18) and SE ( n=24) coast of Peninsular Malaysia (submitted), PD ( n=21), PK ( n=19) (Lim et al., 2015) and SS ( n=18) (Schmoker et al., 2016). For each station, the trophic state index (TSI) was calculated according to Carlson (1977) by the following equations:

TSI=60-14.41 ln(Secchi depth, m),

TSI=9.81 ln(chl a, μg/L)+30.6.

The average TSI values were then used to characterize these habitats into oligotrophic (TSI<40), mesotrophic (TSI=40 to 50), and eutrophic (TSI>50) waters.

2.4 Statistical analysis

All values were reported as mean±standard deviation unless otherwise stated. Least-square linear regression analysis was carried out for each dilution experiment whereas analysis of variance and Tukey’s test were used to compare values. Correlation analysis was used to show the relationships between the different parameters measured. For statistical analyses, the software PAST Version 2.17c (Hammer et al., 2001) was used.

3 RESULT

3.1 Environmental parameters

The environmental parameters measured at BMRS are shown in Table 1. The surface seawater temperature, salinity, DO and pH observed were relatively similar throughout the sampling period ( P >0.05), and ranged from 28.6 to 31.1 °C, 27 to 30, 214 to 329 μmol/L and 8.04 to 8.11, respectively. For TSS, the concentration varied from 103 to 156 mg/L, whereas for POM, the concentration fluctuated from 18 to 25 mg/L and represented 14% to 19% of TSS. POM was also higher on the first day of sampling relative to the second and third day ( F=9.756, df=5, P <0.05).

In term of nutrients, dissolved inorganic nitrogen (DIN=NH4+NO3+NO2) ranged from 1.89 to 2.75 μmol/L throughout the monitoring period at BMRS. Among the nitrogen species measured, NH4and NO3+NO2fluctuated from 0.21 to 1.38 μmol/L and 1.24 to 2.02 μmol/L (Table 1), respectively with NO3+NO2contributing 50% to 83% of the fluctuation in DIN. NH4was significantly different between the first and third day ( F=11.38, df=5, P <0.01), whereas NO3+NO2was higher in the first than the second and third day of sampling ( F=11.61, df=4, P <0.05). On the other hand, PO4and SiO4at BMRS varied from 0.05 to 0.71 μmol/L and 1.50 to 4.39 μmol/L, respectively. On average, PO4was significant higher on the second day of sampling ( F=7.534, df=6, P <0.05), whereas SiO4was not different throughout the observation period ( P> 0.05).

In this study, total chl- a concentration varied over two-fold at BRMS from 2.90 to 15.78 μg/L. Theaverage total chl- a concentration was however not significantly different throughout the observation period ( P >0.05) (Table 1). From the size-fractionation of chl- a experiment, we found that the dominant fraction was the 2 to 20 μm (50±11)%, followed by the >20-μm fraction (48±12)% and the <2-μm (2±2)% fractions. Through correlation analysis, the >20-μm fraction could have explicated (74±2)% ( F=998, df=11, P <0.001) of the fluctuation in total chl a, whereas the 2-20-μm fraction could have elucidated only (26±2)% ( F=128, df=11, P <0.001) of the variation in total chl a. The <2-μm fraction was a minor component and independent of total chl a ( P >0.05).

Table 1 Environmental parameters measured at BMRS (mean±standard deviation)

3.2 Phytoplankton growth and grazing loss rates

On all occasions, statistically significant regression slopes were obtained but one point (the undiluted fraction at 6 PM, Nov. 3, 2017) was removed (Fig.2). At BMRS, μ fluctuated daily from 1.02 to 1.58/d (coeき cient of variation, CV=14%), whereas g varied over a wider range from 0.07 to 0.88/d (CV=62%). Both μ ( P= 0.73) and g ( P=0.25) were not significantly different throughout the short-term monitoring at BMRS. Among the different chl- a size fractions, changes were apparent after the 24-h incubation. The >20-μm size fraction increased in proportion whereas the <2-μm size fraction decreased to near undetectable levels.

3.3 Systematic review of waters around Peninsular Malaysia

The environmental parameters for the different habitats are shown in Table 2. The average temperature at BMRS (29.6±0.7 °C) was lower than at NW (30.7±0.7 °C) ( F= 6.288, df=45, P <0.001), whereas salinity was higher at BMRS (28.9±1.0) than at PK (28.9±1.0) ( F= 20, df=46, P <0.001). On average, Secchi depth (0.3±0.1 m) was also lower at BMRSthan at NW (9.41±3.37 m), SE (6.83±2.37 m) and SS (3.50±0.71 m) ( F=78.34, df=9.63, P <0.001). In term of dissolved inorganic nutrients, NH4, PO4and SiO4were generally lower at BMRS than at PK ( F >9.23, df>40, P <0.001), except for NO3where the concentration was lower at BMRS (1.31±0.18 μmol/L) than PK (5.69±3.01 μmol/L) and SS (6.09± 9.10 μmol/L) ( F=21.79, df=48, P <0.001). Likewise, the estimated total phosphorus was about three times lower at BMRS (1.49±0.52 μmol/L) than at PK (4.07±3.70 μmol/L) ( F=10.68, df=54, P <0.001).

Table 2 The range of environmental parameters compiled from different habitats across Peninsular Malaysia

Fig.2 Linear regression analysis for all the incubation experiments carried out at BMRS

In term of total chl a, the concentration across different habitats typically ranged from 0.20 to 33.56 μg/L. On average, chl a was higher at BMRS (7.17±3.94 μg/L) than at PD (2.31±0.70 μg/L), NW (0.75±0.57 μg/L) and SE (0.72±0.32 μg/L) stations ( F=26.73, df=41, P <0.001) (Fig.3). From the TSI value, stations at the NW (27±6) and SE (30±4) coast of Peninsular Malaysia represented oligotrophic habitats, PD (49±4) and SS (57±5) were mesotrophic habitats whereas BMRS (63±3) and PK (52±6) were eutrophic habitats. Average μ was higher at BMRS (1.34±0.19/d) than other coastal water stations i.e. NW (0.46±0.20/d), SE (0.55±0.26/d), PD (0.52±0.25/d), SS (0.67±0.53/d), and PK (0.61±0.19/d) ( F (5, 106)=34.4, P <0.001) (Fig.3). In contrast, g at BMRS (0.44±0.27/d) was not different relative to other stations at NW, SE, PD, SS, and PK that averaged 0.41±0.22/d, 0.37±0.22/d, 0.53±0.23/d, 0.40±0.21/d, and 0.54±0.18/d, respectively (Fig.3).

4 DISCUSSION

4.1 Environmental parameters

Fig.3 Box-and-whisker plots showing the range and median of in-situ chl a (μg/L), phytoplankton growth rate ( μ, /d), grazing loss rate ( g, /d) and grazing pressure measured at BMRS, northwest (NW) and southeast (SE) coast of Peninsular Malaysia, Port Dickson (PD), stations at the Strait of Singapore (SS) and Port Klang (PK)

The temperature, salinity and pH observed at BMRS were high and constant (Table 2), typical of coastal waters around Peninsular Malaysia (You et al., 2016; Lim et al., 2018). For DO, the concentration observed at BMRS was generally at healthy levels (>200 μmol/L), and not below the 63 μmol/L consensus concentration for hypoxia (Rabalais et al., 2002). TSS measured at BMRS was also similar to other nearshore coastal waters (≥100 mg/L; Bong and Lee, 2008), and was mostly non-biogenic as POM was ≤19% of TSS. BMRS was characterized as eutrophic system, primarily due to the lower Secchi depth and higher total chl- a concentration observed (Fig.3) (Lim et al., 2015, 2018). As the sampling was not carried out during the possible upwelling period (June to September) (Syafrina et al., 2015; Suratman et al., 2017), the variation in Secchi depth and chl a were more likely influenced by nearby river output and nearshore physical processes (Shaari and Mustapha, 2017; Suratman et al., 2017; Lim et al., 2018). In addition, size-fractionation of chl a also concurred with other tropical studies where ≥52% of primary producers were >2 μm (Jagadeesan et al., 2017; Madhu et al., 2017).

4.2 Phytoplankton growth and grazing in BMRS

At BMRS, statistically significant linear regression slopes were obtained on all occasions, suggesting that nutrient was generally not limiting phytoplankton growth (Ayukai, 1996; Paterson et al., 2008). Through correlation analysis, we found that μ at BMRS was independent of temperature whereas g correlated with temperature ( R2=0.464, df=11, P <0.01). Although temperature observed at BMRS varied within a narrow range (≤2.5 °C), the dependency of g could be due to the intrinsic control of grazers through their metabolic rate or extrinsic control that alter the predator-prey encounter rates (Chen et al., 2012; Caron and Hutchins, 2013; López-Urrutia and Morán, 2015; Zhang et al., 2018). We also have evidence that temperature played an important role in regulating different microbial processes in this region (Lee et al., 2009, 2013). In addition, we also found that g at BMRS was significantly correlated with the initial chl- a concentration, specifically the 2-20-μm and >20-μm fractions (Fig.4). When we compared the regression slopes, g was clearly more dependent on the 2-20-μm fraction of chl a ( F=12.27, df=22, P <0.01). This observation was common as the 2-20-μm fraction of phytoplankton are the basic prey for the main grazers i.e. microzooplankton (>20 μm), accounting for 44% to 79% of g (Schmoker et al., 2013).

In this study, g correlated with μ at BMRS ( R2=0.349, df=11, P <0.05), and reflected the predatorprey coupling that is often observed (Calbet and Landry, 2004; Schmoker et al., 2013; Lim et al., 2015). The relationship between g and μ at BMRS could be expressed via the following linear regression equation: g=0.837 μ-0.685 ( F=5.350, df=11, P <0.05). From the equation, g could have accounted for 35% of the fluctuation in μ. The grazing pressure ( g/ μ) which estimates the fraction of primary production being grazed upon ranged from 0.07 to 0.63 (0.31±0.18) and was always below unity (<1), suggesting that primary production at BMRS was able to sustain the marine food web. The high chl- a concentration at BMRS was more likely influenced by external sources rather than localized production as the difference in intrinsic phytoplankton growth and grazing loss rate did not reflect in-situ chl a.

4.3 Systematic review of waters around Peninsular Malaysia

In this study, an extensive range of μ and g by the dilution method (Landry and Hassett, 1982) werecompiled from different coastal waters around Peninsular Malaysia. From the rate measurement, we showed that μ was significantly higher at BMRS whereas g was relatively constant throughout the waters of Peninsular Malaysia (Fig.5). In order to determine the possible environmental factors that drive μ, these data were clustered according to their TSI value (Table 2).

Fig.4 Correlation analysis between phytoplankton growth rate ( μ, /d) and grazing loss rate ( g, /d) to size-fractionation of chl a at BMRS

Fig.5 Correlation between chl a (μg/L) and phytoplankton growth rate ( μ, /d) for oligotrophic, mesotrophic and eutrophic waters around Peninsular Malaysia

For oligotrophic waters, μ was inversely correlated to temperature and suggested that an increasing temperature could have a negative impact on μ ( R2=0.240, df=41, P <0.001). For mesotrophic and eutrophic waters, we found that μ was significantly correlated to chl a (Fig.5). However, the relationship was only significant at <9-μg/L chl a ( R2=0.337, df=32, P <0.001) for mesotrophic waters and <16-μg/L chl a ( R2=0.384, df=28, P <0.001) for eutrophic waters (Fig.5). The concentration of chl a beyond this threshold was excluded as they are considered outliers (mean±2×standard deviation). Although the abundant nutrients in mesotrophic and eutrophic waters may have complicated this relationship (Yoshiyama and Sharp, 2006; Cloern, 2017; Zhang et al., 2018), the dissolved inorganic nutrients were all independent of μ even in the oligotrophic waters of NW and SE stations ( P >0.05). Our observations seemed to suggest that there was no nutrient limitation in these waters (Lim et al., 2015; Shaari and Mustapha, 2017). The correlation between chl- a concentration and μ could either reflect chl adriven phytoplankton growth or phytoplankton growth stimulating chl a (Laws et al., 1984; Ning et al., 1988; Tada et al., 2001) in these waters. In this study, we were unable to determine the effect of light on μ because of incomplete data, nonetheless the impact would likely be minimal as all samples were collected from surface (<1 m in depth) where light intensity is naturally higher (Chavez, 1989).

Regardless of the trophic state, we observed significant coupling between g and μ at every station. Their relationships are as follows: NW: g=0.704 μ+ 0.084 ( F=11.8, df=17, P <0.01); SE: g=0.661 μ+0.006 ( F=40.6, df=23, P <0.001); PD: g=0.679 μ+0.182 ( F=25.2, df=20, P <0.001); SS: g=0.191 μ+0.274 ( F=4.86, df=17, P <0.05); PK: g=0.499 μ+0.237 ( F=6.18, df=18, P <0.05) (Fig.6). The g/ μ slope at NW, SE, PD, SS, and PK accounted for 70%, 66%, 68%, 19%, and 50% of phytoplankton production, respectively. These observations suggested that g was generally an important loss process for u, except at SS when the dependency was lower than the 44% to 79% cross-compilation study by Schmoker et al. (2013). Despite the slope was significantly different ( F (5, 105)=3.94, P <0.01), the average grazing pressure ( g/ μ) was only lower at BMRS than NW, PD and PK ( F (5, 106)=19.4, P <0.01) (Fig.3). As there were no long-term measurement of μ and g in this region (≥5 years), all the data in this study were clustered together into four different seasons consisted of northeast (Nov. to Feb.), southwest (May to Aug.) and two shorter inter-monsoons (Mar. to Apr. and Sep. to Oct.) (Syafrina et al., 2015) to determine the possible seasonal effect on g/ μ. We found no significant difference for g/ μ ( P >0.05), and thus concluded the seasonal impact on g/ μ would likely be minimal.

Although the composition of organism was not observed, it may play a significant role in regulating the relationship between g and μ across different habitats. For example, the lower dependency of g on μ at SS could be due to the long residence time (about one year) as it was an enclosed bay (Schmoker et al., 2016). As a result, it could have influenced the abundance, food quality, species composition and encounter rates that in turn affected the habitat specific predator-prey relationship (Peters, 1994; Poulin and Franks, 2010). We also observed that all linear regression equations from the g- μ relationship had a positive y-intercept except for BMRS. This observation seemed to suggest alternative prey such as heterotrophic nanoflagellates or bacterial were available at other stations (Jürgens et al., 1996; Aytan et al., 2018). At BMRS, grazing was most probably solely dependent on phytoplankton growth due to the higher concentration of chl a observed. More studies are required to investigate this further.

Fig.6 Correlation analyses between grazing loss rate ( g, /d) and phytoplankton growth rate ( μ, /d) across different habitats around Peninsular Malaysia

5 CONCLUSION

From this study, we reported phytoplankton growth ( μ) and grazing loss ( g) rates at BMRS which revealed a relatively high μ. Systematic review of available data around Peninsular Malaysia showed that μ was sensitive to an increase in temperature in oligotrophic waters, whereas chl a distribution could drive μ in mesotrophic and eutrophic waters when the concentration was <9 μg/L and <16 μg/L, respectively. Dissolved inorganic nutrients were also not limiting across different trophic levels of waters in Peninsular Malaysia. We also found that g was significantly coupled to μ albeit at different degrees, and the availability of alternative prey for the grazers was also different among the stations, except at BMRS where g was solely dependent on μ.

6 DATA AVAILABILITY STATEMENT

All data generated and analyzed during this study are included in this published article.

7 ACKNOWLEDGMENT

We are grateful to the University of Malaya (RU006D-2014 and RU009C-2015) and the Minister of Higher Education (HiCoE Phase II Fund, grant IOES-2014D) for the grants that supported this work. We would like to thank two anonymous reviewers for their constructive comments to further improve this manuscript.

8 AUTHOR CONTRIBUTION

LIM J H, carried out the experiment and wrote the manuscript with data analysis support from LEE C W and BONG C W. LEE CW supervised the project and BONG CW aided in the collection of samples and interpretation of the results. All authors discussed the results and contributed to the final version of the manuscript.