Assessment of Genetic Variation in COI Gene of Horseshoe Crab (Tachypleus Gigas) Population Along the Coast of the Bay of Bengal in Odisha, India

Horseshoe crab (Tachypleus gigas) is an archaic group of marine creature which plays a vital role in the saline ecosystem. Many researchers emphasize and enhance the knowledge about the horseshoe crab's basic biology, morphology, and ecology, whereas very little information is available about its population genetics. We attempted to develop a baseline database about the ecology, phylogeography, and genetic variation among the horseshoe crab population from Odisha, India. We collected 152 samples of horseshoe crab from the coastal area of the Bay of Bengal. The generated Cytochrome C Oxidase Subunit I gene (COI) sequences of T. gigas were compared with the sequences of T. gigas obtained from GenBank. The GenBank sequences were of two populations from South China and Malaysia. A total of 26 unique haplotypes were observed in three populations of T. gigas. Pairwise F-statistic distance (F ST ) between South China-India was 0.708; Malaysia-India was 0.608, and South China-Malaysia was 0.136. It indicated that the South China population was closely related to the Malaysian population and the Indian population was appeared to be genetically distinct from the other two populations. It signies the ecological importance of the Indian population. Furthermore, the migrant per generation (Nm) was 0.16, which indicated a low gene ow among T. gigas populations. The haplotype diversity (Hd) and nucleotide diversity (π) were 0.58826 and 0.00476, respectively. This study would help lay future strategy and conservation of horseshoe crab across the Bay of Bengal. and three out-group sequences of horseshoe crabs were derived from GenBank. Two different methods, maximum likelihood (ML) and Bayesian inference (BI) were adopted for phylogenetic analysis. In RAxML-VI v.7.0.0. Software analysis ML tree estimation was performed [40]. A Bayesian phylogenetic tree generated beast. Based on the average haplotype diversity, genetic diversity (π) was estimated. Population pairwise nucleotide differences and their statistical signicance were checked using Arlequin v.3.5.1.2 with 10,000 permutations. Geographical distribution and genetic differentiation of the haplotype population within the species were calculated by using PopArt software. respectively. We compared mtDNA genetic variation in the Indian horseshoe crab population in the Bay of Bengal with other Tachypleus gigas (South Chinese and Malaysian) populations across its distribution range. We found that Indian horseshoe crab has a comparatively lower level of genetic variation (h=4, s=3, mean pairwise difference =1.53) than the well-connected populations of South China and Malaysia. subdivision was the probable reason for lower genetic variations. It suggested signicant population genetic differentiation throughout the range of T. gigas along the coast of the Bay of Bengal. Among the 26 haplotypes, the major two H2 and H3 dominating subpopulations were observed on the Bengal coast. The haplotypes H13, H26, H27, H16 and H21 were one mutation away from the H2 haplotype and the rest H18, H22, H17, H20, H19, H23, H25, H15 and H24 were also one mutation away of H3 haplotype that means other haplotypes are just ancestor of both haplotypes. This analysis revealed that Tajima’s D was − 2.27032 (P < 0.01), Fu and Li’s was–4.736415 (P < 0.02) while Fu’s Fs value was − 22.25092 (P < 0.05). Our data suggested that the T. gigas has experienced a sudden demographic expansion. In the mismatch distribution, a smooth and unimodal curve was observed among the 14 geographic populations. About 65 − 23 Mya during the Paleogene era, the Asian horseshoe crabs were diversied and split. Our results strongly supported the previous studies in which it was mentioned that the American species L. polyphemus was the sister taxon of three Asian species [21, 45, 48] and the Asian species was diversied around 55 Mya. As a rst Asian species, C. rotundicauda (South Korean) was split approximately 50 (80 − 20) Mya. T. gigas showed a closer genetic distance with T. tridentatus than the C. rotundicauda. The Indian and South China species (T. gigas) divergence occurs during15-2 Mya. T. tridentatus has diverged in around 65 − 15 Mya. C. T. crab era − 23 Mya). two subpopulations


