The Distribution Pattern of Marine Bivalve Death Assemblage From the Western Margin of Bay of Bengal and Its Oceanographic Determinants

Chattopadhyay, Devapriya; Sarkar, Deepjay; Bhattacherjee, Madhura

doi:10.3389/fmars.2021.675344

ORIGINAL RESEARCH article

Front. Mar. Sci., 04 June 2021

Sec. Marine Evolutionary Biology, Biogeography and Species Diversity

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.675344

This article is part of the Research Topic Benthic Biodiversity of the Indian Ocean View all 14 articles

The Distribution Pattern of Marine Bivalve Death Assemblage From the Western Margin of Bay of Bengal and Its Oceanographic Determinants

$\r\nDevapriya Chattopadhyay*$ Devapriya Chattopadhyay^1*

Deepjay Sarkar^2,3

Madhura Bhattacherjee¹

¹Department of Earth and Climate Science, Indian Institute of Science Education and Research (IISER), Pune, India
²Department of Earth Sciences, Indian Institute of Science Education and Research (IISER), Kolkata, India
³Department of Geology, Triveni Devi Bhalotia College, Raniganj, India

The global pattern of shallow marine biodiversity is constructed primarily using the data from extra-tropical sites. A severe knowledge gap in the shallow benthic diversity exists for the tropical Indian Ocean, especially along the coastline of peninsular India. Latitudinal biodiversity gradient (LBG)—a poleward decrease in diversity, even though accepted as a pervasive global pattern, often differs from regional trends. Although several oceanographic variables are known to influence regional patterns, their relative effect in shaping the shallow benthic community in tropical seas remains unclear. The east coast of India bordering the Bay of Bengal (BoB) presents a 2,500 km stretch (8–22°N) of tropical coastline with a spatial variation in oceanographic parameters including freshwater mixing, primary productivity, temperature, and shelf area. Here, we documented the marine bivalve distribution using spatially-temporally averaged beach samples and evaluated their relationship with the oceanographic variables. Our data reveal the existence of a highly diverse fauna, comparable to other tropical shallow marine sites. Overall species composition reflects a typical assemblage of the Indian Ocean, dominated by Veneridae but shows an uncharacteristically low proportion of Tellinidae and Lucinidae. The latitudinal variation in diversity shows a mid-latitude drop at around 14°N—a pattern inconsistent with the prediction of latitudinal biodiversity gradient (LBG). The functional groups are dominated by infauna (65%), unattached groups (69%), and suspension feeders (87%). There is only a slight difference in species composition between southern and the northern sites pointing to a predominantly continuous circulation and considerable mixing within the BoB. Productivity range, shelf area, and salinity emerge as best predictors of the species richness. All environmental variables together explain the species composition across the latitudinal bins satisfactorily. The species composition of the east coast shows no distinct nature in comparison to the Indo-Malayan biodiversity hotspot; the proximity to this hotspot and biological exchange with it may have contributed to the high diversity of the east coast fauna. Our study highlights the complex interplay between multiple oceanographic variables in determining the distribution and diversity of tropical shallow marine benthos at a regional scale generating biodiversity patterns that are at odds with global trends such as LBG.

Introduction

Recognizing the determinants of diversity and distribution of fauna at the global and regional scales is one of the important themes of ecological research, especially in the context of recent changes in the climate (Sankaran and McNaughton, 1999; Willig et al., 2003). The study of biodiversity is heavily dominated by the terrestrial ecosystem although the number of higher taxa is substantially greater in the marine realm than on land (Grassle et al., 1991). The lack of documentation of marine biodiversity is not homogeneously distributed and such spatial gaps often lead to substantial underestimation of the global biodiversity. Insufficiency in data is most pronounced in the tropical seas of Asia, especially along the coast of India even though this landmass harbors close to 6,100 km of tropical coastline (May, 1994). This area has been studied extensively for fishes and other macro-invertebrate groups of economic importance (Appukuttan, 1996; Kripa and Appukuttan, 2003) including molluscs. The primary goal of such studies was to identify and describe new species (Melvill, 1909; Hornell, 1922; Gravely, 1941; Ray, 1949; Satyamurti, 1956; Ganapati and Nagabhushanam, 1958; Kundu, 1965; Appukuttan, 1972; Subba Rao, 2017, ZSI, State fauna Series). The regional diversity pattern of molluscs is largely ignored except for studies focusing on specific species (Bharti and Shanker, 2021) or fauna of the west coast of India (Ingole et al., 2009, 2010). Present understanding of the factors controlling the molluscan distribution along the east coast of India is insufficient. The east coast of India bordering the Bay of Bengal (BoB) presents an example of a tropical coastline with a large latitudinal span (8–22°N) and a very high influx of terrestrial sediments through river runoffs influencing salinity profile along the coast (Shetye, 1993; Shetye et al., 1996). The coast also shows spatial variation in temperature, sediment input, and productivity (Bharathi et al., 2018). Molluscan assemblages have long been recognized to serve as a reliable proxy for diversity and ecosystem health in shallow marine systems (Kidwell, 2001). Hence the study of molluscan diversity is critical to appreciate the regional diversity pattern of shallow marine fauna and to recognize its drivers along the east coast of India.

