Life on a beach leads to phenotypic divergence despite gene flow for an island lizard

The study investigates whether strong habitat differences within small islands can drive phenotypic divergence despite ongoing gene flow. While gene flow generally homogenizes allele frequencies and can impede divergence, examples in model systems show divergence and even speciation can occur with gene flow when selection is strong. Oceanic island lizards occupy diverse habitats over short distances, offering a natural system to test habitat-associated divergence without spatial isolation. The authors focus on Teira dugesii on Madeira, contrasting shingle beach versus inland habitats situated less than 1 km apart. They aim to: (1) test for replicated divergence in dorsal color and head morphology across four beach–inland pairs; and (2) assess whether such divergence has arisen despite ongoing gene flow using population genomic data and coalescent-based modeling.

Prior genomic studies in diverse taxa (e.g., sticklebacks, periwinkle snails, stick insects) demonstrate that divergence can occur with gene flow, particularly under strong divergent selection. On islands, many cases of within-island divergence involve historic or current interruption to gene flow (e.g., volcanism-driven fragmentation). Parallel adaptive divergence across microhabitats within single islands has fewer clear examples, though incipient divergence in Anolis dewlaps over short distances has been observed. In small vertebrates, dorsal color often matches substrate to enhance crypsis (e.g., White Sands lizards; Peromyscus mice). The genetic basis of melanism often involves allelic variation in genes like MC1R, though phenotypic plasticity may contribute in early stages. Head morphology has heritable components in lizards, but genomic underpinnings are less defined. Madeira’s T. dugesii had a previously noted melanic beach population at Caniço, suggesting a candidate system for replicated tests of divergence and gene flow.

Design and sampling: A matched-pairs design compared beach (B) and adjacent inland (I) habitats at four Madeiran localities (1: Caniço; 2: Porto da Cruz; 3: Paul do Mar; 4: São Vicente). B and I sites were 0.2–0.8 km apart; localities were 13–39 km apart. Beach sites comprised grey shingle/cobbles/boulders; inland sites were disused, vegetated agricultural areas with stone walls. Adults were trapped with baited containers: 216 males and 118 females photographed; 93 individuals (9–14 per site) tissue-sampled for genomics. Ethics approvals and permits were obtained. Color (dorsal luminance): Standardized dorsal photographs (Nikon D3300; ColorChecker reference) were normalized using grey standards and analyzed with MicaToolbox in ImageJ to extract mean RGB luminances for six dorsal/head regions. RGB values were log10-transformed. Two-way MANOVAs (factors: habitat, locality) tested effects separately by sex using Pillai’s trace; discriminant function analyses (DFA) visualized group separation. Morphology (head size and shape): Dorsal head images (60 mm macro lens, fixed setup) were landmarked with 35 type-1 landmarks (tpsDig). Geometric morphometric workflow (MorphoJ): Generalized Procrustes Analysis, separation of symmetric component, PCA to obtain PCs explaining 95% of variation. Head size as log-transformed centroid size (CS) analyzed by two-way ANOVA; head shape PCs analyzed by two-way MANOVA (habitat, locality, interaction). DFA on Procrustes coordinates assessed shape divergence. Habitat characterization: At each site, 9–10 quadrat photographs (0.25 m²) captured substrate RGB and vegetation cover. Two-way MANOVA tested RGB differences; vegetation cover quantified (%). Genomics (GBS): DNA from tail tips digested (ApeKI, PstI), libraries size-selected (350–450 bp) and sequenced (Illumina NovaSeq 6000). Reads were trimmed (AdaptorRemoval v2), QC’d (FastQC), and SNPs called with GBS-SNP-CROP v4.1. A mock reference was built from the highest-read individual due to lack of a reference genome. Filters: phred ≥30, depth 3–30, call rate ≥0.90, removal of SNPs with heterozygote excess (VCFtools; HWE 1%). Resulting ALLSNP dataset: 19,311 SNPs in 4,135 tags for 93 individuals; thinned dataset: 4,131 SNPs (1 per tag) with putative selected SNPs removed. Selection scans: pcadapt (R) identified outlier SNPs (MAF >5%, Bonferroni-adjusted p < 0.1); bayenv2 tested association of outliers with habitat (binary B/I) using a covariance matrix estimated from the thinned dataset; 10 independent MCMC runs (1,000,000 steps). Spatial and population structure: Pairwise FST (PopGenome). DAPC (adegenet) on ALLSNPs; number of PCs chosen by cross-validation (RMSE, MSAR). sPCA (adespatial) on thinned data with neighbor network linking B and I within localities; eigenvalue tests (9,999 randomizations) for local and global structure. Phylogenetic clustering: Treemix on thinned data (no outgroup), with 1,000 bootstrap trees to assess support for splits. Demographic inference (gene flow models): fastsimcoal2 compared three scenarios per locality using joint folded SFS: NOGFLOW (no migration), ONEGFLOW (constant migration post-divergence), TWOGFLOW (two periods of migration). INDSNP datasets (thinned, excluding outliers) were used for model selection via AIC; ALLSNP datasets for parameter estimation, with monomorphic sites estimated and fixed mutation rate 1e-8 per generation. Each scenario: 100 optimization cycles with 2×10^5 coalescent simulations per cycle; 100 replicates; best run retained. Parametric bootstrap (ALLSNP) generated 100 SFS replicates per locality (300 bp contigs) to obtain 95% CIs; re-optimization with 1×10^5 simulations, 50 cycles.

