Barrier to Gene Flow of Grey Mangrove Avicennia marina Populations in the Malay Peninsula as Revealed From Nuclear Microsatellites and Chloroplast Haplotypes

Introduction

Coastal communities living adjacent to the mangrove ecosystems in tropical, subtropical and warm temperate regions of the world are receiving manifold ecological and socio-economic benefits since the ancient times (Walters et al., 2008; Spalding et al., 2010; Lee et al., 2014; Seddon et al., 2020; Dahdouh-Guebas et al., 2021). The mangrove species occupying different environmental gradients in a habitat serve collectively for land building, coastal protection, water quality improvement, phytoremediation, carbon sequestration, breeding or nursery ground to various aquatic and terrestrial fauna (Donato et al., 2011; Cohen et al., 2013; Analuddin et al., 2017; del Valle et al., 2020; Osland et al., 2020; Kathiresan et al., 2021). Due to the increased public awareness and extensive mangrove plantation schemes in recent years, the (average) annual loss of global mangrove cover declined to 2–4% and opened the doors for conservation optimism (Dahdouh-Guebas et al., 2020; Friess et al., 2020). Contemporary mangrove forest areas took shape historically and their genetic connectivity depends on sea-faring propagules, subsequent establishment, and persistence in suitable environments (Van der Stocken et al., 2019a). Mangrove seed/propagules experience both short to long-distance dispersal by tidal action (hydrochory) and grow either close to the mother tree or in other suitable locations. In this context, the physical land barriers and water current or circulation patterns were found to be crucial in bringing changes to the species' abundance and distribution of plants (Van der Stocken et al., 2019b) and animals (Fairuz-Fozi et al., 2021).

