Migration-Driven Phylogeographic Homogeneity in the Red Macroalga Halymenia durvillei (Halymeniales, Rhodophyta) Along the Andaman Coast of Thailand

1. Introduction

Southeast Asia, located at the crossroads of the Pacific and Indian Oceans, is characterized by numerous islands and complex geographical features. Despite covering only about 4% of the world’s land area, this region has been recognized as one of the most important marine biodiversity hotspot globally [1]. Central to this biodiversity richness is the well-known “Coral Triangle”, often referred as the “Amazon of the ocean”, which hosts about 76% of the world’s reef-building corals and 37% of coral reef fish species [2]. The combination of biodiversity-rich ecosystems, dynamic climatological processes, and a unique geological history has made Southeast Asia an ideal model region for marine phylogeographic studies over the past decades [3,4,5,6,7,8].

The Thai-Malay Peninsula (TMP) is one of the key regions for biodiversity monitoring and management in Southeast Asia [9]. Geographically, this peninsula extends from north to south, with the northeastern coast bordering the Gulf of Thailand and the southern part connected to the southern South China Sea. On the western side, the northern section is adjacent to the Andaman Sea while the southern section links with the Strait of Malacca (Figure 1). The hydrological and ecological environments of the Pacific and Indian Oceans converge along the peninsula, forming a typical land–sea interaction landscape that supports diverse habitats (e.g. tropical coral reefs, seaweed beds, and mangroves) and has attracted growing research interest [3,8,10]. However, marine phylogeographic studies in this region have primarily focused on reef-building corals and coral reef fishes (e.g. [11,12,13]), mangroves (e.g. [5,6]) and arthropods (e.g. [7]). In contrast, limited phylogeographic efforts have been devoted to macroalgae across the TMP. Available studies include the brown alga Padina boryana Thivy [14], the green alga Halimeda J.V.Lamouroux [15,16], and the red algae Bostrychia tenella (J.V.Lamouroux) J.Agardh/B. binderi Harvey 1849 [17,18] and Gracilaria salicornia (C.Agardh) E.Y.Dawson [19]. Notably, these studies have reported dissimilar and even contrasting genetic structuring patterns, underscoring the need of additional phylogeographic evidence from co-occurring macroalgae. Such evidence would help fill knowledge gaps in understanding the patterns and processes underlying speciation and biogeography, while also informing conservation and management prioritization of foundation species in benthic communities.

Figure 1. Sampling distribution and haplotype networks of Halymenia durvillei along the Andaman coast of Thailand, constructed based on mitochondrial cox1 (a) and chloroplast rbcL (b) sequences, respectively. The circle size is proportional to the occurrence frequency of haplotype, and numbers along the branches indicate the number of mutational steps. “mv” denotes a hypothetical unsampled or ancestral haplotype. Sampling sites are divided into two marine ecoregions, the Myanmar Shelf and the Strait of Malacca Shelf, and all site codes are the same as in Table 1.

Figure 1. Sampling distribution and haplotype networks of Halymenia durvillei along the Andaman coast of Thailand, constructed based on mitochondrial cox1 (a) and chloroplast rbcL (b) sequences, respectively. The circle size is proportional to the occurrence frequency of haplotype, and numbers along the branches indicate the number of mutational steps. “mv” denotes a hypothetical unsampled or ancestral haplotype. Sampling sites are divided into two marine ecoregions, the Myanmar Shelf and the Strait of Malacca Shelf, and all site codes are the same as in Table 1.

It is noteworthy that the Andaman Sea, located along the west coast of Thailand, is a marginal sea at the boundary between the Indian and Pacific Oceans, characterized by a unique ecological background and geological history. The basin was formed during the late Eocene to Oligocene as a result of crustal extension related to the collision between the Indian and Eurasian plates [20]. During the late Miocene (c. 10–11 Mya), the trench-arc-basin system began to form and gradually stabilized throughout the Pliocene and Pleistocene [21]. Continued tectonic activity during the Quaternary further shaped the regional landscape and sedimentary environment [22]. During the transition from the late Pleistocene to the Holocene, dramatic sea level fluctuations strongly influenced land connections, species diversity and biogeographic patterns [1]. These dynamic tectonic and eustatic processes divided the Andaman Sea into the northern Myanmar Shelf and the southern Malacca Strait Shelf Provinces ([23], Figure 1). Geologically and geographically, the Myanmar Shelf features a wide and fragmented coastline dominated by muddy sand and mangrove/seagrass habitats. In contrast, the Malacca Strait has a narrow and steep continental shelf, with coarse sand, gravel coral reefs and rocky habitats [24]. Such long-term geomorphological, hydrodynamic, and ecological differences between these two marine ecoregions are likely to have shaped the phylogeographic patterns of benthic organisms, including sessile macroalgae.

