INTRODUCTION
Genetic variability is an essential component of biodiversity. It allows for populations to persist and increase their adaptation to environmental changes. Recent climatic changes and coastal habitat degradation increase habitat disturbance that is driving a decline in biodiversity (Leigh et al. 2019). In particular, loss of genetic variability is unavoidable in closed brackish water communities that are confronted with drastic environmental changes and habitat fragmentation (Sromek et al. 2019). However, genetic variation patterns of red algae that occur in closed brackish waters remain understudied.
Oceanographic currents such as the Kuroshio Current have been reported one of the major driving forces for species distribution and population structure of seaweeds in Southeast Asian waters (Hu et al. 2017, Liang et al. 2022, Muangmai et al. 2023, Fontana et al. 2024). Winter monsoon currents flow from northeast to southwest direction, while summer monsoon current flow in the reverse direction (Fang et al. 2009). These currents have facilitated both northeastern and southwestern dispersals of marine animals and brown algae and shaped their population structures (e.g., Wang et al. 2016, Hu et al. 2017, Liang et al. 2022). The sea-level changes caused divergence and population structure of brown algae (Chan et al. 2014, Hu et al. 2017, Liang et al. 2022). Biogeographical barriers such as Malay-Thai Peninsula have also formed population structure of brown and red algae (Wichachucherd et al. 2014, Bulan et al. 2022, Muangmai et al. 2023).
Gracilaria tenuistipitata C. F. Chang & B.-M. Xia (Rhodophyta) was selected for deciphering genetic variability in Southeast Asia because of its occurrence in brackish waters, as it has been used for agar production, and because our results could be directly compared to previous studies (Song et al. 2015, Yang and Kim 2015, Wang et al. 2023). It was first named for specimens growing in the sublittoral region of low salinity, Guangdong province, China (Chang and Xia 1976). Later, Zhang and Xia (1988) described var. liui for specimens cultured in a pond in Haikou, Hainan Island, and accordingly, var. tenuistipitata was assigned to the original species. Morphologically, var. tenuistipitata displays 1–2 orders of elongate branches, and var. liui has percurrent axes bearing numerous, delicate, and short to long flagelliform lateral branchlets (Zhang and Xia 1988). Previous molecular studies accepted G. tenuistipitata alone (Song et al. 2014, 2015, Yang and Kim 2015, Gurgel et al. 2018, Iha et al. 2018, Lyra et al. 2021), while a few molecular papers studied on G. tenuistipitata var. liui (Hagopian et al. 2004, Wang et al. 2023). Hereafter, we merged publicly available sequences under the name of G. tenuistipitata.
Gracilaria tenuistipitata occurs in brackish waters of nontidal lagoons, ponds, and water canals (Chang and Xia 1976, Lewmanomont 1994, Titlyanov et al. 2011, Das et al. 2022). It is a euryhaline and eurythermal species having a wide range of salinity, temperature, and heavy metal tolerance (e.g., Haglund and Pedersén 1992, Lee and Chang 1999, Tonon et al. 2018). As an important species for food and agar (Montaño et al. 1999, Yarnpakdee et al. 2015), G. tenuistipitata has been cultivated for biomass or co-cultivated with fishes and shrimps in Vietnam and Bangladesh (An and Anh 2020, Ullah et al. 2023). The first complete chloroplast genome of G. tenuistipitata was published among florideophycean red algae (Hagopian et al. 2004), and complete mitochondrial genomes were later analyzed in separate laboratories (Iha et al. 2018, Liu et al. 2018). Interestingly, the Southeast Asian G. tenuistipitata shared the most recent common ancestor with the South American G. chilensis C. J. Bird, McLachlan & E. C. Oliveira (Gurgel and Fredericq 2004, Yang and Kim 2015, Iha et al. 2018, Wang et al. 2023). Previous genetic variability studies on G. tenuistipitata raised a number of questions that remain unresolved. Song et al. (2015) reported a low cox1 (aligned to 1,240 bp) diversity within the species in Southeast Asian waters (about 0.5%). Yang and Kim (2015) further reported a low pairwise divergence (about 0.7%) of 5′ region of cytochrome c oxidase subunit I (COI-5P; aligned to 616 bp), being a lower value compared to those of other Gracilaria species (Wang et al. 2023). However, range-wide sampling studies are required in order to better understand patterns of genetic variability diversity of this economically important species.
The DNA barcoding marker COI-5P has been used for investigating genetic variability within populations as well as identification of species in Gracilaria and other red algae (e.g., Yang et al. 2008, Song et al. 2015, Yang and Kim 2015, Boo et al. 2019, 2023, Muangmai et al. 2023). Boo et al. (2020) reported similar haplotype numbers and nucleotide diversity for cox1 (1,602 bp) and COI-5P (664 bp) from 10 populations of Gelidiella fanii S. M. Lin, supporting the utility of the barcoding marker in population studies. Maternal inheritance also makes mitochondrial DNA more sensitive to demographic events of populations due to non-recombination (Ludt et al. 2012). In the present study, we generated new COI-5P sequences from populations from Vietnam (n = 91), and combined these sequences with publicly available data from GenBank (n = 70), to provide a broader phylogeographic understanding of G. tenuistipitata. The objective of this study was to examine genetic variability of Gracilaria tenuistipitata and estimate its demographic histories and range expansion.
MATERIALS AND METHODS
Sampling and habitat survey
Samples were collected at six locations in 2022 and 2023 based on previous records of G. tenuistipitata (as G. vermiculophylla) (Pham 1969, Nguyen 1992) along the coastline of Vietnam (Fig. 1, Supplementary Table S1). At each location, 4 to 27 individuals were collected at 5 m intervals. Habitats included nontidal lagoons, estuaries, shrimp ponds, and water canals (Fig. 2). At each location, the salinity concentrations were measured using a salinity refractometer Kruss HR27-100 (Kruss, Hamburg, Germany). Water depths were recorded using a Geotech Interface Meter (Field Environmental Instruments Inc., Pittsburgh, PA, USA) and surface water temperatures were measured with the Multiparameter Hi982 (Hana Instruments, Smithfield, RI, USA). Annual surface water temperature data were retrieved from the NCHMF Vietnam (https://nchmf.gov.vn/). Representative specimens were mounted on herbarium sheets as vouchers and deposited in the Institute of Tropical Biology, Ho Chi Minh City, Vietnam. Fragments of each specimen were cleaned to remove epiphytes and then stored in silica gel for DNA extraction.
DNA extraction and COI-5P sequencing
DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing procedures followed Boo et al. (2019). Genomic DNA was extracted from ~5 mg dried thalli ground in liquid nitrogen using the Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer’s protocol. The COI-5P primers used for amplifying and sequencing were GazF2 and GazR2 (Saunders 2005). The PCR was performed using the Gotaq Green Master Mix (Promega) with 5 μL of DNA template. The PCR thermocycling protocol consisted of an initialization step of 94°C for 4 min followed by 35 cycles of 94°C for 30 s (denaturation), 45°C for 30 s (annealing), and 72°C for 1 min (extension), and a final extension of 72°C for 10 min. PCR products were purified using the LaboPass Gel & PCR purification Kit (Cosmo Genetech, Seoul, Korea) and sequenced commercially by the 1st Base Company (Selangor, Malaysia). Newly generated sequences were deposited in GenBank (Supplementary Table S1). All sequences were aligned using the MUSCLE algorithm in MEGA7 (Kumar et al. 2016) with default parameters, followed by manual adjustments and trimming to match the length of the shortest sequence (465 bp).
Genetic variability and phylogeographic analyses
The genetic variability indices were calculated using DnaSP v.6 (Rozas et al. 2017) at population, country, and ecoregion levels following Spalding et al. (2007): the number of haplotypes (h), the number of polymorphic sites (S), haplotype diversity (Hd), and nucleotide diversity (π). The haplotype network was constructed with PopART v.1.7 (Leigh and Bryant 2015) using the median-joining network (Bandelt et al. 1999). The phylogeny of COI-5P haplotypes was reconstructed using the Maximum likelihood (ML) optimality criterion on the W-IQ-tree webserver (Trifinopoulos et al. 2016). The best-fitting substitution model was determined with the model test option (auto), followed by the ML tree search, and 1,000 ultrafast bootstrap replicates. Gracilaria chilensis and G. vermiculophylla (Ohmi) Papenfuss were included as outgroups based on their previously determined close evolutionary relationship to G. tenuistipitata (Lyra et al. 2021, Wang et al. 2023). Non-hierarchical and hierarchical analyses of molecular variance (AMOVA) were performed using Arlequin v.3.5 (Excoffier and Lischer 2010) with Φ-statistics to quantify the proportion of total genetic variance, with significance of fixation indices tested using 10,000 permutations.
Mismatch distribution analysis (MDA) was conducted to test the null hypothesis of spatial expansion using Arlequin v.3.5. For the expansion model, goodness-of-fit was tested with the sum of squared deviations (SSD) and Harpending’s raggedness index (HRag) using 1,000 parametric bootstrap replicates (Schneider and Excoffier 1999). The neutrality tests (Tajima’s D and Fu’s FS) were also performed to infer potential population growth and expansion (Tajima 1989, Fu 1997).
BEAST v.2.7.3 (Bouckaert et al. 2019) was used to infer demographic histories by constructing Bayesian skyline plots (BSPs) of effective population size through time (Heled and Drummond 2009). The substitution rate for COI-5P was set at 7.6 × 10−9 substitutions site−1 y−1, a value previously calculated for the red algae (Bringloe and Saunders 2019). The Markov chain Monte Carlo was run for 1 × 107 generations with trees sampled every 1,000 generations and the first 10% of the samples discarded as burn-in. The result was visualized by Tracer v.1.7 (Rambaut et al. 2018). Three replicate runs using different random seeds were conducted to confirm convergence.
To compare the probabilities of dispersal routes, migration rates were estimated using Migrate-n v.3.2.7 (Beerli 2006, Beerli and Palczewski 2010). Populations from India and the Philippines were excluded from the analysis due to having fewer than five individuals. Twenty populations were pooled based on geographic patterns and the haplotype network as follows: (1) China, (2) Vietnam, and (3) the Malay-Thai Peninsula (Malaysia, Singapore, and Thailand). The coalescent analysis focused on three potential migration routes: China vs Vietnam, Vietnam vs the Malay-Thai Peninsula, and China vs the Malay-Thai Peninsula.
An unconstrained migration model was applied to estimate migration rates and population sizes using a uniform prior divided into 1,500 bins. Bayesian inference was performed with two long chains at different heating levels (1, 1.5, 3, and 10,000), sampling 4 million genealogies at an increment of 100. The first 25% of the samples were discarded as burn-in. Convergence criteria included an effective sample size greater than 1,000 for all parameters and consistency in parameter estimates across three independent runs with different starting points.
RESULTS
In total, 161 sequences of Gracilaria tenuistipitata including 90 newly generated in this study were used for the analyses (Supplementary Table S1). G. tenuistipitata had high haplotype diversity (Hd = 0.725 ± 0.030) but low nucleotide diversity (π = 0.00243 ± 0.00020) (Table 1). The Chinese populations (Hd = 0.564–0.833, π = 0.00133–0.00645) showed higher haplotype and nucleotide diversity than the Vietnamese (Hd = 0.186–0.556, π = 0.00047–0.00123). At the population level, the haplotype and nucleotide diversities were highest at Qingdao (Hd = 0.833 ± 0.222, π = 0.00645 ± 0.00280). At the country level, haplotype diversity was highest in Malaysia (Hd = 0.649 ± 0.126) and lowest in Vietnam (Hd = 0.424 ± 0.055), while nucleotide diversity was slightly higher in China (π = 0.00182 ± 0.00058) and Malaysia (π = 0.00155 ± 0.00041) compared to Thailand (π = 0.00115 ± 0.00026) and Vietnam (π = 0.00111 ± 0.00020). At the ecoregion level, the Yellow Sea (Hd = 0.833 ± 0.222, π = 0.00645 ± 0.00280), and Southern China (Hd = 0.638 ± 0.129, π = 0.00164 ± 0.00043) showed higher haplotype and nucleotide diversity than those from other ecoregions (Table 1).
Haplotype network and distribution roughly grouped haplotypes into three geographical groups; the northern comprised H12, H13, H14, H15, and H16 from China, the central included H1, H2, H3, H4, H5, H6 from Vietnam and also H1 from China, India and Thailand, and the southern included H7, H8, and H9 from Malaysia, Singapore, and Thailand (Fig. 3). Haplotypes between and within these groups were connected by single mutational steps, except for H10 and H11, which were separated from H1 by 3–4 mutational steps.
The ML tree of 16 G. tenuistipitata haplotypes was monophyletic (Supplementary Fig. S1). The pairwise distances within species ranged from 0.02 to 1.5% (Supplementary Table S2). Sixteen haplotypes were identified in G. tenuistipitata, with 10 private haplotypes found in single populations (Fig. 3, Supplementary Table S1). Haplotype H1 (n = 76) was the most dominant and widely distributed haplotype in China, India, Thailand, and Vietnam. Within China, H12 was common (n = 24, 66.7% of specimens analyzed), along with individuals from the Philippines. Additionally, four other haplotypes (H13–H16) from China were linked by a single mutation to H12. Haplotype H7 (n = 19) was found in Singapore, Malaysia, and Thailand, and two other haplotypes (H8 and H9) from Malaysia were connected by a single step.
Non-hierarchical AMOVA showed that 67.97% (p < 0.001) of the genetic variation was found among populations (Table 2). The hierarchical AMOVA showed that 56.37% (p < 0.001) of variation occurred among countries, followed by within population variation (25.52%, p < 0.001). At the ecoregion level, the genetic variation was similar across different sources of variation (among ecoregions: 34.93%, among populations within ecoregions: 35.20%, within populations: 29.87%).
The values of neutrality test (Tajima’s D and Fu’s FS) at the species level of G. tenuistipitata were significantly negative (D = −1.80477, p < 0.05; FS = −9.353, p < 0.01), as well as in China (D = −1.82446, p < 0.05; FS = −3.066, p < 0.05) and Vietnam (D = −1.46632, p < 0.05; FS = −3.696, p < 0.05) (Table 1). At the ecoregion level, only Fu’s FS value was statistically significant in the Gulf of Tonkin (FS = −4.454, p < 0.01), which may indicate population expansion or a recent selective sweep.
The mismatch distribution of all 22 populations exhibited a smooth unimodal with low and non-significant SSD and Harpending’s raggedness (HRag) values, indicating a good fit between the observed and expected values, and supporting population expansion (Table 3, Fig. 4A). The estimated time of expansion for G. tenuistipitata was about 235,000 years ago with a range of 154,000–1,432,000 years ago (Table 3). Vietnamese populations (about 193,000 years ago with a range of 108,698–312,259) underwent earlier population expansion than Chinese populations (178,000 years ago with a range of 86,549–780,044). The BSP showed that the divergence time from the most recent common ancestor was about 300,000 years ago, and indicated steady population growth up to 50,000 years ago, followed by stable population sizes until the present (Fig. 4B).
The estimates of migration rates (M) and effective population size (Ɵ) are shown in Table 4. Migration from Vietnam to China was found to be more extensive than migration in the opposite direction (mean rate: 336.2 vs. 135.3). Similarly, southward migration from Vietnam to the Malay-Thai Peninsula was estimated to be higher than migration in the opposite direction (mean rate: 513.9 vs. 167.3). Theta (Ɵ) estimates indicated that Ne was similar in China and Vietnam but lowest in the Malay-Thai Peninsula (Table 4).
DISCUSSION
Our analysis of mitochondrial DNA data in Gracilaria tenuistipitata provides fundamental insights into the evolutionary history and genetic adaptation of the brackish water red algae in Southeast Asia. Although haplotype and nucleotide diversity were similar with those of other marine congeners (e.g., Muangmai et al. 2023), the distribution and relationship of haplotypes appear to suggest a northern refugium between Vietnam and China. Migration analysis indicated that G. tenuistipitata populations spread both northward and southward from Vietnam. Three latitudinal haplogroups (northern, central, and southern) might have been shaped by a combination of climate and currents. Both mismatch analysis and BSP revealed a gradual demographic or range expansion in the middle Pleistocene without signature of bottlenecks. Similarly, Sargassum polycystum C. Agardh from Southeast Asia showed flat curve without bottleneck signature (Chan et al. 2013), however, most marine animals and seaweeds in temperature regions exhibit bottlenecks associated with low sea-level stands or other environmental stresses (Ludt and Rocha 2015, Hu et al. 2017, Boo et al. 2019). Further studies on diverse species are required to conclude whether the exception is species-specific or a more general characteristic of tropical red algae.
Haplotype and nucleotide diversity of G. tenuistipitata (Hd = 0.725 ± 0.030, π = 0.00243 ± 0.00020) are similar to or lower than those of G. salicornia from Thailand and G. vermicullophylla from Northeast Asia (Kim et al. 2010, Zhong et al. 2020, Muangmai et al. 2023). In contrast, Gelidiophycus divaricatus (G. Martens) G. H. Boo, J. K. Park & S. M. Boo, having a similar distribution in Southeast Asia, revealed much high genetic variability (Hd = 0.920, π = 0.01716) (Boo et al. 2019). Interestingly, the genetic diversity values of G. vermiculophylla decreased when introduced specimens were included (Hd = 0.517, π = 0.00178) (Kim et al. 2010). In addition, with increased sampling (n = 611), the haplotype diversity decreased (Hd = 0.6632), while the nucleotide diversity increased (π = 0.003766) (Zhong et al. 2020). Comparison of haplotype and nucleotide diversity therefore should cautiously be compared between populations or species.
Sixteen haplotypes comprise the Southeast Asian intraspecific variability in G. tenuistipitata. Haplotype H1 is considered as the ancestral haplotype because its common occurrence in countries surveyed and the central position in the star-shaped haplotype network. Migration analysis revealed that G. tenuistipitata populations dispersed southwards and northwards from Vietnam, considered a genetic center of the species. This result suggests that regional monsoon currents promoted gene flow and share genetic variations between distantly located populations, as reported in other seaweeds in previous studies (Chan et al. 2014, Hu et al. 2017, Liang et al. 2022, Fontana et al. 2024). Considering the small number of haplotypes in Thailand and Malaysia, as reported in previous studies (Song et al. 2015, Yang and Kim 2015), the reverse dispersal by summer monsoon currents may have less work. However, Sargassum polycystum is reported to have dispersed southward from a refugium, Hainan Island, to the Gulf of Thailand and eventually to the Strait of Malacca (Hu et al. 2017, Liang et al. 2022).
Together with H1, haplotypes H7 and H12 shared in different countries indicate the genetic connectivity as well as possible long dispersal of this brackish water species. It is hypothesized that fragmented branches or holdfasts of the ancestral haplotype, attached on shells or gravels, likely drifted by monsoon currents or other local currents (Winston 2012, Boo et al. 2023). The genetic connectivity in Southeast Asia has been reported in various marine algae (Chan et al. 2013, 2014, Hu et al. 2017, Ng et al. 2017, Liang et al. 2022, Muangmai et al. 2023). Further mutations occurred from H12 to H13, H14, H15, and H16 in the Chinese region (except one individual in H12 from the Philippines) and also from H7 to H8 and H9 from Malaysia, Singapore, and Thailand (Fig. 3A). Private haplotypes derived from H1, H7 and H12 revealed local adaptation that is typically a recent phenomenon at least on an evolutionary time-scale (Sjöstrand et al. 2014). Various pools with different salinities formed due to sea-level changes in Southeast Asia and subsequent environmental changes (Voris 2000) might have induced such genetic variability.
Both MDA and BSP, common tools for reconstructing historical demography (Grant 2015), revealed that demographic or range expansion of G. tenuistipitata populations started to occur during the Pleistocene. Interestingly, populations from China and Vietnam underwent expansions during a similar time period (Table 3). Our data is likely supported by recent molecular clock analysis that G. tenuistipitata and G. chilensis diverged during the early Cenozoic (Lyra et al. 2021). This may be the characteristics of tropical benthic species that have persisted in waters of less climatic changes during the Pleistocene than temperate species. However, Pleistocene population expansion is common in marine brown algae, seahorses, and intertidal limpets (e.g., Zhang et al. 2014, Wang et al. 2016, Liang et al. 2022). The lack of bottleneck signatures may be the result of habitat isolation of G. tenuistipitata during low sea level in the Pleistocene that has thereafter come back into contact, as reported in oceanic lagoon species (Ludt et al. 2012).
The habitats of G. tenuistipitata, mostly nontidal lagoons, ponds, and water canals along the coast of Vietnam showed a salinity range of 12–25 parts per thousand (ppt) during the surveys (Supplementary Table S3); however, they were sometimes full of freshwater during rainy seasons in a very short period. G. tenuistipitata is likely able to occur in less saline seawaters as well as brackish water. Range-wide sampling is needed to confirm whether G. tenuistipitata occurs in marine waters.
The range of G. tenuistipitata in China, Indonesia, Malaysia, Singapore, Thailand, and Vietnam, has been demonstrated using COI-5P and or rbcL in previous studies as well as the present study (Song et al. 2014, 2015, Yang and Kim 2015, Wang et al. 2023). The finding of haplotype H1 from India and Hawaii (Yang et al. 2008) requires analysis of additional specimens to decipher whether Indian specimens are homogeneous due to aquaculture or they have high diversity. The reports of G. tenuistipitata in Rhode Island and Virginia, USA (Gurgel and Fredericq 2004, Garcia-Rodriguez et al. 2013) resulted from a misidentification of G. vermiculophylla (Gurgel et al. 2018). Taken together, G. tenuistipitata is assumed to be native in Southeast Asia.
It is remarkable for G. tenuistiupitata to show a slight signal of genetic structuring among three geographical group (the northern Chinese, the central Vietnamese, and the southern in Malaysia, Singapore and Thailand). This result suggests that the phylogeography of G. tenuistipitata has likely been shaped by climate and environmental factors and reflects potential structuring by marine ecoregion, as shown in biogeographic patterns of marine organisms (Spalding et al. 2007). If this will be supported by further study, geographical structure by climate may be characteristics of brackish water species. Mueller et al. (2018) reported that Hormosira banksia (Turner) Decaisne, a brown alga in marine and estuarine waters in Australia, revealed similar patterns of genetic structuring to other marine species in the region. However, Coleman et al. (2019), using five microsatellite markers, reported complex demographic processes of the same species, with strong isolation by distance, including clonality in estuarine populations. Generalizing of phylogeographic pattern of G. tenuistipitata as a brackish water species needs further understanding of other brackish water marine algae.
Our analysis of mitochondrial COI-5P sequences covering the Southeast Asia formed a monophyletic clade, indicating a single genetic species. The highest pairwise divergence of COI-5P was 1.5% (Supplementary Table S2), including sequences in previous studies (Song et al. 2015, Yang and Kim 2015, Wang et al. 2023). This value is lower than or similar to the intraspecific range of Gracilaria and other red algae (Yang and Kim 2015, Ng et al. 2017, Fumo and Sherwood 2023). These molecular data don’t support previous reports on the occurrence of var. tenuistipitata an var. liui from Thailand and Vietnam (Lewmanomont 1994, Tsutsui et al. 2005). Wang et al. (2023) reported 0.2% pairwise divergence in both COI-5P and rbcL within var. liui from China. Moreover, Haplotype H12 (MH396021) (Iha et al. 2018) from Haikou city, Hainan (the type locality of var. liui) differed just by a single nucleotide substitution from H1 (OP669464) (Wang et al. 2023) from Diancheng, Maoming city, near Bohe, Guangdong Province (the type locality of G. tenuistipitata). It is proposed that G. tenuistipitata is likely a single species including morphotypes with variable branching patterns. DNA analysis of the type specimen of G. tenuistipitata is needed to confirm this proposal, but it is beyond the scope of the present study.
In conclusion, our study reveals that Gracilaria tenuistipitata comprised 16 haplotypes and revealed a slight signal of three geographical groups in Southeast Asia. This pattern has most likely been driven by a combination of habitat availability, currents, sea level, and climate. Analyses using more variable markers such as single nucleotide polymorphism and microsatellite could detect fine-scale genetic breaks of G. tenuistipitata by ecological barriers such as salinity, temperature, habitat, and others. Recent development of coastal areas in Southeast Asia may deprive habitats of brackish water species and lead to extirpation of G. tenuistipitata populations therein. Minimizing the change of brackish water habitats is a key factor for conserving populations with endemic haplotypes as well as this agar-yielding species. Continued study is needed to detect changes in population size and range shift of G. tenuistipitata due to coastal habitat degradation and climatic change.
