Impact of insular landscape features on the population genetics of a threatened climbing palm, Korthalsia rogersii Becc., endemic to the Andaman Islands

Sample collection and DNA extraction

All seven known populations of K. rogersii, including the newly identified ones, were surveyed during field expeditions in 2022 and 2024 in collaboration with the Andaman and Nicobar Forest Department (Field permit No.F.10 (G-I)/39/Vol.X/652 & No.F.10 (G-I)/39/Vol.XI/462). Leaf samples were collected from the encountered individuals across the populations. To minimize the risk of sampling clonal individuals, especially in small populations (Radhanagar, Betapur, Bakultala and Baratang), sampling from closely stand individuals was avoided. In larger populations (Interview islands, Havelock and ChidiyaTapu), a minimum distance of five m was maintained between sampled individuals. Sample sizes per population ranged from nine to 22, proportional to the population size. Species identification and selection of individuals were guided by detailed morphological descriptions of the species (Renuka & Vijayakumaran, 1995; Mathew et al., 2007). A total of 105 individuals were sampled (Table 1), and collected leaves were silica-dried in the field and subsequently transported to the Kerala Forest Research Institute for genetic analysis.

DNA was extracted using a modified cetyl trimethyl ammonium bromide (CTAB) protocol, which included 2% polyvinylpyrrolidone (PVP) in the extraction buffer, repeated washes with chloroform:isoamyl alcohol (24:1), and extended the incubation with isopropanol up to 30 min to enhance DNA yield and purity (Doyle & Doyle, 1987).

Population genetic analysis of K. rogersii

Microsatellite loci, originally developed and standardised for K. laciniosa (Dasgupta et al., 2021), were screened for cross-amplification in K. rogersii. Of these, 16 loci were cross-amplified; however, only seven loci exhibited clear, unambiguous amplification and polymorphism during an initial screen of 14 individuals (Table S1). These seven loci were selected for SSR genotyping across 105 individuals. Forward primers were labelled with 6-Carboxyfluorescein (6-FAM)/ Hexachlorofluorescein (HEX) and used for PCR amplification. Amplicons were analysed using an ABI 3730XL sequencer (Applied Biosystems) with Gene Scan 500LIZ (Applied Biosystems) as the internal size standard. Detailed PCR conditions, including reaction components and annealing temperatures, for each locus are provided in Table S1. Allele sizing and scoring were performed using GENEMAPPER software v4.0 (Applied Biosystems).

Genetic diversity analysis

The genotype data set was filtered for loci and individuals with >20% missing data. Duplicate genotypes were removed, and unique genotypes were retained for further analysis. MICRO-CHECKER was used to estimate the presence of null alleles and large allele dropout (Van Oosterhout et al., 2004). Chi-squared and exact tests were used to test Hardy–Weinberg equilibrium (HWE), as implemented in the Pegas package in R (Paradis, 2010). Summary statistics, including the number of alleles per locus per population, allelic richness (Ar), number of private alleles (Np), observed heterozygosity (Ho), expected heterozygosity (He) and fixation index (Fis) were estimated. Data filtering and summary statistics were carried out using the Adegenet package in R (Jombart, 2008).

Since the Fis is also influenced by null alleles or other deviations from HW, we also explored different inbreeding coefficients, (Li & Horvitz, 1953), as modified by Ritland (1996) (LH), and Lynch & Ritland (1999) (LR) using the related R package (Pew et al., 2015). We used the grouprel() function from the related R package to estimate the average pairwise relatedness within each population. To evaluate whether the observed within-group relatedness was greater than expected by chance, the function performed random permutations of individual group assignments, maintaining the original group sizes. Within-group relatedness was recalculated for each permutation, repeated 10,000 times. This generated a null distribution of expected relatedness values under random group membership. The observed within-group relatedness values are summarized in Table 1.

Population genetic structure

