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).

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.

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).

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).

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.

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.

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.

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).