Genetic Diversity, Population Structure and Demographic History of Orange Mud Crab Scylla Olivacea from the Bay of Bengal, Bangladesh and Adjacent Seas in the Northern Indian Ocean Based on Mitochondrial COI Gene Sequences

In this study, we analyzed mitochondrial COI gene sequences to reveal genetic diversity, population structure and demographic history of two Bangladeshi (BD) populations (SB and CK) of the orange mud crab Scylla olivacea of the Northern Bay of Bengal (BoB), and compared these two with other four populations in the Northern Indian Ocean region (Arabian Sea, Andaman Sea and Malacca strait) and South China Sea. For all of the populations, nucleotide diversities were low (0.005–0.01) while the haplotype diversities were as high as 0.70–0.96, indicating that the S. olivacea has undergone a recent population expansion after experiencing bottleneck. The pairwise population statistics (FST) revealed that no genetic variation was made between SB and CK populations of BD in BoB. However, these two BoB populations showed separate genetic structure with each of the Andaman Sea (Myanmar coast, MM) and Malacca strait (West coast of Malaysia, MS) populations. On the other hand, two BoB populations did not form separate genetic structure from the population of Arabian Sea (AS). Larval dispersal-based migration by the East and West India coastal currents probably caused this genetic homogeneity between BoB and AS populations.The MM population had separate genetic structure from all of the populations studied in the present study. The Hierarchical analysis of molecular variance (AMOVA) revealed four large population groups of S. olivacea within its distribution range in the Indo-west Pacific region namely, Arabian Sea, Bay of Bengal, Andaman Sea and South China Sea groups. Some geographical barriers (e.g. Indian peninsula, Andaman and Nicobar Islands) along with seasonally formed marine gyres in the Andaman Sea are responsible for separate genetic structure among different populations and also for establishing four population groups. Star-shaped patterns of haplotype network and neutrality test corroboratethe recent population expansion of all populations except MM and CK. Mismatch distribution analysis reveals that the demographic expansion of the species started during the late Pleistocene period approximately 125,000 to 365,000 years ago. These results will help to establish the conservation and management strategy for orange mud crab in the Northern Indian Ocean region including the Bay of Bengal.