Several researchers have observed the land-barrier effect of Malay Peninsula and found genetically different animal populations (e.g., Tachypleus gigas, T. tridentatus, Carcinoscorpius rotundicauda, Barbonymus schwanenfeldii, Varuna litterata) between East and West coasts (Kamarudin and Esa, 2009; Yang et al., 2009; Ismail et al., 2011; Adibah et al., 2015; Liew et al., 2015; Suppapan et al., 2017; Fairuz-Fozi et al., 2021). In case of mangrove, species such as Avicennia alba Blume, Bruguiera gymnorrhiza (L.) Lamk., Ceriops decandra (Griff.) Ding Hou, Lumnitzera racemosa Willd. and Sonneratia alba J. Smith also revealed such genetic separation (Su et al., 2006; Liao et al., 2007; Huang et al., 2008; Wee et al., 2020). However, Wee et al. (2020) elucidated that the Malay Peninsula is acting like a filter, rather than a barrier, to the gene flow of mangroves by considering dispersal potential of its seeds or propagules. This East-West divide is therefore an increasingly documented phenomenon explained from the Sunda shelf flooding during the Holocene and likewise re-population of East and West Malay coastlines from different source populations since the Last Glacial Maximum. During the global de-glaciation period, the flooding of the Sunda Shelf facilitated mangrove species to undergo range shift toward the Malay Peninsula East and West coast (Wee et al., 2014, 2015). Peninsular Malaysia has ca. 110,952 ha of mangrove cover (18% of the country's mangrove extent), of which nearly 17,570 ha (3%) are distributed along the East coast and 93,382 ha (15%) on the West coast (Hamdan and Misman, 2020). The less mangrove cover on the East coast is due to strong waves and current from the South China Sea, especially during the northeast monsoon, whereas the West coast facing the Strait of Malacca and sheltered by Indonesia receives a weaker current (Akhir et al., 2015; Zainol et al., 2021). Such contrasting ocean dynamics may also result in a different spatial genetic structure of mangroves, both at regional and local scales (Triest, 2008).

Mangrove propagules do not have a dormant stage and propagule dispersal is affected by factors such as buoyancy, propagule viability and timely establishment (Rabinowitz, 1978). Thus, the persistence of the population depends solely on the formation, release, distribution and establishment of propagules as clonal growth and vegetative dispersal is absent in mangroves. Once established, and before reaching sexual maturity, propagules and seedlings are subject to predation (Dahdouh-Guebas et al., 2011), environmental factors affecting early growth (Krauss et al., 2008) and anthropogenic pressure (Cannicci et al., 2008). Tidal influence, ocean currents and wind action predict possible distribution patterns and colonization of mangrove species over long oceanic distances (Clarke, 1993; Van der Stocken et al., 2015, 2019a). Many studies illustrated a long-range connection of mangrove trees in relation to ocean currents and direction, especially along same coastlines (Mori et al., 2015a; De Ryck et al., 2016; Ngeve et al., 2017; Hodel et al., 2018; Van der Stocken et al., 2019b). Among others, a stepping-stone model of migration between estuaries was portrayed for Avicennia species (Do et al., 2019; Wee et al., 2020; Triest et al., 2021a,b,c). This contributes to the importance of coastal connectivity through dispersal of propagules, thus gene flow, which is the only natural cohesive force between longer-term estuaries for a species to maintain its evolutionary units. However, as explained previously, the barriers to genetic connectivity may come from land masses and different migration histories (Triest, 2008; Hodel et al., 2018; Wee et al., 2020; Triest et al., 2021b), opposing ocean currents (Mori et al., 2015b; Ngeve et al., 2017) or very large rivers (Triest et al., 2018). In general, the expectation that migration routes followed major ocean and coastal currents is largely confirmed by genetic diversity and population-level genetic structure approaches.

However, connectivity patterns of mangrove settlement between estuaries of the same coastline and different habitats thereof are more complicated. Estuarine landscapes are highly diverse and unique in their complexity such that the establishment of mangrove propagules is dependent on suitable habitats due to sedimentation patterns of coastal and major river systems, channels or creeks and sandbar dunes (Triest and Van der Stocken, 2021). Spontaneous processes of propagule dispersal and mangrove vegetation formation are expected to result in different neighborhood sizes of individual trees in exposed seaward vs. sheltered landward populations (Triest and Van der Stocken, 2021), distance from open sea (Triest et al., 2021a), river flow intensity (Ngeve et al., 2017; Chablé Iuit et al., 2020), channel structures (Triest et al., 2020) or degree of fragmentation (Hasan et al., 2018; Bryan-Brown et al., 2020). Persistence of Avicennia trees is influenced by the coastal landform (Triest and Van der Stocken, 2021). The positioning of mangroves along the coast can show gradients of oceanic influences, from highly exposed coastal areas—projecting into open seas—to far inland estuaries along rivers (Ngeve et al., 2017; Chablé Iuit et al., 2020), or even no tidal influence at all (Triest et al., 2021a). Dispersal and settlement that have occurred over the past few and overlapping generations along a gradient of different tidal and ocean currents left traces in the locally recorded amount of genetic diversity (Triest et al., 2021a) and very often in their neighborhood associated fine-scale spatial genetic structure (FSGS) (Triest et al., 2020; Triest and Van der Stocken, 2021). Sufficient prior knowledge of polymorphic genetic markers is available for the geographically widespread gray mangrove Avicennia marina (Forsk.) Vierh. and related species (Triest, 2008) to allow resolution within an estuary and even of individual trees at the local fine-scale level (Hasan et al., 2018; Do et al., 2019; Triest et al., 2020, 2021a,c; Triest and Van der Stocken, 2021). Microsatellite loci of Avicennia species partially cross-amplify or show unique alleles at taxon level, so cases of hybridization are detectable (Mori et al., 2015a). Hybridization is not uncommon in mangroves (Ragavan et al., 2017) and admixtures or genetic introgression may represent cryptic evolutionary units or even be subject to “cryptic ecological degradation” (Koedam and Dahdouh-Guebas, 2008), a concept of conservation relevance.

The Malay Peninsula is expected to be a land barrier hindering gene flow between the East and West coast (Duke et al., 2002). Nevertheless, the difference in genetic structure across mangrove species suggest that the Malay Peninsula act as a “filter” rather than a strict geographical barrier, depending on the different dispersal capacities of mangrove species and on the ocean currents (Wee et al., 2020). Avicennia spp. having a low long distance dispersal potential are more likely to be restricted than other mangrove representatives. In this study, we verify for an East-West genetic break of A. marina in Malay Peninsula and more profoundly, at regional scale, focused on the connectivity pattern along each coastline. We hypothesize that connectivity between locations would be more pronounced along the East than to the West and that directionality of historically accumulated gene flow would follow main contemporary ocean currents. We specifically aim to (1) analyze the genetic diversity and structure of A. marina populations on both East and West coasts; (2) estimate the likelihood of different migration models between populations along a same coastline; and (3) compare the extent of a FSGS of A. marina at different distances from the open sea or tidal influence. To achieve this, we considered a combination of nuclear microsatellites and a chloroplast intron sequence (trnH-psbA). Understanding the processes and traces of spontaneous dispersal, establishment and persistence of mangrove areas in a regional context is extremely useful for future research, conservation and management priority settings.

Materials and Methods

Study Area

A total of 13 sampling sites were selected from nine mangrove locations in Peninsular Malaysia (Figure 1). Whereas, six sites from three locations namely, Tumpat (E1-A and E1-B), Tok Bali (E2-A and E2-B), and Setiu Wetlands (E3-A and E3-B) are located along an 80 km stretch of the East coast, seven sites from six locations namely, Penang (W1), Matang (W2-A and W2-B), Kuala Selangor (W3), Malacca (W4), Sungai Sebatu (W5), and Pontian (Johor) (SW) were along an 586 km stretch of the West coast (Figure 1; Supplementary Figure 1). On the East, mangroves in Tumpat were patchy, separated by intermittent water canals and largely exposed to the South China Sea (Satyanarayana et al., 2010), unlike they are river adjoining and sheltered from open sea by a sandbar at Tok Bali and Setiu Wetlands. For the West, all locations are directly exposed to the Strait of Malacca and only Matang mangroves were patchy and separated by major river channels. The pioneer A. marina is widely distributed in all locations of the present study. In fact, this species is common throughout its range and found across the Indo-Pacific (Tomlinson, 2016).

Figure 1. Overview of the studied populations of Avicennia marina along the West coast (W1 to W5 and SW) and of introgressed Avicennia marina on the East coast (E1 to E3) of the Malay Peninsula. Population codes are denoted as in Table 1 and in Supplementary Figure 1 with images of all locations.

Sample Collection

Fieldwork was conducted during the period of May-July 2017 and a total number of 486 samples were collected (179 from East coast and 307 from West coast). Young healthy leaves of A. marina, not damaged by insects or diseases, were handpicked from the tip of the branches of adult as well as juvenile vegetation. We considered ≥2.5 cm stem diameter (D₁₃₀) for the adult trees and saplings with three leaf pairs or less for the juveniles (Brokaw and Thompson, 2000; Goessens et al., 2014). Mature leaves were not considered due to the high content of polysaccharides, polyphenols and tannins (Ibrahim, 2011), which may disturb the amplification of DNA during the polymerase chain reaction. The sample collection was done through a transect approach for pure stands and random methods for mixed distant patches. The GPS (Garmin 45, USA) location of each sampled tree or juvenile was recorded to estimate the length of the transect line. The sampled leaves were collected into numbered paper envelopes, allowed to air or sun dry until no moisture traces, and further preserved with silica gel for transportation and handling within 1 month.

While observing transects at a given distance to the open sea (Table 1; Supplementary Figure 1), the interval between each sampled individual along the transect was also measured. The average number of sampled trees was 37 per transect though ranged exceptionally from 16 to 60 for a few sites (Table 2). Shortest distance intervals between the neighboring trees, mostly within 10 m, were considered. The gaps between mangrove patches were included into the total distance as fine-scaled analysis focuses on pairs of individual trees within a distance class below 100 m, regardless of the patch location.

Table 1. Location details of Avicennia populations studied along the Malay Peninsula.

Table 2. Population genetic variables of Avicennia marina transects in mangrove areas along Malacca Strait on the West Malay Peninsula and of introgressed Avicennia marina (with captured cpDNA of A. alba) on the East Malay Peninsula.

For comparative reasons we additionally considered samples of Avicennia species from outside the study area (Supplementary Table 2), namely A. marina from China, Vietnam, Bangladesh, Kenya and South Africa; A. alba, from Vietnam and The Philippines; A. rumphiana Hallier f. from The Philippines and Avicennia officinalis L. from Bangladesh that were available in the BRVU (part of BR herbarium) collection at the Vrije Universiteit Brussel.

DNA Extraction, Microsatellite Genotyping, and Chloroplast DNA Sequencing

Genomic DNA was extracted from approximately 20 mg of each dried leaf tissue using the E.Z.N.A. SP plant DNA Mini kit (Omega bio-tek, Norcross, GA, USA). The concentration of individual samples ranged from 10 to 200 ng/μl. The multiplexed PCR reactions consisted of 13 microsatellite markers: Avma1, Avma02, Avma6, Avma8, Avma10, Avma14, Avma17 (Geng et al., 2007); Am3, Am81, (Maguire et al., 2000a); Aa22, Aa23, Aa67 (Teixeira et al., 2003); and AMK6 (Triest et al., 2020). Primers were fluorescence-labeled with four different dye-labels (6FAM/VIC/NED/PET), and a mixture of 0.2 μM of each primer. About 6.25 μl master mix (Qiagen Multiplex PCR kit), 1.25 μl primer mix, 2.5 μl H₂O and 2.5 μl of genomic DNA were used for multiplex PCR reactions. PCR was performed in a thermal cycler (Bio-Rad MyCycler) under the conditions of initial denaturation at 95°C for 15 min followed by 35 cycles of 30 s denaturation at 95°C, 90 s annealing at 57°C, 80 s elongation at 72°C and a final extension of 30 min at 60°C. All PCR products were separated on an ABI3730XL sequencer (Macrogen, Seoul, Korea) and allele sizes were determined with GeneMarker v.2.60 (SoftGenetics LLC, State College, USA).

We used a non-coding region (trnH-psbA) with primer pairs aF-aR and cF-cR (Sang et al., 1997). PCR amplification was carried out in 25 μL of reaction mixture containing 2.5 μL of genomic DNA, 2.5 μL 10 × PCR buffer, 0.2 mM of each dNTP, 1.6 mM MgCl2, 200 nM of the forward and reverse primer, 80 μg mL⁻¹ bovine serum albumin (BSA) and 1 U Taq DNA polymerase (Promega GoTaq DNA Polymerase; Promega, Madison, USA). PCR reactions were performed in a thermal cycler (Bio-Rad MyCycler) and started with 2 min at 94°C, followed by 30 cycles of 45 s at 92°C, 1 min at 50°C for the reactions with trnH-psbA primers, 2 min at 72°C, and a final extension at 72°C for 5 min. All PCR products of trnH-psbA were sequenced on an ABI3730XL sequencer (Macrogen, Seoul, Korea).

Chloroplast Sequence Data Analysis

With an assumption of potential different gene pools, we tested the samples for chloroplast identity. The trnH-psbA sequences of the chloroplast genome in 3 to 7 samples from each site (43 in total) were identified and compared with Avicennia species from Asian (A. marina, A. alba, A. rumphiana and Avicennia officinalis) and African (A. marina) regions (Supplementary Table 2). Sequences of trnH-psbA were aligned (492 bp) and analyzed in Geneious Prime v.2019.2.1 software using MAFFT alignment and MrBayes Maximum Likelihood phylogenetic analysis (GTR substitution model with gamma rate variation; 4 heated chains of 1,100,000 chain length and 100,000 burn-in; Avicennia officinalis from Bangladesh was taken as outgroup).

Microsatellite Data Analysis

Prior to population and individual sample-based analysis, we tested data for genotypic disequilibrium, potential null alleles and overall resolution of the selected microsatellite markers. A linkage test between all pairs of loci (1,000 permutations) identified no genotypic disequilibrium at the 0.05 level (FSTAT v.2.9.3) (Goudet, 2001). MICRO-CHECKER indicated no scoring errors or large allele dropouts (Van Oosterhout et al., 2004). However, null alleles were present in all populations of the East coast for two loci (Avma14 and Avma2) and in all populations of the West coast (Avma14 and Aa22). We decided to omit these loci for profound data analysis of each regional group separately and thus considered only 11 out of 13 loci for analysis at either East or at West coast level (Supplementary Table 1).

The probability of identity (PI), namely whether two individuals could share an identical multilocus genotype by chance, gave a cumulative probability of identity for all polymorphic microsatellite loci in each site of 6.0 × 10⁻³ on average, thereby providing ample resolution (done with GenAlEx v.6.5, Peakall and Smouse, 2012). Basic population genetic variables were measured separately for each coastline and site: total number of alleles (A_T), mean number of alleles (A_M), effective number of alleles (A_E), allelic richness (A_R), observed heterozygosity (H_O), unbiased expected heterozygosity (uH_E), population inbreeding coefficient (F_IS)—with 1,000 permutations test—using FSTAT and GenAlEx. Pairwise genotypic differentiation at individual and population level was used to produce a PCoA in GenAlEx.

On basis of co-dominant alleles (F-statistics) and microsatellite repeat lengths (R-statistics), the genetic structure among sites (F_ST, R_ST), inbreeding within sites (F_IS, R_IS), overall inbreeding (F_IT, R_IT), and a pairwise genotypic differentiation matrix (F_ST) of each species were calculated via AMOVA-F_ST and AMOVA- R_ST at 999 random permutations (GenAlEx). This allowed further to test for an evolutionary signal (when R_ST > F_ST) and to roughly estimate overall connectivity levels as Nm = F_ST/(1-4F_ST) under the assumption of an island migration model, very likely to be violated. Therefore, specific hypotheses to estimate gene flow were tested with Migrate-n (Beerli, 2006; Beerli and Palczewski, 2010) from the mutation-scaled population sizes (Theta) and immigration rates (M). The Brownian model was tested locus by locus along with the product of all distributions of all loci. Uni- and bidirectional historical migration/expansion models were tested. Uniform prior distribution settings (min, max, delta) were as follows for Theta = 0.0, 20.0, 2 and for M = 0.0, 20, 2. The number of recorded steps was 10⁶ at a sampling frequency of 10³ after an initial burn-in. The effective number of immigrants per generation (N_em) was calculated as [Theta × M]/4. Specific hypotheses testing on directionality were considered in panmixia, source-sink and stepping-stone models for the migration between mangrove locations situated along the same coastline. More precisely, we considered six A. marina populations of the West coast to test the hypothesis of a northward or southward coastal current direction in the Strait of Malacca. Similarly, we considered three Avicennia populations from the East coast to test the hypothesis of prevailing southward coastal current direction of the South China Sea. The Brownian motion mutation model was adopted for randomly generated subsamples of 20 individuals in a transect, following the abovementioned settings, computing two replicate chains (with different seed), and using the Bezier thermodynamic integration (Beerli and Palczewski, 2010) for calculation of the Bayes factors from marginal likelihoods giving the model probabilities.

A Bayesian clustering analysis at individual level for three data sets was carried out in STRUCTURE v.2.3.4 (Pritchard et al., 2000) using an admixture model with correlated allele frequencies. We considered all 13 Avicennia populations (13 polymorphic microsatellite loci) and in further detail only the Eastern (11 polymorphic loci, excluding Avma14 and Avma2) and Western (11 polymorphic loci, excluding Avma14 and Aa22) populations separately. The model ran 10 iterations for each K value from 1 to 13; the burn-in period was 50,000 with 500,000 Markov chain Monte Carlo (MCMC) repeats. The optimal K was inferred with the ΔK statistic (Evanno et al., 2005) and LnPK using Structure Harvester (Earl and von Holdt, 2012) calculated with StructureSelector (Li and Liu, 2018).

The BARRIER v.2.2 software (Manni et al., 2004) was used to detect the location of sharp genetic changes between neighboring populations on basis of 10 pairwise F_ST matrices of every microsatellite locus, allowing a maximum of one barrier per matrix. We opted to calculate from superposition of raw data from F_ST matrices at locus level. The thickness of barrier lines thus will be based on the additivity of matrices accounting for the independent microsatellite loci that we consider as a preferred informative and valid method over bootstrapping a single mean F_ST matrix. A more detailed test of the West coast populations using pairwise F_ST and R_ST was performed at four distance classes of 127, 205, 388, 586 km (automatically generated when considering the option of a similar amount of pairwise comparisons) and 1,000 permutations with SPAGeDi 1.5a (Hardy and Vekemans, 2002), thereby allowing to test for an evolutionary signal (when R_ST > F_ST) specifically within each distance class. These distance classes represent threshold values instead of a full regression of a Mantel test. Linear-regression and log-regression over the full distance range were calculated with b-slope significance level (one-sided test; 1,000 permutations) for both F_ST and R_ST approaches.

The overall F_IJ kinship coefficient (Loiselle et al., 1995) for all pairs of individuals of within-site comparisons was tested for equal distance intervals in four classes (up to 25, 50, 75, and 100 m) by SPAGeDi 1.5a (Hardy and Vekemans, 2002). The FSGS i.e., the spatial autocorrelation of individual populations, were all tested for significance with 1,000 permutations. We computed the log-slope (-b) of linear regressions between pairwise kinship coefficients and geographical distance over restricted distance with 1,000 permutations. The potential influence of the open sea on transect FSGS was verified by plotting F_IJ kinship values of the first distance class (25 m) against the km distance from starting point of transect to open sea. We also checked the possible effect of resolution to detect FSGS, caused by the lower total number of alleles by plotting kinship values over 25 m distance against the total number of alleles within each population.

Results

Chloroplast Sequences and Haplotype Differentiation

The sequences of trnH-psbA from 43 samples resulted in three distinct haplotypes along Malay Peninsula (Figure 2; Supplementary Table 2). All populations from the East coast revealed a particular A. alba haplotype (451 bp) with a sequence similarly to A. alba from Vietnam and closely related to A. alba from The Philippines. This taxon, from the East Malay Peninsula, has identical morphological features resembling A. marina and thus it could be hypothesized as an introgressive A. marina with a captured chloroplast of A. alba following an earlier hybridization event. In contrast, populations on the West coast harbored two different trnH-psbA variants of A. marina that occurred geographically well-separated. Populations from the northern W1 to W3 sites contained a single haplotype (485 bp) that appeared to be more related to A. marina from Bangladesh than to those from Vietnam and China, whereas populations from the southern W4, W5, and SW sites had a unique haplotype (487 bp) clustered separately and closer to African A. marina than to the Asian samples. The two chloroplast sequences of A. marina along the West coast differed consistently in three mononucleotide repeats (T7 or T8; C7 or C10; T8 or T6), one transversion (A/C), and one transition (A/G). The captured A. alba haplotypes along the East coast showed 21 substitutions, three indels and two mononucleotide repeats when compared to A. marina. We could exclude the presence of A. rumphiana (458 bp) or A. officinalis (445 bp) haplotypes in the studied samples.

Figure 2. Maximum Likelihood tree of trnH-psbA intron with chloroplast identity of samples from East and West Malay Peninsula when compared to Avicennia marina from other countries and to other Avicennia species. Note the captured chloroplast of A. alba on the East Malay Peninsula and the two different haplotypes of A. marina along the West Malay Peninsula showing a break between W3 and W4.

Genetic Diversity of Nuclear Microsatellites

There are 54 alleles (on average 33 and 26-44 per transect) in 11 loci from the seven sampling sites along the West coast of Malay Peninsula (586 km stretch), with a mean (A_M) number of 2.6–4.4, effective (A_E) number of 1.6–2.0 and an adjusted richness (A_R) of 2.3–3.5 (Table 2). The overall observed heterozygosity (H_O = 0.343) was slightly lower than the expected heterozygosity (uH_E = 0.383). The within-population inbreeding (mean F_IS = 0.082) ranging from −0.061 to 0.290 was significant (p < 0.05) for four sites (Table 2).

AMOVA- F_ST results from the West coast revealed 37% genetic variation among A. marina populations, 5% among individuals and 58% within individuals, giving an estimate of F_ST = 0.369, F_IS = 0.082 and F_IT = 0.421 (Table 3). AMOVA- R_ST gave higher estimates of genetic variance among the populations than F_ST with R_ST = 0.572 (55%), hence an evolutionary signal was indicated by larger differentiation due to allele repeat length of the microsatellites (mutational steps) instead of allele identity only (Table 3). Considering the cpDNA differentiation between the northern (W1, W2, W3) and the southern (W4, W5, SW) sampling sites, a pairwise genetic differentiation of the nuclear microsatellites was highest for all populations between the northern (W1, W2, W3) and the southern (W4, W5, SW) sampling sites (F_ST = 0.356–0.541), whereas differences between the locations were more distinct for northern sites (F_ST = 0.120–0.301) than to the southern sites (F_ST = 0.010–0.071) (Supplementary Table 3). This indicates a strong connectivity among the southern A. marina populations with a very low or no connectivity with the northern ones. The PCoA fully separated these two geographical groups at the population as well as at the individual levels (Supplementary Figure 2).

Table 3. Summary of AMOVA with F-statistics and R-statistics of Avicennia marina (Western populations) and introgressive Avicennia marina (Eastern populations with captured cpDNA of Avicennia alba) mangrove fragments along coasts of the Malay Peninsula.

On the other hand, six sampling sites along the East coast of Malay Peninsula (80 km stretch) indicated 24 alleles (on average 20 and 18–22 per transect) in 11 loci, with a mean number (A_M) of 1.8–2.2, effective number (A_E) of 1.3–1.4 and an adjusted richness (A_R) of 1.8–2.1 (Table 2). The observed heterozygosity (H_O = 0.172) was slightly lower than the expected heterozygosity (uH_E = 0.204). The within population inbreeding (mean F_IS = 0.161) ranging from −0.088 to 0.232 was significant (p < 0.05) for four sites (Table 2).

AMOVA- F_ST results explained 9% genetic variation among the introgressive A. marina populations on the East coast, 15% among individuals and 76% within individuals, giving an estimate of F_ST = 0.089, F_IS = 0.161, and F_IT = 0.236 (Table 3). AMOVA- R_ST provided much lower estimates of genetic variance among these populations (4%) than F_ST, hence differences in allele repeat length of the microsatellites were non-explanatory (Table 3). Pairwise genetic differentiation was highest between transects E1 and E2 (F_ST = 0.151–0.207), whereas F_ST values between locations E2 and E3 remained low though significant (0.037–0.089) (Supplementary Table 4). All within-site comparisons of transects were non-significant (Supplementary Table 4). The PCoA showed a single cluster at individual level, but clustered transects separately at site level (Supplementary Figure 2).

Genetic Structure, Coastal Connectivity, and Gene Flow Models

The Bayesian clustering analysis with 486 A. marina samples indicated well-separated gene pools and assigned every individual to a single gene pool without any admixture (Figure 3). Delta K was high for K = 2 and K = 3. The latter outcome of K = 3 was supported by each iteration (low standard deviation), approached the LnP(K) plateau and fully corresponds to the previously defined haplotype groups. Also, the STRUCTURE analysis of each coastline revealed K = 1 for the Eastern (N = 179) and K = 2 (N = 307) for the Western populations (namely comprising a northern vs. a southern cluster on the West coast), hence we considered an overall structure of K = 3 clusters. The BARRIER analysis showed two major genetic breaks, one between East and West coasts of Malay Peninsula and another between the populations of W3 and W4 in the Strait of Malacca on the West coast (Figure 4). Implicit spatial analysis of the populations along the West coast using four distance classes revealed a significant pairwise lower F_ST and R_ST in the first distance class up to 127 km and a marginal significant b-log regression over the full distance of 586 km (one-sided test: F_ST r² = 0.47, p = 0.01; R_ST r² = 0.40, p = 0.02) (Supplementary Table 5).

Figure 3. STRUCTURE results of Avicennia marina populations showing a break between W3 and W4 along the West coast and for introgressive Avicennia marina populations reflecting a single gene pool on the East coastline of the Malay Peninsula.

Figure 4. Genetic barriers (red lines) computed from F_ST distance matrices of ten independent microsatellite loci using BARRIER analysis. The thickness of the red line reflects the importance of the barriers. Note the barrier between W3 and W4 in the Strait of Malacca in addition to a previously well-documented East-West barrier.

Specific testing with migrate-n on the directionality for A. marina across six mangrove locations of the West coast indicated that panmixia or bidirectional stepping-stone models (from South to North as well as from North to South) are less likely than a customized unidirectional South to North migration as well as a local bidirectionality between W3 and W4 sites (Table 4; Figure 5). Highest estimated gene flow values were observed from W3 toward W2 (N_em = 1.4) and from W5 toward W4 (N_em = 1.1). The lowest gene flow estimate for the bidirectional scenario of W4–W3 (N_em = 0.5) and an absence of gene flow for the reciprocal (N_em = 0.01) were all in agreement with the genetic break evidenced from BARRIER and STRUCTURE analyses. Migration of the introgressive A. marina across three mangrove locations on the East coast were supported by a unidirectional stepping-stone model (with gene flow estimates of N_em = 0.3–0.5), following the main ocean current directionality of the South China Sea (Table 4; Figure 5).

Table 4. Comparison of migration models on the directionality of (a) Avicennia marina populations for West Malay Peninsula sites and hypothesis of Northward Coastal current direction in the narrowest and shallow part of the Strait of Malacca; and (b) Introgressive Avicennia marina (with captured cpDNA of Avicennia alba) populations for East Malay Peninsula sites and hypothesis of prevailing Southward Coastal current direction in the South China Sea; with Nem included for each scenario.

Figure 5. Migration model with directions and gene flow estimates between populations (black rhombs). Strength of principal ocean currents of the study area are indicated in a gradient from blue and green (strong) to yellow and red (weak). Oceanic currents around Peninsula Malaysia in July 2017. Data from ESR (2009), OSCAR third degree resolution ocean surface currents. Ver. 1. PO.DAAC, CA, USA. Dataset accessed (2018-07-28) at http://dx.doi.org/10.5067/OSCAR-03D01.

Fine-Scale Genetic Structure

The overall F_IJ kinship coefficient for within-site comparisons revealed positive values for only three transects with the shorter distance classes of 25 or 50 m (Table 2), while the other transects lacked a fine-scale structure. The (F_IJ)₁ value of the first distance class (< 25 m) indicated a potential relationship with the distance to the open sea, because of an absence of FSGS (zero values) for mangroves along the West coast when situated closer to the coast than about 250 m (Supplementary Figure 3). However, absence of a detectable and significant FSGS could also be explained from lower allelic diversity and smaller sample size in transects on the East coast (Supplementary Figure 3).

Discussion

Evolutionary Significant Units and Hybridization

Avicennia marina, with a least concern (LC) status of the IUCN, has limited focus for conservation and management in many parts of the world including Malaysia. In Peninsular Malaysia, Rhizophora apiculata Blume and R. mucronata Lam. are more popular than Avicennia spp. due to their abundant seed bank, highest growth rate and commercial importance (Goessens et al., 2014; Tangah et al., 2020). However, the forestry measures such as planting, sustainable use, legal protection of areas, etc., are still applicable to A. marina in all countries wherever is needed (IUCN, 2021). Inbreeding and genetic diversity that are common to mangrove populations were also observed for A. marina (Maguire et al., 2000b; Dodd et al., 2002; Yan et al., 2016). Therefore, continued monitoring and in-depth research are necessary to provide further scientific insights on this population diversity, especially from the areas like Malay Peninsula where land blockage effect is protuberant (Azman et al., 2020; Wee et al., 2020; Fairuz-Fozi et al., 2021).

Though all sampled populations showed identical morphological features of A. marina, its actual abundance and distribution seems different along the East and West coasts of Malay Peninsula. Three evolutionary significant units could be confirmed from both nuclear and cytoplasmic markers. Avicennia marina along the 586 km stretch of West coast differentiated strongly from an introgressive A. marina featuring a captured chloroplast of A. alba along the 80 km stretch of East coast. Over and above this expected East-West division, an intra-regional subdivision was detected among A. marina populations in the narrowest region of the Strait of Malacca. This genetic break was clearly supported by the AMOVA, STRUCTURE, and BARRIER analyses whereas R_ST > F_ST indicates an evolutionary signal of this long-lasting divergence. Two different haplotypes along the West coast showed phylogeographic relationship with either a northern or a putative southern lineage, thereby assuming two Avicennia sources facing each other during the Holocene occupation with a prolonged separation in the Strait of Malacca.

The Malay Peninsula is expected to be a land barrier hindering gene flow between the East and West coast (Duke et al., 2002). The East-West genetic divergence across the Malay Peninsula has already been highlighted for several mangrove species such as S. alba (Wee et al., 2017; Yang et al., 2017) and A. alba (Wee et al., 2020), and to a lesser extent for B. gymnorhiza (Minobe et al., 2010; Wee et al., 2020). However, it has been shown that Malay Peninsula does not act the same depending on the species or genera studied. Indeed, Rhizophora mucronata and R. apiculata did not show the same population structures across Malay Peninsula: while R. apiculata showed some divergence between the two coasts, this was not the case for R. mucronata which showed higher dispersal potential (Inomata et al., 2009; Wee et al., 2020). The difference in genetic structure across species suggest that the Malay Peninsula act as a “filter” rather than a strict geographical barrier, depending on dispersal capacities of mangrove species and ocean currents (Wee et al., 2020). Avicennia spp. having a low long distance dispersal potential are more likely to be restricted to within their oceanic basin (Maguire et al., 2000b), as supported by our results and thus potentially indicating correlation between dispersal potential and genetic admixture after the LGM for Avicennia species.

Along the East coast an introgressive A. marina featuring a captured chloroplast of A. alba could be demonstrated because we could compare with trnH-psbA of pure A. alba samples from different locations in Vietnam and different islands in The Philippines. We mentioned one from each country in this article because the sequences were similar for the Bayes ML and the focus here was on Malay Peninsula diversity. Therefore, the cpDNA of the East coast samples belongs to A. alba and is not a misidentification of trnH-psbA. Since for the East coast samples, their nuclear genome largely amplifies the A. marina microsatellite loci and their cpDNA fully matches A. alba, then the chloroplast must be derived from an A. alba egg-cell in the case of maternal inheritance. Although pollination can go either way, the recipient mother plant in this case was definitely A. alba. Hybridization is not uncommon in mangroves. According to Ragavan et al. (2017), the natural hybrids in mangroves are predominantly found in seven genera namely, Acrostichum, Avicennia, Bruguiera, Ceriops, Lumnitzera, Rhizophora, and Sonneratia. In the case of Avicennia, two hybrids - one from Thailand with the parental species of A. marina and A. rumphiana and another from Brazil with the parental species of A. germinans (L.) Stearn and A. bicolor Standley were reported (Huang et al., 2014; Mori et al., 2015b). For Peninsular Malaysia, both East and West coasts are known to have the mangrove hybrids like Bruguiera x hainesii C. G. Rogers, Bruguiera x rhynchopetala (W. C. Ko) X.-J. Ge and N. C. Duke, Rhizophora x annamalayana Kathir., and Sonneratia x hainanensis W. C. Ko, E. Y. Chen and W. Y. Chen (Wan Juliana et al., 2014; Jamilah et al., 2015; Razali et al., 2016). In addition, presence of Rhizophora x lamarckii Montr., Rhizophora x mohanii P. Ragavan, Sonneratia x gulngai N. C. Duke, in the Malaysian mangroves was informed by Ragavan et al. (2017). We assume that field reports on possible cases of Avicennia hybridization were not available, because the leaf features of those trees that were found to be introgessed hybrids in our study resembled A. marina morphologically. Therefore, it is necessary to supplement nuclear DNA markers with chloroplast sequences in biogeographical regions where more than one Avicennia species occurs.

Connectivity Along Same Coastline

Migrate-n model testing supported northward unidirectional stepping-stone along the West coast, though directionality at the genetic break point (narrowest region of the Strait of Malacca) was unclear most likely due to weaker oceanic currents. Long-distance dispersal of A. marina propagules can be attributed to its buoyancy and viability characteristics as well as the hydrodynamics i.e., tidal inundation (Breitfuss et al., 2003), current speed and direction (Steinke and Ward, 2003). Wind may play a limited role in the dispersal of Avicennia propagules as they show low surface water contact and less influence by the drag force (Van der Stocken et al., 2015). Successful establishment and growth of the propagules after short to long-distance dispersal patterns (Rabinowitz, 1978; Clarke, 1993) becomes a responsible source for historically accumulated gene flow between populations (Duke et al., 1998). In this post-establishment phase, propagule predators are known to play a key role in the Americas, Africa and Asia (Smith, 1987; McGuinness, 1997; Dahdouh-Guebas et al., 1998; Langston et al., 2017). Once the crypto-viviparous propagules of A. marina released into seawater, they are expected to undergo different phases like floating with unshed pericarp, floating with shed pericarp, sunken with unshed pericarp, sunken with shed pericarp, etc., over a period of time (B. Satyanarayana, pers. comm.). If propagules do not reach to a suitable substratum immediately, then they may remain viable in seawater for several days to weeks (Steinke, 1986; Clarke and Myerscough, 1991; Clarke, 1993; Clarke et al., 2001) and would be able to disperse between geographically disjunct estuaries. Extended floating periods of several months have been observed for other Avicennia species, e.g., A. germinans (Rabinowitz, 1978; Alleman and Hester, 2011).

Our estimation of genetic connectivity among A. marina populations along the West coast of Malay Peninsula through comparison of migration models supported a putative unidirectional dispersal route that was congruent with prevailing ocean currents (Rizal et al., 2010, 2012; Haditiar et al., 2020). This further emphasizes the relevance of coastal connectivity to mangrove persistence (Van der Stocken et al., 2019a,b) as well as the importance of apparently discrete estuaries. The connectivity between adjacent mangroves was supported by the stepping-stone migration model. The dispersal of Avicennia propagules is likely restricted to a few tens of kilometers (e.g., Clarke, 1993; Duke et al., 1998; Melville and Burchett, 2002). Recently, Binks et al. (2018) found a long-distance dispersal up to 100 km, which is comparable to the 80 km stretch between northern and southern populations (E1–E3) of the East Malay Peninsula to demonstrate their genetic connectivity in a stepping-stone manner. The present findings also add to the emerging evidence of Avicennia species in general following an adjacent migration, such as the unidirectional way obtained for A. alba on the West coast (Wee et al., 2020), A. marina in eastern Africa (Triest et al., 2021c) or along Leyte island in The Philippines (Triest et al., 2021a) and bidirectional ways for A. germinans along each of the Caribbean and Pacific coasts of central America (Ochoa-Zavala et al., 2019). Although repeated bottlenecks or founder effects of the pioneering Avicennia species may have caused the differentiation of populations (Maguire et al., 2000a; Arnaud-Haond et al., 2006), we noticed only limited evidence of such conditions in the present study region.

The weaker currents slow down propagule dispersal and increase the risk of their entrapment among mangrove roots or settlement on unsuitable substratum such as a sand bar or beach. Low or absent tidal currents in the mangrove patches are responsible for leaving a trace of elevated kinship values. Such FSGS traces were also estimated within the spatial extents of a few meters to several hundreds of meters (Mori et al., 2015b; Do et al., 2019; Chablé Iuit et al., 2020; Triest et al., 2020; Triest and Van der Stocken, 2021). When the populations were sheltered (Triest and Van der Stocken, 2021) or severely fragmented and confined to artificial dikes (Hasan et al., 2018), these kinship values became enhanced. For instance, R. mangle showed a FSGS up to 90 m in different estuarine conditions of the Caribbean mangroves (Yucatan, Mexico), and even further up to 240 m when along a river (Chablé Iuit et al., 2020). In a high rainfall area of the Cameroon Estuary Complex, R. racemosa showed no or only limited autocorrelation within 25 m due to strong hydrodynamic situations (Ngeve et al., 2017). Nonetheless, one may not forget that the FSGS and mangrove populations diversity are also determined by the cumulative effect of insect, wind and bird pollination (Hermansen et al., 2014; Wee et al., 2015) and not solely by propagule dispersal patterns. The flowers of Avicennia are visited by numerous pollinators that largely comprise insects (i.e., honeybees – one of the frequent visitors), bats and birds (Clarke and Myerscough, 1991) that shows the possibility of cross-pollination as well as mating of siblings. The elevated levels of inbreeding in many sites should be explained also from a lack of pollen flow and from non-random mating, which requires a different design of study.

Conservation Relevance

Avicennia marina is largely distributed at the river mouths as a pioneer group of vegetation in Peninsular Malaysia (WIM, 2015). Besides the less human aided propagations for Avicennia spp., they are still considered as protective forest along the coastlines (Roslan and Nik Mohd Shah, 2014). Hence, these populations are mostly of natural origin and strongly support the genetic diversity observations in the present study. The clear-cut distinction of A. marina and introgressive A. marina populations on the East and the West coasts aid for species-level conservation and benefit researchers for future ecological genetic studies, without taxonomic errors. We call for the precautionary principle to not just treat the introgressive A. marina populations that captured a chloroplast of A. alba as Avicennia marina stands sensu stricto. The ecological justification for this resides in the concept of ‘cryptic ecological degradation', defined as 'functional ecological degradation that involves a qualitative decline of typical, stenotopic, vulnerable, valuable and functional species that is masked by a quantitative increase of less typical, eurytopic, disturbance-resistant, less valuable and less functional species. In a more general context, it is a qualitative ecological and socio-economic degradation of one ecosystem component that is masked by an easily detectable quantitative status quo or even increase of another component (Dahdouh-Guebas et al., 2005a,b; Dahdouh-Guebas et al., 2021). Cryptic ecological degradation has been reported in plant ecology (Koedam and Dahdouh-Guebas, 2008) and animal ecology (Bartolini et al., 2011). In the context of the present study the genetic introgression may imply yet another form of the concept of cryptic ecological degradation.

The present study also explains local FSGS with the farthest distance to the open sea in sheltered mangroves, although its absence in some areas could come from lack of distinctive power. The latter can be relevant for gene pools such as found in the introgressive A. marina on the East coast with an intrinsically low allelic and genetic diversity. Micro-evolutionary processes that involved introgressive hybridization with chloroplast capture of another species might leads to severe bottlenecks and pose constraints on the nuclear genome, allowing only a limited amount of variation to be transmitted for future generations. Low levels of diversity and a southward connectivity were also detected for Avicennia populations on this East coast. Irrespective of the morphological entity of A. marina in Malay Peninsula, there are three evolutionary significant units characterized by nuclear and chloroplast markers that merit conservation attention. The high likelihood of historical unidirectional stepping-stone migration between adjacent mangroves and the persistence of sharp genetic breaks, both add to the importance of considering coastal connectivity and conservation strategies of designing marine parks that should aim to preserve the natural cohesive forces on longer term.

Data Availability Statement

The datasets presented in this study can be found in online repositories. Genbank accession numbers of microsatellite loci and of trnH-psbA sequences are provided in Supplementary Tables 1, 2 of this article.

Author Contributions

LT: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (lead), resources (lead), visualization (supporting), writing-original draft (lead), and writing-review and editing (lead). BS: conceptualization (lead), data curation (supporting), formal analysis (supporting), funding acquisition (lead), investigation (lead), methodology (supporting), resources (lead), visualization (supporting), writing-original draft (lead), and writing-review and editing (supporting). OD and KS: conceptualization (supporting), data curation (supporting), formal analysis (supporting), funding acquisition (supporting), investigation (supporting), methodology (supporting), resources (supporting), visualization (supporting), writing-original draft (supporting), and writing-review and editing (supporting). TS: conceptualization (supporting), data curation (lead), formal analysis (supporting), funding acquisition (supporting), investigation (lead), methodology (lead), resources (supporting), visualization (supporting), writing-original draft (supporting), and writing-review and editing (supporting). FD-G: conceptualization (supporting), data curation (supporting), formal analysis (supporting), funding acquisition (lead), investigation (supporting), methodology (supporting), resources (supporting), visualization (supporting), writing-original draft (supporting), and writing-review and editing (supporting). All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Vrije Universiteit Brussel (VUB; grant number BAS42), and by the European Commission-funded Erasmus Mundus Masters Course in Tropical Biodiversity and Ecosystems (TROPIMUNDO). Also, the Higher Institution Centre of Excellence (HICoE) grant (vote #66928) awarded to the INOS, UMT has extended partial support to complete the fieldwork.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors sincerely thank EMMC TROPIMUNDO for facilitating scholarship holder's research. Authors are grateful to the Forestry Department for their kind permission to conduct this research. Mr. Mohammad Nazrin Ishak (INOS) extended full support to organize the fieldwork most successfully and collect the samples from each place.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcosc.2021.727819/full#supplementary-material

References