Aim: To evaluate whether Pleistocene sea-level changes can explain the timing and demographic processes of population divergence between continental and insular plant populations on a land-bridge island.

Location: Alcatrazes Island and adjacent mainland areas along the southeastern coast of Brazil.

Taxon: Three bird-dispersed plant species (Chiococca alba, Eugenia sulcata, and Schinus terebinthifolia)

Methods: We used multi-species population genomic data (RAD-Seq) to test alternative demographic scenarios of population divergence. Approximate Bayesian Computation and Supervised Machine Learning were applied to evaluate divergence timing relative to the Last Glacial Maximum (LGM) and to assess the presence or absence of gene flow between mainland and insular populations. Genetic structure was inferred to identify population clustering within each species.

Results: Genetic structure analyses consistently identified three populations per species, corresponding to one insular and two mainland populations. For all species, demographic models supporting post-LGM isolation were strongly favoured, and no detectable genomic signal of recent gene flow was recovered under the selected models. These results are consistent with postglacial sea-level rise contributing to population isolation and divergence.

Main conclusions: Our comparative genomic analyses across multiple species support the hypothesis that Quaternary sea-level rise contributed to population divergence on land-bridge islands. Recurrent cycles of connectivity and isolation between islands and the mainland may have left a strong imprint on the demographic history of bird-dispersed forest plants.

Population sampling and data collection

To enable comparative investigation of island–mainland divergence, we prioritized unrelated taxa that (1) occur in both insular and adjacent mainland environments, (2) are sufficiently abundant for population-level genomic sampling, and (3) share bird-mediated seed dispersal, thereby minimizing major differences in dispersal syndrome among taxa. All three focal species are widely distributed along the coastal Atlantic Forest of southeastern Brazil, occurring in restinga and adjacent forest formations across a broad regional range. Based on these criteria, we selected Chiococca alba (L.) Hitchc. (Rubiaceae), Eugenia sulcata Spring ex Mart. (Myrtaceae), and Schinus terebinthifolia Raddi (Anacardiaceae). While C. alba and S. terebinthifolia are found predominantly on restinga (Brazilian coastal vegetation) slopes, near dunes on the Brazilian coast, E. sulcata is associated with restinga forests. Regarding life forms, C. alba occurs as a shrub or subshrub, either erect or climbing (Barbosa, 2020; Motley et al., 2005). E. sulcata consists mostly of trees (Zamith & Scarano, 2006), while S. terebinthifolia ranges from shrubs to trees (Nickerson & Flory, 2015).

Because previous information on phylogeographic structure and population relationships was unavailable for the focal species, we prioritized sampling mainland populations geographically closest to Alcatrazes Island to investigate island–mainland demographic differentiation and connectivity, following approaches commonly adopted in comparative insular systems (e.g., Biddick et al., 2019; Benítez-López et al., 2021). This sampling design targets mainland populations proximal to the island and is not intended to capture the full geographic genetic structure of each species across its broader distribution in the Atlantic Forest. In addition to collections conducted on Alcatrazes Island, mainland populations were sampled in the coastal municipalities of Bertioga, Ubatuba, and São Sebastião, São Paulo State. We collected eight to ten individuals per species at each locality. Leaf material was preserved in RNAlater™ and stored at −20 °C until DNA extraction.

DNA extraction and RAD sequencing

Genomic DNA was extracted from leaf material of a total of 86 individuals of the three species following the Sorbitol-CTAB protocol combined with the DNeasy® Plant Mini Kit Qiagen protocol (Qiagen, Finland). The extracted DNA was sent to Floragenex Inc. (Eugene, Oregon) for RAD-seq library preparation using the restriction enzyme SbfI (Baird et al., 2008; Etter et al., 2012). Fragments were ligated to specific barcodes and adapters and subsequently sequenced in a single lane on the Illumina NovaSeq-6000 Sp300 platform, 2 × 150 bp.

Data quality control, filtering, and SNP Calling

