Influence of different types of sessile epibionts on the community structure of mobile invertebrates in an eelgrass bed

Study area

The Akkeshi-ko estuary (locally called Akkeshi Lake) and Akkeshi Bay are located in Northeastern Hokkaido, Japan (Fig. 1) and are connected to each other through a narrow channel (width: approximately 500 m, depth: approximately 10 m). The Akkeshi-ko is a brackish estuary, shallow water (depth range in most of the lake: 0.8–1.7 m with the maximum difference in tide levels of up to approximately ±0.6 m), with an area of approximately 32 km². Most bottom areas of the Akkeshi-ko estuary are muddy and covered with eelgrass (Zostera marina) except for the aquaculture farms of the clam Venerupis philippinarum in the intertidal zone near the channel (Kasim & Mukai, 2006; Hasegawa, Hori & Mukai, 2007; Yamada et al., 2007a; Yamada et al., 2007b). Here, freshwater input from the Bekanbeushi River, which accounts for 98.8% of all of the flow volume (Iizumi et al., 1996), and tidal seawater input from the Akkeshi Bay cause steep physical and chemical environmental gradients (Iizumi et al., 1996; Yamada et al., 2007a).

Akkeshi Bay has an area of approximately 110 km² and opens to the Pacific Ocean at the south end. Two seagrass species Z. marina and Z. asiatica are present from the intertidal zone to the subtidal zone (5 m below mean low water); the former occurs at depths shallower than 2 m and the latter dominates in deeper water (Watanabe, Nakaoka & Mukai, 2005). The influence of the freshwater discharge on species composition of seagrass community is observed near the channel connecting the bay to the Akkeshi-ko estuary (Yamada et al., 2007a).

In this study, we established stations in the Akkeshi-ko estuary (BK: river mouth of the Bekanbeushi River, HN: Horonitai, TB: Toubai, SL: the southern lake, CL: the central lake and CK: Chikarakotan) and one station in Akkeshi Bay (SR: Shinryu) (Fig. 1). BK (mean sea level, MSL hereafter: 0.9 m) is located at the mouth of the Bekanbeushi River and is strongly affected by freshwater inflow. The vegetation is dense with small-sized Z. marina (average shoot length of 1.0 m in August). HN (MSL: 1.1 m) is in a location with a high water temperature and medium salinity relative to the other stations. In addition to Z. marina, Ruppia maritima, a seagrass species that is more tolerant to low-saline water, occurs at HN. The eelgrass beds at HN are mostly continuous but have some gaps, and the average shoot length in August is 1.3 m. TB (MSL: 1.1 m) and SL (MSL: 1.0 m) have a relatively low salinity compared to that of the other stations and are the furthest stations from the Akkeshi Bay. Although these two stations are in a similar environment, the water is often more turbid and the eelgrass bed is patchier at TB than SL. SL has a higher seagrass biomass and shoot density than TB. The average shoot length of eelgrass is approximately 1.3 m in August at both of these stations. CL (MSL: 1.4 m) and CK (MSL: 1.5 m) are deeper stations with a higher salinity and are dominated by longer eelgrass (shoot length: 1.5–3.5 m at the peak season). The eelgrass at SR (MSL: 1.5 m) in the Akkeshi Bay, has a similar shoot size to that in CL and CK. Here, the dominant seagrass species changes from Z. marina bed to Z. asiatica at a depth of approximately 2 m, as mentioned above.

According to Yamada et al. (2007a), salinity varies significantly among stations but does not vary greatly among seasons. During the summer (from July to August), eelgrass biomass, microalgal biomass and mobile invertebrates reach their peak (Hasegawa, Hori & Mukai, 2007; Yamada et al., 2007b). Seasonal changes in the mobile invertebrate species richness are not clearly understood (Yamada et al., 2007b).

Field sampling

We conducted a field survey in August 2012. Sample collection was performed when the tidal current was slow. We collected mobile invertebrates on eelgrass blades when the water level was deeper than the average sheath length of the eelgrass at each station (BK: 20 cm; HN, TB, SL: 30 cm; CL, CK, SR: 40 cm). Because the eelgrass at our study stations is tall (>1 m) compared to the average water depth of each station, the canopy usually reaches the water surface (except for at extremely high tides). All samplings were performed under these conditions. We targeted mobile invertebrates but excluded some species with remarkably high mobility and low dependency on eelgrass habitat, such as mysids and decapods (Yamada et al., 2007b), which were not quantitatively collected by our method (see below).

We measured water temperature and salinity once at each station using a memory sensor (AAQ-175 RINKO; JFE Advantech Co. Ltd., Japan). To obtain the representative values, the sensor was carefully placed approximately 50 cm from the bottom to accurately reflect the environment inside of the seagrass meadow.

We collected three replicate samples (a total of 21 samples from all stations) of mobile invertebrates, spirorbid polychaetes and epiphytic macroalgae together with the entire above-ground parts of the eelgrass using a mesh bag (bore diameter: 20 cm, mesh size: 0.1 mm) based on the mouth area of the mesh bag (314 cm²). Upon collection, we counted the number of eelgrass shoots to determine shoot density. For microalgae, five replicate samples were collected per station, using separate plastic zip bags for each eelgrass shoot, because microalgae easily fell off from eelgrass blades when collected with the mesh bag.