Pairwise population differentiation (Fst; Weir & Cockerham, 1984) was estimated using the diffCalc function in the diveRsity R package with 1,000 bootstrap replicates (Keenan et al., 2013). Individuals were resampled with replacement within populations to generate confidence intervals, and differentiation was considered significant when the lower bound of the 95% interval exceeded zero.

Population genetic structure was further analysed using discriminant analysis of principal components (DAPC) implemented in the Adegenet package in R (Jombart, 2008). To determine the optimal number of principal components (PCs) to retain, cross-validation was carried out in DAPC using xvalDapc function. During cross-validation, the dataset was partitioned into a training set (90%) and a validation set (10%) using stratified random sampling to ensure representation from all populations in both datasets. The DAPC analysis was performed on the training set using varying numbers of retained PCs, and the optimum number of PCs was determined based on the assignment success of individuals in the validation set. The number of PCs that resulted in the lowest mean squared error (MSE) was considered optimal. This optimal configuration was then selected to generate a scatterplot based on the first and second linear discriminants of DAPC, providing a visual representation of genetic clustering among populations.

Bayesian analysis of population genetic structure and admixture patterns was performed using STRUCTURE with a burn-in period of 10⁵ iterations followed by 10⁶ Markov Chain Monte Carlo (MCMC) replications (Pritchard, Stephens & Donnelly, 2000; Hubisz et al., 2009). For each K value (number of genetic clusters), ranging from 1 to 10, 25 replications were performed to ensure consistency of results. STRUCTURE analyses were carried out with and without the LOCPRIOR model. In the LOCPRIOR model, sampling locations were included as prior information to assist clustering. This approach is particularly useful in the case of data sets that contain relatively few markers, small sample sizes, or very weak signals of population structure that may not be captured by the standard STRUCTURE model (Hubisz et al., 2009). To determine the optimal number of genetic clusters (K), we used Structure Selector, which integrates multiple estimation methods (Li & Liu, 2018). Given the tendency of the widely used ΔK statistic to favour K = 2, we also considered alternative metrics such as Ln Pr(X—K), MedMedK, MedMeanK, MaxMedK, and MaxMeanK, as implemented in Structure Selector (Janes et al., 2017).

Isolation by distance (IBD) and environment (IBE)

To better understand the factors shaping genetic differentiation in Korthalsia rogersii, we evaluated the influence of geographic and environmental distances among populations. Among the Great Andaman populations (Radhanagar, Betapur, Bakultala, Baratang, and Chidiya Tapu), both geographic distance and environmental variation, particularly latitudinal climatic gradients, are likely contributing to the observed patterns of genetic structure.

To assess the influence of geographical distance and environmental variation on genetic differentiation among individuals in spatially separated populations, Mantel and partial Mantel tests were performed. Pairwise genetic distances among individuals were calculated using dist.gene function in the ape R package, while a pairwise geographical distances matrix was constructed using geodist. The environmental layer of nineteen bioclimatic variables was sourced from CHELSA Version 2.1 (Karger et al., 2017). Environmental values at each sampling point were extracted using raster package in R, and Euclidean distances between sites were calculated using dist () function. The Mantel test was used to assess isolation by distance (IBD) and isolation by environment (IBE) using geographical distance and environmental distance as predictors, respectively, for IBD and IBE, and genetic diversity as the response variable. For the partial Mantel test, one predictor variable (either geographic or environmental distance) was kept as a covariate to evaluate the independent effect of the other on genetic distance. Both Mantel and partial Mantel tests were carried out using the vegan R package, with statistical significance assessed through 10,000 permutations (Oksanen et al., 2015).

Bottleneck test

Recent bottleneck events were evaluated using BOTTLENECK v 1.2.02 (Piry, Luikart & Cornuet, 1999). An increase in the heterozygosity (He) over the heterozygosity expected (Heq) at the mutation drift equilibrium was considered an indication of a recent bottleneck event. Given that fewer than 20 loci were used in the analysis, Wilcoxon’s signed-rank test was applied, which is recommended for studies with a limited number of loci. The test was run with 10⁴ iterations under both the stepwise mutation model (SMM) and the two-phase model (TPM; and 95% SMM) for the analysis.