Quality verification of the raw reads generated by RAD-seq was performed using FastQC v.0.12.1 (Andrews, 2010). The raw data obtained for the three species were assembled de novo separately using Ipyrad 0.9.93 (Eaton & Overcast, 2020), with the default configuration, except for max_Indels_locus, set to a maximum of 2 indels per locus. At this stage, low-quality reads were filtered out, and specific sequencing adapters were removed. We used VCFtools v.0.1.17 (Danecek et al., 2011) for biallelic loci filtering, indel removal, and missing data configuration. To recover a larger number of SNPs for C. alba, which had lower sequencing quality compared to the other two species, loci with a maximum of 20% missing data were filtered, while for E. sulcata and S. terebinthifolia, which performed better in sequencing, loci with a maximum of 10% missing data were retained. Since most analyses require unlinked SNPs, we used a custom Python script to randomly select one SNP per locus (vcf_single_snp.py; https://github.com/pimbongaerts/radseq).

Population structure and genetic diversity

To detect patterns of genetic structuring within each species, we implemented three analyses. The first was a Discriminant Analysis of Principal Components (DAPC; Jombart et al. (2010)), implemented in R with the package ‘poppr’ (Kamvar et al., 2014). This analysis uses the k-means clustering algorithm and the BIC (i.e., Bayesian Information Criterion) to identify the optimal number of clusters (k). The second analysis was estimated using a program based on sparse non-negative matrix factorization (sNMF; Frichot et al. (2014)), implemented in R through the package ‘LEA’ (Frichot & François, 2015). Ancestry coefficients were estimated for 1–10 ancestral populations (K) using 100 replicates for each K. To determine the optimal number of ancestral clusters, the cross-entropy criterion was used. The results were then visualized through a biplot, showing, for each individual, the relative contribution of each ancestral cluster. The third structure analysis was based on a Bayesian clustering model, implemented in STRUCTURE 2.3.4 (Pritchard et al., 2000), a method that allows detecting the most likely number of genetic clusters (K) present in mainland and insular samples. The ‘admixture’ model was implemented and applied to 10 independent runs, varying K from 2–4. Each run consisted of 500,000 iterations with the Markov Chain Monte Carlo (MCMC) method and a burn-in phase of 100,000 iterations. STRUCTURE HARVESTER 0.7 (Earl & vonHoldt, 2012) was then used to infer the ΔK statistic suggested by Evanno et al. (2005). The results of the 10 runs were summarized using CLUMPP 1.1 (Jakobsson & Rosenberg, 2007), and the visualization of the results was processed using DISTRUCT 1.1 (Rosenberg, 2004).

With the concordance among the results of the structure analyses, indicating three populations for each species, we calculated allelic richness (AR), the proportion of alleles (%A), observed heterozygosity (HO) and expected heterozygosity (HE), and Wright’s F-statistics, such as the inbreeding coefficient (FIS), together with their respective confidence intervals, for each population inferred from 1,000 bootstrap resamplings, as well as genetic divergence (FST). These calculations were performed using the R packages ‘diveRsity’ (Keenan et al., 2013), ‘PopGenKit’ (Rioux Paquette, 2011), and ‘hierfstat’ (Goudet, 2005). These parameters were subsequently compared among populations within each species and across different species.

Inference of species tree topology

To construct a coalescent tree showing the relationship among populations of each species, we used SNAPP v.1.2.5 (Bryant et al., 2012), implemented in BEAST2 v2.7.6, which allows the inference of topology directly from unlinked biallelic markers (Bouckaert et al., 2014; Bryant et al., 2012). The filtered VCF files were converted into Nexus format using a custom Python script (vcf2snapp.py; https://github.com/CoBiG2/RAD_Tools), and the XML input files for BEAST2 were created using BEAUti v.2.7.6 (Bouckaert et al., 2014). We used a chain length of 1 million generations, sampling parameters every 1,000 generations. Convergence was monitored using ESS values from the log files in Tracer v.1.6 (Rambaut et al., 2018), removing the first 10% of trees as burn-in. To construct the “Maximum clade credibility” trees and annotate them with their posterior probabilities, we used TreeAnnotator v1.10.4 (Drummond et al., 2012), making it possible to visualize the results with DensiTree (Bouckaert & Heled, 2014).

Demographic models

We tested three models to explore the demographic history of mainland and insular populations of the three plant species from the southeastern coast of Brazil. The hypotheses for model construction followed the scenarios proposed by Gehara et al. (2017), with adaptations for models including three populations and with adjustments in population size to ensure a better fit between simulated and observed data. For all models, a prior of population size change was set between 0–25,000 years, covering the entire period during and after the LGM. However, no specific condition of population size change was imposed, allowing the program to estimate changes based on summary statistics, without assuming a priori expansion, contraction, or constancy in population size. For all scenarios, the population size prior ranged from 1,000 to 200,000. Moreover, a condition for divergence times was established considering prior knowledge of the topology generated by SNAPP. For C. alba and E. sulcata, we assumed that the divergence between mainland populations occurred more recently than the separation between insular and mainland populations. For S. terebinthifolia, on the other hand, we assumed that the  nsular population ddivergedmore recently from the Bertioga population (the closest to the island), while the Ubatuba population shows an older divergence.

The PreLGM_Is model represents a divergence occurring before the LGM, without migration. The PreLGM_IM model shares similar characteristics but includes constant migration until the LGM. After this period, migration between insular and mainland populations ceases, remaining only between the two mainland populations, with the possibility of changes in the migration rate between them after the LGM. The prior migration between populations was set to range from 0.001–1 migrant per generation. These models were built based on the geological formation of the island, described as an inselberg-like environment (Sabino et al., 2025a). The rationale of these models stems from the hypothesis that insular populations were already genetically separated from mainland populations before sea-level fluctuations in the Late Pleistocene. Because they were isolated in an inselberg-like environment, divergence is thought to have occurred at some point during the Pleistocene, before the LGM. Both models assume a divergence prior to the LGM, with a minimum interval of 20,000 years and a maximum of 2,500,000 years.

The PostLGM_Is model represents a divergence occurring after the LGM, consistent with sea-level variation at the end of the Late Pleistocene and the beginning of the Holocene. In this scenario, divergence priors were set to range from 0–20,000 years. Additionally, population size changes were expected to occur close to the divergence between populations. The rationale of this model derives from the hypothesis that insular and mainland populations had land connections on the continental shelf of southeastern Brazil and, therefore, variations prior to the LGM would not be detected by the models.

Model testing using ABC and SML

To evaluate the alternative models for the three species, we used the ‘.alleles’ file with phased loci, one of the outputs from iPyrad. These files were converted into FASTA format using the function iPyrad.alleles.loci2fasta, implemented in the package PipeMaster v.0.2.3 (Gehara et al., 2020; https://github.com/gehara/PipeMaster) in R. Using the same package, we calculated the following summary statistics for the observed data for the three species: segregating sites, nucleotide diversity (π), Watterson’s theta, Tajima’s D, fixation index (FST), pairwise FST, % shared polymorphisms for subpopulations, % private polymorphisms for subpopulations, and % fixed polymorphisms for subpopulations. The mean and variance of each summary statistic were calculated across all loci for each population.

Genomic data were simulated under the three demographic models (PreLGM_Is, PreLGM_IM, and PostLGM_Is) using PipeMaster, which relies on a modified version of msABC (Hudson, 2002; Pavlidis et al., 2010). For each species, a total of 10,000 simulations were performed for each demographic scenario in order to carry out model selection analyses that best fit the observed data, using both ABC (Approximate Bayesian Computation) and SML (Supervised Machine Learning). Model fit was assessed by summarizing the simulated datasets with PCA and comparing them to the observed data using the function PlotPCs implemented in PipeMaster. With the same package, posterior probabilities of each model were calculated using the function postpr, which applies a rejection algorithm with a tolerance of 0.0001, retaining only the 100 simulations closest to the observed data. Model identifiability and accuracy were further evaluated through cross-validation using the function cv4postpr (Csilléry et al., 2012), also based on the rejection algorithm.

For model selection using SML, we used the Keras interface (Chollet, 2015) implemented in R. We split the data into 80% for training and 20% for testing an artificial neural network. The network architecture consisted of three hidden layers with 32 nodes each, followed by an output layer with 3 units. The hidden layers used the ReLU activation function, and the output layer used softmax. The model was trained with the Adam optimizer and the loss function sparse_categorical_crossentropy. Training was performed for 1,000 epochs with a batch size of 10,000 and a validation split of 0.1.

Parameter estimation using ABC and SML

To estimate the posterior distribution of model parameters for each species, a total of 100,000 simulations were performed for the model with the highest posterior probability. For this purpose, we applied ABC, assessing the reliability of the results with cv4abc using the rejection algorithm with a tolerance of 0.0001, and employing the same summary statistics used for model selection. In addition to ABC, we also estimated the posterior distribution of parameters using SML. For this, we implemented another artificial neural network with Keras (Chollet, 2015), consisting of three hidden layers with 32 nodes each and ReLU activation functions, similar to the architecture used for model selection. In this case, however, the output layer contained a single node with a linear activation function. RMSprop was used as the optimizer, with the loss function set to mean squared error (MSE). We used 88% of the data for training and 12% for testing. The neural network was trained for 1,000 epochs with a batch size of 10,000 and a validation split of 0.1. This approach allowed us to generate density plots showing the posterior probabilities of demographic parameter estimates in comparison with the simulated models.

Synchronicity of demographic events

To evaluate whether divergence events between insular and mainland populations occurred synchronously or independently among the three species, we employed Ecoevolity v.1.0.0 (Oaks, 2019; Oaks et al., 2019, 2020; https://github.com/phyletica/ecoevolity), which implements a full-likelihood Bayesian framework. It directly estimates the likelihood of population history from genomic data (Bocalini et al., 2021) and samples the posterior probabilities of competing models using MCMC (Oaks, 2019).

For this analysis, we generated biallelic matrices with mainland populations of C. alba and E. sulcata grouped, while for S. terebinthifolia only the divergence between Alcatrazes Island and Bertioga was retained, following the topology obtained with SNAPP. This resulted in three pairwise comparisons between insular and mainland populations. Three values were adopted for the Dirichlet process concentration parameter (DPprobs): 1) α = 0.554, corresponding to a prior favoring a single divergence event; 2) α = 1.5, corresponding to a prior favoring two divergence events; and 3) α = 3.564, corresponding to a prior maximizing the probability of three independent divergence events. The mutation rate was rescaled according to the generation time of each species, assuming a generation time of 2 years for C. alba and S. terebinthifolia, and 5 years for E. sulcata (Carvalho, 2003). The estimate for S. terebinthifolia was based on published cultivation and life-history information, whereas generation times for C. alba and E. sulcata were inferred from closely related taxa with similar growth form and reproductive characteristics within Rubiaceae and Myrtaceae, respectively. To evaluate sensitivity to divergence-time priors, three values were tested: 1) t ~ Exp(mean = 0.001); 2) t ~ Exp(mean = 0.005); and 3) t ~ Exp(mean = 0.00001).

Three independent MCMC chains of 50,000 generations were run, sampling parameters every 50 generations for each prior. Convergence was assessed using pyco-sumchains and visualized with Tracer v1.6 (Rambaut et al., 2018). To calculate posterior probabilities and summarize the results obtained with Ecoevolity, we used sumcoevolity, performing simulations to compute the Bayes Factor for all possible numbers of divergence events (Oaks, 2019), discarding the first 101 samples of each chain as burn-in. The tools pyco-sumevents and pyco-sumtimes were used to generate R scripts plotting, respectively, the posterior probabilities for the number of divergence events and the posterior distributions of divergence times for each species comparison.

To detect whether there was concordance between divergence time and sea-level variation, we used the PleistoDist package in R (Tan et al., 2023), which quantifies the historical effects of eustatic sea-level change in the Pleistocene on the continental shelf. In addition, it was possible to visualize sea-level variation up to 1 mya, at intervals of one thousand years, based on the eustatic Pleistocene levels proposed by Bintanja & van de Wal (2008), detecting moments of connection and isolation between Alcatrazes Island and the mainland. For this, we used bathymetry data obtained from the General Bathymetric Chart of the Oceans (GEBCO; https://www.gebco.net).