Species richness of marine bivalves has been documented globally as part of a few large studies (Stehli et al., 1967; Crame, 2000a, b). Although measuring species richness is important for global-scale assessment, it cannot be taken as an absolute index, since it varies with habitat type, complexity and also does not reflect functional diversity. Furthermore, such indices do not reflect changes in species composition or even changes at higher taxonomic levels. Therefore, species richness, although useful for global studies, should be accompanied by more detailed information on species composition and ecology to capture the true essence of regional biodiversity. The only regional-scale studies to understand the drivers of molluscan distribution along the east coast of India are based on literature-based occurrence data (Sivadas and Ingole, 2016; Sarkar et al., 2019). Sarkar et al. (2019) attempted to also evaluate the functional diversity, community composition but recognized that the conclusions may have been affected by the sampling inconsistencies across the literature. A detailed sampling effort following a consistent protocol is necessary to evaluate the true pattern. It has been documented that “time-averaged” (temporally condensed record of multiple generations of fauna living in the same place) and “spatially-averaged” (spatially accumulated by random post-mortem transportation from different habitats at the same time) molluscan death assemblages could yield a near-perfect reflection of the true regional biodiversity despite the potentially different taphonomic history of the specimens (Kidwell and Bosence, 1991; Kidwell and Flessa, 1995; Kidwell, 2001; Warwick and Light, 2002; Mondal et al., 2021). A detailed within-habitat live sampling of marine molluscs along a ∼2500 km stretch of the east coast of India, while desired, proved to be logistically challenging. Sampling “spatially-temporally averaged” death assemblages as a surrogate for regional biodiversity is a technique suitable for the scale as vast as the present one.

In this study, we attempted to understand the regional pattern of bivalve diversity of the east coast of India by employing a detailed sampling of time-averaged death assemblages to minimize methodological inconsistencies. Using the specimens collected from beach samples from 17 localities spanning over 14 latitudinal bins along the east coast, we addressed the following questions:

1. What is the nature of spatial variation in diversity, ecology, and species composition?

2. How do the oceanographic variables influence shaping the patterns?

Materials and Methods

Study Area

The BoB is a tropical semi-enclosed basin, bordered by India and Sri Lanka on the west, Bangladesh on the north, and Myanmar and Thailand on the east. A stretch of ∼2500 km coast of the Indian subcontinent marks the western margin of BoB (Figure 1A). Western BoB experiences seasonally reversing monsoons, severe cyclonic storms (SCS), a large amount of rainfall, and river run-off (Mohanty et al., 2008). The northwestern coastal BoB receives discharge mainly from glacial rivers (Ganges, Brahmaputra). Peninsular rivers (such as Godavari, Krishna) discharge in the southwestern coastal BoB. The high amount of annual freshwater discharge (∼102 m³) from the Ganges-Brahmaputra rivers system at the northwestern BoB especially during the monsoon (June-October), results in the development of an N-S gradient in salinity and temperature along the east coast of India; the northern Bay is characterized by low salinity and cooler temperature compared to the southern regions (Shetye, 1993).

FIGURE 1

Figure 1. Map of India (A) showing the sampling locations (B). The map is based on Survey of India Outline map (1996).

Field Sampling and Identification

We selected 17 sampling localities that span from Kanyakumari (8.08°N) at the south to Sundarbans (21.94°N) at the north (Table 1) spread across 14 latitudinal bins (1° each, Figure 1B). Each bin is represented by at least one sampling locality. The minimum distance between two sampling localities is 5 km. For collecting bivalve specimens from each locality, we took a traverse of ∼1 km along the seashore; all visible molluscan shells were collected. Specimens that were damaged beyond recognition were ignored. We did not collect sediment samples during our collection. The sampling exercise was conducted during a period of five and half years, from July 2010 to December 2015.

TABLE 1

Table 1. Details of the sampling locations used for the present study.

All bivalve specimens were separated from the collected molluscan sample, washed, sun-dried. The individual number was calculated by counting the articulated valves and adjusting for the disarticulated valves. Each latitudinal bin was represented by a minimum of 200 individuals. Taxonomic nomenclature was primarily based on the published work by Subba Rao (2017) and the World Register of Marine Species (WoRMS Editorial Board, 2020). Because of the unavailability of authoritative monographs on Indian bivalves and controversy on species nomenclature, we used an internally consistent species nomenclature and used it as the operational taxonomic unit. Additional information on identification and ecology was derived from Apte (1998), Zuschin and Oliver (2003), Oliver (1992), ZSI, State fauna Series, and NMiTA¹. The species occurrence of Indo-Malayan Archipelago (south-east Asian biodiversity hotspot) was collected from published literature (Dijkstra, 1991; Dolorosa et al., 2015; Lutaenko, 2016) and compared with the data along the east coast of India from the present study and literature-based occurrence compiled in Sarkar et al. (2019).

