Potential impacts of marine urbanization on benthic macrofaunal diversity

Study area

We established field survey sites in four different habitat types (BW: breakwater wall; EB: eelgrass bed; IF: intertidal flat; SB: subtidal bottom) within Matsunaga Bay, located in the eastern Seto Inland Sea, Japan (Fig. 1). Matsunaga Bay is a small semi-closed bay with an area of approximately 12 km². It is connected with other water bodies through the Tozaki-Seto Strait (width: approx. 400 m) and the Onomichi Strait (width: approx. 200 m)^29,30. Water depths are mostly less than 20 m throughout the bay. The water depths at our four survey sites were approximately 4.5 m, but part of site SB located near a shipping channel reached depths of 10–13 m. Intertidal flats cover 35% (4.3 km²) of the bay area, whereas eelgrass beds cover 1.7% (0.2 km²)³¹. The bottom sediment type is mainly muddy throughout the bay, although some parts of EB and SB have sandy and muddy bottoms (see Supplementary Table S1)³⁰.

Although human activities along the coast of Matsunaga Bay appear to be associated with artificial structures (e.g., industrial plants, marinas, and lumber yards), natural habitats are still relatively well-preserved compared to areas in the eastern Seto Inland Sea²⁹. The total population of towns within 5 km of the coastline of the bay can be expected to exceed 100,000 people, which ranks within the top 20% of administrative districts in Japan, including prefectures, towns, villages, and the 23 wards of Tokyo³².

Data collection

We conducted one field survey in September 2016 (summer) and another in January 2017 (winter) to collect data on benthic community structures and environmental conditions in each habitat. We established five sampling sites within each habitat to obtain replicated samples. To reduce biases due to tidal cycle, we performed all field sampling and measurements at high tide, when all habitats were underwater.

We used a standard sample area (approximately 0.15 m² per sample) at each sampling site irrespective of sampling method to obtain comparable data on benthic communities. Some of our sampling methods involved Smith–McIntyre grab samplers and quadrats that could not cover the standard sample area in a single sample; for these methods, we combined the data from three samples to make up a single sampling site. At BW, benthic macrofauna (hereafter referred to as “benthos”) samples were collected by SCUBA divers. We established 15 sampling positions in a 5 × 3 grid (i.e., five depths [sampling sites] and three replicates) based on distance from the seafloor at the breakwater wall (Supplementary Fig. S1). The SCUBA divers used scraper blades, 0.1-mm mesh bags, and 22.5 cm × 22.5 cm quadrats because the benthic communities were mainly composed of sessile organisms. At EB, we employed different sampling methods for the above- and belowground components. SCUBA divers collected aboveground samples of eelgrass-associated benthos and eelgrass shoots using a mesh bag (mesh size: 0.1 mm; bag diameter: 45 cm). They then cut away the eelgrass shoots near each aboveground site and collected belowground samples of the benthos on top of and within the sediment by using the bucket part of a Smith–McIntyre grab sampler (sampling area: 22.5 cm × 22.5 cm). At IF and SB, benthos and bottom sediment samples were collected from a ship using a Smith–McIntyre grab sampler (sampling area: 22.5 cm × 22.5 cm).