Ecological niche modelling

GPS coordinates and altitudes were documented for each sampled individual during field surveys using a handheld GPS device (eTrex 30; Garmin GmbH, Garching, Germany). In total, 69 occurrence points were recorded. Additional occurrence points were obtained from herbarium collections at the Kerala Forest Research Institute (international acronym ‘KFRI’). Nineteen bio-climatic variables for the period 1981–2010 were retrieved from the CHELSA Version 2.1 with spatial resolutions of 30 arcsec (Karger et al., 2017). Bioclimate layers were trimmed to the distribution area of each species using the drawExtent() and corp () functions in the raster R package. To reduce spatial autocorrelation, occurrence points sharing the same grid (one km²) were thinned using R packages spThin, and 15 spatially distinct points retained were used for modelling (Aiello-Lammens et al., 2015). Multicollinearity among the 19 environmental layers was assessed using the Variance Inflation Factor (VIF), and highly correlated layers were excluded from further analysis. Remaining independent layers were retained, and Ecological Niche modelling was carried out using the SDM R package (Naimi & Araújo, 2016). Initially, we employed the presence-background algorithm (Maxent) to model the potential distribution of K. rogersii. However, the resulting model exhibited poor prediction performance, with a low area under the curve (AUC = 0.5). Consequently, alternative presence-absence algorithms, including generalised linear models (GLM), boosted regression trees (BRT), and random forests (RF), were applied using imputed pseudo-absence points. Despite this, none of the models achieved satisfactory predictive accuracy (AUC < 0.8). Among them, the Maxent model performed relatively better, although still below the threshold for reliable interpretation. Therefore, ENM outputs were not used for further inference but are included as a (Fig. S1) for reference.

The lower prediction accuracy is likely due to a combination of factors: the species’ highly restricted distribution, small number of spatially independent occurrence points (n = 15), and potential microhabitat specificity that is not captured by coarse-scale bioclimatic variables. Rare or patchily distributed species often produce limited occurrence data, which reduces model reliability and increases uncertainty in predictions (Moudrý et al., 2024).

Population genetics

Microchecker analysis identified a heterozygote deficit in two microsatellite loci, Kle23 and Kle18, indicating null alleles. However, since the null alleles were not consistently detected across all populations, these loci were retained for further analyses. The observed heterozygote deficit appeared to be driven more by population genetic phenomena than by the presence of null alleles. Hardy-Weinberg Equilibrium (HWE) tests, using both Chi-squared and Monte Carlo permutation methods, revealed that five of the seven loci deviated from HWE, indicating population subdivision and/or restricted gene flow (Table S2). However, in both Chi-squared and Monte Carlo permutation tests, no loci consistently deviated from HWE across all populations, nor did any population consistently deviate across all loci (TableS S2A & S2B). Genetic diversity parameters, including the number of private alleles (Np), allelic richness (Ar), observed heterozygosity (Ho), expected heterozygosity (He) and inbreeding coefficient (Fis), are provided in Table 1. Expected heterozygosity and allelic richness were generally comparable across populations, irrespective of sample sizes, except in Bakultala, which showed notably reduced allelic richness (Table 1). In contrast, the number of observed alleles was positively correlated with sample size, likely reflecting underlying differences in population sizes (Fig. S2). The highest number of private alleles (alleles secluded to a single population) was observed in ChidiyaTapu (13), followed by Bakultala (six). None of the populations exhibited significant Fis values, indicating the absence of high inbreeding or deviation from HWE.