Raw species richness is defined as the direct count of observed species. Raw and rarefied species richness was calculated using the species abundance in each latitudinal bin. Additional indices of diversity and evenness were evaluated by Shannon–Wiener index (H), Pielou’s evenness index (J), respectively, to account for sample size. The species were classified into three groups based on substrate relationship (epifaunal, infaunal, others to include borer and nestler on or within hard substrates), three finely subdivided substrate guilds (byssally attached, cemented, and unattached), and four trophic guilds (suspension feeders, deposit feeders, and chemosymbionts).

Oceanographic Variables

We have included oceanographic variables (productivity, sea surface temperature, and salinity) that are known to influence molluscan diversity (Roy et al., 1998; Valentine and Jablonski, 2015). We have collected pre-existing data on the annual mean and range of net primary productivity (NPP), sea surface temperature (SST), salinity from the Ocean Productivity database² for each latitudinal bin. Species diversity of shallow marine fauna is also known to vary with the available area of habitat (Smith and Benson, 2013) and we used shelf area as a proxy for available habitat. The shelf area for a latitudinal bin is calculated as a product of coastal length and the average distance from the coast to the shelf-slope break for specific locations. The coastal length and shelf width are obtained from GEBCO Compilation Group (2020). Cyclone frequency is known to affect the molluscan death assemblages especially in tropical siliciclastic settings (Bhattacherjee et al., 2021); we included it in our analyses to rule out the possibility of transportation-related influence on the observed species richness and to establish the reliability of our data. Cyclone frequency data were collected from the global-tropical-extratropical cyclone climatic atlas from the United States Navy National Climate Data Center cyclone records. The details of the processing are discussed at Bhattacherjee et al. (2021).

Statistical Analysis

We used the Spearman rank-order coefficient to measure the correlation of diversity with latitude. A significant negative correlation between diversity and latitude is consistent with predictions from LBG. To evaluate the relationship between various oceanographic variables and diversity, functional groups, we used the Spearman rank-order correlation test. We also estimated the effect of these variables on diversity and proportion of dominant functional group with multiple generalized linear models (GLMs) to analyze all parameters simultaneously and evaluating their partial contributions to the total variation in diversity (Quinn and Keough, 2002). We employed the Akaike information criterion (AICc) to rank alternative GLM models based on the trade-off between the model fit and its complexity (Burnham and Anderson, 2002; Grueber et al., 2011). We generated a set of models with all possible combinations of environmental predictors and scored them according to their relative support as measured by AICc. To account for the uncertainty in model selection, we performed model averaging, applying a cut-off criterion of ΔAICc ≤ 2 to choose the best set of models.

To recognize the compositional distinctness along the coast and between regional occurrences within the BoB, we used ordination analysis. The compositional similarities among latitudinal bins were calculated from an abundance matrix of species in bins using the Bray–Curtis similarity index. For biogeographic comparison between India and Indo-Malayan Archipelago, we used an occurrence matrix of species using the Sørensen similarity index. The similarity matrices were clustered by the unweighted pair group method using arithmetic averages (UPGMA) and visualized as a dendrogram. Two-dimensional ordination assembles were created with non-metric multidimensional scaling (NMDS) using the Bray–Curtis similarity indices for comparison among latitudinal bins and Sørensen similarity indices between sites from India and Indo-Malayan Archipelago. To assess the relative importance of environmental variables on the distribution of species along the coast, we used Redundancy Analysis (RDA) on the Hellinger distance-transformed species data (Legendre and Legendre, 1998).

All univariate and multivariate analyses were performed in R (R Core Development Team, 2012). The ecological analyses were done using the “Vegan” package in the R platform.

Results

We have collected a total of 11,458 individual bivalves from 14 latitudinal bins (Table 1) that represent 287 species, 129 genera, and 35 families. The raw species richness is lower than the observed richness in Indo-Malayan Archipelago (518). There is no significant correlation between location-specific abundance and species richness (Spearman’s rho = −0.06, p = 0.82) (Figure 2). The average raw species richness is 44 and rarefied richness is 36 (Table 2). The most common families are Veneridae, Cardiidae, Arciidae, Mactridae, and Donacidae (Table 2). The relative proportion of infauna (substrate relationship), unattached (attachment type), and suspension feeders (feeding style) seems to dominate various functional groups (Figure 3) (Table 2 and Supplementary Table 1).

FIGURE 2

Figure 2. (A) Relationship between species richness and abundance. (B) Rarefied species richness at different latitude bins.

TABLE 2

Table 2. Summary of the diversity of marine bivalve families and functional groups along with their correlation with latitude.

FIGURE 3

Figure 3. Relative proportion of various ecological groups classified based on (A) substrate relationship, (B) attachment type, and (C) feeding type. The boxes are defined by 25th and 75th quantiles. The thick line represents median value within the box.

There is no significant correlation between species richness and latitude (Spearman’s rho = 0.158, p = 0.588) (Figure 4 and Table 2). However, the diversity appears to drop around 14°N. A similar trend is observed in rarefied richness, evenness, and Simpson diversity indices (Figures 4A–D). The species richness of the five most abundant families reveals no significant correlation with latitude, except for Mactridae (Figures 4E–I and Table 2). The latitudinal distribution of life modes also does not show a significant correlation (Figures 4J–R and Table 2).

FIGURE 4

