Lack of mitochondrial genetic structure in the endangered giant clam populations of Tridacna maxima (Bivalvia: Cardiidae: Tridacninae) across the Saudi Arabian coast

1. Introduction

The evolutionary potential of species depends mainly on the genetic variation of their populations. The amount of genetic diversity within and among natural populations is mainly explained by several historical and contemporary processes, including genetic drift, effective migration, natural selection, fragmentation and range expansion (Slatkin, 1985). Genetic polymorphism, resulting from the potential balance between evolutionary and demographic processes, may generate either heterogeneity or homogeneity among local populations.

In the marine environment, genetic variability and structure of benthic populations are potentially shaped by both historical climatic fluctuations (during glacial and interglacial periods of the Pleistocene) and contemporary processes (including larval dispersal potential, oceanographic circulation patterns, substrate availability and recruitment success) (Nikula and Väinölä, 2003; Patarnello et al., 2007; Maggs et al., 2008; Weersing and Toonen, 2009; Hunter and Halanych, 2010; Reuschel et al., 2010; Deli et al., 2016, 2018).

The Red Sea has recently attracted considerable attention in phylogeographic studies due to its particular topography as well as high levels of genetic diversity and endemism characterizing its biota (DiBattista et al., 2013, 2016; Hui et al., 2016). This semi-enclosed tropical basin is restrictly connected to the Indian Ocean through the Strait of Bab el Mandeb. Furthermore, it was repeatedly isolated during sea-level low stands of the Pleistocene glacial cycles (DiBattista et al., 2013, 2016). Hence, it is an ideal natural laboratory to investigate mechanisms of genetic divergence and endemism of marine fauna due to the interplay between impacts of Pleistocene climate shifts and contemporary environmental barriers. Patterns of restricted gene flow between the Red Sea populations and their Indian Ocean counterparts were detected in various marine species (i.e., the mud crab Scylla serrata (Fratini and Vannini, 2002), the damselfish Chromis viridis (Froukh and Kochzius, 2008), the sponge Leucetta chagosensis (Wörheide et al., 2008), and the spiny lobster Panulirus penicillatus (Iacchei et al., 2016)). These findings, associated with the emerging discovery of new endemic fauna to the Red Sea (Richter et al., 2008; DiBattista et al., 2016), allowed considering this marine basin as an important and distinct centre of evolution (Klausewitz, 1989; DiBattista et al., 2016). Furthermore, genetic isolation was unveiled in the Red Sea marine species, such as in the anemone fish Amphiprion bicinctus (Nanninga et al., 2014) and the reef sponge Stylissa carteri (Giles et al., 2015). Such pattern of population genetic structure was attributed to the particular circulation pattern as well as to the strong latitudinal temperature and salinity gradients characterizing the Red Sea (Sofianos et al., 2002; Sofianos and Johns, 2007; Raitsos et al., 2013; Kürten et al., 2014). These important insights, along with the scarcity of information on patterns of genetic structure and gene flow within the Red Sea, trigger the necessity of detailed genetic investigation in order to disentangle the potential impact of both historical and contemporary factors on the genetic makeup of marine fauna.

The giant clam Tridacna maxima (Röding, 1798) represents a good model to address the above mentioned issues. This bivalve species is one of the most ecologically important coral reef species (Soo and Todd, 2014). It has a wide geographic distribution, stretching from the Red Sea and Indian Ocean across the Indo-Malay Archipelago to the central and western Pacific (Gilbert et al., 2007; Andréfouët et al., 2014). Adults of T. maxima are sedentary. They reproduce by broadcast spawning with high fecundity (releasing up to millions of eggs per female at each spawning (Lucas, 1988; Soo and Todd, 2014)). Besides, they have a short pelagic larval dispersal of about nine days (Lucas, 1988). Reproduction within these simultaneous hermaphrodites is sporadic, and natural recruitment of juveniles is often very low (Shau-Hwai and Yasin, 1998). Taking into account all of these eco-biological aspects, considerable genetic differentiation is highly expected among populations of T. maxima. Various population genetic studies, reporting on nuclear and mitochondrial DNA variation in T. maxima, supported this assumption (Laurent et al., 2002; Juinio-Meñez et al., 2003; DeBoer et al., 2008; Kochzius and Nuryanto, 2008; Nuryanto and Kochzius, 2009; Ahmed Mohamed et al., 2016; Hui et al., 2016). Indeed, restricted gene flow was revealed within and among various Indo-Pacific giant clam populations. It was postulated to be caused by either ocean current patterns (Benzie and Williams, 1992a; Macaranas et al., 1992) or geographic isolation (Benzie and Williams, 1992b). Recent large-scale investigations allowed delineation of deep evolutionary lineages and strong population genetic structure in T. maxima. In particular, limited gene flow was unveiled between the Western Indian Ocean, the Red Sea, the Eastern Indian Ocean, the central Indo-Malay Archipelago, the Western Pacific and the Central Pacific (DeBoer et al., 2008, 2014a, 2014b; Kochzius and Nuryanto, 2008; Nuryanto and Kochzius, 2009; Hui et al., 2016). The outcomes of these latter genetic investigations suggest the involvement of Quaternary sea-level oscillations and oceanographic barriers in shaping the contemporary population genetic structure of T. maxima. Despite the existence of such numerous population genetic surveys, little is still known about patterns of gene flow within peripheral areas of the distribution zone of the giant clam, i.e., the Red Sea. Furthermore, for comparative phylogeographic purpose, previous studies did only focus on few specimens collected from one location across this geographic spectrum (Nuryanto and Kochzius, 2009; Hui et al., 2016).