Estimates of individual inbreeding coefficients (LH and LR) revealed considerable variation, with LH values ranging from −0.43 to 0.79 and LR from −0.38 to 0.63. Within-population relatedness analysis revealed significantly higher kinship among individuals in the Bakultala population compared to expectations under a random mating scenario (p < 0.0158), suggesting restricted mating and potential inbreeding within this group. In contrast, most other populations did not show significant deviations from the expected within-group relatedness, indicating more random mating patterns and a lower likelihood of close kin associations. Overall relatedness across all populations was also significantly high (p < 0.0046), indicating the presence of population structure. Interestingly, the ChidiyaTapu population exhibited a near-zero within-group relatedness (−0.0069329) and a relatively high number of private alleles.

Genetic Bottleneck

BOTTLENECK analysis using one-tailed Wilcoxon test (p < 0.05) assesses heterozygous excess as an indicator of recent bottleneck events. In this study, none of the populations showed significant evidence of a recent bottleneck event (Table S3). However, the ChidiyaTapu population showed significant deviations from mutation-drift equilibrium in the two-tailed Wilcoxon test, indicating the possible influence of other demographic events.

Population genetic structuring

Pairwise Fst values ranged from 0.005 (between Betapur and ChidiyaTapu) to 0.165 (between Havelock and Radhanagar) (Table 2). Among the studied populations, Radhanagar, Havelock and Interview Island consistently exhibited moderate to high Fst values in comparisons with other populations, indicating a relatively high level of genetic differentiation and suggesting their genetic distinctiveness. In contrast, Betapur exhibited weak or non-significant differentiation with most populations, except Interview Island. Likewise, Chidiya Tapu showed largely non-significant differentiation with central and southern populations, except with Havelock, where significant differentiation was observed.

Discriminant Analysis of Principal Components (DAPC) revealed clear genetic distinctiveness of the interview Island and Radhanagr populations. The Havelock population, though overlapping with the Betapur population, also showed genetic distinctness to a certain extent. In contrast, populations from the Great Andaman Island (Betapur, Bakultala, Baratang and Chidiyatapu) clustered closely together, except for Radhanagar, which remained genetically distinct (Fig. 2).

STRUCTURE analysis without LOCPRIOR model failed to detect distinct genetic clusters, suggesting weak overall genetic structure and a largely panmictic population across the archipelago. This lack of distinct clustering identified may be partly attributed to the limited number of microsatellite loci employed (Hubisz et al., 2009). However, when the sampling location was incorporated as prior information using the LOCPRIOR model, a more defined genetic structure emerged. Based on the ΔK method (Evanno, Regnaut & Goudet, 2005), K = 2 was identified as the most likely number of clusters (Fig. 2). This analysis revealed distinct admixture patterns, particularly in the Interview Island and Radhanagar populations, consistent with the findings from the DAPC analysis.

At higher K values (K = 5), supported by the mean log-likelihood (LnP(K)) and additional statistical estimators (MedMed, MaxMed, MaxMean), structure analysis revealed finer-scale population structure consistent with the geographic separation among populations. At K = 5, Radhanagar, Interview Island, Havelock, Bakultala, and Chidiyatapu populations exhibited distinct genetic structure with unique admixture patterns (Fig. 2). In contrast, the Baratang and Beetapur populations showed varying degrees of admixture from multiple clusters, suggesting shared ancestry with genetically distinct groups. This pattern of mixed ancestry is likely influenced by sample size or increased model complexity at higher K.

Among all populations, Interview Island emerged as the most genetically distinct population, showing separation as early as K = 2. Havelock Island and Bakultala populations also demonstrated pronounced genetic structuring in the admixture plot at K = 4, indicating their genetic discreteness. While most populations from Great Andaman Island (Radhanagar, Betapur, Baratang & Chidiya Tapu), except the Bakulta, appeared admixed at lower K values, a clear, distinct cluster was evident at K = 5, particularly in the Radhanagar and Chidiya Tapu populations.

Isolation by distance and environment

