Parallel recolonizations generate distinct genomic sectors in kelp following high-magnitude earthquake disturbance.

Large-scale disturbance events have the potential to drastically reshape biodiversity patterns. Notably, newly vacant habitat space cleared by disturbance can be colonized by multiple lineages, which can lead to the evolution of distinct spatial "sectors" of genetic diversity within a species. We test for disturbance-driven sectoring of genetic diversity in intertidal southern bull kelp, Durvillaea antarctica (Chamisso) Hariot, following the high-magnitude 1855 Wairarapa earthquake in New Zealand. Specifically, we use genotyping-by-sequencing (GBS) to analyse fine-scale population structure across the uplift zone and apply machine learning to assess the fit of alternative recolonizaton models. Our analysis reveals that specimens from the uplift zone carry distinctive genomic signatures potentially linked to post-earthquake recolonization processes. Specifically, our analysis identifies two parapatric spatial-genomic sectors of D. antarctica at Turakirae Head, which experienced the most dramatic uplift. Based on phylogeographical modelling, we infer that bull kelp in the Wellington region was probably a source for recolonization of the heavily uplifted Turakirae Head coastline, via two parallel, eastward recolonization events. By identifying multiple parapatric genotypic sectors within a recently recolonized coastal region, the current study provides support for the hypothesis that competing lineage expansions can generate striking spatial structuring of genetic diversity, even in highly dispersive taxa.

INTRODUCTION

Large‐scale disturbance events have potential to substantially reshape biodiversity patterns (Dornelas, 2010; Fraser, Davies, et al., 2018; Hewitt, 2004; Newman, 2019; Parvizi et al., 2019). Although the immediate biological impacts of major disturbances are almost inevitably negative, such events also provide opportunities for genetic turnover, demographic change and the biogeographical expansion of surviving lineages (Boessenkool et al., 2009; Fraser et al., 2009; Sheil, 2016), as well as diversification (Beheregaray et al., 2003; Wallis et al., 2016). In recent years, the widespread application of genomic techniques has started to provide detailed insights into the evolutionary, demographic and biogeographical implications of large‐scale disturbance events (Cole et al., 2019; Lescak et al., 2015; Parvizi et al., 2020).

Newly formed or vacant habitat space cleared by recent disturbance can potentially be colonized by numerous lineages. Through selection, one lineage may dominate a newly formed habitat by being better able to colonize or survive in the new environment (e.g., Canestrelli et al., 2016; De Meester et al., 2002). Alternatively, neutral and density‐dependent processes can underpin biogeographical patterns resulting from recolonization events. Under the Founder Takes All (FTA) hypothesis (Waters et al., 2013), for instance, the first lineage that colonizes an area can benefit from a “priority effect” (De Meester et al., 2016) by having more time to establish itself and expand spatially and numerically across a habitat. The resultant high‐density population can subsequently exclude late‐arriving dispersers from the newly colonized habitat via “high‐density blocking” (Ibrahim et al., 1996).

In some cases, parallel recolonization events involving different portions of a newly vacant habitat can potentially lead to the evolution of distinct “sectors” of genetic diversity (Waters et al., 2013). Such density‐dependent sectoring processes have been observed in real time during the growth of multilineage microbial colonies (Hallatschek et al., 2007), and biogeographical studies have identified parapatric sister‐lineages that evolved through parallel recolonization events following deglaciation events (Hewitt, 2004; Santucci et al., 1998; Von Saltzwedel et al., 2016; Waters et al., 2013). When and where competing expansion events of sister lineages come into secondary contact, associated phylogeographical boundaries may be defined by distinctively narrow hybrid zones, maintained via either adaptive or density‐dependent processes (Hewitt, 2001; Wallis et al., 2016). Despite the wealth of evidence for such parapatric sectoring processes involving relatively anciently evolved lineages (see deglaciation examples above), few comparable examples have been documented for more recent disturbances. Here, we capitalize on a major historical earthquake‐driven disturbance event affecting central New Zealand (McSaveney et al., 2006) to investigate recolonization and spatial sectoring processes.

The Cook Strait in central New Zealand straddles an active plate boundary in transition from subduction to strike‐slip (Wallace et al., 2012), and the southern coast of the North Island has been subjected to repeated tectonic disturbance for millennia. The M8+ magnitude 1855 Wairarapa earthquake remains the most powerful recorded earthquake in New Zealand since European colonization, and is among the largest recorded coseismic offsets worldwide (Darby & Beanland, 1992; Rodgers & Little, 2006). The earthquake caused uplift along ~50 km of the North Island's southern coast (Figure 1; Figure S1), causing numerous landslides that discharged debris into the sea (Hancox, 2005), as well as a major tsunami and many liquefaction events (Goff et al., 1999). The epicentre of the earthquake was along the Wairarapa Fault, and the greatest amount of coastal uplift of ~6 m occurred between Turakirae Head and Wharekauhau on the western side of Palliser Bay (Begg & McSaveney, 2005; McSaveney et al., 2006; Figures 1, S1). In the west, coastal areas in Wellington received moderate uplift ≤1.5 m, but Sinclair Head is generally accepted as the end of uplift >1 m (Begg & McSaveney, 2005), and uplift of 0.3 m was recorded for Cape Terawhiti (Figure 1; Darby & Beanland, 1992). Porirua Harbour on the southeastern coast of the North Island (Figure 1) did not experience uplift (Eiby, 1990), and no uplift is evident east of Wharekauhau across Palliser Bay (Figure 1; Darby & Beanland, 1992).