In light of these considerations, we intended in this study to investigate the population genetic structure of this poorly dispersive bivalve species across the Red Sea. We aimed at testing the hypothesis that environmental heterogeneity within this semi-enclosed Basin can generate possible genetic heterogeneity. Furthermore, we tried to disentangle the potential processes susceptible of shaping contemporary genetic variation of T. maxima. In order to achieve these targets, five locations, covering almost the entire Red sea off the Saudi Arabian coast, were sampled and analyzed genetically by means of a fragment of the mitochondrial gene cytochrome c oxidase I (COI). This genetic marker was shown to be variable enough for population studies in marine bivalves (Nikula and Väinölä, 2003; Sanna et al., 2013; Hui et al., 2016). Such population genetic investigation will enrich our knowledge about the larval dispersal, patterns of gene flow, as well as evolutionary history of T. maxima. It may also allow providing initial information for the conservation of this endangered giant clam species, known to be classified in the IUCN Red List of Threatened Animals (Nuryanto and Kochzius, 2009; Hui et al., 2016).

2. Materials and methods

2.1 Sample collection and genomic DNA extraction

A total of 44 specimens of the giant clam Tridacna maxima were collected from five sampling locations along the Saudi Arabian coast (Fig. 1, Table 1). During our sampling mission, we insisted to collect as many specimens as can be available in the field that can allow carrying out statistics. The adopted sampling scheme allowed covering the northern and central coasts of the Red Sea. From each specimen, 100 mg of mantel tissue were placed in a microfuge tube and frozen in liquid nitrogen. The tissues were ground using small plastic pestle, and DNA was extracted using PureLink genomic DNA mini kit (Invitrogen^TM, USA) as described by the manufacturer.