Along with the strong influence of restricted gene flow due to the water barriers, both Mantel and partial Mantel tests indicated that environmental factors also play a role in shaping population genetic differentiation. Among the predictor variables, environmental variation showed a weak but significant correlation with the genetic distance (r = 0.221, P < 0.001), while geographical distance exhibited no significant correlation (r = 0.0376, P = 0.1). These results suggest that environmental (climate) factors, although modest, may play a role in shaping the genetic structure of K. rogersii.

Genetic diversity and population bottlenecks

Populations of Korthalsia rogersii exhibited moderate levels of genetic diversity, with expected heterozygosity (He) ranging from 0.59 to 0.66. Despite its restricted distribution and small population sizes, the species has retained considerable genetic variation. While population size is typically a major determinant of genetic diversity, He in this study showed no strong correlation with population size, a pattern also observed in other long-lived plant species (Kramer et al., 2008). Moreover, He alone is not considered a reliable predictor of population viability or extinction risk (Schmidt et al., 2023). Similar trends have been reported in other threatened palms, where substantial genetic diversity persists despite conservation concerns (Shapcott et al., 2007; Asmussen-Lange, Maunder & Fay, 2011; Shapcott et al., 2012b).

In contrast, the number of observed alleles, a more sensitive indicator of recent demographic changes, was positively correlated with population size (Allendorf et al., 2022b), suggesting that smaller populations of K. rogersii may be undergoing partial genetic erosion. Allelic richness (Ar), which reflects the impacts of population size, fragmentation, and habitat degradation (Browne, Ottewell & Karubian, 2015; Carvalho et al., 2015; da et al., 2017; Soares et al., 2019), ranged from 3.842 in Bakultala to 4.597 in Chidiya Tapu. However, Ar did not consistently align with population size. Notably, Bakultala, the population with the lowest Ar, exhibited significantly elevated within-population relatedness, suggesting reduced genetic variability. In contrast, populations such as Chidiya Tapu, which had the highest Ar, showed no evidence of significant relatedness.

Despite this high relatedness, Bakultala had a low fixation index (Fis), potentially reflecting a recent bottleneck in which reduced numbers of mating individuals increased relatedness without causing deviations from the Hardy-Weinberg equilibrium. Individual-level inbreeding coefficients showed wide variation across populations (LH: −0.43 to 0.79; LR: −0.38 to 0.63), with no consistent pattern of inbreeding. This range indicates substantial heterogeneity in mating behaviour. Given that K. rogersii is monoecious and the presence of self-incompatibility mechanisms is unknown, occasional selfing is likely and may explain the high inbreeding coefficients (>0.25) observed in certain individuals, even though population-level inbreeding was generally weak.

Private alleles can indicate restricted gene flow or local adaptation (Slatkin, 1985). Their distribution varied across populations and did not correlate with population size. Chidiya Tapu had the highest number of private alleles, although at low frequencies. Despite low pairwise genetic distances with nearby populations, the presence of numerous private alleles suggests Chidiya Tapu may be a relic population retaining ancestral genetic variation. This may be linked to the relatively stable evergreen habitat in southern Andaman and reflects patterns observed in glacial refugia and ancestral populations of other species (Comps et al., 2001; Márquez-Márquez et al., 2025). Interestingly, such elevated private allele richness in the southern population was not observed in the Andaman day gecko, a species with different habitat preferences and life history traits (Mohan et al., 2020).

Bakultala, despite its small size and increased relatedness, harboured the second-highest number of private alleles, indicating both isolation and unique genetic variation. Surprisingly, the geographically isolated populations of Interview and Havelock Islands did not exhibit particularly high private allele richness. This may be explained by historical land connections during the Last Glacial Maximum (Smith & Sandwell, 1997), which likely facilitated past gene flow and reduced long-term isolation.

Population genetic differentiation, admixture and influence of IBD and IBE