FIGURE 1

(a) Map of sampling locations for Durvillaea antarctica along the southern coast of the North Island (numbered 1–22). Locations are coloured according to sample region (Table 1). The red arrows and text indicate the estimated, smoothed height of tectonic uplift (in metres) for the 1855 Wairarapa earthquake, based on previous research (Begg & McSaveney, 2005; Darby & Beanland, 1992; Eiby, 1990). An inset map indicates the location of the region in New Zealand. Figure S1 presents further information, including all sample locations across the Cook Strait (including Marlborough), and it specifies the amount of uplift at each location and illustrates the principal tectonic faults. (b) A photo of D. antarctica at Manurewa Point (Location 21)

TABLE 1

Sampling of Durvillaea antarctica from central New Zealand. All samples were taken from attached individuals, except for one sample from Rarangi–marked with an asterisk (*)–that was freshly beachcast. The first column indicates the labels for the sample locations in Figure 1, and each row and sample location is coloured according to sample region: Boom Rock (cyan), Wellington (yellow), Western Turakirae Head (red), Eastern Turakirae Head (magenta), Southern Wairarapa (blue), and Marlborough (green). The same colouration is for all figures

Label on Figure 1	Region	Location	Collection date	N
1	Boom Rock	Boom Rock	29/07/2020	10
2	Boom Rock	Boom Rock Bay	29/07/2020	8
3	Boom Rock	Smiths Gully	29/07/2020	10
4	Wellington	Karori Rock	04/02/2020	14
5	Wellington	Sewer outlet before Karori	04/02/2020	2
6	Wellington	NW of Sinclair Head	04/02/2020	2
7	Wellington	West of Sinclair Head	04/02/2020	5
8	Wellington	Sinclair Head	03/12/2019	9
9	Wellington	Red Rocks	03/12/2019	5
10	Wellington	NE of Red Rocks	03/12/2019	4
11	Wellington	The Sirens Rocks	15/07/2021	5
12	Western Turakirae Head	Orongorongo Beach	02/12/2019	6
13	Western Turakirae Head	Turakirae Head	02/12/2019	12
14	Eastern Turakirae Head	Barneys Wharae	02/12/2019	5
15	Eastern Turakirae Head	Fishermans Rock	02/12/2019	4
16	Eastern Turakirae Head	Windy Point	02/12/2019	10
17	Eastern Turakirae Head	Mukamukaiti Rocks	02/12/2019	5
18	Eastern Turakirae Head	Corner Creek	05/12/2019	10
19	Southern Wairarapa	Rocky Point	04/12/2019	7
20	Southern Wairarapa	Cod Rocks	04/12/2019	5
21	Southern Wairarapa	Manurewa Point	05/12/2019	10
22	Southern Wairarapa	Oterei Rocks	05/12/2019	11
A	Marlborough	Rarangi	13/04/2012	2*
B	Marlborough	Cape Campbell	20/11/2016	5
C	Marlborough	Ward Beach	19/11/2016	13
D	Marlborough	Wharanui	19/11/2016	10
			Total	189

Southern bull kelp, Durvillaea antarctica (Chamisso) Hariot, is a large brown macroalga that dominates many wave‐exposed rocky shores throughout the Southern Hemisphere, including New Zealand (Fraser et al., 2020; Hay, 1977; Thomsen & South, 2019). Durvillaea antarctica is restricted to the intertidal zone and is adapted to handle some desiccation at low tide, but it requires regular wave exposure (Hay, 1977). Earthquake uplift can therefore cause mortality in D. antarctica if it is sufficient to move individuals above the high tide mark. For example, shores uplifted by several metres during the 2010 Chilean and 2016 Kaikōura earthquakes resulted in large‐scale die‐offs of intertidal Durvillaea (Castilla et al., 2010; Clark et al., 2017; Peters et al., 2020; Schiel et al., 2019; Tait et al., 2021; Thomsen et al., 2021). The most severe uplift (~6 m) from the 1855 Wairarapa earthquake is therefore almost certain to have resulted in the local extirpation of D. antarctica, whereas shores subjected to more minor uplift (≤2 m) may have potentially retained a few survivors from the previous low intertidal zone. Refugial patches of original populations, even if small, greatly reduce the likelihood of turnover following disturbances (Fraser, Davies, et al., 2018).

The uplift zone of the 1855 Wairarapa earthquake has since been recolonized by D. antarctica (Figures 1, S1). Gametes and fertilized eggs of D. antarctica have a low dispersal capacity (Dunmore, 2006), but buoyant adults can detach from substrates and raft over vast distances, and they can remain reproductively viable (Fraser, Morrison, et al., 2018; Tala et al., 2017). The 1855 uplift zone could therefore have been recolonized by local survivors or instead by the offspring of individuals rafting from other regions. Although a recent study of Durvillaea along the Akatore Fault in Southern New Zealand identified uplift‐associated population structure from a prehistoric earthquake (~1000 years BP), which indicates that tectonic disturbance can indeed generate distinct spatial–genomic sectors in D. antarctica (Parvizi et al., 2020), such processes have yet to be assessed over more recent time frames.