The red algal genus Halymenia C. Agardh belongs to the family Halymeniaceae and contains at least 85 valid species recognized worldwide [25], with a typical triphasic life history (i.e. gametophyte, carposporophyte, and tetrasporophyte). Halymenia is particularly common in tropical and subtropical marine waters [26], and several species within the genus (e.g. H. durvillei Bory 1828) possess significant economic value due to extractable carrageenan and phycoerythrin [27,28]. Traditional taxonomy of Halymenia relied primarily on morphological features of vegetative organs (Tan et al. 2018). However, these characteristics can be strongly influenced by environmental conditions, leading to morphological convergence among different species and hence taxonomic confusion [29]. In recent years, integrating molecular and morphological evidence has shown particular promise for species identification and phylogenetic reconstruction in Halymenia [29,30]. Halymenia durvillei is a representative species of the genus in the Indo-Pacific, with distribution range extending from Southeast Asia to the eastern coast of Africa, Japan, Australia, and some Pacific islands [25,29]. This species typically grows on coral reefs or rocky substrates and can withstand strong wave action, showing both ecological significance and economic potential as a benthic foundation species.

In this study, we explored the phylogeographic patterns of H. durvillei along the Andaman coast of Thailand using mitochondrial cytochrome c oxidase subunit I (mtDNA cox1) and chloroplast ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit (cpDNA rbcL). Specifically, we aimed to address the following questions: (i) Does H. durvillei exhibit detectable intraspecific lineage differentiation and population genetic structuring between the northern and southern Andaman Sea? (ii) Have ocean currents controlled by the Indochinese monsoon system [23] and sea-level fluctuations during the late Quaternary ice ages contributed to the phylogeographic history of H. durvillei? Answers to these questions not only advance our understanding of the origin and diversification of seaweed in Southeast Asia, but also provide practical guidelines for conservation decision-making of benthic seaweeds under climate change.

2. Materials and Methods

2.1. Collection, DNA Extraction, Amplification, and Sequencing

From 2018 to 2024, H. durvillei samples were collected from 20 sites along the coast of the Andaman Sea in Thailand (Table 1, Figure 1). For the site PHA1 (Lam Kean beach, Phangnga), samples were collected in five years (2019, 2021, 2022, 2023, 2024). For another site PHU25 (Tankhen Bay, Phuket), samples were also collected in six different periods (2018–2024, Table 1). Most populations were located from the eastern and western sides of Phuket Island, corresponding to the Myanmar and the Malacca Strait shelfs (Figure 1). At each site, 3–5 cm of fresh thalli free of epiphytes were cut from each individual. After collection, all samples were cleaned, dried, and preserved using silica gel. Total genomic DNA was extracted using the FastPure Plant DNA Extraction Mini Kit (Vazyme Biotech Co., Ltd., Nanjing, China) following the manufacturer’s instructions.

MtDNA cox1 and cpDNA rbcL have proven effective for distinguishing species and revealing phylogenetic relationships in Halymenia [29,30,31]. Accordingly, these two markers were selected to analyze phylogeographic diversity and population genetic structure of H. durvillei. Two primer sets, cox1-F (5′-GGAACACTTTAYTTAATTTTTGG-3′) and cox1-R (5′-TGRTATARAATTGGATCWCC-3′) [32], rbcL1-F (5′-AGAATGAACGCTATGAG-3′) and rbcL1-R (5′-GCTGTTTGTAGAGGACC-3′) (designed in this study), were used to amplify cox1 and rbcL, respectively. Polymerase chain reaction (PCR) amplification was carried out in a 50 μL reaction volume following Nguyen et al. (2023). PCR programs for amplifying cox1 and rbcL each consisted of an initial denaturation at 94 °C for 1 min, annealing at 50 °C for 1.5 min and extension at 72 °C for 1.5 min, followed by 35 cycles of denaturation at 94 °C for 1 min, annealing at 50 °C for 1.5 min, and extension at 72 °C for 1 min. A final extension was performed at 72 °C for 5 min. PCR products were checked using 1% agarose gel electrophoresis. After purification, amplification products were sent to Sangon Biotech (Shanghai, China) for bidirectional sequencing using the same primers described for the amplification.