Laboratory procedures

Immediately after being transported to the laboratory, the microalgae were scraped from the eelgrass blades using a glass slide; separated from other organisms such as macroalgae and spirorbid polychaetes; and then filtered using glass fiber filters (Whatman GF/F filter φ 47 mm; Whatman International Ltd., Maidstone, UK). If other organisms were present in the microalgal samples, we carefully removed them from the filters with forceps. Other epibiotic organisms collected using mesh bags were separated from the eelgrass by scraping them off with a glass slide; these organisms were classified as red algae, green algae, spirorbid polychaetes and mobile invertebrates. To obtain dry mass, eelgrass shoots, red algae, green algae, spirorbid polychaetes and filtered microalgae were dried at 60 °C for 4 days in small aluminum foil bags, and then weighed. We counted and identified the mobile invertebrates after extraction with a sieve (500 µm) and fixation with 70% ethanol. Identification of mobile invertebrates was made to the lowest taxonomical unit possible (mostly to species) using detailed guides from the literature (Gammarid amphipod: Nishimura, 1995; Carlton, 2007, Caprella amphipod, Isopod, Copepod, Cumacea: Nishimura, 1995; Carlton, 2007, Gastropod: Okutani, 2000, Polychaeta: Nishimura, 1992; Imajima, 1996; Imajima, 2001, Turbellaria: Nishimura, 1992, Carlton, 2007, Hirunoidea: Nishimura, 1992) and the World Register of Marine Species online database (WoRMS: http://www.marinespecies.org).

Statistical analysis

We used, as predictors, two abiotic factors (water temperature and salinity) and six biotic factors (eelgrass biomass (g dry weight per unit area: g DW m⁻²), eelgrass shoot density (shoots m⁻²), microalgal biomass (g DW m⁻²), red algal biomass (g DW m⁻²), green algal biomass (g DW m⁻²) and spirorbid polychaete biomass (g DW m⁻²)). For eelgrass biomass, we used the dry weight data collected using mesh bags. Because microalgal biomass was collected by a different sampling procedure from other biotic variables, we used the mean value of five replicates. In this study, one of our interests was the effects of morphological traits of macroalgae and spirorbid polychaetes. Thus, we separated red and green algae by a morphological trait (red algae: branched, green algae: filamentous). All invertebrate biomass (mg ash-free dry weight per unit area: mg AFDW m⁻²) was estimated from the abundance and the size fraction using the empirical equations in Edgar (1990).

To test which of the eight biotic/abiotic factors was a likely predictor of the variation in the mobile invertebrate community, we fitted linear mixed models (LMMs) with a Gaussian distribution (Bolker et al., 2009). The station was used as a random variable. As response variables, we used the biomass of 11 dominant species for the population-level analyses, and total invertebrate biomass, species richness and species diversity (Shannon-Wiener diversity index; calculated based on biomass data) for community-level analyses. The 11 most dominant species were selected by a threshold whereby the biomass proportion accounted for more than 1% of the total invertebrate biomass (see Table S1). Ostreobdella kakibir (Hirudinea) was omitted from the analysis because it occurred only at one station (SR), even though they satisfied the requirement. R software (version 3.1.3) was used for all of the analyses (R Development Core Team, 2015).

Prior to the LMM fit, all of the variables excluding species diversity were square root transformed to improve homoscedasticity and meet the assumptions of normality of the LMMs after checking for normality with the Shapiro–Wilk test. To test for collinearity between the eight environmental variables, we calculated Pearson’s correlation coefficients for all pairs. If the absolute value of the coefficient (r) was greater than 0.7, the level where collinearity does not affect model predictions (Dormann et al., 2013), we removed the relevant predictor as necessary. Because water temperature and microalgal biomass were highly correlated (Pearson’s r =  − 0.82, P < 0.01), we removed microalgal biomass from the models. After this removal, we tested potential multicollinearity among the remaining predictors using the variance inflation factor (VIF) analysis with a cutoff of 10 (e.g., Dormann et al., 2013). VIF values were calculated using the vif.mer function developed by Frank (https://raw.githubusercontent.com/aufrank/R-hacks/master/mer-utils.R). However, all seven predictors were below the VIF value of 10 and remained. We therefore defined a reduced model with the seven predictors as the full model.

We fitted the LMMs using the lmer function in the lme4 package (Bates et al., 2014). To obtain P-values of the LMMs, we used the lmerTest package (Kuznetsova, Brockhoff & Christensen, 2014). We selected the optimal model comparing the candidate models on all combinations of the predictors by the Akaike information criterion as corrected for the small sample size (AIC_c: Burnham & Anderson, 2002). We obtained AIC_c based on the maximum likelihood (ML) for comparisons among the LMMs because the restricted maximum likelihood (REML) is inappropriate in the case when the fixed structure is different between the candidate models (Zuur et al., 2009), but the parameters were estimated by REML. We used the AIC_c tab function in the bbmle library (Bolker & R Development Core Team, 2013) to compare the AIC_c. After setting the optimal models, we obtained the standardized coefficients as effect sizes by re-fitting using standardized variables that were scaled by the sample standard deviation and centered by sample mean values.

Additionally, when the effect of water temperature was detected, we tested the relationship between mobile invertebrates and microalgal biomass which was omitted from the LMM because of the multi-collinearity with water temperature.

Environmental factors

Water temperature was lower at the four stations (BK, CL, CK and SR) near the channel than at the other three stations in the inner parts of the estuary (HN, TB and SL) (Table 1). Salinity was lower at the lake-side stations (BK, HN, TB, SL and CK) that were influenced by freshwater inputs. For these stations, the inter-annual variation was also higher as shown by data collected by ourselves and other studies (Table 1).

Eelgrass biomass varied between 140 and 278 g DW m⁻² among the stations. It was lowest at TB, followed by HN and BK (Table 1). Eelgrass shoot density ranged between 85 and 234 shoot m⁻². It was highest at BK and second highest at SL. The mean densities were not largely different among other stations. Microalgal biomass varied by more than ten-fold between the lowest station (SL) and the highest station (CL). In the latter, the microalgal biomass exceeded the biomass of the eelgrass. Macroalgae were not present at TB and SR. Branched red algae were dominated by Neosiphonia sp. and Chondria dasyphylla, and filamentous green algae were dominated by Cladophora sp. The mean biomass of red algae was highest at HN and that of green algae was highest at CL, though their biomasses were less than 15% of that of eelgrass. Spirorbid polychaetes were not present at SL and CL. They were highly abundant at SR where their biomass was more than three-fold greater than the eelgrass biomass.

Mobile invertebrate community

A total of 32 mobile invertebrate species were collected in this study (Table S1). At taxonomic levels, polychaete worms made up 32.2% of the total biomass, followed by gastropods (31.3%), gammarid amphipods (23.0%), and isopods (8.8%). At the species level, a polychaete Nereis sp. was most dominant (24.6%), followed by gastropods Lacuna spp. (23.4%) and the gammarid amphipod Ampithoe lacertosa (18.0%). For an additional eight species including two gammarid amphipods (Monocorophium spp. and Pontogeneia rostrata), two isopods (Cymodoce japonica and Paranthura japonica), two gastropods (Ansola angustata and Siphonacmea oblongata) and two polychaetes (Exogone naidina and Syllis sp.), their proportions were less than 5% at most.

The mean value of the total mobile invertebrate biomass was the highest at CK and much lower at stations along the coastline (HN, TB and SL). Species richness was the highest at CL, followed by CK and was approximately the same level at the other stations (Fig. 2). The mean value of species diversity was the highest at CL and the lowest at SR (Fig. 2).

Population level analyses

We found that each of the nine invertebrate populations belonging to gammarida, gastropoda and polychaeta was predicted by a different combination of environmental factors in the optimal models (Table 2). The effect size of three epibionts on dominant invertebrate species was either similar or larger than abiotic and eelgrass factors (Fig. 3). For two isopods, no environmental factor correlated with their biomass.

Water temperature was selected as the responsible factor for the variation of A. lacertosa, Lacuna spp. and all three polychaetes. Among them, only Syllis sp. showed a significant correlation (positive). The significant effect of the salinity gradient was detected for A. angustata (negative) and S. oblongata (positive).

For the two predictors relevant to the characteristics of the eelgrass bed, the above-ground biomass showed a significant positive relationship only with Syllis sp., whereas shoot density was significantly correlated with Monocorophium spp. (positive), P. rostrata (negative) and E. naidina (negative) (Table 2).

The biomasses of sessile epibionts (red algae, green algae and spirorbid polychaetes) on eelgrass blades were correlated with many invertebrate populations in different manners, excluding A. lacertosa, Monocorophium spp., two isopods and E. naidina. Red algal biomass was positively correlated with P. rostrata, Lacuna spp. and Nereis sp., but negatively correlated with Syllis sp. and tended to be negatively correlated with A. angustata. Green algal biomasses were positively correlated with P. rostrata and Lacuna spp. The biomass of spirorbid polychaetes was positively correlated with all three species of gastropods and was negatively correlated with Syllis sp.

Although epiphytic microalgae were removed from our analysis because of the collinearity with water temperature, no significant correlation was found for species that were correlated with water temperature (A. lacertosa: Pearson’s r = 4.05, P = 0.25; Lacuna spp.: r =  − 3.42, P = 0.69; E. naidina: r = 0.31, P = 0.92; Nereis sp.: r = 10.87, P = 0.34; Syllis sp.: r = 0.21, P = 0.86).

Community level analyses

The total invertebrate biomass tended to decrease with water temperature, and significantly increased with increasing eelgrass biomass and red algal biomass (Table 2). Species richness showed a negative correlation with water temperature and spirorbid polychaetes. Species diversity was positively correlated with green algal biomass, but was negatively correlated with spirorbid polychaetes (Table 2). The effect size of red algal biomass on the total invertebrate biomass was similar to that of eelgrass biomass, and that of spirorbid polychaetes on species richness was also similar to that of water temperature (Fig. 3).