The 1855 Wairarapa earthquake presents a natural experiment to test genetic recolonization processes over ecologically relevant (decadal) time frames. Here, we conducted fine‐scale geographical sampling and used genotyping‐by‐sequencing (GBS; Elshire et al., 2011) to analyse population structure across the uplift zone. The recolonization of the 1855 uplift zone may have resulted in several different outcomes, so we used a machine learning approach to test the fit of alternative habitat recolonization models. Specifically, the uplift zone may be dominated by a single lineage, or it may have been recolonized by multiple lineages––potentially resulting in distinct genetic sectors. If present, these sectors may have formed through stepwise recolonization or they may be parapatric, representing independent recolonization events. These lineages may descend from survivors in areas affected by moderate uplift, or they may instead represent dispersal over longer distances between regions.

MATERIALS AND METHODS

Sampling

Frond tissue samples were collected from 189 attached Durvillaea antarctica “NZ North” lineage (as described in Fraser et al., 2020) individuals from 26 sample sites across central New Zealand, with sampling focused on the southern coast of the North Island (Table 1, Figures 1, S1). These 26 locations represent the central distribution of the D. antarctica “NZ North” lineage, although populations occur elsewhere throughout the North Island, southward to Banks Peninsula in the South Island, and on the Chatham Islands (Figure S2). Most North Island samples were collected between December 2019 and July 2020, and most South Island samples were collected in November 2016 (Table 1). Samples from The Sirens Rocks in Wellington were collected under a Department of Conservation permit (93234‐MAR). At each location, the sampled individuals were spread ~5 m apart from one another, as a basic precaution to avoid sampling closely related individuals (i.e., kin).

DNA extraction and GBS

DNA was extracted and purified following the same method described in Peters et al. (2020), except that we used the updated Qiagen DNeasy Plant Pro DNA extraction kit. The kit protocol was followed, except that the lysis step was extended to 24 h, and immediately after lysis, samples were treated with 100 μl isopropanol and incubated at 65°C for 30 min, vortexing every 15 min. Three new GBS libraries including 185 D. antarctica samples were used for this study (Table S1; Appendix S1). We also used GBS data for four samples sequenced within a previous study (Peters et al., 2020), bringing the total number of analysed individuals to 189 (Table S1; Appendix S1). DNA was digested using the PstI‐HF enzyme, following the GBS protocol described by Elshire et al. (2011), with the same modifications described by Peters et al. (2020). The size selection varied between 200–500 and 200–600 bp (Table S1). The four libraries were sequenced on four separate runs using midoutput flow cells on the Illumina NextSeq 500 platform (75‐bp paired‐end; Table S1).

Processing of GBS data

A total of 921,175,594 reads were sequenced from the three new GBS libraries (Table S1). The process_radtags component of stacks 2.53 (Rochette et al., 2019) was used to demultiplex samples into paired forward and reverse reads per individual, using inline barcodes. The process_radtags component removed low‐quality reads and reads with missing barcodes or PstI restriction sites (‐c ‐q), which provided a total of 810,697,051 retained reads (Table S1). This process included the rescue barcode and RADtag parameter (‐r) to retrieve additional reads, and reads were truncated to 68 bp (‐t 68).

Demultiplexed pairs of reads were assembled into loci without a reference genome using the de_novo pipeline in stacks. In ustacks, the minimum depth of coverage used to create a stack was two (‐m 3), the maximum distance (in nucleotides) allowed between stacks was two (‐M 2) and the maximum distance allowed to align secondary reads to primary stacks was four (‐N 4). A bounded single nucleotide polymorphism (SNP) model was applied with the error rate not being allowed to exceed 5% (‐‐bound_high 0.05). In cstacks, the number of mismatches allowed between sample loci when building a catalogue was two (‐n 2).

Three different catalogues of loci were used to analyse genetic variation at various geographical scales relative to the uplift zone (Table 2). The “Cook Strait” data set included all 189 samples to provide broad context, whereas the “North Island” data set used only the 159 samples from the North Island, and the “Recolonization” data set focused on 131 samples from areas immediately adjacent to and within the uplift zone of the 1855 Wairarapa earthquake (Wellington, Western Turakirae, Eastern Turakirae and Southern Wairarapa; Table 1). Preliminary examinations of the GBS data revealed hierarchical variation among the sample regions, and so these three data sets were used to investigate genomic variation at different geographical scales. The relatively focused “Recolonization” data set was used to simplify the computationally demanding phylogeographical modelling approach. For the populations component of stacks, samples were organized into populations matching the six sample regions (Table 1, Figure 1). For the “Cook Strait” data set, each locus was required to be present in 4/6 populations and in 70% of samples within each population. In the North Island data set, each locus was required to be present in 3/5 populations and in 80% of samples within each population. In the Recolonization data set, each locus was required to be present in 2/4 populations and in 80% of samples within each population. A minimum minor allele frequency of 5% was enforced for loci (‐‐min_maf 0.05). Only the first SNP of each locus was used (‐‐write_single_snp) and all SNPs were processed as biallelic and assumed to represent nuclear loci.

TABLE 2