Replicated phenotypic divergence: • Color: Beach lizards were consistently darker (lower luminance) than inland lizards across all four localities and both sexes. MANOVA males: habitat effect Pillai’s Trace = 0.450, F6,201 = 27.44, P < 0.001; effect size partial η² = 0.45 (larger than locality 0.16 and interaction 0.10). Females: habitat Pillai’s Trace = 0.463, F6,105 = 15.08, P < 0.001; partial η² = 0.46 (greater than locality 0.13 and interaction 0.11). DFAs showed clear separation of beach vs inland on DF1 with beach individuals having lower luminance. • Head morphology: Head size varied by locality and interaction; no consistent size pattern across habitats (females: habitat F1,108 = 4.30, P = 0.041; locality F3,108 = 3.66, P = 0.015; interaction F3,108 = 2.78, P = 0.044. Males, after removing 5 outliers: habitat F1,200 = 3.70, P = 0.056; locality F1,200 = 3.92, P = 0.010; interaction F1,200 = 10.09, P < 0.001). In contrast, head shape showed strong, replicated divergence: males MANOVA (23 PCs) habitat Pillai’s Trace = 0.477, F23,183 = 7.26, P < 0.001; females (21 PCs) habitat Pillai’s Trace = 0.488, F21,88 = 3.99, P < 0.001. DFAs revealed beach lizards have broader snouts at all localities in both sexes (minor deviation in females at locality 2 with very small beach sample). Habitat differences: Substrate RGB differed strongly between habitats (MANOVA habitat Pillai’s Trace = 0.871, F3,69 = 154.70, P < 0.001); vegetation cover was 0% at beaches vs median >60% inland. Genomic data and selection: GBS yielded 19,311 SNPs (4,135 tags) from 93 individuals; thinned dataset 4,131 SNPs. pcadapt detected 52 outlier SNPs; only four showed significant habitat association in bayenv2 (only one supported in 9/10 runs), and none co-located on the same tag—indicating weak, dispersed signals of selection at the SNP panel scale. Structure and lineage tests: DAPC showed some regional separation (e.g., south vs north/east coasts) but no consistent beach–inland clustering; first two DFs explained 70.0% of variation. sPCA indicated significant local (obs = 34.39, P = 0.0002) and global (obs = 29.02, P = 0.0010) structure. Treemix found no evidence that beach and inland populations form distinct lineages. Demographic inference and gene flow: Across all four localities, TWOGFLOW (two periods of migration) was strongly favored over ONEGFLOW and NOGFLOW (AIC differences extremely large; NOGFLOW consistently worst). In the favored model, initial post-divergence migration rates were lower, followed by a more recent period of higher migration. Estimated split times ranged from ~1.75×10^5 (locality 1) to ~1.71×10^6 generations (locality 4), with overlapping 95% CIs. Recent migration rates were generally higher than ancient ones and often asymmetric, commonly greater from inland to beach (e.g., recent Inland→Beach: ~7.26×10^-4 to 1.64×10^-2; Beach→Inland: ~2.34×10^-4 to 9.81×10^-4), consistent across localities.