Introduction
Horseshoe crabs are one of the archaic groups of marine organisms with signi cant economic value and scienti c signi cance. These animals play an important role in both marine and estuarine ecosystems. Four existing species of horseshoe crabs are reported throughout the world. Of these, Limulus polyphemus [1] is found on North America Atlantic coast, whereas other three species of horseshoe crabs like Tachypleus gigas [2], Tachypleus tridentatus [3] and Carcinoscorpius rotundicauda [4] are found in the Asian countries ranging from 6 o S to 31 o N and between 90 o E and 118 o E [5]. Out of the four horseshoe crab species, Tachypleus gigas and Carcinoscorpius rotundicauda have been reported with overlapping distribution in many regions in the Bay of Bengal, India and South China Sea [6,7]. Tachypleus gigas prefer sandy areas and C. rotundicauda is found in muddy regions [8][9][10][11][12]. However, trispine Asian horseshoe crab Tachypleus tridentatus have a preference to occur in the coastal zone of shallow water in southern China ranges from Vietnam to Japan. In the Indian sub-continent, the maximum numbers of horseshoe crabs are distributed along the Bay of Bengal coast of India, ranging from Sundarban mangrove forest of West Bengal state to Dhamara of Odisha state [13,14]. In Odisha state, the coastal region of the Balasore district is richly populated with horseshoe crabs. Horseshoe crab inhabits three regions like the coastal bay, intertidal mud zone or sand at zone, and deepwater regions [15]. Balasore has quali ed as a preferable place for breeding due to the intertidal zone, which helps for their spawning activity [10]. It has been reported that the T. gigas and C. rotundicauda come to the estuary for spawning purposes in pairs during the tidal movement of spring seasons in India [16,17]. The population of horseshoe crabs is declining rapidly due to over shing, water pollution, habitat destruction, resulting in reduced abundance and quality of prey available as a food source to horseshoe crab [6]. Hence, it is well accepted and considered an endangered marine organism in India and other regions [18][19][20]. Several controversies have been raised among the evolutionary biologist about the horseshoe crab phylogeny from the last several years. They are often considered the foremost example of long-term survival in marine. These marine organisms still retain their body morphology and stable genetic structure over millions of years [21][22][23], hence called as 'Living fossil' [24], 'Phylogenetic relict' [25] and becoming endangered marine animals. Generally, in the marine ecosystem, organisms exhibit strong gene ow for their dissemination activity across the spacious geographic ranges [26]. Some habitually display a low level of genetic differentiation for their incapability of long-distance dispersal over the large geographic scale [27]. Horseshoe crabs are lean migratory animals that are probable to have subdivisions population on a small geographic scale [26, 28]. They reveal a low level of genetic diversity throughout their entire evolutionary history. The allozyme [25,29], restriction fragment length polymorphism (RFLP) of mtDNA [30], and microsatellite analysis [31] revealed that the population of American horseshoe crabs (Limulus polyphemus) was genetically distinct in the Gulf of Mexico and the Atlantic coast. In a small geographic range of North America, limited gene ow was observed through COI gene analysis [28]. Two genetic groups were detected based on the distribution of haplotypes and the signi cant genetic structure of Japanese horseshoe crabs (T. tridentatus) populations by the mtDNA AT-rich region [32]. A few research studies have been conducted to explore the phylogeny and phylogeography of four extant species using various molecular methods. However, very little knowledge is available about the horseshoe crab's population genetics diversity along the Bay of Bengal coast in India. A molecular study is a powerful tool for investigating genetic differentiation and it has widely been used to determine population genetics, phylogeography and genetic structure. It possesses a rapid rate of evolution, lack of recombination, and maternal inheritance [33]. Mitochondrial Cytochrome C Oxidase subunit I (COI) is often considered a conserved region and used as a robust marker for the invertibrates and marine species due to a high degree of polymorphism leading to population subdivisions (Hu et al. 2009) [34]. The objective of this study was to examine molecular variations, genetic diversity and historical demography of Indian horseshoe crab in the select area of the Bay of Bengal coast.

Sample collection and DNA extraction
We collected 152 tissue samples from the dead remains of horseshoe crabs from fourteen sampling sites of the coastal areas of the Bay of Bengal (Table  1). Samples were stored at −20 °C in 70% ethanol for molecular analysis. Genomic DNA was extracted from the tissue sample using the standard Proteinase K digestion/phenol-chloroform method [35], and a standardized protocol was adopted to extract DNA [36]. The extracted DNA was quanti ed in a QIAxpert spectrophotometer. Isolated DNAs were diluted in a nal concentration of 40-50 ng/μl for PCR ampli cation and stored at 4°C.