2.2. Genetic Diversity and Population Structure

The obtained cox1 and rbcL sequences were manually trimmed using BioEdit 7.1.9.0 [33] and edited and aligned using MEGA 11 [34]. MtDNA cox1 and cpDNA rbcL haplotypes were identified using DnaSP 5.10 [35]. Genetic diversity indices, including the number of haplotypes (N_h), haplotype diversity (h), and nucleotide diversity (π), were calculated using Arlequin 3.5 [36]. Haplotype networks for cox1 and rbcL were constructed using Network 10.2.0.0 [37] to illustrate their evolutionary relationships.

Analysis of molecular variance (AMOVA) was conducted to partition the levels of genetic variation among groups, among populations within groups, and within populations. To minimize the bias caused by uneven sampling sizes, the eight populations with only one individual (PHU22, PHU23, PHU37, PHU14, PHU54, PHU19, PHU49, and PHU9) were excluded from the analyses. Two grouping criteria were applied for AMOVA partitioning. Criterion 1: sampling sites were divided into five groups according to geographic proximity – Group 1 (RNG1a and RNG2), Group 2 (PHA1, PHA4 and PHA8), Group 3 (TRA2 and TRA5), Group 4 (PHU31, PHU38, PHU50 and PHU51), and Group 5 (PHU25). Criterion 2: sampling sites were divided into two groups according to the marine ecoregion classification of the Myanmar Shelf and the Malacca Strait Shelf Provinces (Figure 1) – Group 1 (RNG1a, RNG2, PHA1, PHA4, PHA8, PHU51, PHU50 and PHU31), and Group 2 (TRA5, TRA2, PHU38 and PHU25).

