Genetic legacies of mega-landslides: cycles of isolation and contact across flank collapses in an oceanic island

Catastrophic flank collapses are recognised as important drivers of insular biodiversity dynamics, through the disruption of species ranges and subsequent allopatric divergence. However, little empirical data supports this conjecture, with their evolutionary consequences remaining poorly understood. Using genome-wide data within a population genomics and phylogenomics framework, we evaluate how mega-landslides have impacted evolutionary and demographic history within a species complex of weevils (Curculionidae) within the Canary island of Tenerife. We reveal a complex genomic landscape, within which individuals of single ancestry were sampled in areas characterised by long-term geological stability, relative to the timing of flank collapses. In contrast, individuals of admixed ancestry were almost exclusively sampled within the boundaries of flank collapses. Estimated divergence times among ancestral populations aligned with the timings of mega-landslide events. Our results provide first evidence for a cyclical dynamic of range fragmentation and secondary contact across flank collapse landscapes, with support for a model where this dynamic is mediated by Quaternary climate oscillations. The context within which we reveal climate and topography to interact cyclically through time to shape the geographic structure of genetic variation, together with related recent work, highlights the importance of topoclimatic phenomena as an agent of diversification within insular invertebrates.

geographic isolation (Goodman et al., 2012;Machado, 2022;Salces-Castellano et al., 2021).In this vein, it has been recently argued for the need to focus attention on geomorphological dynamics within islands, highlighting the role of volcanic eruptions and other major landform-changing events, such as mega-landslides, in the evolutionary process (Otto et al., 2016).During the typical developmental life cycle of oceanic islands, geological activity acts to remodel existing landscapes.This geological dynamic is represented by three major events: eruptive volcanic activity, catastrophic flank collapse and millennial-scale erosion.Of these three phenomena, the immediate consequences of eruptive events and flank collapses may directly impact island biotas by provoking local extinction in affected zones (Borges & Hortal, 2009).Such extinctions driven by volcanism (Bloor et al., 2008;Carson et al., 1990;Vandergast et al., 2004) or by landslides (Brown et al., 2006;Juan et al., 2000;Macías-Hernández et al., 2013) may be followed by long-term habitat discontinuities generating genetic differentiation among populations (Goodman et al., 2012).As conditions within an ecological discontinuity generated by a catastrophic geological event become suitable, recolonization from adjacent zones could promote episodes of secondary contact among previously isolated populations.Evidence for such a complex dynamic is limited, but several phylogeographic studies within the Canary Island of Tenerife provide some support.For instance, Brown et al. (2006) revealed the potential role of the Güímar flank collapse on cladogenesis through population fragmentation and isolation within the lizard Gallotia galloti, based on genetic and morphological data.Population differentiation within the spider species Dysdera verneaui coincides with the western and eastern divisions of the Anaga peninsula, separated by the scarps of a megalandslide (Macías-Hernández et al., 2013).An estimated landslide age of 0.5-1.0 million years (Ma) (Watts & Masson, 2001) was interpreted as support for a potential causal relationship, as it coincides with the inferred divergence time between the western and eastern D. verneaui lineages (Macías-Hernández et al., 2013).Similar east-west genetic discontinuities coinciding with the limits of the flank collapse within the Anaga peninsula have been found across 13 species of beetle, with genomic evidence for secondary contact and admixture for 12 of the 13 species, and estimated divergence times ranging from ≤ 0.1 to 4.6 Ma (Salces-Castellano et al., 2020).
These idiosyncratic divergence times and contemporary admixture, together with more recent evidence for a much older origin for the landslide between 4.2 and 4.7 Ma (Walter et al., 2005) argue against the generality of a causal relationship between mega-landslides and intraspecific divergence.Rather, they reveal the synergistic role of landslide-sculpted topography and climatic oscillations throughout the Quaternary that have forced distributional shifts across flank collapse limits, promoting cycles of isolation and secondary contact (Salces-Castellano et al., 2020, 2021).
It is yet to be understood whether the dynamic of isolation and secondary contact observed by Salces-Castellano et al. (2020, 2021) is context-dependent, or if it may apply to other flank collapse systems.More studies are needed that integrate detailed populationlevel sampling with genetic markers that provide sufficient resolution to allow fundamental predictions from such complex dynamics to be tested.Catastrophic eruptive and erosional activity are likely to be consistent features throughout much of the life cycle of an oceanic island (Jackson, 2013).This in turn suggests a consequential evolutionary impact, and the potential for hypothesis testing when such events are clearly documented in the geological record.One of the best characterized archipelagos, from a geological point of view, is that of the Canary Islands (Carracedo & Troll, 2016).Within this archipelago, the island of Tenerife presents a complex geological history, in which three older volcanic shields are believed to have been originally isolated and then merged within the last 3.5 Ma due to successive volcanic activity (Ancochea et al., 1990;Cantagrel et al., 1999), or alternatively that two more recent shields formed at the margins of an older and larger central shield, much of which was subsequently overlain with more recent subaerial volcanic activity (Carracedo & Pérez-Torrado, 2013;Carracedo & Troll, 2016;Guillou et al., 2004;Figure S1).During the last 2 Ma, Tenerife has suffered several major eruptions (Ancochea et al., 1990(Ancochea et al., , 1999;;Huertas et al., 2002), and has been the subject of many flank collapses (Hunt et al., 2014;Figure S1), including some of the largest recorded megalandslides within the archipelago.The well-documented geological history of Tenerife provides a suitable framework to investigate the impact of geological events on diversification within oceanic islands.
Taxa with evidence for recent and ongoing diversification, in turn, constitute fertile ground for investigating the mechanisms promoting divergence among populations.
The dispersal-limited Laparocerus tessellatus complex of weevil species has been demonstrated to be a suitable model to assess the influence of landscape history on geographic patterns of individual relatedness (García-Olivares et al., 2017, 2019).The joint action of topographic complexity and Quaternary climate oscillations has been revealed to have shaped patterns of genomic divergence and admixture on the island of Gran Canaria, through climate-induced elevational shifts in distribution limits (García-Olivares et al., 2019).
Beyond topoclimatic variation as an engine for diversification within geologically quiescent islands, geologically active islands are also likely to structure genetic variation within species (Gübitz et al., 2000;Juan et al., 2000;Thorpe et al., 1996).In contrast to the relative geological dormancy of Gran Canaria over the last 3 Ma, Tenerife represents an island characterized by explosive volcanic activity and numerous gravitational flank collapses over the same period of time.The integration of a fine-scale sampling of the L. tessellatus complex within the geologically active island of Tenerife thus provides an excellent framework to investigate the role of catastrophic geological activity, in particular mega-landslides, on evolutionary dynamics within islands.
In the present study, we evaluate the within-island evolutionary consequences of mega-landslides using the L. tessellatus species complex on Tenerife, which is composed of two taxonomically described species, L. tessellatus Brullé, 1839 and L. freyi, Uyttenboogaart, 1940, the latter with four recognized subspecies: L. f. freyi, L. f. vicarius Machado, 2022, L. f. punctiger Machado 2016, and L. f. canescens Machado, 2016.Two testable predictions can be made to evaluate the role of mega-landslides on the demographic and evolutionary history within the L. tessellatus species complex.The first prediction is that geologically more stable areas within an island may act as reservoirs for population persistence, unlike areas that have suffered catastrophic flank collapses.The second prediction is that areas derived from flank collapses are likely to promote geographical isolation followed by admixture among genomically divergent populations colonizing from areas peripheral to the landslide, a process that in turn may result in an increase of genetic diversity under recent or ongoing admixture events (Boca et al., 2020).To test these predictions, the L. tessellatus complex on Tenerife was sampled representatively across its range, and genome-wide data was integrated into a population genetics and phylogenomics framework to examine spatial patterns of genetic variation across the island.Specifically, we firstly quantify geographic patterns of genetic diversity and population structure to evaluate their correspondence with contrasting areas of geological stability and flank collapse.Second, we use a simulation-based approach to evaluate competing scenarios of population isolation and estimate timeframes of divergence and gene flow among the inferred genetic groups.These analyses are in turn used to test for temporal congruence with range fragmentation linked to well-documented flank collapses or, alternately, explained by more recent dynamics of population splitting.Finally, we estimate changes in the effective population size through time in order to identify to what extent populations may have responded in parallel to common processes promoting persistence or driving divergence.

| Sample collection
Representative geographical sampling from the Tenerife species of the L. tessellatus complex was achieved by complementing previous sampling from Faria et al. (2016) and García-Olivares et al. (2017) with 74 specimens from 61 new localities.This additional sampling effort gave rise to a total of 126 individuals from 102 sites in Tenerife (Table S1).We also included five individuals from the L. tessellatus complex, belonging to a sister clade from the nearby island of Gran Canaria (García-Olivares et al., 2019), as an outgroup.

| ddRAD-seq library preparation
We extracted DNA using the Qiagen DNeasy Blood & Tissue kit following the manufacturer's instructions.DNA was processed using the double-digestion restriction-site associated DNA sequencing protocol (ddRADseq, Peterson et al., 2012) as described in Mastretta-Yanes et al. (2015) and García-Olivares et al. (2019).In brief, DNA was digested with the restriction enzymes MseI and EcoRI (New England Biolabs, Ipswich, MA, USA).Genomic libraries were pooled at equimolar ratios and size selected for fragments between 200 and 250 base pairs (bp) and, then, sequenced in a single-end 100-bp lane on an Illumina HiSeq2500 platform (Lausanne Genomic Technologies Facility, University of Lausanne, Switzerland).
Methods S2 provides all details on sequence assembly and data filtering.Unless otherwise indicated (see raxml and bpp analyses), we performed all downstream analyses using datasets of unlinked SNPs (i.e. a single SNP per RAD locus) obtained with ipyrad considering a clustering threshold of sequence similarity of 0.85 (clust_threshold) and discarding loci that were not present in at least 80% of individuals (min_samples_locus).Optimal parameter values in ipyrad were identified according to the sensitivity analyses conducted by García-Olivares et al. ( 2019) for the L. tessellatus species complex.To exclude the possibility that we had sampled close relatives, we calculated the relatedness between all pairs of genotyped individuals using the relatedness2 function in vcftools version 0.1.16(Danecek et al., 2011).

| Genetic clustering analyses
Population genetic structure was assessed using two complementary approaches.First, we used the Bayesian Markov Chain Monte Carlo (MCMC) clustering method implemented in the program structure version 2.3.3 (Pritchard et al., 2000).We ran structure with 200,000 MCMC cycles after a burn-in step of 100,000 iterations, assuming correlated allele frequencies and admixture (Pritchard et al., 2000) and performing 40 independent runs for each value of K ancestral populations (from K = 1 to K = 10).The most likely number of ancestral populations was estimated after retaining the 10 runs per each K-value with the highest likelihood estimates.Convergence across runs was assessed by checking that the 10 retained replicates per K-value provided a similar solution in terms of individual probabilities of assignment to a given ancestral population (q-values; Gilbert et al., 2012).As recommended by Gilbert et al. (2012) and Janes et al. (2017), we used two statistics to interpret the number range of ancestral populations (K) that best describes our data: log probabilities of Pr(X|K) (Pritchard et al., 2000) and ΔK (Evanno et al., 2005), both calculated in structure harvester (Earl & vonHoldt, 2012).Finally, we used the Greedy algorithm in clumpp version 1.1.2to align replicated runs of structure for the same K-value (Jakobsson & Rosenberg, 2007).
We also visualized the major axis of genomic variation by performing a Principal Component Analysis (PCA) as implemented in r version 4.0.3(R Core Team, 2021) using the package adegenet (Jombart, 2008).Missing data were replaced by the mean frequency of the corresponding allele estimated across all samples using the 'scaleGen' function (Jombart, 2008).

| Genetic differentiation among individuals
We estimated individual-based genetic distances using Nei's distance (Nei, 1972) as implemented in the r package StAMPP (Pembleton et al., 2013).To assess the relationship among inter-individual genetic distances and geographic distances (isolation-by-distance scenario, IBD), we calculated pairwise weighted topographic distances between each pair of individuals based on a digital elevation model (DEM) at 90-metre resolution using the r package topoDistance (Wang, 2020).Methods S3 provides details for the calculation of weighted topographic distances.The matrices were analysed using multiple matrix regressions with randomization (MMRR; Wang, 2013).

| Analysis of genetic diversity and admixture
Heterozygosity is expected to be higher in recently admixed individuals owing to the recombination of source ancestral genomes carrying different genetic variants (Boca et al., 2020;Kolbe et al., 2008;Witt et al., 2023).To test for a positive relationship between heterozygosity and admixture, we calculated observed heterozygosity (H O ) per individual using vcftools.Individuals were assigned to admixed and non-admixed population groups, according to the probability of assignment (q-value) to each ancestral population.Individuals showing a high q-value (>95%) to a given ancestral population were assigned as single ancestry, while those showing intermediate q-values (5% < q-value >95%) were considered to be of admixed origin.Sensitivity analyses were conducted with a more conservative threshold for admixed ancestry (10% < qvalue >90%), and with admixed individuals divided into two groups to account for the potential effect of varying admixture proportions on genetic diversity (Boca et al., 2020).The first group represents more even representation of ancestral population genomes (q-values: 30%-70%), and the second group more uneven representation (q-values: between 5%-30% and 70%-95%; and, alternatively, between 10%-30% and 70%-90%).ANOVAs were used to test for significant differences for H O between admixed and non-admixed groups of individuals.To take into account possible geographical clines of genetic diversity (Guo, 2012), relationships between H O and spatial variables (longitude and latitude) were tested with linear regressions.

| Phylogenomic inference
Phylogenetic relationships among all individuals were reconstructed using a maximum-likelihood (ML) approach as implemented in raxml version 8.2.12 (Stamatakis, 2014).Analyses were conducted on a matrix of concatenated SNPs (all SNPs per locus), applying an ascertainment bias correction using the conditional likelihood method (Lewis, 2001).We used a GTR-GAMMA model of nucleotide evolution, performed 100 rapid bootstrap replicates and searched for the best-scoring maximum likelihood tree.Individuals from Gran Canaria were used as an outgroup (Table S1).
We reconstructed phylogenetic relationships among the main genetic groups inferred with clustering analyses, using two coalescent-based methods for species tree estimation: svdquartets (Chifman & Kubatko, 2014) and snapp (Bryant et al., 2012).Individuals were assigned to three genetic groups, hereafter referred to as West, North and South, according to the most optimal clustering scheme as inferred in structure (K = 3, see Results section).Only individuals showing the highest probability of assignment (q-value >95%) to a given ancestral population were used for these analyses, resulting in 36 individuals evenly distributed across the three demes.Additional details on svdquartets and snapp analyses are described in Methods S4.
Finally, we used bpp version 4.4.1 (Flouri et al., 2018) to estimate the timing of divergence among the three main genetic groups identified by structure.Branch length (τ) estimation was performed by fitting the phylogenetic tree inferred in svdquartets and snapp as the fixed topology (option A00 in bpp).These analyses were conducted using the same subset of 36 individuals analysed for group-level phylogenetic inferences.We estimated divergence times using the equation τ = 2μt, where τ is the divergence in substitutions per site estimated by bpp, μ is the mutation rate per site per generation and t is the absolute divergence time in years (Walsh, 2001).We assumed a mutation rate of 2.8 × 10 −9 per site per generation calculated for Drosophila melanogaster (Keightley et al., 2014) which has been estimated to be similar to the spontaneous mutation rate calculated for the butterfly Heliconius melponeme (Keightley et al., 2015).Details on τ estimation in bpp can be found in Methods S4.

| Inference of past demographic history
We inferred the demographic history of each genetic group using stairwayplot2 version 2.1.1,which implements a flexible multi-epoch demographic model on the basis of the site frequency spectrum (SFS) to estimate changes in effective population size (N E ) over time (Liu & Fu, 2020).These analyses were conducted using the subset of 36 individuals with single ancestry assignment.Given a reference genome is not available, we calculated the folded SFS of each genetic group using the script easySFS.py(I.Overcast, https:// github.com/ isaac overc ast/ easySFS).We considered a single SNP per locus to avoid the effects of linkage disequilibrium.Each genetic group was downsampled to ∼66% of individuals (i.e. 8 individuals per deme) to remove all missing data for the calculation of the SFS, minimize errors with allele frequency estimates, and maximize the number of variable SNPs retained.Final SFS for West, North and South contained 13,930, 11,774 and 11,270 variable unlinked SNPs, respectively.Analyses in stairwayplot2 were run fitting a mutation rate of 2.8 × 10 −9 per site per generation (Keightley et al., 2014) and considering a one-year generation time (Machado & Aguiar, 2019).
We performed 200 bootstrap replicates to estimate 95% confidence intervals.

| Testing alternative models of divergence and gene flow
We used a simulation-based approach as implemented in fastsimcoal2 version 2.5.2.21 (Excoffier et al., 2013) to statistically evaluate the fit of our observed data to alternative scenarios of gene flow between genetic groups and estimate the timing of divergence and gene flow between them (Figure S2).These analyses were conducted using the same subset of 36 individuals analysed for group-level phylogenetic inferences.According to phylogenomic analyses, all scenarios considered an early split between the West genetic group and the ancestor of the North and South genetic groups (T DIV1 ), followed by their divergence (T DIV2 ).Using this topology, we evaluated scenarios assuming (i) divergence in strict isolation, (ii) post-divergence gene flow among the three populations and (iii) post-divergence gene flow only between West and North, and North and South.We modelled the timing of gene flow as a time interval, represented by a parameter for the time gene flow was initiated (T MIG1 ) and another parameter for the time that it ended (T MIG2 ).In the two isolationwith-migration scenarios, the timing of gene flow (T MIG1 and T MIG2 ) was modelled to be fixed across population pairs (models iia and iiia) or, conversely, to vary independently across population pairs (models iib and iiib).We additionally tested the aforementioned isolation-with-migration scenarios without estimating gene flow timing parameters (models iic and iiic).Finally, we replicated these seven scenarios, allowing for exponential demographic changes in effective population size (N E ) within each deme since the last historical event (divergence or migration).This approach yielded a total of 14 alternative models (Figure S2).Details on composite likelihood estimation, model selection approach and calculation of confidence intervals for parameter estimates under the most-supported model are described in Methods S5.

| Genomic data
Illumina sequencing provided a total of 379.54 M sequence reads, with an average of 2.90 M sequence reads per individual (SD = 1.74 M) (Figure S3).After the different filtering and assembly steps, each specimen retained on average 50,388 clusters (SD = 12,583), with a mean depth per locus of 35.79 (SD = 13.58)across individuals.All pairs of genotyped individuals had negative relatedness values (ranging from −3.01 to −0.01), which excludes the possibility that close relatives were included in the analysis (Manichaikul et al., 2010).

| Genetic clustering analyses
structure analyses identified the most likely number of ancestral populations to be K = 2 under the ΔK criterion (Evanno et al., 2005; Figure S4).One of the two ancestral populations was sampled in the westernmost part of the island, the Teno massif, with the other distributed across southern, northern and eastern areas of the island.Individuals with strong signatures of genetic admixture were identified within the northern flank, geographically coinciding with the landslides of Roques García (0.6-1.3 Ma), La Orotava (0.54-0.69 Ma) and Icod (0.15-0.17Ma) (Figure S5).However, LnPr(X|K) steadily increased up to an asymptote of K = 4. Within this range of likely K values, the largest probability gain was from K = 2 to K = 3 (Figure S4), suggesting that genetic variation is hierarchically structured within ancestral populations inferred by K = 2 (Janes et al., 2017).Accordingly, when considering K = 3 genetic variation was further structured into geographically coherent populations hereafter referred to as West, North and South (Figure 1a).In a spatial context, individuals inferred to be of single ancestry almost exclusively clustered together geographically, largely within areas that have not been affected by flank collapses, including scarps flanking areas of collapse.Consistent with inferences for K = 2, individuals of admixed ancestry between West and North were distributed across the aforementioned areas of flank collapse.
Additionally within K = 3, individuals inferred to be of admixed ancestry between East and South were sampled within the geographic limits of the Güímar landslide (0.83-0.85 Ma; Figure 1a).Accordingly, ancestry coefficients (q-value) across individuals sampled within landslide limits were statistically higher than those observed in individuals distributed in areas that have not been affected by flank collapses (ANOVA: F 1,124 = 119.20;p < .001; Figure S6).These geographic patterns of clustering and admixture persisted when assuming higher K values.For K = 4 and K = 5, structure analyses identified hierarchical population subdivision within the northern population, with individuals assigned with high probability to two ancestral populations, one distributed in the easternmost area of the Anaga peninsula (K = 4), and the other on the scarps that define the western limits of the Orotava flank collapse (K = 5; Figure S5).For K ≥ 3, individuals of admixed ancestry were consistently found within the La Orotava and Güímar valleys, with a third more limited area of admixture within the western area of the Anaga peninsula with K = 4 and K = 5 (Figure S5).Consistent with inferences from structure, principal component analysis (PCA) also supported genetic variation to be largely organized into three geographically concordant genetic groups, differentiated across the two first components (Figure 1b).Two gradients of admixture were identified within the PCA plot (Figure 1b F I G U R E 2 Phylogenetic relationships inferred with snapp (panel a) and svdquartets (panel b) among the three main genetic groups according to structure results, assuming three ancestral populations (K = 3) as the most likely clustering scheme.Panel (b) shows amonggroup divergence times estimated using bpp with a subset of 5000 randomly chosen loci.The topology was fixed using the phylogenetic tree inferred using snapp and svdquartets.Bars on nodes indicate the 95% highest posterior densities (HPD) of the estimated divergence times.Asterisks denote fully supported nodes in both panels.described along PC1, while differences among individuals from East and South were structured along PC2.Thus, the distributions of individuals from the three ancestral populations within the PCA are congruent with their geographical distributions within the island of Tenerife (Figure 1a,b).

| Individual-based genetic differentiation, diversity and admixture
Individual-based genetic differentiation estimated with Nei's distances were significantly correlated with weighted distances (R 2 = 0.593, β = 0.039, p = .001,Figure S7), consistent with a scenario of isolation-by-distance (IBD).Analyses of average observed heterozygosity (H O ) among individuals of single ancestry (q-value >95%) and admixed individuals (5% < q-value >95%) revealed significant differences (ANOVA: F 4,120 = 21.81;p < .001;This result remained significant across all sensitivity analyses (all ANOVAs: F 4,120 > 11.60; p < .001).Non-parametric Kruskal-Wallis rank sum tests yielded similar results.Finally, Ho was not significantly associated with latitude or longitude (p > .075),indicating that admixture provides a better explanation for geographic variation in genetic diversity over geographic gradients.

| Phylogenomic inference
The phylogenetic tree reconstructed in raxml including all individuals revealed the existence of three principal clades corresponding to West, North and South groups, consistent with results from structure (Figure S9).The raxml tree supported a sister relationship between South and North, within a lineage derived from an earlier divergence that gave rise to West.Individuals of admixed ancestry were phylogenetically placed basally within their respective clades (Figure S9).
Analyses in svdquartets and snapp focused on individuals of single ancestry according to structure yielded similar inferences.Topologies inferred by both svdquartets and snapp showed fully supported phylogenetic relationships among the structure-derived groups and supported an early divergence giving rise to West, followed by a more recent split between North and South (Figure 2), in concordance with the individual-based raxml analyses.Alternative runs in snapp assuming different prior distributions provided similar topologies and branch lengths.According to bpp analyses, the aforementioned splits were estimated to be approximately 0.0033 and 0.0023 τ units (Figure 2).Assuming a mutation rate of 2.8 × 10 −9 per site per generation (Keightley et al., 2014) and a 1-year generation time (Machado & Aguiar, 2019), both diversification events are estimated to have taken place during the Pleistocene, in the Chibanian age, about 0.59 and 0.41 Ma, respectively.Inferences of divergence timing from bpp were consistent across runs based on datasets considering different subsets of 5000 loci.

| Inference of past demographic history
stairwayplot2 analyses indicate that the three genetic groups have experienced parallel demographic responses, undergoing severe demographic declines since the end of the last glacial maximum (LGM, ∼19-21 ka; Figure 3).Consistent with analyses of genetic diversity, West and North presented higher historical estimates of N E than South (Figure 3).

| Testing alternative models of divergence and gene flow
fastsimcoal2 analyses identified the most supported model as being that considering simultaneous gene flow between all population pairs (model iia, Table 1; Figure S2).Considering a 1-year generation time, fastsimcoal2 estimated that all populations diverged from a common ancestor during the Pleistocene, in the Calabrian age, about 1.13 Ma (95% CI: 1.02-1.26Ma; Figure 4).
The posterior divergence between North and South was inferred to take place in the Chibanian age, about 0.50 Ma (95% CI: 0.45-0.56Ma), consistent with estimates using bpp.Historical gene flow F I G U R E 3 Demographic history of the three main genetic groups as estimated in stairwayplot2.Genetic groups correspond to structure results assuming three ancestral populations (K = 3) as the most likely clustering scheme.Continuous and dashed lines represent the median estimate and 95% confidence intervals of the effective population size (N E ), respectively, as inferred in stairwayplot2.Both axes are logarithmically scaled, with the X axis representing thousands of years (ka) and the Y axis representing N E .Highlighted period on the X-axis shows the extent of the last glacial maximum (LGM: ∼19-21 ka).
was estimated to be higher both between West and North, and between North and South, compared to between West and South (Figure 4).This coincides with contemporary patterns, with admixture between West and North and between North and South individuals in the Orotava and Güímar valleys, respectively, but no evidence for admixture between West and South (Figure S5).
Gene flow was estimated to have initiated approximately 154 ka (95% CI: 98-293 ka), considerably after the most recent diversification event (0.50 Ma), and to have ceased approximately 80 ka (95% CI: 26-120 ka; Figure 4).These estimates for the timing of TA B L E 1 Comparison of alternative models tested using fastsimcoal2 (Figure 4; Figure S2).

Note:
The best-supported model is highlighted in bold.Models were built both not considering (i) and considering migration among demes (ii and iii).The timing of gene flow was modelled to be fixed across population pairs (iia and iiia) or, conversely, to vary independently across population pairs (iib and iiib).Gene flow was modelled among all demes (ii) or only among West and North, and North and South genetic groups (iii).We also tested scenarios of gene flow without estimating timing parameters (iic and iiic).For each model, the statistics of its alternative scenario (−gr) accommodating exponential demographic changes in effective population size (N E ) within each deme are also indicated.
Abbreviations: N E , effective population size; lnL, maximum likelihood estimate of the model; k, number of parameters in the model; AIC, Akaike's information criterion value; ΔAIC, value from that of the strongest model; ω i , AIC weight.

F I G U R E 4
Parameters inferred from coalescent simulations with fastsimcoal2 under the best-supported demographic model.For each parameter, we show its point estimate and lower and upper 95% confidence intervals.Model parameters include ancestral (θ ANC , θ ANC-S-N ) and contemporary (θ W , θ S , θ N ) effective population sizes, timing of divergence (T DIV1 , T DIV2 ), timing of gene flow (T MIG1 , T MIG2 ) and migration rates per generation (m).Note that the effective population size of the West genetic group (θ W ) was fixed in fastsimcoal2 analyses to enable the estimation of other parameters.
gene flow largely coincide with the time interval where N E is estimated to have been higher, according to stairwayplot2 analyses (Figure 3).

| DISCUSS ION
It has previously been revealed that mega-landslides can act as drivers for inter-island dispersal of species (García-Olivares et al., 2017).However, the evolutionary and demographic consequences of mega-landslides within islands, while having received more interest (e.g. Brown et al., 2006;Juan et al., 2000;Machado, 2022;Macías-Hernández et al., 2013), remain less clear.

| Geological stability, population persistence and differentiation
The geographic distribution of individuals assigned to single ancestral populations within the L. tessellatus species complex of Tenerife largely corresponds to areas characterized by long-term geological stability (Figure 1).Hierarchically, three main ancestral populations were consistently inferred in structure, of which two were distributed within the northwest and northeast of Tenerife, respectively.These two regions broadly correspond to the Teno massif and the Anaga peninsula (Figure S1), with both regions having remained relatively geologically stable since the end of the Miocene (Carracedo & Pérez-Torrado, 2013), a time interval that encompasses the estimated origin and subsequent diversification of the L. tessellatus complex within Tenerife, as inferred in bpp and fastsimcoal2 (see also Faria et al., 2016;Machado, 2022;Machado et al., 2017).Individuals assigned uniquely to the third ancestral population are almost exclusively associated with areas outside of, but proximate to, scarps defining the Orotava, Roques García and Güímar flank collapse limits.These terrains predate the flank collapses with which they are associated (Figure 1a; Figure S1, prior to the onset of the current interglacial (Berger et al., 2016;Petit et al., 1999).Geographic gradients of admixture between ancestral populations are consistent with more recent and ongoing gene flow (Verdu & Rosenberg, 2011).
Species within the L. tessellatus complex have limited dispersal ability, contributing to the geographic structuring of their genetic variation over small spatial scales (García-Olivares et al., 2019).While limited dispersal would favour a narrow area of admixture, limited genomic incompatibility among populations and time would favour geographically more extensive gene flow (McEntee et al., 2020).
Across the combined northern flank collapses of Roques de García, Icod and La Orotava, a gradient of admixture (ancestry assignment to a single population < 90%) spans a geographic distance of 25 km.
Across the single southern flank collapse of Güímar, a second gradient spans a geographic distance of 10 km.While it remains uncertain when the secondary contact for these two gradients of admixture was initiated, we speculate that it is likely to have been some time after the LGM.Within this temporal window, time since initial secondary contact has been sufficiently long to give rise to geographically extensive admixture within the constraints of limited dispersal.
However, higher observed heterozygosity among individuals of admixed origin, compared to individuals of single ancestry, reveals that secondary contact has been sufficiently recent such that genetic diversity within populations remains above equilibrium expectations (Alcala et al., 2013).Globally, climate transition from the LGM until the present is typically characterized by upslope shifts for both the lower and upper elevation limits of species (Davis & Shaw, 2001;Rahbek et al., 2019).However, understanding how the distribution limits of the L. tessellatus complex may have changed from the LGM until now is complicated by evidence that the lower elevation limits of the orographic cloud bank have been forced downslope since the LGM (Salces-Castellano et al., 2021).Given the consistent structuring of genomic variation for the L. tessellatus complex in Anaga

| CON CLUS IONS
The paradigm view of the potential evolutionary consequences of mega-landslides on oceanic islands is one of instantaneous range disjunction for species with upper elevational limits that fall below maximum scarp heights, potentially followed by secondary con-

ACK N OWLED G EM ENTS
We wish to thank Loren Rieseberg, Melisa Olave and one anonymous referee for their constructive and valuable comments on an earlier version of the manuscript.We also thank Centro de Supercomputación ), both corresponding to the two main geographic clines identified in structure, each delimited within an area of flank collapse.Differences in allele frequencies between West and all other individuals were largely F I G U R E 1 Panel (a) depicts the geographical location of the sampled individuals and their ancestry coefficients (pie charts) as inferred by structure, assuming three ancestral populations (K = 3).Panel (b) shows a principal component analysis (PCA) of genetic variation for the sampled individuals.Pie charts represent the position of each individual along the two first principal components (PCs) and their respective ancestry coefficients based on structure results, assuming K = 3. Size of pie charts in both panels represents individual-based genetic diversity, estimated as observed heterozygosity (H O ).

Figure
Figure S8).Individuals derived from both West-North and North-South admixture presented significantly higher H O than corresponding single-ancestry populations (post hoc Tukey's tests: p < .001 in all comparisons involving admixed and single-ancestry populations).
Sampling across a landscape encompassing a sequence of geographically proximate flank collapses, we reveal a dynamic of population isolation and secondary contact that coincides with the landscape features of past flank collapses.In support of our first prediction, we found that individuals with ancestry assignment to a single population were characteristic of relatively geologically stable areas that have neither suffered recent volcanic activity nor flank collapses.The estimated timings of divergence among ancestral populations fall within the geological age estimates for flank collapses within the intervening landscape between areas characterized by ancestral genotypes.In support of our second prediction, individuals with signatures of mixed ancestry were typically sampled within areas of flank collapse.Overall, our study provides a conceptual framework for evaluating the effects of complex geological dynamics in generating novel genetic variation within islands over short spatial scales, through geographic isolation, population persistence and posterior admixture.

Figure
FigureS5), thus favouring the persistence of genomic variation that was contemporaneous to variation extirpated within the areas of flank collapse.Despite inherent uncertainties in divergence time estimates, the timing of divergence among the three ancestral populations (bpp and fastsimcoal2) aligns with the timeframe of the northern flank collapses of Roques de García (0.6-1.3 Ma), Orotava (0.54-0.69 Ma) and Icod (0.15-0.17Ma), while the split between the South and North populations is estimated to have initiated subsequent to the southern mega-landslide of Güímar (0.83-0.85 Ma;Hunt et al., 2014).The geological events across northern Tenerife are suggestive of a cumulative effect on divergence across the geographically proximate, and in part overlapping, flank collapses of Roques de García, Orotava and Icod.

4. 3 |
Quaternary climate and species range within flank collapse topography Analyses with structure under increasing values of K revealed finerscale geographic structuring of North with K = 4 (Figure S5), where genomic variation is organized into western and eastern ancestral populations across the northeastern Anaga peninsula of Tenerife, with geographically intermediate admixed individuals.This pattern was further corroborated by running structure with only the 11 individuals sampled in the Anaga Peninsula (Figure S10).This structure coincides with that observed for 13 co-distributed beetle species within the cloud forest of the Anaga peninsula (Salces-Castellano et al., 2020), including three related species of Laparocerus.Salces-Castellano et al. (2021) have revealed that this shared structure across species is best explained by a dynamic of isolation and secondary contact driven by climatic oscillations of the Quaternary.Quaternary climate oscillations within a topographically complex landscape have also been found to explain isolation and secondary contact within the L. tessellatus complex on Gran Canaria (García-Olivares et al., 2019).Given the findings ofSalces-Castellano et al. (2020, 2021)  andGarcía-Olivares et al. (2019), patterns of admixture across areas of flank collapse within Tenerife are plausibly mediated by range fragmentation and isolation above scarps during glacial climate conditions, with subsequent range expansion and secondary contact during interglacial periods.Further support for such a dynamic comes from demographic reconstructions (Figure3).In areas occupied by individuals assigned uniquely to one ancestral population, effective population sizes (N E ) are estimated by stairwayplot2 to have decreased substantially since the end of the last glacial maximum (LGM: ∼19-21 ka).This generalized response to warming temperatures across all three populations highlights the sensitivity of the focal taxa to climatic variation under a scenario of niche conservatism(Wiens et al., 2010).

(
Figure S10) with that observed by Salces-Castellano et al. (2020), it can reasonably be assumed that lower elevation limits for the L. tessellatus complex across the northern slopes of Tenerife have shifted downslope since the LGM.One possibility is that the reduced elevation gradients imposed by scarps when lower elevation limits shift downslope may have facilitated establishment and expansion across the more gradual slopes of flank collapse valley floors.Although the orographic cloud formations are a dominant influence across the northern slopes of Tenerife, elements of laurel forest formations within the scarps of the valley of Güímar (del Arco-Aguilar & Rodríguez-Delgado, 2018) also highlight their potential influence across southern slopes.However, more understanding is needed about local variation for the influence of the orographic cloud layer through time.
tact.Here we have found support for a model within which the orographic features left by a flank collapse have a more lasting influence on evolutionary processes within species.We reveal that Quaternary climate oscillations can give rise to a cyclical dynamic of range fragmentation and secondary contact across flank collapse landscapes.This dynamic is most likely to be influenced by the sharp elevation gradients associated with scarp height and should be most consequential for species with limited dispersal ability that occupy higher elevations.More generally, our results highlight the role of climate and topography in enhancing genetic diversity within insular species distributions, through both the establishment of divergent populations and the recombination of their genomes across areas of secondary contact.AUTH O R CO NTR I B UTI O N S BCE, VG-O, JP and VN conceived the original idea.VN and BCE led the study.VG-O, HL and BCE performed fieldwork.VG-O and YA performed laboratory work with assistance of JP and conducted exploratory analyses.VN conceived the methodological approach and conducted formal analyses.BCE and VN wrote the manuscript.All authors contributed critically to the draft and gave final approval for publication.
de Galicia (CESGA) and Teide High-Performance Computing facility (TeideHPC) provided by the Instituto Tecnológico y de Energías Renovables (ITER), S.A. for access to computer resources.Fieldwork was supported by the Cabildo of Gran Canaria (No. Exp.: 167/15) and Cabildo of Tenerife (No. Sigma: 015-00218).We acknowledge support for the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).This work was supported by the Ministry of Economy and Competitiveness (MINECO) through grants CGL2013-42589-P, CGL2017-85718-P and PID2020-116788GB-I00, co-financed by FEDER.VN was supported by a Juan de la Cierva-Formación postdoctoral fellowship (FJC2018-035611-I) funded by MCIN/AEI/10.13039/501100011033.VG-O was funded by a FPI pre-doctoral fellowship (BES-2014-067868) from MINECO.JP was supported by the Juan de la Cierva Program-Incorporation