PCR ampli cation and sequencing
We used conserved primers for PCR ampli cation and sequencing of the mtDNA COI gene [37]. The total 20 μl reaction volumes were prepared for PCR reactions using 1 × PCR buffer (10 mMTris-HCl, pH 8.3,and 50 mM KCl), 1.5 mM MgCl 2 ,1 × BSA,0.2 mM of each dNTPs,3 pmol of each primer, 0.5 units of DreamTaq DNA Polymerase (Thermo Fisher) and 1 μl (~ 30 ng) of template DNA. The PCR conditions were initial denaturation at 95 °C for 5 minutes and 95°C for the 40s, annealing temperature at 45°Cfor 50s,72°C for 50s followed by 5 cycles and second denaturation at 95 °C for 40 s, annealing at 51 °C for 50 s and extension at 72 °C for 50 s. The nal extension was at 72 °C for 10 min, followed by 35 cycles, and maintained at 4℃ until further use. For avoiding any cross-contamination, e ciency, and reliability, control and the negative reaction were taken during DNA extraction and PCR reactions. The PCR ampli cation was con rmed by electrophoresis on 2% Agarose gel stained with safe stain dye (0.5 mg/ml) and visualized under a UV transilluminator. The exonuclease-I and shrimp alkaline phosphatase (Thermo Scienti c) was treated at 37°C for 20 min, followed by 85 °C for 15 min with the ampli ed PCR product to obtained intact DNA. The sequencing was performed from forward and reverse direction using BigDye v3.1 Kit in 3500 XL Genetic Analyzer (Applied Biosystems).

Data alignment and molecular analysis
The obtained DNA sequences were scrutinized for similarity and retrieved through BLAST sequences accessible from NCBI GenBank. Later the sequences were aligned by using Clustal X available in the BioEdit packages [38]. The IUPAC nucleotide was clean in BioEdit software. The numbers of variable sites and population diversity were analyzed using DNASP. A percentage similarity matrix and pairwise distance matrix were generated using MEGA X [39]. Evolutionary analyses were conducted in MEGA X. The haplotype median-joining network was generated by Network v.5.0.0.0 as well as PopArt software. In ARLEQUIN v3.11 software analysis of molecular variance (AMOVA), Tajima's D parameter estimations, gene ow, and mismatch distributions were estimated. Neighbour Joining, Maximum Parsimony and Maximum Likelihood tree were constructed based on 581 bp of the mitochondrial DNA COI gene from T. gigas, and three out-group sequences of horseshoe crabs were derived from GenBank. Two different methods, maximum likelihood (ML) and Bayesian inference (BI) were adopted for phylogenetic analysis. In RAxML-VI v.7.0.0. Software analysis ML tree estimation was performed [40]. A Bayesian phylogenetic tree was generated in the beast. Based on the average haplotype diversity, genetic diversity (π) was estimated. Population pairwise nucleotide differences and their statistical signi cance were checked using Arlequin v.3.5.1.2 with 10,000 permutations. Geographical distribution and genetic differentiation of the haplotype population within the species were calculated by using PopArt software.

Divergence time estimation
The updated version BEAST v.1.7.5. was used for estimation of divergence among the horseshoe crabs species [41]. The single COI gene sequence of each species is considered for analysis. We took T. gigas from the Bay of Bengal coast and another from the South China Sea, expecting a genetic gap in both the population for the large geographic distance. The remaining three species, namely L. polyphemus was from Mexico, T. tridentatus from the China Sea, C. rotundicauda from South Korea, for divergence estimates. The best model HKY was selected for this analysis from the JModel test and the base frequencies were not considered for estimation. We also took equal to empirical values. It was expected to evolve a phylogenetic tree with Yule speciation tree priors and a log-normal uncorrelated tranquil clock [42]. For all these analyses, 25 million generations made for a run. With the help of Tracer v.1.5, the Estimated Sample size (ESS-500) was checked and calculate the value in between the 200-600 ranges of all parameters [40]. We calibrated our tree with the divergence times described by Obst et al. (2012) [43]. According to Fisher [21]   between South China-India was 0.708; Malaysia-India was 0.608, and South China-Malaysia was 0.136. We also con rmed the intraspecies genetic distance, ranging from 0.00 to 0.05, while inter-species distance varied from 0.06 to 0.08 for the species under study. The highest distance was observed between T. gigas and L. polyphemus.