Although there is no universally accepted threshold for interpreting Fst values, values <0.05 indicate low, 0.05–0.15 as moderate, 0.15–0.25 as high, and >0.25 as very high genetic differentiation (Hartl & Clark, 1997). Based on these, Korthalsia rogersii populations exhibited a range from low to high pairwise genetic differentiation. Notably, Interview Island and Radhanagar consistently displayed moderate to high Fst values in comparison with all other populations, suggesting substantial genetic distinctiveness. In contrast, the central population (Betapur) exhibited weak or non-significant differentiation from most other sites, suggesting higher levels of gene flow. However, Betapur showed significant differentiation from Interview Island, despite its geographic proximity, likely due to isolation by water. Similarly, Chidiya Tapu was largely undifferentiated from nearby central and southern populations, except for Havelock, where significant differentiation was observed. These patterns underscore the strong influence of intervening water channels in shaping the population structure of K. rogersii, a conclusion further supported by both DAPC and STRUCTURE analyses.

Although Chidiya Tapu harbours the highest number of private alleles, it consistently clustered with Middle and South Andaman populations in both DAPC and lower K values of STRUCTURE. This suggests a largely shared gene pool of the Chidiya Tapu population. The presence of rare alleles in the absence of marked genetic divergence supports the interpretation that Chidiya Tapu may represent a relic or ancestral population, retain unique lost genetic variants, while remaining part of a broader gene pool (Slatkin, 1985; Excoffier & Ray, 2008).

Most South Andaman populations clustered closely together, although Havelock Island appeared somewhat distinct in both DAPC and STRUCTURE analyses, likely a result of its geographic isolation. Interestingly, the Bakultala population exhibited a unique admixture profile despite its low pairwise Fst with the nearby Betapur population. This distinct pattern may be attributed more to high within-population relatedness than to true genetic divergence. Similarly, Chidiya Tapu only formed a distinct cluster at higher K values in STRUCTURE, reflecting shallow differentiation likely driven by its high frequency of private alleles (Slatkin, 1985).

Overall, the genetic structure of K. rogersii appears to be shaped primarily by restricted gene flow, particularly the water barriers isolating populations on Interview Island and Havelock. In addition to physical isolation, latitudinal climatic variation also seems to contribute to the observed differentiation, as indicated by Mantel test results. Two major genetic clusters were detected within Great Andaman Island, a northern cluster (e.g., Radhanagar) and a southern one (e.g., Chidiya Tapu), with signs of genetic admixture in the central populations (e.g., Betapur and Baratang). Populations on the Interview and Havelock Islands remained distinct from these mainland clusters, likely due to prolonged geographic isolation.

Climatic gradients may further reinforce this structuring; the northern Andaman experiences more pronounced precipitation seasonality than the south, potentially influencing local adaptation. Comparable patterns of north–south genetic structuring have been documented in other insular species, such as the Andaman day gecko (Phelsuma andamanensis), which forms two natural genetic clusters (Mohan et al., 2020). By contrast, the Andaman keelback (Xenochrophis tytleri) shows stronger longitudinal differentiation, driven mainly by saltwater barriers (Mohan, Swamy & Shanker, 2018). These differences suggest that the direction and strength of genetic structuring depend on species-specific life history traits and dispersal capacities, although a general latitudinal trend is evident across the landscape. Limited knowledge of the pollination and seed dispersal mechanisms in Korthalsia rogersii, however, hinders a more comprehensive interpretation of the observed genetic patterns. Although empirical data on seed dispersal are unavailable, birds are likely the primary seed dispersers, as observed in many other rattan species. Similarly, pollen dispersal is presumably insect-mediated (Renuka, Indira & Muralidharan, 1998). Consequently, gene flow in this species is likely influenced by the movement and behaviour of these biotic vectors.

To comprehensively understand the processes shaping genetic variation in the Andaman archipelago, studies from multiple taxa are essential, along with insights into species ecology, including pollination and seed dispersal mechanisms. Future research should integrate high-resolution genomic data, extensive spatial sampling, and environmental variables to capture fine-scale patterns of phylogeographic structure and landscape-level genetic connectivity in this unique insular ecosystem.