Shared genomic outliers across two divergent population clusters of a highly threatened seagrass

Sample collection

Leaf samples (n = 336) were collected from 12 sites, including nine estuaries/estuarine bays along the South African coast, one bay in Mozambique (Inhaca) and one bay in Kenya (Shimoni; Fig. 1A; Table 1). At each location, with the exception of Inhaca, and Shimoni, three leaf samples were collected 2 m apart from five beds at two sub-sites (>100 m) for a total of 30 leaf samples per location. This sampling approach was designed to minimise the sampling of clones. Samples were preserved with silica crystals before being processed. Sampling permits were granted by SanParks and CapeNature (permit number 0028-AAA008-00159); DAFF and DEA (permit number was RES2014/103). Despite intensive questioning of collaborators and other contacts throughout this study, no samples of Z. capensis were obtained from Tanzania (L Nordlund, pers. comm., 2015).

Laboratory protocols

Accurately estimating genome-wide variation and detecting signals of local adaptation in non-model organisms, such as seagrasses, requires many individuals from many sites to be sequenced, which can be prohibitively expensive despite the advances made by high-throughput sequencing methods such as RADseq (Ellegren, 2014; Andrews et al., 2016). As such, a pooled sequencing (pool-seq) approach was utilised, combining the genomic DNA of multiple individuals before sequencing (Sham et al., 2002; Kofler, Betancourt & Schlötterer, 2012). This approach decreases cost whilst increasing the number of individuals analysed and allowing for a more population focussed analysis (Futschik & Schlötterer, 2010; Schlötterer et al., 2014).

Genomic DNA was extracted from leaf tissue using Qiagen DNeasy plant kit (Qiagen, Valencia, USA) following standard protocols, with the exception of eluting the DNA in nuclease-free water instead of elution buffer. Genomic DNA quality was then assessed using gel electrophoresis and DNA concentrations of each sample determined by Qubit analysis at the Central Analytical Facility of Stellenbosch University (CAF). Genomic DNA was pooled by location, with all individuals having equimolar representation, to create a total of 12 samples for Illumina sequencing. The two sites at Langebaan, Oostewaal and Geelbek (Fig. 1A), were kept as separate pools to allow for comparison between the observed morphotypes; one short and stunted on the muddy tidal flats (Geelbek) which experience prolonged exposure to conditions outside the water and the other is longer with a higher biomass on the sandy permanently submerged area (Oostewaal) (D Pillay, pers. comm., 2014).