All benthos was extracted using a 1-mm sieve and preserved in buffered 10% formalin in the field as soon as possible after sampling. The samples were identified to the lowest possible taxonomic unit and counted in the laboratory. After identification, we organized the dominant benthic species (or taxa) according to their primary feeding types and common life forms with reference to the World Register of Marine Species (http://www.marinespecies.org/) and the literature. No vertebrate species were targeted in this study. We defined the primary life forms of adult benthic species on/in their substrates as “common life forms.”

Although differences in environmental conditions were not the focus of our study, we did assess whether there were considerable water quality differences among sites. The purpose of this assessment was to try to identify sites with similar conditions so that exogenous impacts on biological communities could be discounted as much as possible in the analysis. We measured water and sediment conditions at each sampling site (except at BW, where sediment conditions were not measured due to the absence of sediment). Prior to benthos and sediment sampling, we measured water temperature, salinity, pH, and dissolved oxygen concentration at each site at a depth directly above the seafloor by using a multi-parameter water quality meter (AAQ‐RINKO, JFE Advantech Co. Ltd., Japan). At BW, where the substrate (i.e., the breakwater) is oriented vertically (see Supplementary Fig. S1), we measured water conditions in the middle of the water column. We also measured temperature, pH, oxidation–reduction potential (ORP), water content, and median particle size (D₅₀) in the sediment. Sediment temperature, pH, and ORP were measured by using a portable ion meter (IM-32P, DKK-TOA Co., Japan) immediately after each sample was collected. Sediment water content and D₅₀ were measured in the laboratory once benthic species had been removed from the sample.

Data analysis

First, we identified how many species were shared between all habitat pairs to understand inter-habitat species-sharing relationships. Second, we compared species compositions and abundances among habitats using similarity indices and multivariate analysis (described below).

To detect species sharing in terms of species commonality and endemism among the four habitats, we classified benthic species into the following three categories: common, endemic, and shifting. Common species were defined as species that occurred across all habitats. Endemic species are those that were found in only one habitat. Shifting species were defined as those that occurred in multiple habitats (but not across all habitats) and therefore showed a broad allowable range of habitat types or conditions. To analyse the importance of habitat sets in maintaining local species diversity, we further categorized the shifting species into two- or three-habitat users (i.e., those that occurred in two or three different habitats). Moreover, we calculated the numbers and proportions (i.e., using the Jaccard similarity index) of shared species in each habitat pair to evaluate the potential strength of any inter-habitat relationships. The Jaccard similarity index ((J)) is calculated as follows:

where (S_{alpha }) is the number of species in habitat (alpha), (S_{beta }) is the number of species in habitat (beta), and (S_{alpha beta }) is the number of species that are shared among habitats (alpha) and (beta).

In terms of the functional groups, we analysed abundance matrices of abundant species grouped by primary feeding types and common life forms. Focusing on abundant species is a useful way to reflect the functional characteristics of biological communities¹⁴. Therefore, we identified the most abundant species from each sample before constructing the abundance data matrices. To determine how many species to select for analysis, we calculated the number e of equally-abundant species that would be required to obtain the Simpson diversity index of each community (i.e., the effective number of species³³). We then selected e abundant species from each sample in rank order from most to least abundant. If multiple species of the same rank occupied this cut-off threshold, we selected all of them. This selection method, which is unique to our study, successfully identified dominant species that accounted for over 70% of the total abundance in each community in all habitats.

To compare benthic macrofaunal communities among habitats, we performed non‐metric multidimensional scaling (NMDS)³⁴ based on dissimilarity matrices obtained by using the metaMDS function in the vegan package³⁵. To compare species compositions, we constructed Jaccard dissimilarity matrices based on presence/absence data, and to compare species abundance and functional compositions, we constructed Bray–Curtis dissimilarity matrices based on abundance datasets. All abundance datasets were square-root transformed before calculating the Bray–Curtis dissimilarity matrices to reduce the influence of abundance bias. We accepted the NMDS ordinations if stress values were less than 0.2 to maintain the accuracy of the two‐dimensional ordinations³⁴. Then, we tested the effects of habitat type and sampling time by conducting a two-way permutational multivariate analysis of variance (two-way PERMANOVA)³⁶ using the adonis function in the vegan package. Here, we considered four habitats (Habitat), two sampling times (Time), and their interaction term as explanatory factors. Although our main focus was differences in community compositions among habitats at each sampling time, we also examined the magnitude of variation in each habitat by comparing two stable seasonal extremes (i.e., summer and winter). If the results of the PERMANOVA were significant, we performed post-hoc tests (pairwise PERMANOVA) to identify which pairs of community structures were significantly different by using the pairwise.perm.manova function in the RVAideMemoire package³⁷. We used 9999 permutations for the NMDS ordination, PERMANOVA, and pairwise PERMANOVA. P-values calculated during the pairwise PERMANOVA were corrected using the false discovery rate method³⁸. For the benthic community data collected at each sampling point at EB, above- and belowground datasets were integrated to reflect spatial representativeness (see Figs. S3 and S4). All analyses were performed using R version 3.5.1³⁹.

Source: Ecology - nature.com