Figure 4. Latitudinal variation in various diversity indices (A–D), common families [(E) Veneridae, (F) Arcidae, (G) Cardiidae, (H) Mactridae, and (I) Donacidae] and functional groups based on substrate relationship [(J) infauna, (K) epifauna, and (L) others], attachment [(M) byssally attached, (N) cemented, and (O) unattached] and feeding [(P) deposit, (Q) suspension, and (R) chemosymbiotic].

In the NMDS plot (stress value = 0.155), the sites do not form latitudinally separated clusters, although sites north of 15°N show a slight separation along the NMDS axis 2 (Figure 5A). Some species that are abundant in the north of 15°N such as Wallucina assimilis, Trapezium oblongum are absent in the south. Similarly, a few dominant species of the south such as Vepricardium coronatum and Timoclea costellifera are absent in the northern region. A dendrogram constructed by UPGMA, however, does not show any distinct cluster of latitude-specific species composition (Figure 5B).

FIGURE 5

Figure 5. (A) Ordination plot using non-metric multidimensional scaling, based on Bray–Curtis distance calculated from a matrix of marine bivalve species at different latitudinal bins. (B) Dendrogram constructed by UPGMA method using arithmetic averages.

Productivity (range) shows a significant correlation with species richness (Figure 6A and Table 3). None of the environmental variables significantly correlated with the proportion of infauna, unattached, or suspension feeders (Figures 6B–D and Table 3). Among the environmental variables, productivity (mean and range), temperature (mean and range), shelf area, and cyclone frequency do not show any relationship with any other environmental variables (Supplementary Table 2). After excluding the variables based on autocorrelation, the GLM models show that the combination of productivity (range), salinity (mean), and shelf area were the most important predictors of species richness (Table 4 and Supplementary Table 3). Shelf area and temperature (mean) emerge as important predictors for dictating the proportion of infauna. Shelf area, productivity (range), and mean value of salinity and temperature emerge as strong predictors for the proportion of unattached bivalves (Table 4). Salinity (range, mean), temperature (range), and productivity (mean) predict the proportion of suspension feeders (Table 4). Except for species richness, however, the variables have limited explanatory significance as indicated by the low or negative value of adjusted R².

FIGURE 6

Figure 6. The relationship between oceanographic parameters and species richness (A), proportion of infauna (B), unattached (C), and suspension feeder (D).

TABLE 3

Table 3. Correlation between species richness, proportion of ecological groups with environmental variables for latitudinal bins using Spearmann rank order correlation test.

TABLE 4

Table 4. Highest ranked logistic models (ΔAICc ≤ 2) investigating environmental correlates of the species richness and ecological categories of the marine bivalves.

About 79% of the variation in species composition is explained by all environmental variables considered together using RDA on the Hellinger distance-transformed species data (Figure 7 and Table 5). The constrained variance of the global model with all environmental variables 0.66 and the unconstrained variance 0.34. Forward selection to choose a model with fewer variables, however, increases the unconstrained variance to 0.79 pointing to the limited explanatory power of the models with fewer environmental variables. The RDA plot also shows no significant influence of any single environmental variable in explaining the species composition across latitudinal bins.

FIGURE 7

Figure 7. Redundancy analysis (RDA) biplot of bivalve species composition and environmental parameters.

TABLE 5

Table 5. Biplot scores for constraining variables derived from the redundancy analysis (RDA) to determine the relative contribution of each environmental variable for the compositional variation of bivalves.

In an occurrence-based NMDS plot (stress value = 0.131) comparing the present species composition with literature-based occurrence data of the east coast (Sarkar et al., 2019) and Indo-Malayan Archipelago, the sites do not form geographically separated clusters indicating a high degree of compositional similarity between the two regions (Figure 8).

FIGURE 8

Figure 8. Ordination plot using non-metric multidimensional scaling, based on Bray–Curtis distance calculated from a matrix of marine bivalve species occurrence from Indo-Malay Archipelago, literature-based data of east coast from Sarkar et al. (2019) and the present study.

Discussion

Studies on bivalves from the Indian coast primarily focused on the detailed morphological description and identification at a local scale (Subba Rao, 2017, ZSI, State fauna Series among many others) barring a few (Ansari et al., 1977; Khan et al., 2010; Manokaran et al., 2015). Because local communities are an integral part of larger biogeographic regions, the diversity of local scale (spatial scale of meters to hundreds of meters) must be affected by regional-scale processes (spatial scale of 200 to thousands of kilometers) (Witman et al., 2004). The regional-scale studies on Indian bivalve diversity are primarily conducted on the west coast (Ingole et al., 2009, 2010). Regional analysis of the east coast fauna is based on data from literature compilation (Sivadas and Ingole, 2016; Sarkar et al., 2019) and hence suffers from potential sampling inconsistency. In contrast, the present study aims to understand the regional pattern of bivalve diversity instead of discrete local variations by employing a detailed sampling of time-averaged death assemblages to minimize methodological inconsistencies.

Reliability of Beach Sampling