Genetic variation and haplotype diversity
To estimates the diversity indices, we used only the COI gene for molecular analysis. We detected two major haplotypes (H2 and H3) distributed in all sampling sites, while some local haplotypes were speci c to the particular sites. Molecular diversity index such as haplotype (h) and nucleotide (π) diversity was estimated among the horseshoe crab's population in the Bay of Bengal. Nucleotide (π) and haplotype (h) diversity were found 0.00 to 0.00059 and 0.00 to 0.534, respectively. Overall nucleotide diversity and haplotype diversity in the Indian horseshoe crab population in the Bay of Bengal was 0.00265 and 0.578, respectively. We compared mtDNA genetic variation in the Indian horseshoe crab population in the Bay of Bengal with other Tachypleus gigas (South Chinese and Malaysian) populations across its distribution range. We found that Indian horseshoe crab has a comparatively lower level of genetic variation (h=4, s=3, mean pairwise difference =1.53) than the well-connected populations of South China and Malaysia.

Historical demography
analysis. This analysis reveals that Tajima's D was -2.27032 (P <0.01), Fu and Li's was-4.736415 (P <0.02), while Fu's Fs value was -22.25092(P <0.05). Our data suggest that the T. gigas has experienced a sudden demographic expansion. In the mismatch distribution, a smooth and unimodal curve (Fig. 2) was observed among the fourteen geographic populations. The analysis of mismatch distribution was consistent with Fu's Fs statistical test for the population expansion model.

Phylogenetic analysis
Based on the partial sequence mtDNA COI gene, the phylogenetic tree was established. The best model was also considered as GTR+I from the JModel test.
Only T. gigas species were included in the time of phylogenetic analysis. 17 unique haplotypes were obtained from 152 samples. These 17 haplotypes and ten sequences of T. gigas from NCBI GenBank were included for phylogenetic analysis. Interestingly our sample data categorized into two major subpopulations. 11 haplotypes fall under one subpopulation and the remaining six were in another group. In contrast, these were believed a sister group relationship with each other. Each phylogenetic partition was signi cantly different from one another (P<0.01). Its overall analysis yielded constant phylogeny. In our study, open and highly supported the connection between Tachypleus gigas and Tachypleus tridentatus was observed than the C. rotundicauda, which showed sister taxon relationship, i.e., monophyletic genus Tachypleus. It was also strongly supported for Tachypleinae [44] between Tachypleus and C. rotundicauda to establishing phylogenetic clade as a monophyletic relationship.

Divergence time estimates
We obtained 150-130 Mya divergence between Limulidae-Tachypleinae. The diversi cation of Asian horseshoe crab Tachypleinae occurred during 50 Mya