The program Structure v2.3.4 [38] was used to infer potential genetic clustering of H. durvillei populations. To avoid potential bias caused by uneven population sizes among sampling sites, the eight populations with only one individual (PHU22, PHU23, PHU37, PHU14, PHU54, PHU19, PHU49, and PHU9) were removed from the clustering analysis. Genetic clustering was performed based on the identified cox1 and rbcL haplotypes. The number of clusters (K) was set from 1 to 5, with five independent runs for each K value. The Markov chain Monte Carlo (MCMC) simulation was run for 5×10⁶ iterations, with the first 5×10⁵ iterations discarded as burn-in. LnP(K) and ΔK [39] were computed using the online tool StructureSelector (https://lmme.ac.cn/StructureSelector/) to identify the optimal K value and generate clustering diagrams.

2.3. Gene Flow and Population Demographic History

Since mtDNA cox1 contained relatively richer phylogeographic information compared with cpDNA rbcL (see Results below), the cox1 dataset was used to estimate gene flow among populations using Migrate-n v5.0.4 [40]. To minimize the bias caused by uneven sample sizes, the eight populations with only one individual (PHU22, PHU23, PHU37, PHU14, PHU54, PHU19, PHU49, and PHU9) were removed from the analysis. The relatively stable genetic composition detected among a few consecutive years (see Table 1) allowed us to choose PHA1-2022 and PHU25-2023-07 that harboured the maximum number of cox1 haplotypes (see Table 1) as population representatives for estimation of migrants. We also specifically estimated potential ocean-driven migrants among adjacent sampling sites during the northeast monsoon season (November-April). The effective population size (θ = xN_eμ, where N_e represents the effective population size, μ is the mutation rate per generation, and x equals 1 for organellar data) and the migration rate (M = m/μ, where m denotes the migration rate per generation) were substituted into the formula (Nm =θM/x) to calculate the effective number of migrants. The MCMC was run following the method described by [8].

Bayesian skyline plot (BSP) analysis was performed using BEAST v1.10.4 [41] to evaluate historical changes in effective population size of H. durvillei. The best-fit evolutionary model (GTR+F+G4) was selected using ModelFinder implemented in PhyloSuite v1.2.3 [42] based on the Akaike Information Criterion (AIC). The divergence time of Halymenia within the family Halymeniaceae (c. 170 million years ago (Mya)) was inferred from phylogenetic dating of red algal lineages [43] and then used as a calibration reference. A strict molecular clock was applied and the tree prior was set to the coalescent Bayesian skyline. The MCMC analysis was run for a total of 1×10⁷ generations, with the first 1×10⁶ generations discarded as burn-in. Sampling was performed every 1000 generations. The convergence was assessed in Tracer v1.7.2 [44]. When the effective sample size (ESS) values were greater than 200, the results were considered to have reached convergence.

3. Results

3.1. Population Genetic Diversity

A total of 333 mtDNA cox1 sequences were obtained from 20 sampling sites. After alignment and editing, the final cox1 dataset comprised 580 bp with 90 variable sites, yielding 11 haplotypes (C1–C11, GenBank accession numbers PZ816517–PZ816527) (Table 1). Haplotype C1 was the most abundant and geographically widespread across the Andaman coast of Thailand (Figure 1a). Six haplotypes (C6–C11) were private, with C6-C8 restricted to Phuket Island, while three haplotypes C9–C11 were endemic to Phangnga (Pakarang Cape: PHA4). Haplotypes C2–C4 were only found in the Tankhen Bay, Phuket (PHU25) and the Lam Kaen beach, Phangnga (PHA1), and C6–C8 also solely occurred in Nakalay beach, Racha Noi (southeast coast Bay) and Koh He (South), Phuket, respectively, while C9 was only found in Pakarang Cape, Phangnga (PHA4) (Figure 1a, Table 1).

For cpDNA rbcL, 333 sequences (1242 bp after alignment) contained 98 variable sites and resolved four haplotypes (R1–R4, GenBank accession numbers PZ816513–PZ816516) in total. Among them, R1 was the most abundant haplotype, occurring in 19 sites (Figure 1b). Haplotype R2 was found at Lam Kean beach, Phangnga (PHA1-2022), Koh He (South), Phuket (PHU14) and Tankhen Bay, Phuket (PHU25), while R3 and R4 were private to Phangnga (Pakarang Cape: PHA4 and Lam Kaen beach: PHA1-2024, respectively) (Figure 1b). Haplotype network revealed cox1 haplotypes separated by 1–45 bp and rbcL haplotypes by 1–63 bp (Figure 1). Excluding singleton sites (e.g. PHU23, PHU51), the highest cox1 haplotype (h=0.6667) and nucleotide (π=0.056322×10^-2) diversity occurred in the population RNG2 from Koh Kham, Koh Phayam, Ranong, followed by the populations PHU25-2018 and PHU25-2023-07 from Tankhen Bay, Phuket. In comparison, rbcL revealed the highest haplotype (h=0.5000) and nucleotide (π=0.033011×10^-2) diversity in PHU25-2018 (Table 1). In the Myanmar Shelf, the site PHA1 harboured high numbers of haplotypes (C1–C4, C10–C11). In the Strait of Malacca Shelf, the site PHU25 from Tankhen Bay, Phuket exhibited relatively rich haplotype (C1–C4) diversity (Table 1).

3.2. Population Genetic Structure

AMOVA for both cox1 and rbcL indicated no significant genetic variation among groups. Instead, most genetic variation (90.28%–98.63%) resided within populations, with 3.23%–12.57% occurring among populations within groups (Table 2). Bayesian Structure analysis (optimal K=2) partitioned H. durvillei populations into two genetic clusters (K=2, Figure 2a–2b). The cox1-based analysis revealed low admixture in populations from Koh Phayam, Ranong (RNG2), Phangnga (i.e. PHA1, PHA4, PHA8), and Tankhen Bay, Phuket (PHU25), suggesting limited gene flow rather than genetic recombination or introgression (Figure 2c–2d).

3.3. Gene Flow and Population Demographic History

MtDNA cox1-derived migration analysis revealed clear asymmetric gene flows among 12 populations along the Andaman coast of Thailand (Figure 3a). Most populations (RNG1a, RNG2, PHA4, PHU50, PHU51, PHU38, and TRA5) served as both sources and sinks, with bidirectional genetic migrations. In the Myanmar Shelf, the highest effective number of migrants (Nm=89.14) was detected from PHU31 to PHU51, far exceeding the reverse direction (Nm=0.27). Similarly, effective migrants from PHU31 in the Myanmar Shelf to PHU38 in the Strait of Malacca Shelf (Nm=64.57) were nearly 250 times higher than the reverse (Nm=0.2591). In contrast, the populations PHA1, PHU31, PHU25, and TRA2 mainly acted as sources rather than sinks of migration, with very limited migrants received (Figure 3a). Genetic exchange of numbers of migrants and directions among adjacent populations during the northeast monsoon season (Figure 3b) were broadly similar to the patterns detected in Figure 3a.

BSP analysis indicated a stable effective population size for H. durvillei from 125 Mya to 15 Mya (Figure 4). A sharp decline began during the early middle Miocene (c. 15 Mya), followed by gradual demographic expansion starting in the early Pleistocene (c. 2 Mya) (Figure 4).

4. Discussion

4.1. Population Genetic Diversity Along the Andaman Sea

Both mtDNA cox1 and cpDNA rbcL revealed no genetic structuring and low population differentiation in H. durvillei along the Andaman coast of Thailand (Figure 1, Table 1). Such an intraspecific genetic diversity is notably lower than that of its congener H. malaysiana from Malaysia and the Philippines (cox1: h=0.574±0.085; rbcL: h=0.551±0.090) and nucleotide diversity (cox1: π=0.00395±0.00085; rbcL: π=0.00165±0.00027) [29]. However, the observed genetic diversity aligns with patterns reported for other macroalgae in the Indo-Western Pacific [8,16,45]. For example, the red seaweed Gracilaria salicornia harboured low-to-moderate haplotype (h=0.356–0.760) and nucleotide diversity (π=0.00081–0.00261) across the TMP, despite its populations in the Andaman Sea having higher genetic diversity than those in the Gulf of Thailand [19]. Similar population differentiation patterns were also reported in the red macroalgal congeners Bostrychia tenella [17] and B. binderi [18]. In this region, low population genetic diversity has also been documented in the brown macroalga Sargassum plagiophyllum C.Agardh (cox1: h=0.000–0.489, π=0.00000–0.00090) [8] and the green macroalga H. macroloba (tufA: h=0.000–0.343, π=0.000–0.041) [16]. The exception is S. polycystum, which exhibits moderate genetic diversity (cox1: h=0.000–0.583, π=0.000–0.156; cox3: h=0.000–0.655, π=0.000–0.329; cox1+cox3: h=0.314–0.739, π=0.033–0.182) [8], likely due to its broader distribution across the TMP. These series of phylogeographic evidence, together with gene-flow patterns observed in other macroalgae [8,45], imply that ocean-circulation mediated migration interacted with physical barriers such as the biogeographic break north to Phuket province [17], likely played a dominant role in shaping phylogeographic homogeneity along the Andaman Sea.

4.2. Glacial Refugium and Demographic History

Marine glacial refugia, characterized by high genetic diversity and endemism [46,47], are critical for species survival during sea-level fluctuations. During the glacial maximum, extensive exposure of the continental shelf around the Andaman Sea due to around 130-meter of sea-level decline resulted in the compression and fragmentation of large coastal habitats [48]. Such a tectonic re-configuration led to marine species with narrow habitat preferences to be isolated in scattered refugia [4]. In this study, BSP analysis revealed a relatively stable effective population size from 125 Mya to 15 Mya, followed by a sharp demographic contraction from 15 Mya to c. 2 Mya and a progressive expansion during the early Pleistocene (Figure 4). Such a dynamic demography history is consistent with the classic “glacial contraction–interglacial expansion” refugium model [4,8,47]. The presence of endemic cox1 haplotypes (C2–C4) in the population PHU25 (Tankhen Bay) suggests that isolated habitats in the northern Malacca Strait [5,8,49] served as a marine refugium for H. durvillei during the maximum ice extent [6]. Postglacial expansion of H. durvillei outward from glacial refugium was likely driven by ocean currents along the Andaman coast (see discussion below) [8,16,17,45].

However, the demographic estimates for H. durvillei (Figure 4) were much earlier than those reported in other macroalgae in Southeast Asia. For instance, the brackish-water red alga Gracilaria tenuistipitata was documented to expand its distribution range and population size during the middle Pleistocene [50]. The brown algae S. polycystum and S. plagiophyllum experienced demographic expansion from 0.015 Mya to 0.060 Mya, corresponding to the late Pleistocene [8]. These dissimilar results imply that marine macroalgae across the TMP may respond differently in magnitude to paleo-climate change and oceanographic fluctuations. The stable demographic phase (125–15 Mya) may ascribe to extensive, well-connected shallow-water habitats available in the Andaman Sea prior to 15 Mya. During this period, regional tectonic configuration was relatively quiescent, and the Alcock and Sewell Rises were still joined to the continental margin of Sundaland [20]. Such a habitat continuity likely favoured the long-term persistence of H. durvillei populations in the Andaman Sea. The sharp demographic contraction at c. 15 Mya coincides with a major tectonic re-organization in the Andaman Sea. During the early middle Miocene to the earliest late Miocene (c. 15–10 Mya), the previously conjoined Alcock and Sewell Rises began to rift away from the Sundaland margin, forming the east Andaman basin [20,21]. This rifting event fragmented the once continuous shallow-marine habitats, reducing the area and connectivity of suitable substrata for benthic macroalgae such as H. durvillei. The subsequent prolonged population size decline (15–2 Mya) likely reflects persistent habitat instability through the late Miocene and Pliocene. During this interval, regional stress field rotated from east–west extension to NNW–SSE transtension, accompanied by intensified strike-slip faulting along the South Sagaing Fault and other structures [21]. The demographic expansion at c. 2 Mya corresponds to the onset of the early Pleistocene. This period was characterized by high-amplitude glacial–interglacial sea-level fluctuations, which repeatedly exposed and flooded the Sunda Shelf and reshaped coastal marine habitats [1] to facilitate H. durvillei to expand its distribution range.

4.3. Ocean Circulation and Phylogeographic Homogenization

Genetic exchange during the northeast monsoon season show similar patterns to the whole years (Figure 3), indicating that the strong circulations produced driven by this seasonal monsoon considerably shaped the dispersal mode of H. durvillei populations along the Andaman coast of Thailand. The observed phylogeographic homogeneity in H. durvillei likely arises from interactions between ocean circulation and historical biogeographic barriers [8,16,17,45]. For example, the biogeographic break north of Phuket [17] may act as a partial filter, while the Andaman Sea’s dynamic currents (e.g. the North Indian Ocean and Malacca Strait currents) facilitate gene flow. These currents generate coastal upwelling and mixing [51], enabling both short-distance propagule dispersal and long-distance rafting of vegetative fragments [52]. Similar mechanisms have been proposed for Bostrychia tenella [17] and Gracilaria spp. [19], underscoring the role of oceanography in shaping genetic connectivity across the region.

5. Conclusions

The low genetic diversity and phylogeographic homogeneity of H. durvillei highlight vulnerabilities to climate change and anthropogenic pressures. Conservation strategies should prioritize in-situ protection of unique genotypes (e.g. cox1 haplotypes C6–C9) and critical populations (e.g. PHU25, PHA1), which exhibit high diversity and endemism. The gene flow patterns along the Andaman Sea driven by ocean current during the northeast monsoon season also highlight geographic areas where populations were bi-directionally and/or multi-directionally connected, providing additional information for marine spatial management. A network of marine protected areas (MPAs) could leverage the bidirectional gene flow patterns identified here to maintain connectivity. Future studies should integrate dispersal modeling (e.g. Lagrangian particle simulations [53] and multigenerational dispersal probabilities [54] to quantify how ocean currents and habitat continuity shape genetic structure, informing targeted management of this ecologically significant species.