Regional biodiversity assessments are strongly influenced by the sampling protocols. A regional-scale species richness is difficult to measure because the richness increases with the effort of sampling (Warwick and Light, 2002). Grab samples have often been considered to portray a true representation of the biological community (McKinney and Hageman, 2006). However, there are some limitations. A fundamental issue is that the megabenthos—benthic organisms that are large enough to be visible on seabed photographs, cannot be studied reliably from grab samples alone owing to their patchy distribution (Gage and Tyler, 1991). Comparison between grab samples and underwater studies revealed this difference for many regions including the Adriatic Sea (Zuschin et al., 1999; Stachowitsch et al., 2007; Riedel et al., 2008; Zuschin and Stachowitsch, 2009), and Antarctica (Dayton and Oliver, 1977) highlighting the role of sampling technique on the inferred diversity and ecological composition of a region.

Kidwell (2013) claimed that death assemblages are more diverse than the actual living fauna of an area due to time averaging. Owing to temporal accumulation of shifting patches and post-mortem transportation from proximal habitats, time-averaged assemblages tend to capture a more complete picture of biodiversity compared to ecological snapshots. However, post-mortem transport and reworking of older shells are most common in marine coastal areas and could potentially bias the time-averaged data. Processes such as abrasion, fragmentation, bioerosion, and decay could affect the death assemblages on recent marine fauna (Kidwell and Bosence, 1991; Warwick and Light, 2002; Nebelsick et al., 2011). Few marine shells, especially those that are extremely thin, dominated by a high-organic microstructure, with an exceptionally short lifespan, might be less represented in death assemblages due to differential post-mortem durability (Albano and Sabelli, 2011; Kidwell, 2013). High energy events such as tropical storms can also produce habitat mixing of shallow seafloors (Miller et al., 1992; Bhattacherjee et al., 2021). The small original sample size and severe analytical truncation of a data set can magnify the effects of such bias (Kidwell, 2013). Despite all these potential biases, time-averaged death assemblages from beach samples have been documented to substitute a conventional sampling effort of regional diversity (Kidwell, 2013) and may provide a stronger signal for species abundance compared to the original living community (Kidwell, 2001). Consequently, beach sampling is considered a valid sampling technique for diversity studies for Recent marine molluscs (Warwick and Light, 2002; Mondal et al., 2021). It is widely utilized for its logistical ease compared to other techniques for acquiring live community data (Warwick and Light, 2002; Zuschin and Ebner, 2015).

Beach sampling may present a few limitations in reconstructing the community structure of the east coast of India. In a siliciclastic setting frequented by storms in the northern east coast of India, Bhattacherjee et al. (2021) demonstrated a low life-dead fidelity of molluscan species in a local scale. Using the published literature data on bivalve occurrence, they also found that the frequency of cyclones may induce significant mixing at a regional scale. The abundance-based species diversity/composition of the present study, however, does not reveal significant dependence on cyclone frequency pointing to the preservation of the true biodiversity signature of the beach samples (Tables 3, 4). Preservation of original ecological structure in beach samples is another point of worry. All infaunal bivalves have aragonitic shells and major epifaunal shells are calcitic. Due to the higher preservability of calcite, epifaunal shells may dominate a dead assemblage at the expense of the infaunal ones; this bias is prevalent among semi-lithified specimens (Cherns et al., 2011). Our specimens are less likely to suffer from such bias, primarily because of their pristine condition. More importantly, our results show a dominance of infauna—a pattern opposite of what would have been created by this specific preservation bias. This confirms the recent findings of Gupta et al. (2020) that the beach assemblages from the eastern coast of India preserve pristine shells pointing to an assemblage enduring the least taphonomic disturbance.

Nature of Variation in Species Richness and Composition

Regional-scale studies of tropical shallow marine bivalve diversity are limited. The east coast of India is known to harbor a rich assemblage of molluscan fauna (Subba Rao, 2017). Our specimens, collected over 17 sites from the east coast, also show a high species richness of bivalve fauna (∼287), although slightly lower than the occurrence-based reports of 371 species (Sarkar et al., 2019) and 518 species reported from the biodiversity hotspot of Indo-Malayan Archipelago. The sampling protocol of the present study did not include the sediments and small bivalves therein. This exclusion may explain the lower observed diversity in comparison to the occurrence-based reports of BoB. The species richness of this region is comparable to other tropical hotspots such as the Red Sea (Zuschin and Graham Oliver, 2005).

Latitudinal biodiversity gradient is one of the most commonly accepted global patterns of regional faunal distribution including marine bivalve distribution (Roy et al., 1998; Jablonski et al., 2006; Roy and Goldberg, 2007). Global pattern documented for various temperate and polar marine regions often shows varying character; its true nature in species-rich intra-tropical region is still sparsely studied. The present study shows an interesting pattern of latitudinal variation in the tropics that does not follow the classical LBG. Instead of a monotonic increase leading to a peak at the southern point, it shows a low species diversity near 14°N followed by a symmetric increase (Figure 5). The drop in evenness near 14°N along with a change in diversity points to the contribution of multiple oceanographic parameters varying regionally as opposed to the monotonic variation in single variables such as sea surface temperature. In contrast to the claim of SST being the primary predictor of the global distribution of marine biodiversity, our study found a negligible role of temperature in dictating the variation in species richness at a regional scale (Table 4).