The loci data sets generated using stacks. The table lists the number of samples, number of groups (i.e., sample regions), populations settings, and number of loci excluded and retained for each data set after applying filters. The final column lists the mean missing data per sample for each data set, as estimated by vcfr

		Locus representation		Filtering
n	Data set	No. of groups	Populations settings	No. of LD	No. of genotype depth <8	No. of loci depth outliers	No. of highly correlated loci	Total excluded loci	No. of final variant loci (SNPs)	Mean missing data per sample
189	Cook strait	6	‐p 4 ‐r 0.70	28	325	24	N/A	375	8474	17.2%
	Cook strait (no HC)						3232	3607	5242	16.8%
159	North Island	5	‐p 3 ‐r 0.80	30	1526	12	N/A	1552	5426	11.9%
	North Island (no HC)						1383	2935	4043	12.4%
131	Recolonization	4	‐p 2 ‐r 0.80	53	608	33	N/A	680	7119	15.8%

The settings in stacks listed above were selected after iteratively modifying parameters in stacks (Table S2), as recommended by Mastretta‐Yanes et al. (2015). We aimed to maximize the number of variant loci, while paying attention to coverage depth, missing data per sample and per locus, and the risk of erroneously combining too many reads. Overall, these preliminary investigations revealed the data to be highly consistent across parameter changes, and relative to the other runs, the final selected parameter settings yielded a medium number of loci with moderate missing data.

Filtering loci

The initial output loci from stacks were filtered. Loci estimated to be in linkage disequilibrium (LD) within the panmictic population were identified using plink 1.9 with a cut‐off of 0.8 (Purcell et al., 2007). One locus from each pair of loci estimated to be in LD was removed at random. Using vcftools 0.1.16 (Danecek et al., 2011), genotypes for a locus were removed from individuals if they had a coverage depth below 8 reads (‐‐minDP 8), and then after recoding loci, sites with a relatively high proportion of missing data among all individuals (≥40% for the Cook Strait data set, ≥60% for the North Island data set, ≥35% for the Recolonization data set) were removed (‐‐max‐missing 0.4, 0.6, 0.35). The loci estimated to be in LD or to have low genotype coverage depth were organized into a list and excluded (‐B; Catchen et al., 2013) and the populations component of stacks was rerun (same settings as above).

Variation in coverage depth per locus was investigated in the subsequent data set using vcfr (Knaus & Grünwald, 2017), and loci that were outliers for mean coverage depth and the standard deviation of coverage depth were identified. The outlier range was 1.5 times the interquartile range, above the upper quartile and below the lower quartile. These coverage depth outlier loci were added to an updated, second list of excluded loci and the populations component of stacks was rerun to produce the main filtered data set. Alternative versions of the Cook Strait and North Island data sets were also produced to test for the influence of pseudoreplication among loci (Waples et al., 2021). These alternative data sets excluded highly correlated loci (r 2 > .8), which were identified using vcftools (‐‐interchrom‐geno‐r2 ‐‐min‐r2 0.8). These loci were estimated to be highly correlated, but were not identified as being in LD by the earlier, more lenient filtering step in plink. These loci were used to create alternative lists of excluded loci, alongside the coverage depth outliers, and the populations component of stacks was rerun to produce the alternative versions of these data sets.

Missing data per sample and variation in loci coverage depth (same methods as described above) was assessed for the filtered data sets using vcfr. Some output files from stacks were converted to different file formats in pgdspider 2.1.1.3 (Lischer & Excoffier, 2012) for some downstream analyses. For phylogenetic reconstructions, VCF files were converted to phylip format, with loci filtered to require at least four samples per locus (−m 4), using vcf2phylip 2.0 (Ortiz, 2019).

Genotypic analyses

Phylogenetic relationships among samples were inferred by constructing unrooted maximum‐likelihood (ML) and neighbour‐joining (NJ) trees with free rates using iqtree 1.6.12 (Nguyen et al., 2015) and vcf‐kit 0.1.6 (Cook & Andersen, 2017) respectively. The analysis in iqtree was conducted with 10,000 ultrafast bootstrap replicates and the implementation of the modelfinder algorithm. Trees were visualized in figtree 1.4.4 (FigTree, 2018).

Population structure was assessed using principal components analysis (PCA) implemented in adegenet 2.13 (Jombart, 2008; Jombart & Ahmed, 2011). We determined the maximum number of “meaningful” principal components (PCs) to interpret by comparing PC Eigen values. Population structure and admixture among samples was further assessed using lea 2.8 (Frichot & François, 2015), which analysed 10 values of K. The lea analysis was conducted with default settings. Genetic summary statistics including observed (H O) and expected heterozygosity (H E), the inbreeding coefficient (F IS) and allelic richness (AR) were estimated for sample regions in the Cook Strait data set using the R packages adegenet, vcfr, hierfstat (Goudet, 2005) and dartr (Gruber et al., 2018).

Recolonization modelling

We used the Recolonization data set to estimate the most likely pattern of recolonization for D. antarctica on Turakirae Head. Following a similar approach to previous studies on Durvillaea (Parvizi et al., 2021), we performed phylogeographical model selection using the R package delimitr 2.0.2 (Smith & Carstens, 2019).