Figure  1.  Sampling locations of the giant clam Tridacna maxima across the Red Sea coast of Saudi Arabia. Patterns of distribution and proportion of the recorded COI haplogroups, along the examined locations, are exhibited in colored circles. The base map was constructed with the software DIVA-GIS 7.5.0 (http://www.diva-gis.org).

DownLoad: Full-Size Img PowerPoint

Table  1.  Sampling information on the giant clam Tridacna maxima including collection sites, geographic coordinates, and number of sampled specimens (N) per each location

Collection site	Geographic coordinates	N
Wajh	26°17′00″N, 36°25′00″E	10
Umluj	25°02′13″N, 37°17′16″E	6
Yanbu	24°04′56″N, 38°03′56″E	10
Rabigh	22°47′34″N, 39°02′01″E	6
Jeddah	22°17′00″N, 39°06′00″E	12

 | Show Table

DownLoad: CSV

2.2 Mitochondrial DNA amplification and sequencing

A fragment of the mitochondrial cytochrome oxidase subunit I (COI) gene was amplified using a pair of tridacnid-specific primers (COI-Tricro-Frwd: 5′-GGG TGATAATTCGAACAGAA-3′; COI-Tricro-Rev: 5′-TAGTTAAAGCCCCAGCTAAA-3′), previously specified by Kochzius and Nuryanto (2008). PCR reactions were performed in a total volume of 50 µL including 10 pg of DNA template, 1× PCR buffer, 2 mmol/L MgCl₂, 0.2 µmol/L of each primer, 0.2 mmol/L of each dNTPs and 1 unit of AmpliTaq Gold^TM Taq polymerase (Applied Biosystems, USA). PCR amplifications were carried out with an initial denaturation phase of 3 min at 94°C, followed by 35 cycles (each composed of 60 s of denaturation at 94°C, 60 s of annealing at 50°C and 60 s of extension at 72°C). These cycles were followed by 5 min of final extension at 72°C. Amplified PCR products were purified with PureLink PCR purification kit (Invitrogen^TM, USA) and sequenced at King Saud University (Saudi Arabia). The obtained sequences were visually inspected with Chromas Lite, aligned with Clustal W (implemented in BIOEDIT (Hall, 1999)), and then trimmed to a 685 bp fragment for subsequent analyses. All analyzed COI sequences were submitted to GenBank (accession numbers: LC322934-LC322977).

2.3 Statistical analyses

The nucleotide compositions, as well as the number of variable and parsimony-informative nucleotide sites, were assessed with MEGA version 7.0.18 (Kumar et al., 2016). Measurements of DNA polymorphism, including number of polymorphic sites (N_ps), number of haplotypes (N_h), haplotype diversity (h; Nei, 1987), nucleotide diversity (π; Tajima, 1983; Nei, 1987), and mean number of nucleotide differences (K), were calculated for each population as well as for the total dataset using DnaSP version 5.10 (Librado and Rozas, 2009). The number of private haplotypes was determined with ARLEQUIN version 3.1 (Excoffier et al., 2007). Endemism level for each population, as well as for the total dataset, was deduced by assessing the proportion of private haplotypes in relation to the total number of detected haplotypes.

In order to infer evolutionary relationships among the recorded COI haplotypes of T. maxima, a statistical parsimony network was constructed with the software TCS version 1.21 (Clement et al., 2000) under the 95% probability criterion for a parsimonious connection (Hudson, 1990; Templeton et al., 1992). Phylogenetic relationships among COI haplotypes were also inferred from a rooted phylogeny using the Bayesian approach as implemented in BEAST version 1.7.5 (Drummond et al., 2012). Unlike classical methods, this Bayesian inference of phylogeny is more informative about the evolutionary processes. Prior to the analysis, the most appropriate model of sequence evolution for the dataset was selected using the software MODELTEST version 3.7 (Posada and Crandall, 1998). Bayesian analyses were carried out with the HKY model (Hasegawa et al., 1985), found to be the best-fit substitution model for our data. The Markov Chain Monte Carlo (MCMC) simulations were run for 10 million steps and sampled every 1 000 steps. The corresponding outputs were reviewed in TRACER version 1.5 (Rambaut and Drummond, 2009) for robustness. Then, the resultant trees were summarized in TreeAnnotator version 1.7.5 (Drummond et al., 2012). The final results were checked and exhibited in FigTree version 1.4.0 (Rambaut, 2012). The Bayesian posterior probability (PP) was used to assess statistical support of the generated clades. Two published COI sequences of T. maxima (accession numbers: DQ155301 and KX713504), as well as those corresponding to the two sister species (T. crocea (accession number: DQ269479) and T. squamosa (accession number: KP205428)), were retrieved from GenBank and used as outgroup taxa. The particular choice of these COI sequences was based on their total nucleotide length (nearly equal to the base-pairs length generated for the new sequence data from the Red Sea), which allowed minimizing the loss of useful information in the phylogenetic reconstruction.

Population genetic differentiation was assessed by one-level AMOVA (analysis of molecular variance; Excoffier et al., 1992), as implemented in ARLEQUIN version 3.1 (Excoffier et al., 2007). Pairwise comparisons of genetic differentiation were also computed in the software ARLEQUIN. Both nucleotide divergence (Φ_ST: based on the Tamura-Nei model (Tamura and Nei, 1993)) and haplotypic frequency (F_ST) were used to infer these estimations. For both kinds of analyses, the resulting significant values were calculated from 10 000 permutations.

The existence of genetic subdivision within the analyzed dataset of T. maxima was also assessed by comparing levels of the two measurements of population differentiation, G_ST (based solely on haplotype frequencies) and N_ST (taking into account the genetic relationship among haplotypes). These parameters were estimated in the software PERMUT & CPSRR version 2.0 (Pons and Petit, 1996), following the methods described by Pons and Petit (1995, 1996). A significantly higher value of N_ST than that of G_ST usually hints at the presence of phylogeographic structure (Pons and Petit, 1996; Petit et al., 2005). The overall genetic diversity (H_T) and the average diversity within populations (H_S) were also estimated with the program PERMUT and compared.

Signatures of population demographic changes in the giant clam T. maxima were assessed by three neutrality tests (Tajima’s D (Tajima, 1989), Fu’s F_S (Fu, 1997), and Ramos-Onsins and Rozas’s R₂ (Ramos-Onsins and Rozas, 2002)). Both D and F_S indices were estimated in ARLEQUIN, while R₂ statistic was calculated in DnaSP. The examination of deviation from neutrality by all three indices was based on 1 000 coalescent simulations. A scenario of population expansion is likely supported by significantly negative D and F_S values as well as significant R₂ (in small population sizes). The Harpending’s raggedness index rg (Harpending, 1994) was also used to examine demographic changes in T. maxima according to the population expansion model implemented in ARLEQUIN. A total of 10 000 replicates allowed testing the significance of the rg index. The four above mentioned parameters (D, F_S, R₂, and rg) were applied to each examined population, to the overall dataset, as well as to the two identified COI haplogroups.

Other estimates of population size changes, including examination of site mismatch distributions of T. maxima and test of both demographic and spatial expansion models, were provided for the whole dataset as well as for the two detected haplogroups. Contrasting plots of observed and theoretical distributions of site differences yield insight into past population demographics. The expected mismatch distributions under a sudden expansion model were computed in ARLEQUIN. The sum of squared deviations (SSD) between observed and expected distributions were used as a measure of fit. The probability of obtaining a simulated SSD greater than or equal to the expected was computed by 1 000 random permutations. If this probability was greater than 0.05, the expansion model was accepted. Graphical representation was carried out by means of the growth-decline model implemented in DnaSP. Range expansion in T. maxima was investigated by the spatial expansion model, as implemented in ARLEQUIN (Excoffier, 2004).

The Bayesian Skyline plot (BSP) method (Drummond et al., 2005) was also applied to explore the magnitude of historical demographic events. Notably, this Bayesian coalescent approach allows the inference of population size history and also yields accurate estimation of expansion events (Grant, 2015). BSP plot was generated for the whole COI dataset. Analyses were carried out in BEAST version 1.7.5 (Drummond et al., 2012) considering a HKY model of nucleotide substitution (as already calculated by MODELTEST) and a strict molecular clock. The COI gene mutation rate of 2% per million years (Layton et al., 2014), generally adopted and used for marine mollusks, was implemented in the analysis in order to date expansion event. Two independent MCMC (each with 40 million iterations) were carried out. Following the removal of the first 10% iterations (4 millions) as burn-in, the remaining replicates were combined in LogCombiner version 1.7.5 (Drummond et al., 2012) and summarized as BSP after checking their convergence (effective sample sizes (ESS) of all parameters >200) in TRACER version 1.5 (Rambaut and Drummond, 2009).

3. Results

3.1 Genetic diversity

After being proofread and trimmed to an alignment length of 685 base pairs, molecular analysis of 44 sequences of Tridacna maxima (from five locations across the Red Sea) unveiled the existence of 23 variable nucleotide sites. Out of these, five were parsimony-informative. More than 40% of the examined sequences were different and allowed the characterization of 20 haplotypes (Fig. 2, Table 2). The nucleotide composition of the analyzed fragment showed an A-T bias (C=16.80%; T=37.50%; A=24.80%; G=20.90%), which is typical for invertebrate mitochondrial DNA (Simon et al., 1994).

Figure  2.  Evolutionary relationships among the recorded haplotypes of Tridacna maxima inferred from the Bayesian inference of phylogeney (a), based on the alignment of 657 bp of Cox1 gene, and from TCS (the method of Templeton et al. (1992)) parsimony network (b), based on the alignment of 685 bp of COI gene. For the Bayesian phylogenetic tree, the values of the posterior probabilities (Bayesian inference) are highlighted for the main generated nodes. For the TCS haplotype network, small black circles correspond to missing (or hypothetical) haplotypes. Each line between two points represents one mutational step. Circle sizes depict proportions of haplotypes; the smallest corresponds to 1 and the largest to 14 individuals.

DownLoad: Full-Size Img PowerPoint

Table  2.  Genetic diversity (number of polymorphic sites (N_ps), number of haplotype (N_h), number of private haplotype (N_ph), endemism (N_ph/N_h), haplotype (h) and nucleotide (π) diversities, and mean number of nucleotide differences (K)), Neutrality tests (Tajima’s D, Fu’s F_S, and Ramos-Onsins and Rozas’s R₂ tests), and mismatch distribution raggedness index (rg) for each population of Tridacna maxima as well as for the entire sample

Population	N	Nps	Nh	Nph	Endemism (Nph/Nh)	h	π	K	D	FS	R2	rg
Wajh	10	7	6	3	0.500	0.844±0.103	0.003 8±0.000 5	2.600	0.215	−1.164	0.157	0.108
Umluj	6	5	3	1	0.333	0.733±0.155	0.003 9±0.000 7	2.667	1.219	1.574	0.245	0.827
Yanbu	10	9	5	3	0.600	0.756±0.130	0.004 2±0.001 0	2.867	−0.433	0.218	0.137	0.140
Rabigh	6	11	6	5	0.833	1.000±0.096	0.006 7±0.001 1	4.600	−0.273	−2.328	0.121	0.098
Jeddah	12	8	8	5	0.625	0.848±0.104	0.002 8±0.000 6	1.924	−1.093	−4.167	0.092	0.067
Total data	44	23	20	17	0.850	0.845±0.041	0.004 2±0.000 3	2.846	−1.527	−11.647	0.052	0.046
Note: N indicates sample size for each population. Significant values are in bold. Non significant values for the raggedness index accept the null hypothesis of expectation under a sudden demographic expansion model. For haplotype and nucleotide diversities, each value is the mean±standard deviation.

 | Show Table

DownLoad: CSV

Genetic diversity analysis of the total examined COI dataset showed high haplotype diversity (h=0.845±0.041) and low nucleotide diversity (π=0.004 2±0.000 3). The total mean number of nucleotide differences was 2.846. Details on these parameters of genetic variability, for each examined population, are reported in Table 2. Endemism levels, highlighting the proportion of private haplotypes in relation to the total defined halpotypes, ranged from weak values in the populations of Wajh and Umluj to moderate values in those of Yanbu and Jeddah (Table 2). The population of Rabigh was shown to be the most endemic in terms of the COI haplotype composition.

3.2 Phylogenetic relationships among COI haplotypes

The outcome of the Bayesian phylogenetic analysis showed marked genetic divergence between T. maxima and the two closely related species T. crocea and T. squamosa (Fig. 2a). Within T. maxima, a clear separation between two well statistically supported clades was highlighted (posterior probability [PP] for each generated clade is 1) (Fig. 2a). This finding is in accordance with previous ones unveiling genetic distinctiveness of the Red Sea (Nuryanto and Kochzius, 2009; Hui et al., 2016). Furthermore, within the new examined COI dataset from the Red Sea, a clear separation between two sub-clades (or haplogroups) can be noticed (Fig. 2a). The first sub-clade (Haploroup 1) was found to be statistically supported (PP=0.89); while the second sub-clade (Haplogroup 2) was poorly resolved (PP=0.31) (Fig. 2a). This pattern of mitochondrial genealogy was clearly supported by the outcome of the TCS statistical parsimony procedure, exhibiting a double star-like shape (Fig. 2b). For both generated sub-networks (corresponding to Haplogroups 1 and 2), all haplotypes were closely related to the two most common Haplotypes H1 and H2 (Fig. 2b). The most frequent Haplotype H2 (recorded in 14 specimens; Table 3) was found in all populations except Rabigh. The next most common Haplotype H1 (present in 11 individuals; Table 3) prevailed across the entire surveyed region. These common two haplotypes (characterizing the two delineated mitochondrial sub-clades) were separated by four mutational steps (Fig. 2b). Notably, both identified haplogroups were found to be co-distributed across the five sampled locations (Figs 1 and 2b).

Table  3.  Distribution pattern of the recorded haplotypes in the examined populations of Tridacna maxima

Haplotype	Population
	Wajh	Umluj	Yanbu	Rabigh	Jeddah
H1	4	3	2	1	1
H2	2	2	5	0	5
H3	0	1	0	0	0
H4	1	0	0	0	0
H5	1	0	0	0	0
H6	1	0	0	0	0
H7	1	0	0	0	1
H8	0	0	1	0	0
H9	0	0	1	0	0
H10	0	0	1	0	0
H11	0	0	0	1	0
H12	0	0	0	1	0
H13	0	0	0	1	0
H14	0	0	0	1	0
H15	0	0	0	1	0
H16	0	0	0	0	1
H17	0	0	0	0	1
H18	0	0	0	0	1
H19	0	0	0	0	1
H20	0	0	0	0	1
Total	10	6	10	6	12

 | Show Table

DownLoad: CSV

3.3 Population genetic structure

Analysis of molecular variance (one-level AMOVA) revealed lack of genetic differentiation among examined populations of T. maxima (Φ_ST=0.047, df=43, P=0.166, based on Tamura-Nei distance; F_ST=0.020, df=43, P=0.232, based on haplotype frequency) (Table 4). These results were supported by the outcome of pairwise comparisons of genetic differentiation (Table 5). Indeed, all pairwise comparisons, estimated from nucleotide divergence, were not significant, except the one noted between the populations of Wajh and Jeddah (Φ_ST=0.196, P=0.033). Conventional F-statistics, computed from haplotype frequencies, also revealed high genetic connectivity among the examined populations, highlighting non-significant F_ST values (Table 5).

Table  4.  One-level analysis of molecular variance (AMOVA) testing for genetic differentiation among the examined locations, and PERMUT analysis testing for the existence of population subdivision and phylogeographic structure within the examined dataset

Partition tested	Genetic differentiation			Phylogeographic structure
	Tamura and Nei distance	Haplotype frequency		Ordered haplotypes	Unordered haplotypes
Wajh vs. Umluj vs. Yanbu vs. Rabigh vs. Jeddah	ΦST = 0.047 ns	FST = 0.020 ns		NST = 0.017	GST = 0.023
Note: ns means non-significant difference (P>0.05).

 | Show Table

DownLoad: CSV

Table  5.  Pairwise comparisons of genetic differentiation, for Tridacna maxima, estimated from nucleotide divergence (Φ_ST, below the diagonal) and haplotype frequency (F_ST, above the diagonal)

	Wajh	Umluj	Yanbu	Rabigh	Jeddah
Wajh	*	−0.079	0.024	0.018	0.032
Umluj	−0.112	*	−0.016	0.054	0.028
Yanbu	0.076	−0.015	*	0.103	−0.035
Rabigh	−0.063	−0.109	0.034	*	0.071
Jeddah	0.196	0.119	−0.031	0.158	*
Note: Significant values (P<0.05) are reported in bold. Significance level was calculated from 10 000 permutations.

 | Show Table

DownLoad: CSV

Lack of population subdivision and phylogeographic structure, within the studied giant clam species across the surveyed region, was also inferred from the outcome of PERMUT analysis. The overall genetic diversity (H_T=0.856) of T. maxima was not considerably higher than the average diversity within populations (H_S=0.836). Furthermore, the N_ST value was remarkably lower than the G_ST value (Table 4), highlighting lack of a significant relationship between phylogeny and the geographical distribution of haplotypes.

3.4 Demographic history

The applied neutrality tests revealed significant deviations from mutation-drift equilibrium for the two populations of Rabigh and Jeddah (Table 2), suggesting that these populations seem to have experienced a recent expansion event. Evidence of departure from neutrality was also recorded for both identified haplogroups as well as for the whole dataset. Indeed, all analyzed neutrality tests resulted in significant values (with marked negative output for both Tajima’s D and Fu’s F_S), associated with small and non-significant values of the Harpending’s raggedness index, for Haplogroup 1 (D=–2.160, P = 0.003; Fu’s F_S=–9.570, P=0.000; R₂=0.050, P=0.000; rg=0.046, P=0.790), Haplogroup 2 (D=–1.876, P=0.013; Fu’s F_S=–3.893, P=0.005; R₂ = 0.073, P=0.000; rg=0.039, P=1.000) and total data (D=–1.527, P=0.040; Fu’s F_S=–11.647, P=0.000; R₂=0.052, P=0.009; rg=0.046, P=0.408).

Mismatch distribution analyses showed unimodal variation for Haplogroup 1 as well as for the whole dataset (Fig. 3), highlighting a scenario of a sudden demographic expansion. This trend was not observed for Haplogroup 2 (Fig. 3). Statistical analysis of mismatch distributions allowed accepting the model of demographic expansion for both Haplogroup 1 and total data (exhibiting non-significant values for the two demographic indices SSD and rg), but not for Haplogroup 2 (Fig. 3, Table 6). Nevertheless, the spatial expansion model was accepted for all three analyzed genetic entities (Table 6), suggesting that both haplogroups, as well as the whole dataset of T. maxima, have recently expanded their distribution range.

Figure  3.  Mismatch distribution for the two identified haplogroups as well as for the whole dataset of Tridacna maxima. Observed frequencies (dotted line) were compared to the expected frequencies (continuous line), estimated by the growth-decline model implemented in DnaSP. The demographic expansion parameters used, τ and θ_initial, were calculated in ARLEQUIN; and θ_final value was fixed at 1 000 to simulate the infinite. The two demographic indices SSD and rg, provided for the observed mismatch distributions, were calculated under the assumption of a demographic expansion model, implemented in ARLEQUIN.

DownLoad: Full-Size Img PowerPoint

Table  6.  Test of both demographic and spatial expansion models for the two identified haplogroups as well as for the whole dataset of Tridacna maxima

Examined group/tested expansion model	N	Demographic expansion			Spatial expansion
		SSD	P		SSD	P
Haplogroup 1	26	0.004	0.582		0.003	0.594
Haplogroup 2	18	0.291	0.001		0.004	0.798
Whole dataset	44	0.018	0.292		0.019	0.384
Note: N represents sample size, and SSD sum of squared deviations between observed and expected distributions under the tested expansion model. The probability of obtaining a simulated SSD greater than or equal to the expected was computed by 1 000 random permutations. If this probability (P) was greater than 0.05, the expansion model is accepted.

 | Show Table

DownLoad: CSV

Demographic history of the metapopulation of T. maxima, from the Red Sea, was also inferred and detailed from the coalescent approach of the BSP analysis. Pattern of effective population size evolution over time was shown to suddenly and recently increase following a relatively stationary period (Fig. 4). Overall, BSP results aligned with those inferred from neutrality tests and mismatch distribution analyses. Specifically, the outcome of the analysis revealed that the sudden expansion event started roughly at about 36 500 years ago (CI: 22 000–48 000 years ago). Since approximately 20 000 years ago, the effective population size has remained constant until present.

Figure  4.  Bayesian skyline plot (BSP) for the whole dataset of Tridacna maxima, showing changes in effective population size (Ne multiplied per generation time) over time (measured in years before present). The thick solid line depicts the median estimate, and the margins of the blue area represent the highest 95% posterior density intervals.

DownLoad: Full-Size Img PowerPoint

4. Discussion

The results of the present study revealed high genetic connectivity among populations of the giant clam T. maxima across part of the Red Sea off the Saudi Arabian coast. Pattern of genetic homogeneity was also recorded within other marine species across the same surveyed geographic spectrum, such as in the reef-building coral Pocillopora verrucosa (Robitzch et al., 2015). However, our finding clearly contrasts with those of earlier population genetic studies of T. maxima (using nuclear and mitochondrial DNA markers) unveiling significant genetic differentiation and marked patterns of restricted gene flow across other parts of the distribution zone (Laurent et al., 2002; Juinio-Meñez et al., 2003; DeBoer et al., 2008; Kochzius and Nuryanto, 2008; Nuryanto and Kochzius, 2009; Ahmed Mohamed et al., 2016; Hui et al., 2016).

The obtained results may provide two major insights on the impact of biotic (species biological characteristics) and abiotic factors (impact of environmental factors of the surveyed geographic spectrum) on the observed pattern of genetic homogeneity within the analyzed dataset from the Red Sea. First, the short planktonic larval dispersal period, characterizing the examined giant clams (Lucas, 1988), does not seem to have a significant impact on population genetic structuring. Our finding consolidates the notion stating that population genetic structure does not always reflect a species’ dispersal capabilities. Indeed, an increasing number of investigations have unveiled lack of direct relationship between planktonic larval duration and discerned levels of genetic divergence among populations (Weersing and Toonen, 2009; Hunter and Halanych, 2010; Pannacciulli et al., 2017; Deli et al., 2017, 2018). Second, the prevailing strong environmental and physico-chemical gradients in the Red Sea were not shown to potentially induce genetic divergence among populations of T. maxima. Indeed, the water exchange occurring between the Indian Ocean and the Red Sea through the Strait of Bab el Mandeb allows the entrance of cooler and less saline oceanic water to the Red Sea. This water current is subsequently warmed up in the shallow southern Red Sea and exhibits northward gradual decrease in temperature, but increasing levels of salinity due to high evaporation rate (Sofianos et al., 2002). Despite this induced natural environmental heterogeneity within the Red Sea, expected gradual shifts in temperature and salinity across the surveyed geographic spectrum do not seem to have significant impact on population genetic structure within T. maxima, contrary to what has been shown recently in other marine species such as in the coral reef fish Amphiprion bicinctus (Nanninga et al., 2014), and the reef sponge Stylissa carteri (Giles et al., 2015). Noteworthy, the latitudinal range covered by the latter investigations was relatively more complete than the one in the present study. Accordingly, we assume that the observed genetic variability within examined populations of the giant clam does not seem to be achieved by the potential impact of positive selection on mitochondrial DNA haplotypes, as has been already unveiled for other marine invertebrate species such as the walleye pollock Theragra chalcogramma in the North Pacific (Grant et al., 2006) and the purple sea urchin Paracentrotus lividus (Penant et al., 2013) across the Eastern Mediterranean (Levantine Basin). However, with the still persisting lack of samples from the southern Red Sea, the potential impact of natural selection on forging genetic diversity of the giant clam needs to be investigated. Furthermore, analysis of other molecular markers, prone to reflecting the impact of positive selection on genetic variability, is required in order to assess the finding of the present investigation. Variable nuclear microsatellite markers, which have been already isolated and characterized in T. maxima (Grulois et al., 2015), are planned to be analyzed in Red Sea populations of this species in order to obtain a more complete view about the basis of their differentiation.

In light of these considerations, several explanations can be advanced to explain the origin of high genetic connectivity among the examined Red Sea populations of the giant clam. For instance, the lack of mtDNA genetic structure in T. maxima could be linked to incomplete lineage sorting owing to a possible slow COI mutation rate in this species. Noteworthy, despite being highly variable in many marine invertebrate species including bivalves (Nikula and Väinölä, 2003; Reuschel et al., 2010; Sanna et al., 2013; Hui et al., 2016; Deli et al., 2018), the COI gene was shown to exhibit limited genetic diversity in other marine invertebrates. For example, genetic variation level at this mitochondrial marker was shown to be so low that it was unable to discern marked genetic divergence patterns in the poorly dispersive red sea star Echinaster sepositus (Garcia-Cisneros et al., 2016). Nevertheless, recent large-scale investigations have revealed the existence of highly divergent lineages of COI gene in T. maxima (DeBoer et al., 2008, 2014a, b; Kochzius and Nuryanto, 2008; Nuryanto and Kochzius, 2009; Hui et al., 2016). Hence, the above-mentioned hypothesis cannot explain the observed mtDNA homogeneity in the examined giant clam species.

Given that several isolating contemporary (i.e., discontinuous habitat and/or heterogeneous oceanographic features) and historical (sea level regression during glacial maxima of the Pleistocene that have contributed to genetic divergence and isolation in marine fauna) processes have been proposed to explain the strong genetic divergence among COI lineage of T. maxima across macro-geographic scales within the distribution area (Nuryanto and Kochzius, 2009; Hui et al., 2016), we rather propose that the COI gene may have little resolution in yielding genetic divergence within the species only at local scales (lacking known physical and/or historical barriers to gene flow). Accordingly, the observed mitochondrial genetic homogeneity within T. maxima across part of the Red Sea may not be caused by incomplete lineage sorting nor by the limited number of samples, but rather stem from the lack of isolating processes (that can significantly affect spatial distribution of genetic diversity and/or restrict gene flow) across the limited surveyed geographic spectrum. Focusing on more populations of T. maxima (especially from the southern Red Sea coast), and contrasting the newly obtained results with those inferred from the analysis of variable nuclear markers, will shed more light on pattern and origin of genetic structuring in this species.

Knowing that the COI gene is variable enough to allow detecting patterns of genetic divergence within T. maxima, we may think of the potential effect of the particular oceanographic features of the Red Sea in maintaining high levels of gene flow within this species along the studied sampling area. Indeed, the noticeable presence of the two common haplotypes (H1 and H2; Fig. 2b) in all examined populations, as well as the marked co-distribution of both delineated haplogroups across the five surveyed locations (Fig. 1), indicate that gene flow is likely involving in maintaining genetic connectivity among the sampling sites. Noteworthy, given that the short and limited larval dispersal potential within T. maxima may not be efficient enough to trigger widespread gene flow, this explanation is still questionable. However, Hui et al. (2016) claimed that in the presence of a short pelagic larval duration, giant clams can potentially disperse over long distances and exhibit low genetic differentiation. Furthermore, Slatkin and Barton (1989) stated that migration of few individuals per generation is sufficient to oppose the impact of restricted gene flow (possibly stemming from short larval dispersal potential) and hence prevent genetic differentiation. In this context, marine currents can act as powerful transport routes for marine larvae (Pineda et al., 2007) and play a crucial role in assuring gene flow in marine environments. Despite the scarcity of information on current flow patterns in the Red Sea, the main driver of surface currents in this semi-enclosed Basin is presumably a dominant North to South wind (Sofianos et al., 2002; Sofianos and Johns, 2007; Raitsos et al., 2013). Furthermore, permanent anticyclonic and cyclonic eddies can be found in the Red Sea (Sofianos and Johns, 2007; Raitsos et al., 2013). It has been shown that even if larvae dispersal is restricted to the range of the eddy, cyclonic water currents might significantly assure larval transport from one shore to the other owing to the considerably narrow width of the Red Sea (Lobel and Robinson, 1986; Sammarco and Andrews, 1989; Wolanski et al., 1989). Hence, despite the limited and short planktonic larval potential characterizing the examined giant clam species, the combined effect of the wind patterns and permanent eddies in the Red Sea can potentially promote larval dispersal, maintaining gene flow over long distances, and thereby lead to panmixia within T. maxima across the surveyed geographic spectrum.

Alternatively, the noticeable pattern of population genetic structure within T. maxima may likely reflect the impact of historical processes, rather than contemporary ones. Indeed, the observed high genetic diversity along with low nucleotide diversity, altogether with the shallow haplotype genealogy can be considered as residual effect of a recent evolutionary history of the species in the Red Sea. This assumption is likely supported by the significant deviation from neutrality (as inferred from the outcome of the three applied neutrality tests) as well as by the recent demographic expansion event (retrieved from both mismatch distribution analysis and BSP plot). Nevertheless, pattern of COI genealogy may not be strongly consistent with a recent dramatic population expansion from a small population size, as at least two haplotypes (H1 and H2) were shown to be relatively common and geographically widespread. A possible explanation to this finding is that the current geographical distribution of extant haplotypes of T. maxima across the Saudi Arabian coast may have been shaped by the impact of recent range expansion. Hence, the lack of population genetic structure observed across the surveyed distribution range may not necessarily reflect high levels of contemporary gene flow. Instead, it may potentially originate from the impact of recent spatial and demographic expansions, as clearly recorded in our study. The spatial expansion model accepted for the two haplogroups, defining the two common haplotpes (H1 and H2), as well as their remarkable sympatric occurrence in all surveyed locations, strongly favor this possibility, and suggest potential admixture between both genetic entities following their recent range expansion. This scenario has been recently proposed to explain similar pattern of spatial distribution of genetic diversity in other marine invertebrate species such as in the black sea urchin Arbacia lixula across the central African Mediterranean coast (Deli et al., 2019).

While all of the above-mentioned scenarios may provide possible explanations to the high genetic homogeneity within T. maxima, they need to be considered with caution and regarded as tentative owing to the limited samples examined in this study as well as to the restricted sampling area (focusing only on the Saudi Arabian coast). Indeed, increase in the sample size of the giant clam, along with extended sampling covering the entire Red Sea, may unveil new insights into the phylogeographic structure of the species. Earlier population genetic investigation showed that exhaustive increase in sample size allowed unraveling significant genetic structuring (with mtDNA) in marine invertebrate species, previously regarded as panmictic due to geographically and numerically limited datasets (Fratini et al., 2016). Furthermore, the observed genetic pattern in our study, based on the analysis of only one mtDNA marker, may not be discerned when analyzing independent genetic markers (such as nuclear markers) with different evolutionary rates and coalescent histories. Hence, the use of several genetic markers in conjunction with well representative sampling will provide an accurate view on pattern and origin of population genetic diversity and structure of T. maxima. In this context, analysis of more powerful nuclear markers (such as microsatellite loci) is required and planned to be carried out in order to confirm the retrieved genetic signal. Furthermore, analysis of microsatellite markers in more populations from the Red Sea (specifically from the southern part) will shed more light on the impact of environmental heterogeneity on genetic variability of this poorly dispersive bivalve species.

At last, regardless the mechanisms responsible for shaping genetic diversity and structure of examined populations of T. maxima, the preliminary genetic signal, obtained in this study, is of practical importance for understanding the evolutionary history as well as for the build up to the management and conservation of this endangered species in the Red Sea.

Acknowledgements

We thank the University of Taibah (Saudi Arabia), and the University of King Saud (Saudi Arabia) for supporting the current investigation. Special thanks are extended to the latter for help with mitochondrial DNA sequencing. Molecular analysis in this work was funded by the Research Supporting Project (RSP-2019/75), King Saud University, Riyadh, Saudi Arabia. The corresponding author (Temim Deli) carried out the statistical analyses and wrote the manuscript. We also thank two anonymous reviewers for their helpful and constructive comments that greatly contributed to improving the manuscript quality.