The distribution of species is influenced by a combination of species-specific traits (physiological tolerance, nature of reproduction, and dispersal) and the nature of the available habitat (oceanographic variables and geographic barriers) (Hesse et al., 1951). The overall composition of bivalve species of our study shows a characteristic assemblage of the Indian ocean, dominated by Veneridae. Species of Timoclea, Sunetta, Meretrix, Anadara, Mactra, and Donax constitute more than 70% of the overall species composition. The assemblage, however, differs from a typical Indian ocean assemblage by its uncharacteristically low proportion of Tellinidae and Lucinidae. It is important to note that these three families are the most important components of tropical infaunal bivalve communities in both species richness and individual abundance (Jackson, 1974). Species association along the coast neither shows a continuous change in diversity as predicted by LBG nor any distinct cluster reflecting the existence of any oceanographic barriers. Occurrence data on species composition points to the existence of such a strong barrier around 15°N along the west coast separating the species association of the northwestern sites from the rest of the Indian sites (Sarkar et al., 2019). This barrier on the west coast is primarily developed due to differences in salinity and productivity between the northern and southern Arabian Sea (Madhupratap et al., 2001; Sarkar et al., 2019). Based on the same occurrence data, Bhattacherjee et al. (2021) demonstrated a higher level of similarity in species associated with the northern part of the east coast compared to the southern part due to the high intensity of cyclone-induced mixing of the bottom fauna. Considering the methodological inconsistencies in sampling that may have influenced the literature-based occurrence data, we should take the results with caution. The abundance-based data of the east coast of the present study, in contrast, shows a more continuous variation in species composition along the east coast (Figure 6) pointing to a relatively well-mixed circulation of the BoB without any major oceanographic barrier in contrast to the western coast of India (Sarkar et al., 2019). The lack of sharp change in oceanographic variables along the coast (Supplementary Figure 1 and Supplementary Table 4) also reaffirms our claim of predominantly continuous circulation.

The functional composition of the bivalve species is primarily dominated by infauna and bivalves without attachments. Among infaunal genera without any attachment, Cardium, Sunetta, and Mactra are the most common ones. The species of Ostrea dominates the epifauna with attachment. The species of Sunetta, Donax, and Anadara are the most common suspension feeders. Only a very few species of Wallucina represent chemosymbionts in this assemblage. Sarkar et al. (2019) did a preliminary analysis on the diversity and distribution of Indian marine bivalves based on occurrence data. A comparison of the major results between these two studies regarding the east coast, suggests a high degree of congruence despite possible methodological inconsistency and limitation in sampling effort in the first approach. Both the datasets revealed a higher richness at all taxonomic level in the east coast where the two most common families are Veneridae, Arcidae; both the studies showed a dominance of specific functional groups such as infaunal and unattached. The lack of a clear LBG and a non-monotonic increase in species richness along the east coast is demonstrated by both datasets.

Role of Physical Environment Influencing Diversity and Distribution

The eastern coast of India harboring the BoB displays variation in several oceanographic variables (Supplementary Figure 1). The northern region is characterized by high siliciclastic input from large rivers and variable salinity (Subramanian, 1993; Ganesh and Raman, 2007) while the southern region is characterized by relatively low siliciclastic input and presence of reefs (UNEP-WCMC, WorldFish Centre, WRI, TNC, 2010). The observed variation in bivalve distribution along this coast, however, is not controlled by any single physical parameter and points to the emergent influence of multiple oceanographic variables creating a complex biological response.

Productivity is widely believed to be a significant controller of marine benthic diversity which may increase between regions of low to moderate productivity, and then abruptly decline toward regions of higher productivity (Levin et al., 2001; Rex et al., 2005). A strong positive correlation is found between eutrophication and bivalve diversity throughout the Indo-West Pacific tropical province (Vermeij, 1990; Taylor, 1997). At low levels of productivity, food limitation is thought to limit the number of species that can survive. The reason behind the decline in species diversity in regions of highly fluctuating productivity is attributed to demographic stochasticity that may bring about diversity declines (Levin et al., 2001). The highly diverse east coast is characterized by low but generally stable productivity. Despite a high riverine flux bringing nutrients, the productivity remains lower (Nair et al., 1973; Madhupratap et al., 2003; Kumar et al., 2007) because of the loss of nutrients in the deep ocean (Qasim, 1977; Sen Gupta et al., 1977; Radhakrishna et al., 1978). The chlorophyll α concentrations (10–20 mg m^–2) and primary productivity values (40–502 mg C m^–2 d^–1), although lower than the Arabian Sea values for the same season, show significantly lower seasonal fluctuation in the east coast (Madhupratap et al., 2003). Our study identifies the important role of productivity variation in guiding species richness and proportion of dominant functional groups (Table 4). Although mean productivity is not correlated with richness, an increase in productivity range appears to cause a decrease in species richness. Such a decrease is probably caused due to the fluctuating oxygen concentration produced from the degradation of organic flux and related stress produced by the fluctuating productivity (Berger et al., 1989). This phenomenon is not unusual for benthic communities such as bivalves. Global data shows that a highly diverse benthic community is more likely to be supported by a stable primary production than a fluctuating production summing up to a higher value of annual productivity (Valentine and Jablonski, 2015). Moreover, productivity is also found to be correlated with change in dominance of functional groups for bivalves; oligotrophy supports a community of exposed sedentary suspension feeders in contrast to active benthic communities of highly productive regions (McKinney and Hageman, 2006). The dominance of mobile infauna (Figure 4) in our study also corroborates the same finding.