Library preparation and sequencing followed the ezRAD method (Knapp et al., 2016; Toonen et al., 2013; Nielsen et al., 2018), which obtains a reduced representation sequencing library using high frequency restriction enzymes. Pooled genomic DNA was freeze dried before library construction following the protocol of Knapp et al. (2016) and Mi-Seq Illumina sequencing (V3 2x300) conducted at the Genetics Core Facility (GCF) of the Hawaii Institute of Marine Biology (HIMB). These data are stored in the National Center for Biotechnology Information’s (NCBI’s) Sequence Read Archive (SRA; PRJNA503110) and georeferenced at GeOMe (https://geome-db.org) at the project: Zostera capensis pooled RADseq.

Data processing and alignment

Quality of raw reads was analysed using FastQC (Andrews, 2010) and quality filtering carried out using FastQC Toolkit (Andrews, 2010), removing low quality bases (<20 phred score). Additionally, TrimGalore! v 0.4.4 (available at: http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) was used to remove any remaining adapter sequences or reads shorter than 30 bp. BWA-MEM (Li, 2013) was used to map filtered paired-end reads from each pooled sampling site to the genome of sister species, Zostera marina (available from NCBI, BioProject number PRJNA41721, GenBank accession number LFYR00000000 (Olsen et al., 2016), with a minimum mapping quality of 20. Ambiguously mapped reads, PCR duplicate reads, reads with less than 20 mapping quality and less than 20 base quality were filtered out before converting SAM files to BAM files with SAMtools (Li et al., 2009). Number of mapped and unmapped reads were then calculated using the idxstats command in SAMtools. As pools with a higher number of mapped reads may have an artificially inflated number of SNPs, mapped reads were subsampled to median coverage in SAMTools using the view command with the ‘-s’ flag. Although subsampling results in a loss of data, it is nonetheless important for correctly interpreting true differences between the sampling sites, as opposed to differences in data quality or quantity (Schlötterer et al., 2014). To confirm that subsampling removed any possible correlation between the number of mapped reads and the number of SNPs and outlier loci identified downstream, Spearman’s correlation coefficients were calculated using the rcorr function of the ‘Hmisc’ (Harrell Jr & Dupont, 2006) package in R (R Core Development Team, 2008). BAM files were sorted and indexed before being used to create pileup files for each individual sampling site with the mpileup command in SAMTools (Li et al., 2009), using a minimum quality score of 20 and maximum read depth of 10,000. Finally, a pileup file combining all sites was created using the same parameters in SAMtools and converted to a sync file using PoPoolation2 (Kofler, Pandey & Schlötterer, 2011) for downstream use.

Calling SNPs and simulating data

The total number of SNPs and private SNPs were identified using snp-frequency-diff.pl in PoPoolation2 (Kofler, Pandey & Schlötterer, 2011) with genomic sites required to have a minimum minor allele count of four, and coverage between 10 and 500 across all 12 sampling sites. SNPs were then filtered to retain only those present among sampling sites, and not those present due to differences between the reference sequence (Z. marina) and Z. capensis. As many software cannot handle pooled data, requiring individuals to be specified within sampling sites, subsample_sync2GenePop.pl in PoPoolation2 was used to simulate a multi-locus dataset of a subset of SNPs identified by PoPoolation2. Because this programme cannot simulate different numbers of individuals across sites, the median sample size of 30 individuals was selected for every site. The resulting GenePop file was then converted to various formats in PGDspider (Lischer & Excoffier, 2012) for downstream analyses.

Outlier loci identification and functional annotation

Due to the uncertainty surrounding RADseq, Pool-seq and outlier detection methods (Narum & Hess, 2011; Da Fonseca et al., 2016; Mckinney et al., 2016; Lowry et al., 2017; O’Leary et al., 2018), four outlier detection methods were employed, namely BayeScan v2.1 (Foll & Gaggiotti, 2008), Lositan (Antao et al., 2008), BayeScEnv (De Villemereuil & Gaggiotti, 2015) and PCadapt (Luu, Bazin & Blum, 2016), which includes F_ST-based approaches, genotype-environment correlations and principle component analyses (see supplementary materials for details).

Outlier loci identified by two or more methods were considered “candidate outliers” and their functional roles evaluated by subjecting 1,000 base pairs upstream and downstream of each of the identified outlier SNPs to BLASTx searches, with the non-redundant protein sequence database and an E-value cut off of 10⁻⁵ (Altschul et al., 1997) using Blast2Go (Conesa et al., 2005). In addition to BLASTing against the general NCBI database, these searches were also carried out against the Zosteraceae family in general, specifically, Zostera marina (Olsen et al., 2016) and Z. muelleri (Lee et al., 2016). Gene Ontology (GO) mapping, Interproscan (Jones et al., 2014) and annotation were performed with Blast2Go default parameters. The number and proportion of candidate outliers unique to each site and shared between pairwise sites was calculated. The overlap of outliers identified between the different approaches was visualised using the ‘VennDiagram’ package (Chen & Boutros, 2011) in R. The frequency of outlier loci identified by all four approaches (Lositan, BayeScan, BayeScEnv and PCAdapt) was plotted across sampling sites using the ‘ggplot2’ package (Wickham, 2009) in R and listed in Table S3. A one-way analysis of variance (ANOVA) and a post hoc TukeyHSD test, were performed with the ‘aov’ and ‘TukeyHSD’ functions in R, to compare outlier frequencies across sites.

Neutral variation

All identified outlier loci were flagged as being putatively under selection for the purposes of this analysis and were therefore removed from the dataset to isolate neutral drivers of patterns of population structure. The neutral-only multi-locus dataset set was then re-simulated, using subsample_sync2GenePop.pl in PoPoolation2 with 30 individuals per site as described above, and used to calculate measures of neutral variation.

Genome-wide variation and differentiation

To characterise genetic diversity, Tajima’s nucleotide diversity (π), Watterson’s theta (θ) and Tajima’s D were estimated for the complete and neutral-only datasets using a sliding window approach with Variance-sliding.pl in PoPoolation v1.2.2 (Kofler et al., 2011). For these comparisons, filtering stringency was reduced to a minimum minor allele count of two and coverage between 10 and 500 per sampling site. As the estimation of allele frequencies in pooled individuals is highly reliant on sequence coverage, a high sequence coverage and large sliding windows were used in order to increase accuracy (Kofler et al., 2011). Observed and expected heterozygosity and the inbreeding coefficient (F_IS) was estimated from the simulated datasets with the divBasic function of the ‘DiveRsity’ package (Keenan et al., 2013) in R.

To investigate genome-wide levels of differentiation, the fixation index (F_ST) for pairwise comparisons of populations was estimated using a sliding window approach with fst-sliding.pl in PoPoolation2, using a minimum minor allele count of four and a coverage between 10 and 500. Fisher’s exact test was carried out with fisher-test.pl in PoPoolation2 to estimate the significance of allele frequency differences between sites. Patterns of differentiation were visualised on a principle coordinates analysis (PCoA) plot generated in R using the pco function of the ‘labdsv’ package (Roberts, 2007). The PCoA plot was generated both with and without Kenya in order to account for sampling bias. The simulated neutral dataset was used to investigate population clustering by means of Bayesian Analysis of Population Structure (BAPS) software (Corander & Marttinen, 2006; Corander, Marttinen & Mäntyniemi, 2006) testing K = 1–10.

Habitat suitability for Z. capensis in the LGM

In order to understand the influence of historical environmental conditions on the contemporary patterns of genomic variability, the suitable habitat for Z. capensis was hindcast to the Last Glacial Maximum (LGM; 21 kya). Zostera capensis occurrence data was obtained from Adams, Veldkornet & Tabot (2016) and environmental data downloaded from MARSPEC at 5 arcminute resolution for both the present-day (Sbrocco & Barber, 2013) and the LGM (CNRM-CM33 model; Braconnot et al., 2007; Sbrocco, 2014). Following Chefaoui, Duarte & Serrão (2017), only Sea Surface Temperature (SST) of the coldest month (Biogeo14) and warmest month (Biogeo15) were utilised to avoid using strongly correlated variables for Species Distribution Modelling (SDM; Guisan & Thuiller, 2005; Braunisch et al., 2013; Chefaoui, Duarte & Serrão, 2017). Precipitation variables were excluded to decrease model uncertainty (Varela, Lima-Ribeiro & Terribile, 2015). These variables represent relevant present-day and LGM conditions, which are recognised as important determinants of intertidal seagrass habitat suitability (Short Neckles, 1999; Short et al., 2010; Valle et al., 2014) and they are projected along the present-day (Sbrocco & Barber, 2013) and LGM coastlines (Braconnot et al., 2007; Sbrocco, 2014), respectively. Environmental variables such as salinity and oxygen content were not included due to the high natural variability in estuarine systems over tidal and seasonal time-scales. QGIS (QGIS Development Team, 2012) was used to crop raster extents, by means of the buffer and crop tools, to focus on the coastal areas including and surrounding the present-day range of Z. capensis.

Ecological niche modelling was implemented through an ensemble approach with the ‘biomod2’ package (Thuiller et al., 2016) in R. As in Chefaoui, Duarte & Serrão (2017), the following six presence-absence algorithms were included in the ensemble models: generalized additive model (GAM), flexible discriminant analysis (FDA), generalized boosting model (GBM), multiple adaptive regression splines (MARS), generalized linear model (GLM), and random forest (RF). Default parameters were used for all algorithms, except for the GLM which was fitted with a quadratic term, the GBM which was run with 1,000 trees, and the GAM which was executed with the GAM_mgcv function. As the occurrence data (Adams, Veldkornet & Tabot, 2016) included reliable presence and absence records for estuaries along the entire South African coastline, no pseudo-absence selection was required. The data was split into a calibration (80%) and a validation (20%) set and three iterations were performed for each algorithm with three permutations to estimate and weight variable importance, for a total of 18 models. Models were assessed with the true skill statistic (TSS; Allouche, Tsoar & Kadmon, 2006) and the area under the receiver operating characteristic (ROC) curve (AUC; Fielding & Bell, 1997), considering both specificity (true negatives) and sensitivity (true positives). Only models scoring TSS >0.55 and AUC >0.8 were used to produce ensembles. Retained models were ensembled to produce a weighted mean SDM and first used to project the present-day habitat suitability, in terms of SST, along the South African coastline, and then used to hindcast the habitat suitability to the LGM. The present-day and LGM habitat suitability projections, as well as the changes in habitat suitability between the present-day and LGM were plotted in R.

Disentangling contemporary signals of IBD and IBE

A redundancy analysis (RDA) (Legendre & Legendre, 2012) was conducted to evaluate the relative contribution of spatial and environmental variation to genomic variability and patterns of genetic differentiation. RDA can be useful as a multivariate regression technique when running regression analyses with multivariate predictors (space and environment) and multivariate responses (here, minor allele frequencies of SNPs). As spatial distances are not suitable for constrained ordination or regression as implemented in RDA, geographic distances were transformed to Principal Coordinates of Neighbourhood Matrix (PCNM) distances with the pcnm function in the ‘vegan’ package (Oksanen et al., 2015) in R. Environmental distances were calculated within the RDA function from the variables in Table 2 (excluding the macrophyte species measure, which were only available for South Africa). The ordistep function from the ‘vegan’ package was used to select the most informative variables and build the ‘optimal’ model. Four separate RDAs were conducted with minor allele frequency as the response. Predictor variables in the first RDA were transformed geographic distances, and in the second were environmental distances. Lastly, two partial RDAs were performed, partitioning out the effect of transformed geographic distance and environmental variation from the total variation respectively. The anova function of the package ‘vegan’ was performed with 999 permutations (Legendre, Oksanen & Ter Braak, 2011) to test the significance of RDAs.

Sequencing and mapping

54,982,056 paired reads were obtained, with paired reads from each sampling site ranging from 1,368,372 to 7,429,328 (Table 3). After filtering reads for quality and adapters, and subsampling to a median, 7,432,397 reads, ranging from 222,741 to 750,736 per site, were aligned to the Z. marina reference genome (Table 3). The number of filtered subsampled mapped reads had no correlation with the number of SNPs (r = 0.17; p > 0.05) or outlier loci (r =  − 0.05; p > 0.05) identified.

Neutral and outlier loci

The complete simulated dataset consisted of 308 loci (Fig. 3). From this dataset, 101 outlier loci were detected by Lositan, while BayeScan and BayeScEnv detected 25 and five outlier loci respectively. The five outlier loci identified by the ecological association approach in BayeScEnv were correlated with precipitation of the driest quarter and annual mean moisture levels. By analysing allele frequencies of the non-simulated dataset, PCAdapt identified 38 outlier loci. All 169 outlier loci were removed from the complete dataset in order form a neutral-only dataset with which to examine patterns of neutral variation.

Genome-wide variation

The number of SNPs identified by PoPoolation2 in the neutral and complete datasets ranged from 845 to 1,683 and 913 to 1,784 per sampling site respectively. The number of private SNPs per site was generally low and five populations did not have any private SNPs (Tables 3 and 4). The genome-wide average nucleotide diversity (Tajima’s π) of the neutral and complete datasets ranged from 0.023 to 0.041 and 0.023 to 0.035, respectively. Allelic richness did not vary much between sites, ranging from 1.23 to 1.36. θ_W of the neutral and complete datasets were the same, ranging from 0.029 to 0.043 (Tables 3 and 4). The west and south coast sites, with the exception of Oostewaal (L2), exhibited marginally higher nucleotide diversity and θ_W than the east coast sites. Tests for deviations from neutrality produced genome-wide average Tajima’s D that were negative for all sampling sites and ranged from −0.723 to −0.275 and −0.706 to −0.273, for the complete and neutral datasets respectively. Genetic diversity metrics calculated from the simulated datasets included expected heterozygosity (0.04 to 0.06) within each sampling site (Tables 3 and 4), and the inbreeding coefficient, F_IS, which was uniform across sampling sites and equal to 1, for both the complete simulated dataset and the neutral simulated dataset. Although an F_IS of 1 is partly due to the nature of the simulated data, generating multiple individuals from a highly clonal pool, it nonetheless indicates extremely high levels of inbreeding.

Genome-wide differentiation and clustering

F_ST values were estimated from the complete non-simulated dataset for pairwise comparisons of sites (Tables S1 , S2), with Fisher’s exact tests showing no significant differentiation between pairs of sites for either dataset. Similarly, clustering analysis conducted in BAPS on neutral loci revealed no structure across sites, with all sites falling into one cluster (K = 1; p < 0.05; Fig. S1). Although there is no significant population structuring, the PCoA (Fig. 2) of pairwise F_ST values for neutral loci suggests that the west and south coast sites, (except for Oostewaal-L2), are more closely related than the east coast sites. The same pattern was observed for the PCoA generated without Kenya (Fig. 2).

However, when the clustering analysis in BAPS included outlier loci, two clusters were detected (p < 0.05; labelled cluster one and two), with cluster one comprising samples from the west and south coasts, and cluster two including populations from the east coast of South Africa in addition to Mozambique and Kenya (Fig. 1B). Notably, one west coast site in Langebaan, Oostewaal (L2), groups with cluster two rather than cluster one (Fig. 1B).

PCoAs of pairwise F_ST comparisons from the complete dataset and all outlier loci resulted in a similar, but slighter denser pattern than observed for the neutral dataset (Fig. 2). Sites from cluster one formed a tight group, relatively separate from the remaining sites. Sites from cluster two did not group as closely as those from cluster one, with Mozambique most differentiated. Moreover, Mozambique, followed by Kenya, exhibited much higher outlier allele frequencies than other sites (Table S3).

While some outlier loci were identified by more than one method, there was little overlap between outlier loci identified using the four different approaches (Fig. 3), with only three outliers shared between all four methods. However, irrespective of how many outlier loci are included, the frequency at which outlier loci occurred at each site reflects the two clusters identified using BAPS (Fig. 1B). ANOVA and post hoc TukeyHSD tests revealed that outlier frequencies were not significantly different across sampling sites (p > 0.05) but that they were significantly different across the two clusters, with higher frequencies observed in cluster 2 than cluster 1 (F_1,24 = 66.61, p < 0.001; Fig. 3). No private outliers were identified as all outlier loci occurred at two or more sites and all candidate outliers (identified by more than one method) occurred at most of the sites (Table S3).

Functional annotation of candidate outlier loci

Two-thousand base pairs surrounding each of 10 candidate outlier loci were subjected to the Blast2Go pipeline. Although all of the 10 candidate outliers yielded significant hits when BLAST searches were conducted against the general NCBI database, Zosteraceae, Z. marina and Z. muelleri, the majority of these hits did not fall within gene regions of known function. GO terms (GO:0016020-IEA ‘membrane’ and GO:0016021-IEA ‘integral component of membrane’) were assigned to five of the 10 candidate outlier loci with BLAST matches to hypothetical and predicted proteins (Table S3).

Habitat suitability for Z. capensis in the LGM

Multiple models from each algorithm met the TSS > 0.55 and AUC > 0.8 criteria and were retained to produce ensembles. Ensemble models obtained the following average validations scores: TSS = 0.654, AUC = 0.904, sensitivity = 92.11, specificity = 73.29. Predicted distributions of suitable habitat, in terms of SST, differed between present-day and LGM conditions, in terms of geographic location, extent and probability of occurrence (Fig. S2). The highest probability of occurrence can be seen on the south coast (up to ∼25° longitude) and west coast (up to ∼18° latitude) for the present-day projection, and on the western-south coast (up to ∼21° longitude) and west coast (up to ∼18° latitude) for the LGM projection. Ensemble models project an 11.05% loss and a 10.79% gain of suitable habitat from the LGM to present-day, with a 26.1% range shift. These shifts are most evident in the loss of suitable habitat on the south and south-east coasts (∼21–27° longitude), southern-west coast (∼30–35° latitude), and west coast (∼12–18° latitude), as well as the gain of suitable habitat on the northern-east coast of South Africa (∼30–35° latitude), the south coast of Madagascar and the northern-west coast of Africa (∼3–8° latitude; Fig. 4). Further, within a South African context, the western-south coast represents an area of stable temperature regime, where suitable habitat has occurred from at least as far back as the LGM until the present day (Fig. 4). This can also be seen in patches on the west coast and on the east coast of Africa (∼5–25° latitude).

IBD vs IBE

Of the 11 environmental variables, seven were selected by the RDA as the most informative (Table 2). The pure RDA of genetic variation against transformed geographic distance was not significant (P > 0.05), but was significant when carried out against environmental variation, with 70.4% of the variation in the data explained by the retained environmental variables. Unexpectedly, neither partial RDA analyses, the first conditioned on transformed geographic distance and the other on environmental variation, were significant. Although environmental variation explained such a high percentage of the variation observed in the data, partitioning out the effect of geographic distance on environmental variation rendered the association with genetic variation non-significant.

Genomic diversity of a threatened seagrass

Genomic variability did not differ greatly between populations with all sites displaying very low heterozygosity and a high inbreeding coefficient (F_IS = 1; Table 4). However, sites in cluster one did exhibit slightly higher levels of variability than those in cluster two (Table 4). With such high inbreeding coefficients, it is likely that this species does indeed rely heavily, if not solely, on clonal growth and vegetative reproduction, rather than sexual reproduction (Tables 3 and 4). In terms of reproductive strategy, clonality in seagrasses can vary between species with a continuum from monoclonality to meadows with high clonal diversity (Van Dijk & Van Tussenbroek, 2010), and the predominance of certain clonal lineages may indicate long-term selection on phenotypes. Such selection may be in response to environmental variables, where conditions are more favourable for clonal lineages, but may also represent shared ancestry prior to historic sea-level fluctuations reshaping the topography of the South African coastline (Ramsay & Cooper, 2002; Compton, 2011; Ludt & Rocha, 2015). In the case of Z. capensis, it is likely that a combination of historical (see below) and contemporary factors shape the patterns of observed genomic diversity. Importantly, Z. capensis is unlikely to be influenced by contemporary gene flow between its fragmented and isolated meadows, being restricted to sheltered and low wave action environments (Van Niekerk & Turpie, 2012). In addition, the lack of recorded sexual reproduction in this species through flowering (McMillan, 1980; D Pillay, pers. comm., 2014) is likely to contribute to maintaining clonal populations throughout the range, with important implications for potential restoration efforts in the region.

Shared adaptive divergence across two genomic clusters shaped by historical dynamics

Although neutral variation can reveal much about a species demographic history, in many cases, patterns revealed from outlier loci can provide unique insights into evolutionary potential and patterns of resilience (Stapley et al., 2010; Guo et al., 2015; Funk et al., 2016; Gaither et al., 2018). Particularly in marine systems where gene flow is generally presumed to be high, signals of outlier loci can help detect population structure (André et al., 2011; Freamo et al., 2011; Hess et al., 2013; Candy et al., 2015; Araneda et al., 2016; Tigano & Friesen, 2016; Attard et al., 2018). For example, Atlantic herring in the Baltic and North Sea (André et al., 2011), Atlantic salmon in eastern Canada (Freamo et al., 2011), and Chilean blue mussels (Araneda et al., 2016), all exhibit little to no structure in terms of neutral variation but reveal significant population structure for outliers putatively under selection. In Z. capensis, despite the generally low levels of genomic variation detected across its range, clustering analyses revealed differentiation of all populations into two major clusters when outliers that may represent putatively adaptive variation were considered in addition to neutral data (Fig. 1B). Cluster one comprised of sites from the west and south coast, and cluster two sites from the east coast in addition to one west coast site at Langebaan, Oostewaal, also grouping with this cluster (Fig. 2). In addition, PCoAs support the clustering analyses, with both neutral and adaptive variation (Fig. 2). Interestingly, the PCoAs indicate that east coasts sites are more distantly related than the west and south coast sites. Therefore, east coast sites may have had an earlier origin which supports the likelihood of a refugial area for Z. capensis on the east coast.

Temperature-based ensemble models, however, suggest reduced, and more fragmented seagrass habitat along the South African south coast, that is likely to have divided Z. capensis into two clusters with the south-western and east coast possibly representing refugial areas during the LGM, with subsequent dispersal into its present-day distribution. A refugial area on the south-western coast may explain the presence of both clusters in Langebaan Lagoon. Notably, this split between the clusters roughly coincides with the split between described temperate and sub-tropical bioregions (Sink et al., 2012) along which phylogeographic breaks have been recorded for marine coastal species (Von der Heyden, 2009; Teske et al., 2011), including one other saltmarsh plant (Potts, Veldkornet & Adams, 2016). Given lack of gene flow and apparent high levels of clonality in Z. capensis, the structure detected through outlier loci most likely reflects ancestral adaptation during conditions more conducive to gene flow or incomplete lineage sorting during post-LGM expansion.

Interestingly, although our historical models were based on environmental variables, they broadly mirror the findings of changes in topology and composition of the South African coastline during the last 70,000 years (Toms et al., 2014). During the past 45,000 years, lowered sea level stands of up to −120 m caused significant shifts from rocky to sandy/muddy shorelines which isolated populations of obligate rocky shore species. Although Z. capensis is found in present-day estuaries, the latter would also have been affected, although the extent of change is unknown. Our findings, in combination with Toms et al. (2014), do however show that combinations of abiotic changes have the potential to influence the population dynamics of marine and estuarine species in the Atlantic/Indian Ocean transition zone, as they have done in Australia over exceedingly short time scales (Puritz et al., 2012).

We provide evidence for a pattern of shared outlier loci for populations across a distinct environmental gradient and large geographic span, despite the presence of two genomic clusters that appear to represent distinct historical lineages. This is in contrast to our hypothesis of distinct population-level signals of adaptive variation as seen in previous studies (Williams & Oleksiak, 2008; Perrier et al., 2013; Ravinet et al., 2016). All outlier loci were shared among sites with no populations harbouring private outliers (Table S3), suggesting the same genomic basis for each Z. capensis meadow with the same suite of genes under selection across sites in response to the various environmental gradients. However, differences in the frequencies of these outlier loci across sites provides the foundation for the two clusters, with sites from cluster one exhibiting outlier loci at lower frequencies compared to cluster two (Fig. 3). Notably, this pattern of differential outlier allele frequencies could be observed even when only considering three outlier loci (Fig. 3). The same clustering pattern is detected with both non-simulated and simulated datasets, which were necessary in order for the analysis of pooled data with certain software. This similarity demonstrates that simulated datasets can be used to detect biologically significant evolutionary patterns, regardless of the over-simplifications these simulated datasets may introduce, or the number of SNPs one chooses to employ.

There have been numerous other studies that also report high levels of shared adaptive variation across sites, such as that in Atlantic salmon of eastern Canada, where the allele frequencies of shared outlier loci were used in to assign individuals to their region of origin, assisting with stock management (Freamo et al., 2011). Similarly, in Pacific and Atlantic sticklebacks different allele frequencies of shared outlier loci have been used to distinguish marine and freshwater populations (Jones et al., 2012). At a smaller scale in western Canada, most outlier loci in sticklebacks were specific to single watershed regions (Deagle et al., 2012). Likewise, few shared outlier loci were observed in the periwinkle, Littorina saxatilis, in Sweden (Ravinet et al., 2016), suggesting that the shared or private nature of outlier loci might be highly context specific. Despite the potential for high levels of gene flow and similar selective pressures, L. saxatilis populations displayed a considerable amount of unshared genomic divergence, possibly due to complex polygenic traits involved in habitat adaptation.

The contribution of IBD and IBE towards the spatial arrangement of genomic variability in Z. capensis

Despite the low probability of connectivity between sites, due to both the isolated nature of estuaries and the lack of sexual reproduction recorded for this species, geographic distance (IBD) was not a significant driver of the observed genomic variation, with some evidence for IBE in this system. However, because of spatial autocorrelation with environmental variables chosen in this study, there appears to be a large spatial component shaping genomic variation, which cannot be separated from the effect of environmental variables. Our results suggest that IBE plays an important role in shaping genomic variation in this seagrass, in particular dissolved oxygen, annual mean moisture, precipitation and temperature related environmental variables, that were all significantly associated with outlier loci. Although the functional relevance of these outliers is unclear, this may indicate some level of adaptive variation. In addition, it is likely that fine-scale environmental variation, specific to each of the estuaries from which Z. capensis were sampled additionally contributes to IBE. However, given the lack of in-situ measurements of important environmental variables for coastal South and eastern Africa, this is not possible to determine at this stage.