Discussion
The marine biodiversity crisis has increased due to heavy anthropogenic activity and resulted in the depletion of habitats for marine organisms. Some of them have got extinct without our knowledge or ever having registered their existence. To assess the divergence time, genetic differentiation, phylogeographic pattern, and genetic structure of the horseshoe crab populations, we used the mtDNA gene (COI) for molecular analysis. Our nding was congruent with a previous study [45], which displayed the morpho-static structure of horseshoe crabs due to the lack of variation in molecular genetics. Horseshoe crab species showed a slow divergence rate compared to other arthropods [46] due to their long generation time and sexual maturity after ten years [47]. We included three major T. gigas populations such as India, South China and Malaysia. The South China population was closely related to the Malaysian population, and the Indian population was appeared to be genetically distinct from the other two. It can be explained due to its least dispersal capability in the marine ecosystem [26]. We con rmed the intraspecies genetic distance ranged from 0.00 to 0.05 while inter-species distance varied from 0.06 to 0.08. The highest genetic distance was observed between T. gigas and L. polyphemus. The L. polyphemus was restricted in the Gulf of America and the Atlantic coast only. It showed different genetic appearances. Due to this characteristic, the geographical barrier played an important role in dividing the genetic break for which they may not participate in genetic exchange and the species. The changes in sea level are the major factors to acclimatize/adapted the species to the particular environment resulting in genetic diversity. Hence, we can assume that they are pretty far from the genetic distinct of this species. We also conclude that the Indian horseshoe crab has a comparatively lower level of genetic variation (h = 4, s = 3, mean pairwise difference = 1.53) than the well-connected populations of South China and Malaysia. But it was quite similar to the geographically distinct populations. Therefore, it might also be possible that the Indian horseshoe crab population's isolation due to habitat subdivision was the probable reason for lower genetic variations. It suggested signi cant population genetic differentiation throughout the range of T. gigas along the coast of the Bay of Bengal. Among the 26 haplotypes, the major two H2 and H3 dominating subpopulations were observed on the Bengal coast. The haplotypes H13, H26, H27, H16 and H21 were one mutation away from the H2 haplotype and the rest H18, H22, H17, H20, H19, H23, H25, H15 and H24 were also one mutation away of H3 haplotype that means other haplotypes are just ancestor of both haplotypes. This analysis revealed that Tajima's D was − 2.27032 (P < 0.01), Fu and Li's was-4.736415 (P < 0.02) while Fu's Fs value was − 22.25092 (P < 0.05). Our data suggested that the T. gigas has experienced a sudden demographic expansion. In the mismatch distribution, a smooth and unimodal curve was observed among the 14 geographic populations. About 65 − 23 Mya during the Paleogene era, the Asian horseshoe crabs were diversi ed and split. Our results strongly supported the previous studies in which it was mentioned that the American species L. polyphemus was the sister taxon of three Asian species [21,45,48] and the Asian species was diversi ed around 55 Mya. As a rst Asian species, C. rotundicauda (South Korean) was split approximately 50 (80 − 20) Mya. T. gigas showed a closer genetic distance with T. tridentatus than the C. rotundicauda. The Indian and South China species (T. gigas) divergence occurs during15-2 Mya. T. tridentatus has diverged in around 65 − 15 Mya. Although, C. rotundicauda has a contrasting pattern of connectivity with T. gigas, which can reproduce its isolated lineage. The diversi cation of Asian horseshoe crab occurred during the Paleocene era (65 − 23 Mya). We con rmed two distinct horseshoe crab (T. gigas) subpopulations along the Bay of Bengal Bay in India.
The molecular analysis and evolutionary relationship among the Indian horseshoe crab (T. gigas) species signi cantly displayed a more distinct population. Additionally, in the phylogenetic tree, the genetic relationship of Tachypleus gigas is extreme nearer to Tachypleus tridentatus other than the Carcinoscorpius rotundicauda species. Within the T. gigas species of India and Southeast Asian populations, a restricted gene ow was observed. Interestingly, in the horseshoe crab population of India and the Southeast Asian population, a limited gene ow was observed. However, more detail and further scienti c analysis are required to investigate and establish the genetic structure and the evolutionary history of the Indian Horseshoe crab population. An advanced molecular technique should be adapted to conserve valuable marine organisms. The Indian horseshoe crab was geographically distributed in a limited area on the Bay of Bengal coast, but now the current status is under severe exhaustion across all over its range even species is vanishing in many estuaries. Due to anthropogenic activities, seawater pollution, over shing and coastal development caused the loss of habitat vegetation. Accordingly, suitable conservation strategies should be adopted to protect the species and its habitats from extinction.

Availability of data and material
The sequence data is submitted to GenBank under Accession No. MZ262422-MZ262447.
Ethical statement: The study was conducted using a tissue collected from the naturally dead animal and therefore, no Institutional Animal Ethics Committee approval was required.

Consent to Participate (Ethics)
Since the experiment was not conducted on humans, no consent to participate was required.

Consent to Publish (Ethics)
This study has not used any secondary data for publication, hence no consent to publish was required.      C G  C G  A T  G  T  C T  G  G  G  G  G  G  C T  G  T  C T  A  T  C T   H2  Location of 14 sampling sites along the coast of the Bay of Bengal in Balasore District, Odisha. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.  Phylogeographic tree for the COI fragment. The color represents the localities.

Figure 4
Chronogram based on the BEAST divergence time analysis showing the estimated evolutionary radiation of Asian horseshoe crabs. Blue bars indicate 95% con dence levels. The times scale below the chromatogram features all major geological periods in Mya.