Temperature is known to play an important role in controlling the diversity, composition, and geographic distribution of marine species at a global scale (Barry et al., 1995; Walther, 2002; Parmesan and Yohe, 2003; Hiscock et al., 2004; Tittensor et al., 2010). Bivalves are exposed to a wide range of temperatures depending on their substrate relationship. Species that live at shallow depths (1 m or less) experience a greater range of temperature in a single locality than are most deep-water species over their entire geographic range (Jackson, 1972, 1973, 1974). Consequently, the shallow marine bivalves show a relatively high-temperature tolerance. Such inherent physiological tolerance of shallow marine bivalves along with a relatively narrow variation in temperature along the east coast (2.5°C) (Supplementary Figure 1 and Supplementary Table 4) explains why temperature does not emerge as a strong predictor of diversity in our study (Figure 7 and Table 4). This also supports the existing concerns about generalizing the effect of temperature on marine bivalves (Nawrot et al., 2017; Chattopadhyay and Chattopadhyay, 2020).

Another important factor influencing marine biodiversity is salinity. Salinity is found to play a role in guiding the species diversity of estuarine molluscs (Montagna et al., 2008). The relationship between salinity and diversity is more complex for marine ones (Sanders, 1968; Roy et al., 1998). An inverse relationship between bivalve species richness and salinity range is observed for global distribution (Valentine and Jablonski, 2015). The lack of correlation between salinity and overall species composition in our study may be influenced by the fact that marine organisms at high temperatures tend to tolerate the variation in salinity much more easily than in low temperatures (Panikkar, 1940; Sanders, 1968). Our results also show that salinity (mean) is the common factor influencing the richness and proportion of dominant functional groups (Table 4).

The long-held axiom in ecology states that larger areas hold more species than smaller areas (Rosenzweig, 1995) and predicts a positive influence of available area on diversity. Available habitat area plays an important role in shaping the community structure of shallow marine benthos through time and the continental shelf is recognized to be one of the most speciose habitats (Piacenza et al., 2015). Reconstructing the true shelf area is far from challenging and hence, several proxies are used (Holland, 2012). A few studies have used coastal length as a proxy (Cain, 1938; Connor and McCoy, 1979; Tittensor et al., 2010) and others estimated shelf area as a product of coastal length and width (Sanciangco et al., 2013). All of these global studies confirmed a positive correlation between habitat area and biological diversity of marine groups including bivalves (Valdovinos et al., 2003). Literature-based occurrence data shows a correlation between coastal length and marine molluscan diversity along the Indian coast (Sivadas and Ingole, 2016; Sarkar et al., 2019). Our results also show an influence of shelf area on species diversity (Table 4) supporting the validity of the habitat-area hypothesis in a small regional scale.

Among regional environmental parameters, no single environmental variable explains the major canonical axes of species variation satisfactorily (Figure 7 and Table 5) indicating the influence of several environmental variables acting together in shaping the regional nature of species distribution. Occurrence-based data demonstrate a partitioning between northern and southern sites where the southern eco-region is characteristically dominated by borers like Martesia striata, and the northern eco-region is dominated by various species of Donax (Sarkar et al., 2019). Our abundance data, however, shows only a slight difference in species composition between southern and northeastern sites. The lack of strong compositional dissimilarity along the east coast probably owes to a highly mixed circulation pattern of the BoB.

Apart from the selected oceanographic variables considered in the study, the nature and amount of sediment influx may play an important role in guiding the distribution, especially influencing the distribution of functional groups in this region. The eastern coast of peninsular India is characterized by extremely high sediment influx (>1,350 million tons of suspended sediments/year) brought by the Ganges-Brahmaputra River systems together with other rivers to the BoB, especially toward the north (Milliman and Meade, 1983; Subramanian et al., 1985). This high riverine input makes it easier for the infaunal bivalves to thrive and flourish. The southern part of India (8–15°N) receives less than a fourth of the sediment supply brought by smaller rivers (Vamsadhara, Hyadri, Godavari, and Krishna) compared to the larger rivers of the north (Bharathi et al., 2018). The presence of coralline hard substrate in the south (UNEP-WCMC, WorldFish Centre, WRI, TNC, 2010) along with relatively low riverine input make it more habitable for epifaunal forms—a slight increase in epifaunal proportion in the southern latitudinal bins (Figure 3) reflects this change in habitat. In contrast to the reported dominance of epifauna in the west coast (Sarkar et al., 2019), the east coast demonstrates the dominance of the infaunal sediment burrowers (Figure 4). A similar pattern is observed for the dominance of the unattached form in the east coast compared to that of the attached forms in the west coast. The dominance of infaunal, unattached species in the east coast is most likely a reflection of higher siliciclastic sediment input. Even though we did not measure the sediment load and the nature of substrate explicitly in this study, the observed pattern of coastal difference in dominant functional groups underscores the role of substrate in structuring shallow benthic communities.