The study’s core question—can habitat-associated divergence arise within a small island despite ongoing gene flow—was supported by replicated patterns of darker dorsal coloration and broader snouts in beach-dwelling T. dugesii relative to inland counterparts less than 1 km away. Genomic analyses across thousands of SNPs found no evidence that beach populations constitute distinct evolutionary lineages, and coalescent model comparisons strongly favored scenarios with gene flow at all four matched localities. The replication of direction and form of phenotypic divergence across sexes and sites, together with substantial gene flow, supports divergent selection between beach and inland habitats as the driver. Potential selective agents include predator-mediated crypsis (e.g., kestrels) favoring darker coloration on dark shingle beaches, while morphological differences in head shape may relate to foraging or microhabitat use, though mechanisms remain to be tested. The genomic signal of selection was modest (few outliers associated with habitat), consistent with polygenic or subtle genetic architectures, gene flow diluting signals at neutral loci, and/or phenotypic plasticity contributing. Asymmetric gene flow (inland→beach > beach→inland) matches expectations of a large mainland-like source (inland) supplying a peripheral habitat (beach). The surprising trend of higher recent compared to ancient migration could reflect changes in coastal connectivity/topography and sea-level history on Madeira. Overall, the findings demonstrate that strong, parallel phenotypic divergence can be maintained across very short distances in the face of gene flow, highlighting microgeographic adaptation within islands.

Within-island phenotypic divergence in Teira dugesii on Madeira is repeatedly observed between shingle beach and nearby inland habitats, with beach lizards being darker and exhibiting broader snouts. Population genomic and demographic modeling reject the formation of distinct beach lineages and instead support persistent, often asymmetric gene flow, with recent migration higher than ancient migration. These results show that divergent selection across microhabitats can overcome gene flow to generate and maintain phenotypic divergence on small spatial scales. Future work should identify the genetic basis of dorsal coloration and head morphology (e.g., candidate genes like MC1R and genome-wide scans with a reference genome), quantify the roles of phenotypic plasticity versus heritable variation, experimentally test selective mechanisms (predation, habitat use, diet), and refine demographic timelines with improved ecological and paleoenvironmental data.

• Phenotypic overlap between habitats indicates substantial gene flow and/or plasticity; the genetic basis of traits was not directly validated. • Small sample sizes for some groups (e.g., female beach at locality 2) and removal of outlier individuals (five male size outliers) may affect some morphological inferences. • Assumption checks indicated deviations from normality and unequal covariance matrices; robust statistics were applied but residual biases may remain. • GBS without a reference genome required a mock reference; linkage and physical context of SNPs are unknown, limiting mapping of causal variants. • Selection scans detected few habitat-associated outliers and results varied across bayenv runs, potentially reflecting power limitations, polygenic adaptation, or gene flow. • Treemix lacked an outgroup, precluding inference of migration edges. • Demographic inference used estimated monomorphic sites and a fixed mutation rate; generation time for T. dugesii is not empirically established, adding uncertainty to absolute timing estimates. • The mechanism underlying higher recent than ancient migration rates is speculative, and alternative historical scenarios cannot be excluded.