Uplift associated demographic history
We used multidimensional folded site frequency spectrum (mSFS) to infer post-uplift recolonization dynamics using the continuous-time coalescent simulator fastsimcoal v.2.6 (Excoffier et al. 2013) on those species that showed a distinct cluster in the uplifted region. Based on the clustering results we defined three-population models for D. poha and L. segnis (i.e. North, Uplift, and South lineages) and four-population models for D. antarctica and O. neglectus(i.e. Far North, Nearby North, Uplift, and South lineages; see Supplementary Methods and Fig. S7-S8 for details).
Given the absence of pre-quake genetic data for the studied species, we considered a constant population size for all non-uplift lineages and characterized all models by a bottleneck event at the time of the earthquake. For each demographic model we performed 50-100 independent runs of 500,000 genealogical simulations and 40 cycles of Brent algorithm to estimate the expected SFS and composite likelihood for a given set of demographic parameters. We chose the best fit recolonization scenario by calculating the Akaike Information Criterion (AIC) and Akaike’s weights.
To infer the degree of synchrony in size change among the uplift lineages of the host kelps and their epifauna we performed a hierarchical co-demographic modeling using the R packageMulti-DICE (Xue & Hickerson 2017). We first projected down the SFS of each ‘uplift’ lineage to 38 haploid individuals. We simulated 100,000 scenarios with species expanding from a relative low population size by a factor 100 - 10,000. Effective population size before expansion was constrained between 100 and 1000. Species were considered co-expanding in scenarios where their expansions started within the same 50 years (τ buffer prior = 50 years). Additionally, a more relaxed 200 years co-expansion threshold (τ buffer prior = 200 years) was applied in a different set of simulations. In order to infer the most likely scenario, the proportion of co-expanding species, ζ, and the timing of the expansion events, τ, we generated the aggregate site frequency spectrum (aSFS; Xue and Hickerson 2015) for each simulation. We used R to sample the best 5% of the simulations of the aSFS (5000 simulations) using the rejection method implemented in the ABC package (Csilléry et al. 2012). We visualised the fit of the observed aSFS within the retained simulations using a PCA produced with the R package ade4 (Dray & Dufour 2007).