delimitr uses custom demographic model sets defined in fastsimcoal26 (Excoffier & Foll, 2011; Smith & Carstens, 2019) and simulates a folded multidimensional site frequency spectrum (mSFS) under each user‐defined model. Next, it summarizes the data by coarsening the mSFS and creating a binned mSFS (Smith et al., 2017; Smith & Carstens, 2019). It then uses the simulated data to construct a random forest (RF) classifier and calculates out‐of‐bag (oob) error rates. The machine learning RF classifier is then applied to the observed data to find the model producing mSFS most similar to the observed mSFS (Smith et al., 2017; Smith & Carstens, 2019). delimitr calculates a posterior probability for the best fitting model by regressing against oob error rates (Smith & Carstens, 2019), following Pudlo et al. (2016). Under this approach, delimitr only simulates unlinked SNPs and ignores invariant sites, and so it does not require an estimated mutation rate, and it can be used for demographic model selection without parameter estimation in nonmodel organisms (Smith & Carstens, 2019; Lado et al., 2020).

Each sample region within the Recolonization data set (Wellington, Western Turakirae, Eastern Turakirae and Southern Wairarapa) was treated as a separate lineage for the coalescent models. We tested eight post‐uplift recolonization scenarios to determine the potential source population of the Western and Eastern Turakirae Head lineages and to resolve whether the two recolonization events occurred parallel or in a stepwise manner (see Figure 2a for illustration of models). We also tested two pre‐uplift null models, where all four lineages independently diverged from a common ancestor before the earthquake event (Figure 2a). We restricted our analysis to these 10 models because, given the absence of prequake genetic data, we could not make further assumptions regarding additional demographic processes that may have affected the studied populations. The VCF output file from stacks was converted to an observed folded mSFS input file for delimitr, using the easysfs.py wrapper (https://github.com/isaacovercast/easySFS) for dadi 2.1.1 (Gutenkunst et al., 2009).

FIGURE 2

Recolonization modelling for the Recolonization data set (7119 loci) using delimitr. (a) The 10 models tested to estimate the most likely recolonization route for Turakirae Head. Model 3 (highlighted in green and marked with a star) was estimated to be the most likely model. (b) A simplified map of the four geographically adjacent sample regions included in the Recolonization data set. Dashed arrows indicate the recolonization process under Model 3, and the solid arrow indicates the prevailing current of the Cook Strait

Given the absence of prequake genetic data for D. antarctica in the region, we assumed constant effective population sizes (N ) for the Wellington and Southern Wairarapa lineages as these areas received minor uplift (≤2 m) and no uplift respectively. In contrast, we considered the two heavily disturbed lineages on Turakirae Head (2–6 m of uplift) to have evolved through post‐uplift recolonization and subsequent range expansion. These post‐uplift recolonization events were characterized by a bottleneck event that persisted for 20 generations before the Western and Eastern Turakirae lineages returned to their current N . To evaluate how the length of the bottleneck events affected our phylogeographical model selection, we also replicated all 10 recolonization models considering a longer bottleneck duration that persisted for 40 generations. Detailed descriptions of the demographic parameters and prior ranges used for the models are provided in Table S3.

For all demographic models, we used geological information to choose prior search ranges for divergence time parameters. Accordingly, we fixed the divergence time of the heavily uplifted Turakirae lineages to the timing of the Wairarapa earthquake at 1855. For the older divergence time events in our models (Figure 2a), we set a uniform prior between a lower bound of the 1855 Wairarapa earthquake event and an upper bound of the Holocene sea‐level stability in New Zealand (164–7000 years BP; Clement et al., 2016). In all models we incorporated contemporary gene flow, drew current and historical N from a wide uniform prior U (10–200,000) and converted time priors to generations ago using the minimum 2‐year generation time of D. antarctica (Hay, 1977).

Using the fastsimcoalsims function in delimitr we simulated 10,000 data sets under each model. Following the delimitr user manual to choose the number of bins to coarsen and summarize the mSFS, we used 10 classes. We then used 500 decision trees to construct the RF classifier, calculated error rates and approximated the posterior probability of the best model.

RESULTS

After filtering, the Cook Strait, North Island and Recolonization data sets contained 8474, 5426 and 7119 loci respectively (Table 2). For each data set, a small number of loci (<100) were excluded for being in LD or for having low or outlying coverage depth (Table 2). A larger number of loci (325–1526) were excluded for low genotype depth (Table 2). The alternative versions of the Cook Strait and North Island data sets contained 5242 and 4043 loci respectively, after highly correlated sites were excluded in addition to the standard filtering (Table 2). Mean missing data per sample was relatively low (11.9%–17.2%) across all three data sets (Table 2), and missing data was consistent among most samples (Figures S3 and S4). After filtering, the coverage depth for loci was well constrained (Figure S5). Missing data and coverage depth was very similar for the alternative Cook Strait and North Island data sets without highly correlated loci.

Most samples from the same region clustered together in the ML trees for both the Cook Strait and North Island data sets (Figures S6 and S7), which is concordant with strong phylogeographical structuring observed in previous GBS data sets for Durvillaea (Parvizi et al., 2020; Peters et al., 2020). A small number of individuals (~13) clustered unexpectedly among the Wellington, Eastern Turakirae, Southern Wairarapa and Marlborough regions (Figures S6 and S7). In the Cook Strait data set, samples from Southern Wairarapa were the most distinct, but low support values for basal nodes meant that samples from the remaining regions effectively formed a polytomy (Figures S6 and S7). With highly correlated loci removed, relationships among regions changed slightly but support for basal nodes was still limited (Figures S8 and S9). Results for the NJ trees were similar, with and without highly correlated loci, although individuals from Boom Rock formed a nested clade among Wellington samples (Figures S10–S13).

For PCA, the broken‐stick test indicated that between three and six PCs should be retained for each data set (Figure S14), but only the first two PCs for each data set seemed to exhibit meaningful patterns and therefore we focus on plots using those axes. For the Cook Strait data set, PC1 (11.7% of variance) distinguished samples from Southern Wairarapa, and PC2 (7.9%) distinguished many samples from Marlborough. Samples from the remaining regions were clustered closely together across the two axes (Figure 3a). In the North Island data set, PC1 (15.2%) again distinguished samples from Southern Wairarapa, and the remaining sample regions were spread across PC2 (9.4%) with Wellington and Western Turakirae on one side of the axis, Boom Rock in the middle, and Eastern Turakirae at the opposite end (Figure S15). For both data sets, PCA results were almost identical with highly correlated loci removed (Figures S16 and S17), and the plot for the Recolonization data set was almost identical to the North Island data set (Figure S18). Comparing all data sets, PCA indicated that Wellington and Western Turakirae were the most similar sample regions, but that Eastern Turakirae was relatively distinct (Figure 3a; Figures S15–S18). Again, ~13 individuals clustered unexpectedly among the same regions—matching the phylogenetic results (Figure 3a; Figures S15–S18).

FIGURE 3

Genotypic variation in Durvillaea antarctica across central New Zealand. (a) Principal components analysis plot for the Cook Strait data set (8474 loci) PC1 (11.7% of variation) and PC2 (8.4%). No further patterns were observed across other retained PCs. Ancestry matrices generated by lea for the (b) Cook Strait data set and (c) North Island data set (5426 loci), showing two values of K each

Summary statistics for the Cook Strait data set for the six sample regions were mostly similar (Table S4). Observed heterozygosity was lower than expected for all regions, corresponding to positive inbreeding coefficient (F IS) values. Allelic richness was similar among all regions (Table S4).

According to cross‐entropy values, higher values of K (i.e., 6 or 7 onwards) were favoured for all data sets using lea admixture analysis (Figure S19). However, clustering among samples was clearly hierarchical and population structure was highly consistent across all values of K (i.e., 2–10) and among all data sets (Figures S20–S24). For example, the Wellington and Western Turakirae sample regions were always clustered together at lower values of K across all data sets, and the Boom Rock sample region was distinguished with high confidence for almost all values of K for the Cook Strait and Wellington data sets. At higher levels of K (i.e., 8–10), additional clusters identified particular sample locations within sample regions or seemed to identify consistent aspects of statistical uncertainty among samples (Figures S20–S24). In concordance with the phylogenetic and PCA results, ~13 individuals clustered unexpectedly among the sample regions, which probably represents ongoing dispersal (Figures S20–S24).

In the Cook Strait data set, samples from Boom Rock and Southern Wairarapa were readily distinguished (Figures 3b; S20). Although a cluster (from K = 4 onwards) did not align with any sample region and appeared to represent clustering uncertainty, individuals from Eastern Turakirae, Marlborough and Western Turakirae were each in turn distinguished from samples from Wellington (Figures 3b, S20). Lower values of K (i.e., 2–4) indicated an appreciable degree of admixture between Eastern Turakirae and Southern Wairarapa (Figures 3b, S20). Samples from Marlborough were only distinguished from the Wellington and Western Turakirae regions from K = 6 onwards, which potentially indicates low genetic differentiation across the Cook Strait (Figures 3b, S20). Western Turakirae was the last region to be separated from Wellington, indicating a close relationship between these regions and matching the PCA results (Figure 3a). Additional clusters identified at high values of K (i.e., 8–10) distinguished particular sample locations within regions or seemed to reflect further clustering uncertainty (Figure S20). Results with highly correlated loci removed were almost identical (Figure S21). Overall, the results for PCA and lea analysis were more informative than phylogenetic reconstructions, probably because the former analyses are better suited to interpreting admixture among samples (Meisner & Albrechtsen, 2018).

Results for the lea analysis of the North Island data set aligned closely with the Cook Strait data set (Figures 3c, S22). Boom Rock and Southern Wairarapa were readily distinguished from other regions, and at K = 3 the complete uplift zone of the 1855 earthquake (i.e., Wellington, Western Turakirae, Eastern Turakirae) was clearly identified as a single distinct unit of population structure (Figures 3c, S22). Subsequent values of K (i.e., 4, 5) hierarchically distinguished Eastern and Western Turakirae from Wellington (Figures 3c, S22), matching the pattern observed in the Cook Strait data set. However, low values of K (i.e., 2–4) did not indicate admixture between Eastern Turakirae and Southern Wairarapa (Figures 3c, S22), in contrast to the Cook Strait data set (Figure 3b). High values of K in the North Island data set (i.e., 6–10) continued to distinguish the six sample regions and additional clusters probably reflected clustering uncertainty (Figure S22). Results were again consistent with highly correlated loci removed (Figure S23), and results for the Recolonization data set were concordant with both larger data sets (Figure S24).

Using the Recolonization data set (7119 loci), phylogeographical model selection using delimitr broadly indicated that the Western and Eastern Turakirae Head coasts were probably recolonized by D. antarctica from Wellington rather than Southern Wairarapa (Table 3). Of the 10 tested scenarios, Model 3, where the Wellington lineage recolonized Western and Eastern Turakirae through two parallel, independent recolonization events, was the best fitting model (posterior probability = 0.89; Figure 2; Table 3). An alternative set of simulations using longer durations for post‐uplift bottleneck events (40 instead of 20 generations) also found best support for Model 3 (posterior probability = 0.96; Table S5). Model 3 also received the highest number of votes during model selection (Table S6). All six parallel and stepwise models without admixture had a low classification error rate (Tables 3, S5), which indicates a high statistical power for distinguishing different recolonization routes. Models 7 and 8 had a higher classification error rate due to confusion with Models 1 and 3, and Model 5 respectively (Tables 3, S5). This result is intuitive as Models 7 and 8 are derived, admixed versions of these models (Figure 2a). For the pre‐uplift null models, Model 9 had a high classification error rate due to confusion with Models 3 and 4 (Tables 3, S5), whereas Model 10 had a low classification error rate.

TABLE 3

Cross‐validation results for phylogeographical model selection in delimitr for the 10 models (Figure 2a), using 20 generations for the duration of genetic bottlenecks

Model	1	2	3	4	5	6	7	8	9	10	Classification error
1	8281	85	108	193	9	7	783	69	458	7	0.1719
2	64	8925	137	231	8	11	122	7	486	9	0.1075
3	197	296	7507	76	252	1	754	83	820	14	0.2493
4	260	230	68	8196	3	251	168	25	790	9	0.1804
5	12	15	287	24	9421	197	18	73	126	7	0.0759
6	7	11	8	359	197	9172	14	103	125	4	0.0828
7	2752	160	2115	142	39	4	3645	348	785	10	0.6355
8	1728	25	600	126	2161	478	1204	3244	419	15	0.6756
9	278	465	1199	1901	222	329	498	181	4892	35	0.5108
10	7	9	5	2	2	4	4	5	14	9948	0.0052

Model 3 was selected as the best fitting model (highlighted in grey) and had a posterior probability of .89. See Table S5 for alternative results using 40 generations for the duration of bottlenecks, and Table S6 for votes given per model.

DISCUSSION

Disturbance‐driven spatial‐genomic sectors along a continuous rocky coast

Our GBS analysis reveals that Durvillaea antarctica collected within the complete uplift zone of the 1855 Wairarapa earthquake are readily distinguishable from specimens sampled from elsewhere across central New Zealand (Figure 3). Our analysis identified two distinct, parapatric spatial–genomic sectors of D. antarctica on Turakirae Head (Figure 3), and the combined geographical boundaries of these two genetic units closely align with the major uplift zone for the 1855 Wairarapa earthquake (2–6 m; Figure 1). Previous genetic studies focused on the Akatore fault in south‐eastern New Zealand have similarly revealed uplift‐associated genetic structure for intertidal Durvillaea (Parvizi et al., 2019, 2020). However, our current study is the first to reveal multiple parapatric genetic sectors within a disturbed region (see Hewitt, 2001), which were probably generated by separate recolonization events (Figure 2).

The fine‐scale hierarchical phylogeographical variation detected across the complete 1855 uplift zone (Wellington, Western Turakirae and Eastern Turakirae) is notable given the small spatial scales involved. In particular, the Western and Eastern Turakirae spatial–genomic sectors occur across a continuous rocky coastline (~7 and 12 km respectively; Figure 1). This striking genetic structure most plausibly evolved as a consequence of the recolonization process itself (i.e., density‐dependent processes; see Hewitt, 2001, 2004; Ibrahim et al., 1996; Waters et al., 2013), rather than through any contemporary biogeographical barrier. Selection is unlikely to explain the observed genetic structure given the ecologically similar habitats involved, and the limited evolutionary time (maximum of 82 generations) since the earthquake. These data thus provide some of the first nonexperimental evidence of parapatric sectoring evolving over ecologically relevant (human) time frames as a result of natural disturbance.

Sources of recolonization

The Wellington region immediately west of Turakirae Head experienced only moderate tectonic uplift in 1855 (≤1.5 m; Figure 1), making it likely that many D. antarctica specimens survived in this western region. By contrast, the severe uplift across Turakirae Head (2–6 m; Figure 1) makes it highly unlikely that D. antarctica could have survived within this heavily impacted zone. For comparison, following the 2016 Kaikōura earthquake, intertidal Durvillaea populations were heavily reduced at sites uplifted by <2 m, but appeared to have been extirpated at sites of greater uplift (Peters et al., 2020; Tait et al., 2021; Thomsen et al., 2021). The hierarchical clustering observed in the lea analyses, and the strong support for Model 3 in the delimitr analysis, suggests that surviving specimens from the Wellington region were the source for parallel, independent recolonizations along the heavily uplifted coastline of Turakirae Head (Figure 2). Notably, the support for Model 3 over Model 5 in the delimitr analysis indicates that Turakirae Head was recolonized through two founding events, rather than through a stepwise recolonization process from west to east (Figure 2a; Table 3). Given that the null models were not supported over Model 3 or the admixture models (7 and 8), it is also unlikely that the observed pattern of population structure existed prior to the earthquake. The RF model selection approach that we applied here is, however, limited in assessing how well the observed data fit a model, and can only estimate the probability of the best model among those tested (Pudlo et al., 2016; Smith et al., 2017).

Recolonization from Wellington instead of Southern Wairarapa is plausible given the prevailing eastward current through the Cook Strait (Chiswell et al., 2017; Figure 2b). A small number of specimens (~13) exhibited unexpected genotypes associated with other regions (Figure 3), which probably represents ongoing dispersal and gene flow across the Cook Strait subsequent to the recolonization of the 1855 uplift zone. Three of these individuals appear to represent bi‐directional dispersal between Eastern Turakirae and Southern Wairarapa (Figure 3). In the lea analyses, low values of K (i.e., 2–4) for the Cook Strait data set indicated some admixture between Eastern Turakirae and Southern Wairarapa (Figures 3b, S20). However, this signal is absent in the other data sets (Figures S20–S24), and westward recolonization models in delimitr received very low support (Table 3).

All sample regions exhibited a deficit in observed heterozygosity (Table S4), which may indicate population substructure (the Wahlund effect), potentially reflecting the identification of some sample locations at higher levels of K in the lea analyses. The deficit may also have been caused by sampling a higher number of related individuals than expected by chance, which could reflect genuine elevated rates of inbreeding within local areas. Elevated rates of inbreeding are plausible given the limited dispersal ability of Durvillaea gametes (Dunmore, 2006) and the patchy and high‐density distribution of individuals at most locations. If true, the deficit in observed heterozygosity suggests that dispersal is inadequate to prevent elevated rates of inbreeding within most regions. There is almost certainly ongoing dispersal throughout the wider range of the D. antarctica “NZ North” lineage (Figure S2), but overall the results here indicate that wide‐scale extirpation remains a key requirement for the success of immigrant lineages and population genetic turnover in Durvillaea bull kelps (Fraser, Davies, et al., 2018). Future research involving broader‐scale geographical sampling of D. antarctica could reveal further evidence of dispersal, and indicate how the studied populations in the Cook Strait fit within the overall phylogeography of the species. Although recent research on another macroalga found that high relatedness among individuals within sample locations did not significantly impact estimates of population structure (Nepper‐Davidsen et al., 2021), there is clearly a need for future studies of Durvillaea to investigate the relationship between sampling distance and relatedness.

Multispecies impacts of disturbance

Major disturbance events can potentially have parallel effects across large numbers of codistributed species (Hewitt, 2001; Parvizi et al., 2021; Wallis et al., 2016). Along these lines, a recent GBS study of another southern bull kelp species, Durvillaea poha C.I. Fraser, H.G. Spencer & J.M. Waters, 2012, concluded that the 1855 Wairarapa earthquake may have created an ecological opportunity for a recent northward range expansion from its otherwise southern distribution (Vaux et al., 2021). Additionally, the absence of another congener, D. willana Lindauer, 1949, from the southern coast of the North Island has been suggested to reflect disturbance‐related impacts of the 1855 earthquake (Hay, 2020). Future studies could shed light on the biogeographical histories of various codistributed macroalgae across this region.

While the observed fine‐scale intraspecific patterns for D. antarctica are striking, previous phylogeographical studies of other coastal macroalgae (e.g., Buchanan & Zuccarello, 2012; Muangmai et al., 2015) and intertidal invertebrates (e.g., Ross et al., 2009; Veale & Lavery, 2012; Walton et al., 2019) have not identified corresponding biogeographical disjunctions within the 1855 uplift zone or linked to Turakirae Head. However, most of these studies used relatively coarse‐grained geographical sampling across New Zealand and analysed microsatellites or short‐range PCR (polymerase chain reaction) products, which may have provided less insight compared to the fine‐scale sampling and GBS used for this study. As such, there is potential for ongoing research testing for biogeographical anomalies associated with disturbance from the 1855 Wairarapa earthquake that may have been overlooked by previous studies.

Our study demonstrates that machine learning approaches can provide useful insights for phylogeographical and population genomic research. Given the wealth of data provided by reduced representation and whole genome sequencing, software such as delimitr provide an exciting opportunity to formally test biogeographical hypotheses. New machine learning methods continue to be developed (Greener et al., 2021), and so similar and improved analyses could be applied to other organisms potentially affected by large‐scale disturbance events.

AUTHOR CONTRIBUTIONS

J.R.W., C.I.F. and D.C. conceived the study; F.V. led the fieldwork, prepared the GBS libraries and conducted most genetic analyses; E.P. conducted the phylogeographical model selection; F.V. led the writing with the assistance of J.R.W., C.I.F., D.C. and E.P.; F.V. produced the figures.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

OPEN RESEARCH BADGES

This article has earned an Open Data badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.jdfn2z3bt.

Appendix S1

Click here for additional data file.

Figure S1–S24

Table S1–S6

Click here for additional data file.