Introduction
Mud crabs of the genus Scylla de Haan, 1833 (Crustacea: Decapoda: Brachyura: Portunidae) inhabit mangrove areas and estuaries having wide distribution range throughout the Indo-Pacific region, including the Bay of Bengal. Mud crab species are important crustacean resource in commercial coastal fisheries and aquaculture industry of several Asian countries, including Bangladesh (Overton et al.,1997). Mud crab in Bangladesh breeds throughout the year but March to April was identified as the peak breeding season and the second peak was observed in August to September (Ali et al., 2020). For Andaman Sea, the breeding period was reported in June to November (Koolkalya et al., 2006). Reproduction of mud crab takes place in the open sea and the presence of a planktonic larval stage lasts 16 to 75 days depending on environmental conditions (Baylon, 2010). The larvae have high dispersal potential during the planktonic life stages through passive transport by regional and global currents. This enables mud crab to connect with other geographic areas and expand to new habitats, resulting in panmixia of the species (Hill, 1994;Gopurenko et al., 1999Gopurenko et al., & 2003. Identification of different species of Scylla is often difficult due to their high morphological plasticity, overlapping morphological, and morphometric traits. Recently partial fragments of nuclear and mitochondrial DNA of Scylla species collected from Bangladesh coast including Sundarbans confirmed that the mud crab species commonly found in Bangladesh is orange mud crab Scyla olivacea (Sarower et al., 2016;Habib et al. 2016;Akter, 2016).
Mud crab is exported alive from Bangladesh to various countries of the world such as China, Taiwan, Thailand, Singapore etc. These live crabs are harvested from the mangrove area mainly from Sundarbans and also from Chakaria located in the South-west and South-east coast of the country, respectively. Some mud crabs are also cultured in the coastal brackish water ponds, locally called 'gher'. The present culture system mostly involves fattening of small underweight crab to produce market size hard shell crabs with supplementary feeding. Another culture system includes the production of soft-shell crabs harvested just after molting when the exoskeleton is removed which are exported as frozen food. The crab seeds required for aquaculture in the country are fully collected from the wild stock of coastal mangroves. However, indiscriminate catching of mud crab from nature is not sustainable and has a potential risk of damaging the wild population. Mud crab has been selected as one of the target species for stock enhancement programs in different countries. For example, stock enhancement activities, such as broodstock management, induced spawning, crablet production, and acclimatization for release, have been accomplished by Japan, Malaysia, Philippines (Fushimi and Watanabe, 1999;Lebata et al., 2009;Ikhwanuddin et al., 2018;Asmat-Ullah et al., 2021). Recently, the government of Bangladesh has taken initiative to develop production technology for producing crablets in hatchery to serve the growing demand of seed, and initial responsibility for this has been delegated to the Bangladesh Fisheries Research Institute (BFRI). The institute has become successful initially with limited survivability of larvae and is continuing to scale-up the production (GIZ, 2017;Sakib et al., 2022).
Investigation of molecular genetic variation within species has become a common tool for studying patterns of biodiversity and diversification. Genetic variation within and among populations is vital for the management of endangered or commercially important taxa. Generally, a single panmictic population could be recovered through increased recruitment by propagation or short-term conservation initiative (Munro and Bell, 1997). On the other hand, different populations with unique genetic structures should be managed as distinct units, and such units require separate monitoring and management (Salgueiro et al., 2003). Maintenance of genetic variation is also an important concern for the management of hatchery stocks of a species and development of breeding programs (Fuji et al., 2007;Mastrochirico-Filho et al., 2019). Further, an understanding of the processes that lead to observed patterns of genetic diversity and structure in a species of interest can give insights into the mechanisms responsible for divergence, demographic history and speciation (Beheregaray and Sunnucks, 2001;Momigliano et al., 2017;Habib et al., 2022). Such insights can provide guidance for the management and conservation of vulnerable marine populations (Mathews, 2006). Further, it is also important to know the genetic diversity of each of the existing populations of a species. Genetic variability influences both the health and long-term survival of populations because decreased genetic diversity has been associated with reduced fitness, such as high juvenile mortality, diminished population growth, reduced immunity, and ultimately, higher extinction risk. To the best of our knowledge, no research has been conducted in Bangladesh on the intraspecific patterns of genetic diversity and the population genetic structure of S. olivacea using advanced molecular marker like mitochondrial DNA (mtDNA).
In the present study, we compare the population genetic differentiation of S. olivacea in the South-West (Sundarbans mangrove area) and South-East (Chakaria mangrove area) coast of Bangladesh using nucleotide sequence variation at cytochrome oxidase I (COI) gene in mtDNA marker. This represents the first appraisal of population structuring of S. olivacea in Bangladesh. Furthermore, we have also gathered all available COI sequence data from different scientific articles and GenBank to evaluate the levels of variation in genetic diversity and gene flow among different populations along with their historical demography throughout the Indian Ocean region.

Collection of samples and sequence data
A total of 48 orange mud crab individuals were collected from two locations namely, Sundarbans (SB) and Chakaria (CK), located on the South-West coast and the South-East coast of Bangladesh, respectively (  (Table 2). For Myanmar population, we randomly selected a single location (PZ) from 7 locations of the study of Segura-Garcia et al. (2018). In the study of Saher et al. (2019), authors wrongly identified nine sequences (Seq no: 55 to 63 as given in the article) of S. olivacea as S. serrata. We sampled these 9 sequences with more 10 sequences of S. olivacea from GenBank (Table 2).

DNA extraction, PCR amplification and sequencing
Using a DN easy blood and tissue kit (Qiagen, Germany), genomic DNA was extracted from the crab tissue following the manufacturer's protocol. The target region of the mtDNA COI gene was amplified from the genomic DNA through the polymerase chain reaction (PCR) using the forward and reverse primers, Mtd-10 5'-TTGATTTTTTGGTCATCCAGAAGT-3' and C/N 2769 5'-TTAAGTCCTAGAAAATG TTGRGGGA-3', respectively as used by Rosly et al. (2013). The PCR was performed using a 50 μl PCR master mixture containing 5 μl of 10×PCR buffer, 10 mM each of the dNTPs (2.5 μl each), 25 pmoles of each primer, 2.5 unit of Taq DNA polymerase, and 0.5-1.0 μg template DNA. The temperature profile was as follows: 5 min of preheating at 94°C followed by 35 cycles of denaturation at 50°C for 30 seconds, annealing at 72°C for 1 min, extension for 3 min at 94°C, and completion with final extension for 5 min at 72°C. The PCR products were visualized on 1% agarose gels, and by electrophoresis with a standard size marker and purified with a PCR purification kit (TIANGEN-China). The purified PCR products were   The sequence data were edited and aligned using Clustal Wimplemented in MEGA-6 ( Tamura et al., 2013). Molecular diversity indices such as nucleotide diversity (π), haplotype diversity (h), average number of nucleotide differences (k), number of haplotypes (N h ), polymorphic sites (S), transitions (ti) and transversions (tv) for each population were obtained using the program ARLEQUIN (version 3.5, Excoffier et al., 2005). Pairwise genetic variation and structure between sampled populations were assessed using fixation index (F ST ). An analysis of molecular variance (AMOVA; Excoffier et al., 1992)   was used to examine the level of genetic variability among populations that employs an approach to identify groups of populations genetically separated from one another. The statistical significance of the estimates was tested using 10,000 per mutations. The exact test of haplotype differentiation among populations (Raymond and Rousset, 1995) was also used to test the null hypothesis of population panmixia. All the above calculations were performed in ARLEQUIN. Each of the populations including two Bay of Bengal, Bangladesh populations (SB, CK) was treated as an independent unit in the population genetic analyses (Fixation index, Exact test, AMOVA) except the two populations of Pakistan (PK) and India (IN) which were pooled and treated as one population (designated as AS) considering that these localities were in the same sea basin (i.e.Arabian Sea) and the sample size of IN was very small (only 03), and all of the IN haplotypes share with PK haplotypes except a single sequence.

Phylogenetic and haplotype network analysis
The program MEGA 6 (Tamura et al., 2013) was used for reconstructing the phylogenetic relationship among haplotypes using the Neighbor-Joining (NJ) method (Saitou and Nei, 1987), first for only two Bangladeshi populations (SB, CK) and secondly, for all of the eight populations including two sequences (two haplotypes) of Australia (AU) population found in GenBank (Accession no: AY373355, AY373356). Robustness of the phylogenetic relationships was evaluated using bootstrap analysis with 10,000 replications (Felsenstein, 1985). Genealogical relations among haplotypes of all of the studied populations were further assessed by creating the network of haplotypes using a statistical parsimony method with TCS 1.18 (Clement et al., 2000). In the haplotype network, a haplotype is denoted by a circle and a line indicates one mutational step between haplotypes.

Demographic analyses
The demographic histories were examined by two different methods. Firstly, the neutrality statistics Fu's F S of Fu (1997) was used to test if the neutrality holds for the S. olivacea populations. Significant negative F S statistics can be interpreted as indications of demographic expansion. Secondly, historic demography of expansion of the populations was also investigated by examination of frequency distributions of pairwise differences between sequences by mismatch distribution analysis. The distribution is usually multimodal when a population is in demographic equilibrium or is subdivided into several units, but it is unimodal in populations following a recent population demographic expansion and range expansion (Slatkin and Hudson, 1991;Rogers and Harpending, 1992;Ray et al., 2003). In the mismatch distribution analysis, the history of demographic expansion was denoted by three parameters: τ (tau, date of onset of population expansion expressed in units of mutational time, θ 0 and θ 1 (mutational parameters of population size representing the population sizes before and after the expansion, respectively) (Rogers, 1995). The concordance between the observed and the expected frequency under sudden expansion model was further tested by means of the distribution with the Harpending'sreggedness index (Hri) (Harpending, 1994). The mutational time value of τ was transformed to estimate of real time in years with the equation, τ = 2ut, where u is the mutation rate per generation for the whole sequence and t is the time since expansion. The mutation rate was calculated from the divergence rate (D) by the equation, u = D/2. Divergence rate of 1.66-2.33% per million years (MY) for the COI gene was adopted as it was used for other brachyuran crabs' genus: Sesarma and Scylla by Schubart et al. (1998) and He et al. (2010), respectively. The generation time of about 0.5 year (average) was mentioned for S. paramamosain by Ng (1998); Palomares and Pauly (2022) which was used in calculation of the expansion time (years) in the present study.

Genetic diversity and haplotype distribution
The length of COI gene region sequences obtained after removing the ambiguous sequences close the primers at either end from 48 individuals of two Bangladeshi populations (viz.SB and CK) was 512 base pairs (bp). The COI sequences of different populations other than two Bangladeshi populations (SB, CK) collected from different studies as well as GenBank were about 540 to 650 bp in length. These sequences were aligned with the sequences of Bangladeshi populations and only the aligned 512 bp portion of the sequences were used.
In the 48 sequences of two Bangladeshi (SB and CK) samples, a total of 26 haplotypes with 33 polymorphic sites containing 26 transitional and 7 transversional changes were identified. On the other hand, all of the 213 sequences of seven populations comprised a total of 69 haplotypes. Most of the nucleotide substitutions of the sequences were transitional and no indels were detected in the sequences (Table 3).
The nucleotide diversities (π) were very low in each sample, 0.005-0.01 nucleotide differences per site while the haplotype diversities (h) were high, 0.70-0.96 in all populations. Overall, the mean h among the 48 Bangladeshi samples (SB & CK) was estimated to be 0.93, and the mean π was 0.009. On the other hand, these values (i.e. mean h and π) were 0.93 and 0.008, respectively among the 213 samples of seven populations analyzed in the present study (Table 3).
High levels of haplotype diversity in contrast with very low nucleotide diversity (h>0.5 and π<0.5) of Bangladeshi populations (SB and CK) and other five populations of mud crab analyzed in the present study (Table 3) indicates that S. olivacea has experienced rapid population expansion after bottleneck with accumulation of mutations (Grant and Bowen, 1998). In the case of a rapid expansion of a population, the rate of stochastic loss of haplotypes slows down and more haplotypes become retained than lost by random genetic drift (Avise et al., 1984). This kind of genetic make-up has also been observed in other crab species such as shore crab, Pachygrapsus crassipes (h = 0.92, π = 0.009) (Cassone and Boulding, 2006), swimming crab, Portunus trituberculatus (h = 0.66, π = 0.10) (Xu et al., 2009a) and mangrove crab, Perisesarma guttatum (h = 0.77, π = 0.004) (Silva et al., 2010) and mud crab, Scylla paramamosain (h = 0.841, π = 0.002) (Ma et al., 2011). Another factor that maintains high haplotype diversity is environmental heterogeneity (Nei, 1987). S. olivacea occurs in the estuaries and brackish water rivers and channels of mangrove areas throughout the Indo-Pacific region from Southeast Asia to the Bay of Bengal and Arabian Sea, and also from Japan to northern Australia. Such large environmental heterogeneity might have caused the high levels of haplotypic diversity observed for mud crab S. olivacea. High haplotype diversity has also been suggested for large and stable effective population size over time (Stepien, 1999).

Population genetic structure
We calculated pairwise F ST (fixation index) between each pair of populations to test genetic differentiation among populations. Genetic differentiation between MM (Myanmar) and each of the other 6 populations exhibited significant genetic structure (P<0.05). The MM population even reveals separate breeding unit (significant exact P values) from all other populations. Andaman Sea is a nearly enclosed marginal sea bounded by the coastlines of Myanmar, Thailand, and west side of the Malay Peninsula in the North and East side. It is separated from the Bay of Bengal to its west by the Andaman Islands and the Nicobar Islands. Its southern end is at Breueh Island just north of Sumatra. Such land boundaries around the Andaman Sea should limit the movement or dispersal of mud crab between inside and outside which causes the unique population structure in the Andaman Sea. Further, a clockwise gyre during winter monsoon (December to March) and an anticlockwise gyre during summer monsoon (May to October) are created inside sea (G1 & G3 in Fig. 2; Varkey et al., 1996). In the Andaman Sea, S. olivacea spawn all the year round, but peak spawning occurs from the onset of rainy season to the early winter from June to November (Koolkalya et al., 2006) which largely coincide with summer monsoon periods. In the year-round breeding period including peak season the larvae may potentially be transported for considerable distances but, alternatively, these rotational currents (i.e. gyres) prevent dispersal which is also responsible for creating unique population genetic structure within this sea. Though this sea is connected with South China sea through the Strait of Malacca in its southeast side, another anticlockwise gyre (G2) is located in the north Malacca Strait near north Sumatra (Fig. 2) during winter (Rizal et al., 2012) which also substantially prevents the movement of larvae and make separate population genetic structure between these sea areas (MM vs. MS & SC) ( Table 4). The F ST values also showed that there is no population genetic variation among SB, CK and AS (i.e. PK+IN) populations. Larval transportation of S. olivacea through the West India Coastal Current (WICC) and the East India Coastal Current (EICC) of the northern Indian Ocean (Fig. 2) are probably responsible for genetic homogeneity among these populations.
The genetic structure of the S. olivacea was further investigated by hierarchical AMOVA that took into account geographical distribution (different sea environment). The AMOVA result showed significant genetic variation (6%, F CT = 0.06; P= 0.04) when the populations were divided into four groups (Fig. 3) viz. Arabian Sea (AS i.e. PK+IN), Bay of Bengal, BoB (SB+CK), Andaman Sea (MM) and South China Sea (SCS i.e. MS+SC) region including Malacca Strait. Three clear geographic barriers in the Indo-west Pacific sea areas restricted the gene flow and made these four population groups. As shown in Fig. 3, the main island of Indian Peninsula acts as a physical barrier (Barrier-1) dividing the Bay of Bengal and Arabian Sea into two large marine area. This geographical separation of marine habitats and long distance (i.e. isolation by distance) between populations limit gene flows between AS and BoB. In addition, there is a 48 km long land link namely Adam's Bridge, formed by limestone shoals, coral reefs, and tombolo situated between Pamban Island off the southeastern coast of Tamil Nadu state of southern India, and Mannar Island off the Northwestern coast of Sri Lanka (Weerakkody, 1997;Fig. 3). This Bridge together with Pamban Island, Mannar Island, and Sri Lanka formed about 500 km long land barrier which might also act as associated barrier of Barrier-1, and prevent the migrations and dispersal of adult and larvae of mud crab between AS and BoB sea basins and contributes to form separate genetic population groups within each sea. Secondly, the Andaman-Nicobar Ridge and an associated chain of sea mounts along the Indo-Burmese plate boundary in the Indian Ocean which made the Andaman Sea as a semi-closed marginal seaby partly isolating from the Bay of Bengal act as a barrier (Barrier-2) and prevent substantially dispersal of mud crab between BoB and MM populations. The anticlockwise gyre (G2) located at north of Sumatra (Fig. 2), in the north Malacca Strait acts as another barrier (Barrier-3) of gene flow between Andaman Sea and South China Sea region (i.e. between MM and SCS population group).

Phylogeny and haplotype network
The phylogenetic reconstruction of the haplotypes of two Bangladeshi populations (i.e. SB and CK) showed clear divisions making two clades (Clade-1 and Clade-2) (Fig. 4). Clade-1 was shared by the haplotypes of both populations but dominant with the CK population (56%) and Clade-2 was highly dominated by the haplotypes of SB population (67%). The phylogeny also shows that only four haplotypes were shared by two populations. On the other hand, the topology of the neighbor-joining tree of 70 COI haplotypes of seven populations (SB, CK, MM, MS, PP, IN and AU) showed no significant geographical clusters in any clade belonging to any specific sampling locality (Fig. 5). Haplotypes of a population were distributed throughout the tree and well mixed with those from other populations.
Haplotype network (i.e. parsimony networks of haplotypes) showed two haplogroups ( Fig. 6) with star-burst structures within each group. The star-like shape of the network suggests recent population expansion of each group. No geographical association was found for each group. Haplotype frequencies of the COI sequences among different populations and their GenBank accession number are given in Table 5. In general, the most dominant haplotype (H10) was shared with 16.74% (36/215) individuals from six populations (localities), 6 haplotypes (H8, H14, H17, H12, H4, H29) were shared by 4 to 5 populations, 7 haplotypes (H23, H33, H36, H11, H20, H48, H30) were shared by 2 to 3 populations and the rest 56 haplotypes were unique to their own population (  Fig. 6).

Demographic history
Estimates from Fu's F S test were significantly negative (P<0.05 i.e. population expansion) for SB, MS, SC populations ( Table 6). The Fs values for CK, MM, and AS populations were not significant (P>0.05). Therefore, the neutrality tests established that these populations have not experienced demographic expansion; rather these populations have undergone demographic equilibrium. We have also studied demographic history of S. olivacea by another approach namely mismatch distributions (MMD) analysis of the six populations. Though the observed mismatch distributions somewhat looked like multimodal pattern for each of the populations (Fig. 7) indicating the history of demographic equilibrium, the goodness-of-fit tests were consistent with the hypothesis of demographic expansion for SB, MS and SC [Hri(P)>0.05] as found in Fu's F S test (Table 6). On the other hand, raggedness index (Hri) was found significant (P<0.05) for CK and MM suggesting demographic expansion for these two populations which is also consistent with Fu's F S test. Only contradictory result between Fu's Fs test and MMD analysis was found for AS population. Ramos-Onsins and Rozas (2002) showed that the Fu's F S test is more sensitive than the analysis of mismatch distribution. Therefore, interpreting between F S statistics and MMD analysis finally confirmed that demographic expansion occurred in the AS population.Conversion of the τ values of these four expanding populations (SB, AS, MS, SC) into the time in years since expansion resulted approximately 125,000 to 175,000 years ago for SB, 200,000 to 280,000 years ago for AS, 256,000 to 359,000 years ago for MS, and 260,000 to 365,000 years ago for SC population during the late Pleistocene (the past 1 million years). The late Pleistocene period of glacial cycles was episodic by a series of large glacial-interglacial changes (Imbrie et al., 1992). Glaciations-deglaciations and their interface, and associated temperature and salinity shifts, global changes in ocean circulation patterns, and altered productivity regimes of this historical time period probably had great effects on the demographic history that caused population expansion of S. olivacea in the above mentioned four populations (SB, AS, MS, SC) of Indo-west Pacific. Since the Andaman Sea is almost enclosed sea blocking mud crab larval dispersal, outside of sea, demographic history of MM population has been characterized with long stable neutrality. In a number of other crab species in the coastal waters of the Western Pacific, the population expansion occurred in a similar time frame as the orange mud crab showed, such as green mud crab, S. paramamosain (91,900 to 161,200 years, He et al., 2010), swimming crab, Portunus trituberculatus (127,000 and 429,000 years, Xu et al., 2009a), mitten crab Eriocheir sensu (25,000 and 130,000 years, Xu et al., 2009b).
Assessment of genetic structure of the present study will provide a baseline from which mud crab populations can be appraised more in the future because assessment of population genetic structure is a vital tool for maintaining a productive fishery while managing a species as a resource (Seeb et al., 1990). The present mtDNA analysis from historic process and contemporary factors provides no evidence for genetic subdivision of S. olivacea within Bangladesh coast, however, Bangladeshi populations showed genetically different structure from some other sea areas of the Indian Ocean. This study would be useful for genetic stock estimation, sustainable exploitation, conservation and management, and artificial breeding program of this important crab species not only within the country but also in the global conservation aspect.