We classify stratigraphic bounding surfaces defined by sedimentologic, palaeoecologic or taphonomic, and ichnologic criteria. Our case study refers to a particularly well-preserved outcrop of Miocene sedimentary rocks on the Pacific Coast of central Costa Rica. These rocks (mainly fine sandstones) accumulated in a tropical, shallow-marine, strongly subsiding, sand-swamped forearc coastal embayment. Sedimentary bounding surfaces (or lithofacies-bounding surfaces) can be subdivided into erosional and non-erosional boundaries between depositional systems. In our case study erosional boundaries (angular unconformities and erosional disconformities) were produced by block tilting, tectonic reactivation of cliffs, or wave and tidal ravinement during sea-level rises, and are interpreted as sequence boundaries and bases of parasequence sets, respectively. Non-erosional transitions delimit genetically closely related depositional systems. Sharp contacts are considered as bases of parasequence sets, diffuse transitions are interpreted as parasequences.
Palaeoecologically and taphonomically defined fossiliferous bounding surfaces or layers mostly occur within depositional systems; they seldom coincide with erosional lithofacies-bounding surfaces. Two types can be distinguished: autochthonous and allochthonous fossil concentrations. Autochthonous fossil concentrations can be grouped into layers with epifaunal and semi-infaunal molluscs in growth position, layers with infaunal bivalve associations in growth position, and coprolith concentrations. The majority of these layers is developed as concretion horizons. Autochthonous fossil concentrations invariably signal a relative rise of sea level and are interpreted as bases of parasequences or parasequence sets. Allochthonous fossil concentrations can be grouped into fossil concentrations produced by non-turbulent shelf currents, driftwood concentrations, tempestite fossil concentrations, fossil concentrations produced by near-shore turbulent flow regimes, and fossil concentrations produced by tidal currents. They are only seldom developed as concretion horizons. The appearance of allochthonous fossil concentrations in certain intervals of a given section, and disappearance in others, is clearly a signal of shallowing or deepening seas. Fossil concentrations produced by non-turbulent shelf currents and driftwood concentrations can be interpreted as bases of parasequences.
Ichnofabrics can be used as an additional tool for the identification and interpretation of bounding surfaces. Erosional boundaries may be marked by well-known features such as Gastrochaenolites. Sequence boundaries may be preceded by Thalassinoides boxworks strongly thickened by concretionary growth. At non-erosional boundaries, however, trace fossils commonly fail to mark the crucial turning point. Yet they seem to record environmental changes (for instance, a coming transgression) far earlier than any sedimentary phenomenon or body fossil: after prolonged swamping with sand during a relative sea-level lowstand, they ‘explode’ in diversity during the following transgression before this transgression produces its own typical sediment. Equally, after the dumping of organic matter during the transgression trace-making organisms announce the coming highstand by reworking their habitat to the degree of total destratification, before the highstand deposits proper force them to keep up with aggradation, which results in the typical highstand selection of deep tier structures. Thus, the degree of reworking by bioturbation seems to be a very good and easily recognizable indicator of the internal architecture of parasequences and parasequence sets.
Estuaries are located at the boundary between land and sea, and are a particular environment influenced by periodic tides and waves [1,2]. Estuaries provide multiple ecological services, such as nutrient cycling, climate change adaptation, and function as habitats and spaces for recreation . Costanza et al. (1997) estimated that the total value of annual ecosystem services of estuaries is $22,832 ha−1·yr−1, which is among the highest of 21 biomes . In contrast, estuaries were reported to be highly degraded by human activities because of their rich biodiversity, high levels of nutrients in the land, and abundant natural resources [5,6]. Approximately four billion people now inhabit land within 60 km of the world’s coastlines, and they have placed considerable pressure on estuaries . Anthropogenic impacts induce changes of estuary environments, such as habitat loss, deterioration of water quality, and degradation of resources [7,8,9,10,11,12].
Knowledge about the relationship between estuarine biota and habitat characteristics is important to establish appropriate conservation plans. A multidisciplinary approach was adopted to evaluate estuarine ecological integrity in the United States using the indicators of water quality, bed material, and habitats . Ecological conditions were evaluated using the Benthic Quality Index, based on the recognition of regional benthic fauna . In Australia, a conservation plan for each estuary has been designed and conservation actions have been progressing. In the Australian system, estuaries are classified by the strength of impact of the energy of their tide, waves, and rivers [15,16]. In this Australian system, estuaries are categorized into six main types (tide-dominated, tide-dominated estuary, tidal flat, wave-dominated delta, wave-dominated estuary, and strandplain), and potential habitat types associated with each estuary type are suggested . Research on the relationships between fauna and estuary type and reference conditions for each estuary type was also conducted . In European countries, a Water Framework Directive (WFD) is used to define a framework for the protection of European waters in order to reach “good status” objectives for water bodies throughout the EU ; many studies have also evaluated the environmental condition of estuaries to achieve their integrated management . Phytoplankton or benthic fauna have been used as biological indicators [21,22] and physicochemical indicators have also been used to evaluate the condition of estuaries [23,24,25]. In addition, Galván (2010) classified Northeast Atlantic estuaries using hydrological and geological data and revealed that transitional waters with a complex morphology showed the highest values of species diversity, while those with a smaller tidal influence showed lower species diversity .
A conservation plan for estuaries should be based on knowledge of an estuary’s classification results based on biological or physical/environmental factors. Many studies have thus established systems for classifying estuaries; their geography, physical environment, and biota are representative of the indicators that have been used. Pritchard (1952) focused on estuarine geography and established three classification types (coastal plain estuaries, fjords, and bar-built estuaries) . In contrast, Davies (1964) classified estuaries into microtidal, mesotidal, macrotidal, and hypertidal, in accordance with the difference of tidal levels . In addition, Fairbridge (1980) presented a comprehensive classification method using geomorphological history, river discharge of water and sediment, tidal current and waves, and coastal processes . Darlymple et al. (1992) further developed the classification concept of Fairbridge and presented the estuarine habitat type based on river energy, wave energy, and tide energy . On the other hand, fish fauna was used for classifying South African estuaries and the classification results were associated with water temperature or geological formation [30,31,32,33,34]. Colloty et al. (2002) categorized 92 estuaries of Eastern Cape Province using 54 plant species and suggested a relationship between the classification results and topographical factors (permanently open and temporarily open/closed) .
The examples of conservation efforts or research tools developed in the US, Europe, Australia, and South Africa mentioned above are mainly intended for estuaries of large rivers or coastal zones. However, few studies or conservation plans have been established for river estuaries of small and medium-sized rivers. In Japan and Southeast Asia, where the land is composed of many peninsulas or islands, there is an extended and complex coastline, geographically variable tides, and rivers of various sizes. Therefore, there are many types of estuaries according to the complex geography or physical condition of these areas, and the relationship between the ecological system and the physical environment has not been revealed. Furthermore, the anthropogenic impact on estuaries is extremely high in Japan because of the high population density in lowland areas which means that the reference condition is difficult to determine. Therefore, to establish conservation strategies for the great diversity of existing estuaries, it is important to comprehend the environmental condition by drawing comparisons with other estuaries. Furthermore, predicting the habitat or biota occurring in the river estuaries from physical factors enables the establishment of specific targets for environmental restoration.
Against this background, the aims of this study were to reveal the relationship between physical factors and the composition of the biotic community. To achieve this objective, we investigate the relationship between physical indicators and molluscan fauna and compare the classification results of physical indicators and molluscan fauna and consider the cause of non-correspondence between physical indicators and molluscan fauna.
We selected molluscan fauna to evaluate the integrity of the river estuaries. Molluscan fauna respond sensitively to water quality or bottom sediment and include species that inhabit only one particular environment or have a low capacity to thrive in different habitats. Molluscan species at individual locations directly reflect the environmental conditions at these sites . Therefore, molluscan species are ideal for evaluating the environmental conditions or determining the impact of human activities [37,38,39].
2. Materials and Methods
2.1. Study Area
The present study focused on 19 river estuaries in Kyushu, Japan (Figure 1). River estuaries are characterized by large environmental gradients of water quality, riverbed material, and microtopography in the longitudinal and transverse directions . We selected river estuaries for study here by including those with a variety of physical environments.
2.2. Molluscan Fauna Data
Molluscan fauna collections were performed at low tide at the middle and spring tides from 28 April 2015 to 22 January 2016. The molluscan species were sampled at a dried out habitat one reach (about ten times of the river width) from the mouth. We defined the land of low flow channel as a habitat and set one to three sampling points at each habitat according to the habitat area. Eight kinds of habitats were set from the viewpoint of particle size of the sediment, vegetation, and artificial structure (silt, sand, gravel, boulder, bedrock, riprap, concrete construct, and vegetation). We collected molluscan (Bivalvia and Gastropoda) fauna in a region of a 50 cm square and 10 cm deep at each onshore habitat in a stream reach.
2.3. Physical Environment Data
We collected riverbed material at each molluscan sampling site. The grain size accumulation curve was obtained by a sieve analysis test, as specified in JIS Z 8815 , for the grain size range of 3.0–75 mm. For grain sizes greater than 75 mm, we measured the short diameter (ds) and long diameter (dl) of 50 pebbles and calculated the grain diameter as (ds × dl)0.5. When the sampled material did not contain particles larger than 3.0 mm, the grain size accumulation curve was calculated by laser diffraction particle size analysis (SALD, 2000). In addition, the morphology of each habitat and the cross section and longitudinal shape of the river were surveyed by Real Time Kinematic-Global Positioning System (RTK-GPS). The habitat planform and cross-sectional profiles of each river were created by the measurements.
2.4. Calculation of Physical Indicators
In this study, we focused on physical indicators at the watershed scale and the habitat scale. The indicators of river energy, tide energy, and wave energy were adopted for the watershed scale, and the indicators of river bed material and habitat type were adopted for the habitat scale. Specific indicators were as follows: indicator of tide energy: (1) tidal range; indicators of wave energy; (2) wave exposure; (3) direct fetch; indicators of river energy; (4) topographic gradient; (5) form ratio; (6) specific discharge of occurrence probability 1/5, and indicators of habitat scale; (7) silt; (8) sand; (9) gravel; (10) boulder; (11) bedrock; (12) riprap; (13) concrete construct; and (14) vegetation.
(1) Tidal range was calculated by the average of the difference of monthly maximum and minimum tide levels for the period from January–December 2014. The tidal data were obtained from the Japanese Coastguard. To represent the tide at the investigated site, the recorded value from the nearest Japanese Coastguard observation point was adopted; (2) Wave exposure was quantified using the Baardseth Index . To calculate this index, the center of a transparent circular disc with a radius of 15 cm (representing 7.5 km) was placed at the investigation site on a 1:50,000 scale chart. The disc was divided into 40 sectors, with the angle of each sector being 9°. Sectors containing peninsulas, islands, or parts of the mainland shore were ignored ; (3) We calculated direct fetch by integrating the distance to the opposite shore from the investigation point, extending the radius in each direction with a width of ±22.5° from the center point . A maximum distance of 200 km was set as the distance to the opposite shore, in accordance with the concept of wave height saturation ; (4) The topographic gradient was obtained by dividing the altitude of the headstream by the length of the main river; (5) The form ratio was calculated by dividing the basin area by the square of the length of the main river; (6) The specific discharge of occurrence probability 1/5 was calculated by a rational run-off formula using a rainfall intensity formula released by each river’s administrator; For (7) silt; (8) sand; (9) gravel; and (10) boulder; based on the grain size accumulation curve of riverbed materials collected at each habitat, we defined silt as an average grain diameter of 0.005–0.075 mm, sand as an average grain diameter of 0.075–2.0 mm, gravel as an average grain diameter of 2.0–75 mm, and boulder as an average grain diameter greater than 75 mm; For (12) riprap and (13) concrete construct, when the function of a habitat was changed between a concrete construct or stone material, we defined large stones acting as coastal defenses as riprap and bed beaching as a concrete construct. We used presence-absence data of the indicators (7–14).
2.5. Statistical Analysis
To identify the possible groupings of similar molluscan fauna between river estuaries, we conducted non-metric multidimensional scaling analysis (nMDS) to summarize the composition of molluscan fauna using the Bray–Curtis similarity. nMDS was conducted for 14 physical indicators and species that appeared in over three rivers, to exclude the influence of species appearing at a low frequency. To approach a normal distribution, the number of individual molluscs was used after logarithmic conversion [log(e + 1)]. As a result of the permutation test, physical indicators (p < 0.05) were shown as vectors. Additionally, indicator species of each group were obtained using the indicator value method (IndVal) . Secondly, molluscan fauna and physical indicators were classified by hierarchical cluster analysis (Ward method). Euclidean distance was used for calculating the distance between objectives. To investigate the difference of a physical indicator among groups, the average values of physical indicators of each group were analyzed with the Kruskal–Wallis test and Steel Dwass test. These analyses were conducted using the statistical analysis software R.
3.1. Molluscan Fauna Survey
In total, 27 families, 55 species, and 6003 individuals were collected in 19 river estuaries (Table 1). The species for which the highest number of individuals were collected was Batillaria multiformis, for which a total of 1105 individuals were confirmed in 12 rivers. At the other end of the spectrum, species belonging to Bivalvia, such as Nuttallia commode, Moerella iridescens, and Laternula boschasina, were confirmed in a few rivers. Regarding the numbers of species and individuals in each river, the Hai River was associated with the highest numbers of species and individuals (17 species and 1249 individuals). In contrast, no molluscan species were found in the Obaru River.
3.2. Relationship between Molluscan Fauna and Physical Indicators
Figure 2 indicates the plotted molluscan species and physical indicators on nMDS dimensions. As a result of the permutation test, tidal range, direct fetch, silt, sand, boulder, concrete construct, and vegetation were selected. Littoraria intermedia, Cerithidea rhizophorarum, Assiminea sp., or Cerithidea cingulate were plotted in the first quadrant, characterized by a large tidal range and abundant silt and vegetation habitat. The second quadrant represented habitats with large direct fetch and an abundant concrete construct. Mytilus galloprovincialis and Littorina brevicula were plotted in this quadrant. Ruditapes philippinarum, Monodonta labio form confusa, or Nipponacmea radula belonged to the third quadrant, characterized by habitats with large direct fetch and sand habitat. Pattelloida pygmaea form heroldi, Nipponacmea nigrans, and Turbo cornatus coreensis were plotted in the fourth quadrant, characterized by abundant sand and boulder habitat.
In the nMDS dimensions, the Cerithidea genus, which inhabits muddy bottoms or reed fields in tidal wetlands, was plotted in the first quadrant. The first quadrant in this analysis represented a high tide and vegetation-rich environment. Additionally, Mytilus galloprovincialis, which inhabits wave-dominated rock reefs or concrete constructs via byssus attachment, was plotted into the second quadrant. The second quadrant in this analysis represented high direct fetch and a gravel-abundant environment. Therefore, the nMDS results appropriately reflect the habitat environment of these species.
3.3. River Estuary Classification Using Molluscan Fauna
From the results of the nMDS analysis, the molluscan fauna of the 19 rivers were divided into four groups (Figure 3). Group A comprised fauna from three rivers flowing into the Genkai Sea and the Kurosaki River flowing into the Sumou Sea. Group B comprised five rivers flowing into the Genkai Sea and Imari Bay. Group C comprised four rivers flowing into the Imari Sea and the Kusami River. Group D comprised six rivers flowing into the Imari Sea and the Buzen Sea. The Obaru River was not classified into any groups because no molluscan species were found.
Table 2 shows the IndVal index of molluscan fauna in the four groups.
In group A, the IndVal score of all species was lower than 20. The IndVal indices of Littorina brevicula and Mytilus galloprovincialis were relatively higher than those of the other species. The IndVal indices of Batillaria multiformis and Batillaria attramentaria exceeded 40 in group B. In group C, the IndVal indices of eight species exceeded 40, and especially, the value of Cerithidea ornata was the highest. The IndVal indices of Littoraria intermedia and Assiminea sp. exceeded 50 in group D.
3.4. River Estuary Classification Using Physical Indicators
Figure 4 shows the results of the cluster analysis using physical indicators. In this case, the 19 rivers were divided into three groups (groups α–γ). Group α was composed of the rivers flowing into the Ariake Sea and the Genkai Sea. All of the rivers flowing into the Imari Sea belonged to group β. Group γ was mainly composed of the rivers flowing into the Buzen Sea and the Ariake Sea.
The differences of six quantitative physical indicators (tidal range, wave exposure, direct fetch, topographic gradient, form ratio, specific discharge) among the three groups are indicated in a boxplot (Figure 5). Group γ exhibited the highest tidal range (388.4 cm), followed by group α (238.5 cm) and group β (220.0 cm). The results of the Kruskal-Wallis test indicated a significant difference in this regard among the three groups (p < 0.01). The Steel-Dwass multiple comparison test also showed a significant difference between group γ and the other two groups (p < 0.05). Direct fetch varied in the order of group α (583.0), group γ (110.7), and group β (11.7). The Kruskal-Wallis test showed a significant difference among the three groups (p < 0.001), as did the Steel-Dwass multiple comparison test. Wave exposure was the highest (9.4) in group α. The results of the Kruskal-Wallis test also indicated a significant difference among the three groups for this variable (p < 0.001), and the Steel-Dwass multiple comparison test indicated a significant difference between group α and the other two groups (p < 0.05). In contrast, topographic gradient, form ratio, and specific discharge were not confirmed to differ significantly among the three groups in one-way analysis of variance. The average values of eight qualitative physical indicators in each group are shown in Table 3