Another important regional parameter, imparting a likely control over the species richness and composition, might be the location of the nearest biodiversity hotspot. Crame (2000a, b) has documented the presence of a bivalve biodiversity hotspot near the Indo-Malayan Archipelago and claimed that species are radiating from there. He put this as a mechanism to explain the clines of species richness decreasing radially from this area in a north-south latitudinal pattern and east-west longitudinal pattern. The higher species richness of the east coast may also be influenced by this proximity. It has long been recognized that the BoB faunal assemblage largely represents the Indo-Malayan archipelago and is considered a contiguous part of it. However, it has not been documented for bivalve fauna. Both occurrence data (Sarkar et al., 2019) and the present study reaffirm the lack of distinctness of bivalve fauna of the east coast from the Indo-Malayan species association (Hocutt, 1987). The present study suggests that the high-diversity regions in the east coast may have developed due to a combination of several oceanographic variables and influenced by larval transport with biodiversity hotspots like the Indo-Malayan Archipelago.

No single environmental predictor emerged to explain the diversity in the east coast. All the environmental variables together explain the species composition satisfactorily (Figure 7). This underscores the complex interplay between multiple oceanographic variables in determining the distribution and diversity of tropical shallow marine benthos in a regional scale—a pattern that is noted for benthic biodiversity in large marine ecosystems such as the western margin of the United States (Piacenza et al., 2015).

Conclusion

This study using “time-averaged” and “spatially-averaged” bivalve death assemblages from beach sampling is one of the first attempts to document regional species diversity of marine bivalves along the east coast of India and to investigate its oceanographic correlates. By following a rigorous and consistent sampling protocol, this study reveals the regional pattern of bivalve diversity and distribution in this largely unexplored region. The overall diversity is comparable to other tropical shallow marine sites. Bivalve assemblages along the east coast are characterized by a higher proportion of infauna and groups without any attachments than in the west coast. This difference points to the inherent difference in the physiography of the coasts; while the east coast receives a higher amount of riverine sediments that is conducive for infaunal bivalves, a low influx of sediments makes the west coast favorable for epifaunal groups. When compared to classical trends like LBG, the diversity change with latitude does not show a consistent pattern pointing toward a more regional nature of diversity variation. There is no strong distinction in species composition along the east coast indicating a continuous circulation and considerable mixing of water within the BoB. Productivity range, shelf area, and salinity appear to influence the species richness. Along with these oceanographic variables, salinity contributes significantly to determine the proportion of the dominant functional group (unattached, suspension feeders). All the environmental variables together explained the species composition across latitudinal bin pointing toward a complex interaction between multiple environmental variables and their contribution to species composition. Lack of distinctness in bivalve species composition in comparison to those of the Indo-Malayan archipelago suggests the influence of this biodiversity hotspot in the development and maintenance of the species-rich fauna of the east coast of India. Apart from demonstrating that regional patterns in bivalve distribution are influenced by a complex interplay between several environmental factors, our study also establishes the temporally-spatially averaged molluscan death assemblages as a reliable proxy for regional biodiversity studies.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

DC designed the research. DS collected the data. DS and MB analyzed the data. DC wrote the manuscript with input from DS and MB. All authors contributed to the article and approved the submitted version.

Funding

The research was funded by IISER Kolkata research grant (ARF-2012, ARF-2013), IISER Pune Startup Grant, SERB core research grant (CRG/2018/002604). DS was supported by IISER Kolkata Ph.D. Fellowship. MB was supported by IISER Kolkata IPh.D. fellowship and IISER Pune Ph.D. fellowship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Chattopadhyay et al. Bivalve Distribution Along Indian Coast S. R. Prasanjit and Saswata Roy for their help in collecting specimens. We would also like to thank the help of Sudipta Sarkar in generating the data on shelf area and preparing the map. Comments from two reviewers significantly improved the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.675344/full#supplementary-material

Supplementary Figure 1 | Latitudinal variation in the oceanographic parameters along the east coast of India.

Supplementary Table 1 | Summary of Wilcoxon test comparing the proportion of various ecological types. The significant results are in bold.

Supplementary Table 2 | Significance of Spearmann correlation test between environmental variables. The significant results are in bold.

Supplementary Table 3 | Results of multiple GLM analyses to assess contribution of environmental variables in determining species richness, substrate relationship (infaunal proportion), attachment type (unattached proportion), feeding (suspension proportion). The significant results are in bold.

Supplementary Table 4 | Summary of the environmental variables used in the present study.

Footnotes

References

Albano, P. G., and Sabelli, B. (2011). Comparison between death and living molluscs assemblages in a Mediterranean infralittoral off-shore reef. Palaeogeogr. Palaeoclimatol. Palaeoecol. 310, 206–215. doi: 10.1016/j.palaeo.2011.07.012

CrossRef Full Text | Google Scholar

Ansari, Z. A., Parulekar, A. H., Harkantra, S. N., and Nair, A. (1977). Shallow water macrobenthos along the central west coast of India. Mahasagar 10, 123–127.