Phylogenetic analysis
We inferred the phylogenetic relationship between individuals, using the permissive and stringent datasets (“69inds_40MD”, “64inds_20MD”).
First, we estimated a maximum likelihood (ML) tree, concatenating all loci. We used fasta files to estimate the most likely evolution model in ModelFinder (Kalyaanamoorthy et al., 2017), and computed a ML tree in IQ-Tree v.2.1.2 (Minh et al., 2020). To evaluate the support values of the inferred topology, we used the ultrafast bootstrap approach (Hoang et al., 2018) with 1,000 replicates. We ran the analysis five times, chose the tree with the highest likelihood, and visualized it with ggtree v. 2.4.2. (Yu et al., 2017).
Second, we inferred a coalescent tree that incorporates incomplete lineage sorting, using the SVDquartets (Chifman & Kubatko, 2014) and gQMC (Avni, Cohen & Snir, 2015) algorithms implemented in tetrad, within ipyrad. This approach uses one randomly chosen SNP per locus for inferring the topology between any combination of four individuals and joins them into a super tree. In addition to the permissive and stringent datasets, we also ran tetrad on the “full” dataset, using 100 non-parametric bootstraps, and constructed a majority rule consensus tree.

Population summary statistics

To assess if coastal and island lineages of O. insularis differ in their levels of genetic diversity, we considered both the permissive and stringent datasets (“69inds_40MD”, “64inds_20MD”).
First, we measured the following individual-level statistics: 1) expected (He) and observed (Ho) heterozygosity, representing the average and observed number of polymorphic sites in each individual across ddRAD loci; 2) inbreeding coefficient (FIT), representing the heterozygosity of an individual relative to that observed in the total sampling; and 3) nucleotide difference per SNP (π*), which represents the average number of nucleotide differences between the two haplotypes from one individual. We converted the datasets into fasta format using PGDSpider (Lischer & Excoffier, 2012), calculated He, Ho and FIT in vcftools (Danecek et al., 2011), and π* in DnaSP6 (Rozas et al., 2017). We tested if there was a significant correlation between the percentage of missing data and FIT or π* in each individual, using a linear regression model with R (R Core Team, 2017). To test if population-wide π* and FIT were different between lineages, we used a pairwise Student’s t-test with Bonferroni-Holm p-value correction in R (R Core Team, 2017).
Second, to quantify how genetic diversity is partitioned between lineages, we calculated the pairwise fixation index (FST). We converted the vcf-files to the Arlequin-input format using PGDSpider (Lischer & Excoffier, 2012), and used Arlequin ver. 3.5.2.2 (Excoffier & Lischer, 2010) to compute pairwise FST, considering the six lineages determined previously. We performed 5,000 boostraps to estimate significance from the null hypothesis of no population structure. We assessed how much genetic variation is explained by four levels of population structure (between 2 major regions, within 2 major regions, between 6 minor regions and between individuals), by performing a locus-by-locus Analysis of Molecular Variance (AMOVA, Excoffier et al., 1992). We used 10,000 permutations and excluded loci missing in at least one of the lineages.
Lastly, to understand if the genetic diversity observed within each lineage is consistent with changes in effective population size, we computed Tajima’s D using DnaSP (Rozas et al., 2017). We expect lineages that experienced a recent range expansion to have negative values of Tajima’s D, significantly different from zero (alpha= 0.05).

Demographic analysis

To estimate the past demographic history of O. insularis , we applied demographic modelling with diffusion approximation methods as implemented in the program δaδi (Gutenkunst et al. 2009). First, to assess changes in effective population size within each lineage, we fit one-population models for the lineages containing more than 3 individuals (i.e. N-Coastal, S-Coastal, S-Oceanic). We compared four demographic models that differ in the number of population size changes (none, one, two) and in the mode of change (none, discrete, exponential) (details in Table S13). Second, to assess gene flow between lineages, we performed two-population models for the two pairs of adjacent lineages with a larger number of individuals (i.e. N-Coastal vs S-Coastal and S-Oceanic vs S-Coastal), which provide some power of estimating demographic parameters. We fitted 15 different demographic models (details in Table S14, following Portik et al., 2017) that incorporate different scenarios of presence and directionality of migration (m), time of a split between lineages (T), and change of effective population size (nu).
For both classes of demographic modeling, in order to retain a maximum number of segregating sites, we filtered the raw vcf output file (“full”) to only include individuals from the relevant lineages, as in the more stringent dataset (“64ind_20MD dataset”). We excluded loci with more than 40 % missing data, and retained the first SNP per ddRAD locus, minimizing physical linkage. We converted the filtered vcf tiles to the dadi file format, projected down to the number of variants maximizing the number of segregating sites present (Table S15), and calculated the observed one- or two-dimensional folded site-frequency spectra (SFS). We optimized these models using the pipeline developed by Portik et al. (2017) with the following settings (one population modelling/ two population modelling): three/four rounds of optimization with 10, 20, 30/ 15, 10, 5, 5 replicates in each round; with maximum iterations of 5, 10, 50 / 15, 15, 25, 50 per replicate in each round; and a parameter fold of 3, 2, 1 / 3, 2, 1, 1 using the default Nelder-Mead method. To guarantee model convergence, we ran each optimization five times. We selected the model with the lowest Akaike information criterion (AIC) score (Akaike 1974), which penalizes the likelihood by the number of parameters, and estimated dAIC and wAIC (Burnham & Anderson, 2002). For the two-population models, we carried out a goodness-of-fit test of the most likely demgraphic model (Barrat et al., 2018) to test general model fit of our inferred parameters. The test is considered passed if the empirical log-likelihood value falls within the distribution of log-likelihood values fitted to simulated SFS.

Results

Assembly of RAD-loci and filtering

Assemblies with alternative similarity thresholds have equivalent summary statistics (Table S3; supplementary Fig. S1), and thus we chose a similarity threshold of 0.9 to avoid over splitting of loci. Assuming this threshold, the “full” dataset consisted of 299,304 homologous loci, with 572,012 SNPs, with 66,34% of missing data (Table S3); the more permissive “69inds_40MD” dataset consisted of 34,455 loci with 99,915 SNPs, and the more stringent “64inds_20MD” consisted of 25,702 loci with 72,816 SNPs. We found a significant correlation between the number of raw reads and percentage of missing data per individual (adjusted R2: 0.9128, p-value: 2.2 x 10-16).

Population structure

In the EMU-PCA, for the “69inds_40MD” dataset, PC1 and 2 explain 19 and 14% of the genetic variance, respectively. Individuals are clustered into six groups that recapitulate their geographic location (Fig. 1B, Fig. S2A-B): 1) South-coastal (including AL and BA), 2) South-oceanic (the oceanic islands TM), 3) North-SPS (the oceanic islands SPS), 4) North-Caribbean (the coastal OIC), 5) North-Atlantic (the oceanic islands STH, ASC), and 6) North-coastal (including CE, RN, FN, and AR). Results for the more stringent dataset were qualitatively similar (Fig. S2A).
In the TESS3R analysis, using the more permissive dataset (Fig. 2A-B, Fig. S3), the first population split separates samples from the southern coast of Brazil from all the others, with Bahia (BA) and Alagoas (AL) bordering these clusters and containing individuals with admixed ancestry. At K=3, the northern oceanic islands (São Pedro / São Paolo SPS) form a cluster without showing shared ancestry. At K = 4, the southern populations are divided into a coastal (AL and BA) and an oceanic island cluster (TM), with the individuals from Bahia (BA) showing admixture between these two clusters. At K=5, the Caribbean individuals are assigned to a new cluster, showing admixture with the northern-coastal cluster. At K = 6, individuals from the most remote oceanic islands (ASC and STH) form their own unmixed cluster. The exception is the individual from Saint Helena, which contains the largest amount of missing data and therefore is assigned to every cluster with some probability. Assuming further number of ancestral clusters, we find that previous clusters split equally into two clusters, as thus we refrain from interpreting these biologically. Using the more stringent dataset, we obtained overall concordant results (Fig. S4).
For the fineRADstructure analysis (Fig. 2C) with the “full” dataset, the highest degree of co-ancestry was observed in individuals from São Pedro / São Paulo (N-SPS), followed by individuals from Ascension Island that cluster together with Saint Helena (N-Atlantic), and individuals from the Caribbean (N-Carribean). These three groups are nested within a northern group that also includes individuals from N-Coastal localities (Ceara, Rio Grande do Norte) and nearby islands (Fernando de Noronha, and Rocas Atoll). All individuals from the archipelago of Trindade and Martim Vaz (S-Oceanic) have a relatively high co-ancestry and are nested in a southern group encompassing individuals from the S-Coastal localities (Alagoas and Bahia), reflecting hierarchical population structure also in the southern group. Results from runs with the more permissive (“69inds_40MD”) and more stringent (“64inds_20MD”) datasets reflected the same relative levels of co-ancestry and hierarchical population structure (Fig. S5, S6).

Phylogenetic analysis

The ML trees estimated with the most likely model (TVM+F+R3) were largely consistent between datasets. Individuals belonging to the island lineages (N-Atlantic, N-SPS, and S-Oceanic) and to the northernmost coastal lineage (N-Caribbean) form well supported (ultrafast bootstrap > 90) monophyletic clades (Fig. 3A). In contrast, individuals from the geographically widespread coastal lineages (N-Costal, and S-Coastal) form paraphyletic clades, with the respective oceanic clades nested within. The N-Caribbean lineage is the sister of the N-Coastal lineage. The clades leading to the two northern oceanic islands (N-Atlantic and N-SPS) show remarkably long branches, reflecting high genetic divergence. While in the most restrictive dataset these two clusters are sister taxa (Fig. S7, S9), this relationship is not supported in the most permissive dataset (Fig. 3A, S8).
The topology of the coalescent trees was largely congruent with that from the ML trees. Differences include the monophyly of the S-Coastal lineage with high bootstrap support (>90), the paraphyly of the N-Caribbean lineage and the monophyly of the N-Coastal lineage with bootstrap support of >50 (Fig. 3B). Comparing topologies based on the three datasets, we observe that the N-SPS and N-Atlantic lineages form a monophyletic clade in the more permissive dataset (Fig. 3B, Fig. S10), while this monophyly is not recovered in the “full” and more stringent datasets (Fig. S11, S12). Similar to our ML trees, the branch length of lineages from the northern oceanic islands (N-Atlantic and N-SPS) are substantially longer than those of coastal lineages (Fig. S10).

Population summary statistics

Using both datasets, we find that π* and observed (Ho) heterozygosity are lower in the oceanic islands relative to the coastal lineages (Tables S7, S8). Conversely, the inbreeding coefficient FIT is lower in the coastal lineages and highest in all oceanic islands, with N-SPS showing 4 times higher inbreeding relative to the coast (Fig. 4, Fig. S13, Tables S7, S8). For both π* and FIT, all comparisons between coastal and oceanic lineages yielded significant p-values (Tables S5, S6), except the comparisons between S-Oceanic and N-Caribbean (FIT: 0.111, π*: 0.182). We did not find a significant correlation between the percentage of missing data and individual measures of diversity (p-values: 0.324 for FIT, 0.145 for π*), confirming that the observed patterns are not driven by missing data.
For both datasets, FST was significant between all six lineages (Tables S9, S10), except for pairwise comparisons where lineages containing three or less individuals (N-SPS, N-Caribbean, N-Atlantic), reflecting the limited sampling. FST values ranged from 0.056 observed between N-Coastal and N-Atlantic, to 0.703 observed between the two oceanic islands (N-SPS and N-Atlantic) of the northern region. The AMOVA of the more permissive dataset (Table S11) showed that variation was highest within individuals (67.88 %), followed by variation among the lineages within groups (14.62 %), variation explained by the broader regions (North-South, 8.93 %), and finally variation within lineages (8.56%). The AMOVA of the more stringent dataset was qualitatively similar (Table S12).
For both datasets, Tajima’s D (Tables S7, S8) was negative for the three coastal lineages, however it was only significantly different from demographic stability in the N-Coastal lineage. The three lineages of oceanic islands had Tajima’s D values closer to zero.

Demographic history

For the single population modelling, in all cases the neutral model showed the highest AIC values (Table S16), showing that any model accounting for size changes in effective population size (Ne) explains the observed SFS significantly better. For the S-Coastal and S-Oceanic lineages, the simpler “two_epoch” model showed the lowest AIC score (Fig. S14) and similar residuals to other expansion models (Fig. S15). Assuming this model, parameter estimates indicate an increase in Ne of 1.34 (S-Oceanic) and 4.17 (S-Coastal) times, relative to the Ne of the ancestral population “Na” (Table S16, Fig. S16). For the N-Coastal lineage, AIC scores favored the “three_epoch” model incorporating an additional change of Ne. Here, the population first increased to 3.16 times and subsequently to 11.89 times the Na (Table S16, Fig. S16).
For two population modelling, in both cases we obtained the highest AIC values for the neutral model of no population split, rejecting this scenario (Table S17, Fig. S17). In the case of N-Coastal vs S-Coastal, the model integrating a period without gene flow, followed by secondary contact and asymmetric migration (sec_contact_asym_mig_size) shows the lowest AIC value (Fig. S17). Assuming this model, we find that after the initial population split, the S-Coastal population shrinks to 0.0124 times the Na, while the N-Coastal population increases to 2.33 times of the Na. After secondary contact, both effective population sizes increase to 1.3 and 8.46 times of the Na, respectively (Fig. 4A, Table S17). After secondary contact, migration rates are highly asymmetric, being 2.6401 from N-Coastal to S-Coastal and 0.0628 in the opposite direction. For S-Coastal vs S-Oceanic, the model incorporating a population split followed by continuous asymmetric migration (asym_mig, Table S14, S17, Fig S17) showed the lowest AIC values. Assuming this model, S-Coastal increases to an effective population size of 14.1 times the Na after the split, while S-Oceanic increases to 1.19 times the Na. Migration rates are moderately asymmetric at 1.01 from S-Oceanic to S-Coastal and 1.52 in the other direction (Fig. 4B). Residuals for the of the modelled JSFS are higher in the Coastal vs S-Coastal comparison, relative to the S-Coastal vs S-Oceanic comparison (Fig. 4A). Yet, both models passed the goodness-of-fit test (Fig. S18), showing that these are idealized models are fair representations of the evolutionary history acting in these populations.

Discussion

Oceanic currents and depth cause cryptic diversification in O. insularis

Whereas diversification processes in terrestrial habitats have often been linked to geographic barriers inhibiting gene flow, oceanic barriers driving diversification in marine organisms are less well understood, particularly in cephalopods. Our genomic methods inOctopus insularis recovered 299,304 ddRAD loci, containing 572,012 linked SNPs, offering the first insights into the genome-wide patterns of genetic differentiation in this species, and in the oceanographic barriers associated with it.
At a deeper phylogenetic scale, we find that populations of O. insularis are structured into two widely distributed northern and southern groups (Fig. 1B, 2B-C, Fig. 3). This deeper division between the northern and southern groups coincides with the South Equatorial Current (SEC) (Fig. 1A). This finding is in agreement with the distribution of the two major haplogroups in mitochondrial DNA (Lima et al., under revision). Given that O insularis produced up to ~ 95,000 eggs (Lima et al., 2014) and that planktonic paralarvae disperse with oceanic currents (Lima et al., 2017), it is perhaps not surprising to find such a strong role of the SEC in the genetic divergence of this species. This current-mediated North-South division has been reported for other co-distributed species showing a pelagic propagule or larval dispersal, such as corals (Peluso et al., 2018) and mangrove trees (Francisco et al., 2018). Accordingly, a recent review of marine barriers to gene flow, Martins et al. (under review) found the SEC to compose the largest value of phylogeographic concordance among Brazilian coastal organisms, suggesting that this current imposes a major biogeographic barrier across Atlantic species. Studies in other species of octopuses (Octopus vulgaris , Melis et al. 2018; Macroctopus maorum , Higgins et al. 2013) have shown that population structure coincides with oceanic currents in the Mediterranean Sea and the southern Pacific, suggesting that oceanic currents might pose a strong oceanographic barrier for sedentary species with pelagic paralarvae.
At a shallower phylogenetic scale, we find six evolutionarily independent lineages (Fig. 2). The broader northern group is substructured into four lineages: a more widespread coastal lineage encompassing four localities near the coast of Brazil (N-Coastal), a Caribbean lineage (N-Caribbean), a first oceanic lineage encompassing the archipelago of São Pedro and São Paulo (N-SPS), and a second oceanic lineage encompassing the two islands off the coast of Africa (N-Atlantic), Saint Helena and Ascension. In turn, the broader southern group is substructured into two lineages: a more widespread coastal lineage (S-Coastal), and an oceanic lineage representing the archipelago Trindade and Martim Vaz (S-Oceanic). Importantly, these six lineages explain almost double of the genetic variation explained by the two broader groups (AMOVA, Table S11), underscoring their evolutionary significance. This finer population structure is coincident with previous ecological models of habitat suitability for this species that are largely driven by ocean depth (Lima et al., 2020). Oceanic lineages are coincident with volcanic islands that provide suitable habitat highly isolated from the continuous habitat along the coastal shelf of the American continent, with the exception of the S-Coastal lineage, which is connected to the southern coast of Brazil through a chain of seamounts. Genetic differentiation between populations from oceanic islands and the Brazilian coast have also been reported for dolphins (Oliveira et al., 2019), rockpool blennies (Neves et al., 2016), posobranch gastropods (Barroso et al., 2016) and corals (Peluso et al., 2018), confirming that such breaks in habitat constitute important drivers of divergence for coastal taxa with very different dispersal rates.

Island colonization is associated with reduction of genetic diversity

A well-supported observation in terrestrial island biogeography is the lower genetic diversity of recently colonized island populations compared to their mainland source populations (White & Searle, 2007; Boessenkool et al., 2007). Decreased diversity in insular populations leads to increasing inbreeding, decreasing their adaptive capability (Spielman, Brook & Frankham, 2004), and further exposing these to a higher risk of extinction (Frankham, 1997). Yet, it is less clear how genetic diversity changes during island colonization in marine systems, such as in O. insularis.
Our phylogenetic analyses shed some light on the evolutionary relationships between the six lineages (Fig. 3). We show that within the southern group, the S-Oceanic clade is nested within the broadly distributed S-Coastal clade, consistent with a colonization of this oceanic archipelago from the southern coast of Brazil. This direction of colonization is concordant with phylogenetic studies in co-distributed fish species (Macieira et al., 2015; Simon et al., 2021; Pinheiro et al., 2017), which suggest that the colonization of the oceanic archipelago Trindade and Martim Vaz occurred during the LGM (Pinheiro et al., 2017), when seamounts now submerged likely formed an ecological corridor for coastal taxa (Mazzei et al., 2021; Simon et al., 2021). Within the northern group, while the coalescent tree is consistent with a single colonization of the islands from the coast (i.e. N-SPS and N-Atlantic are sister clades; Fig 3B), the ML tree is consistent with two independent colonization events (i.e. N-SPS and N-Atlantic are not sister clades; Fig. 3A). Yet, the ML analysis using the more stringent dataset is again more consistent with a single colonization event (Fig. S7, S9), corroborating our findings with the coalescent methods, which are more appropriate to study recent radiations due to the expected large amounts of incomplete lineage sorting (Giarla & Esselstyn, 2015). Given our finding that these two northern oceanic lineages show the highest genetic differentiation observed in the species (FST: 0.703) and show the longest phylogenetic branches (Fig. 3), more samples would be needed in these islands in order to establish if colonization occurred once or twice within the northern group. Nevertheless, our results conclusively demonstrate that the island populations represent at least two independent colorization events from the coast – one from the South-Coast and one or two from the North-Coast – allowing us to test evolutionary consequences of island colonization for patterns of genetic diversity.
Regarding genetic diversity segregating within the six lineages ofO. insularis , we observe that values of nucleotide diversity (π*) are significantly lower in the three oceanic lineages (N-SPS, N-Atlantic, and S-Oceanic) compared to their respective source of colonization (N-Coastal or S-Coastal; Fig. 4B, Fig. S13, Table S5). Conversely, individual inbreeding coefficients (FIT) are significantly higher in the oceanic lineages (Fig. 4A, Table S6), reflecting the same pattern observed in the co-ancestry matrix (Fig. 2A). These results suggest that island colonization is associated with strong decreases in genetic diversity, both due to the demographic founder effect reducing effective population size, and to the smaller availability of suitable habitat in the islands relative to the coast (Lima et al., 2020); e.g. the available habitat up to 50 m depth in N-SPS is of 1.1 km2 (Ávila et al., 2018). In accordance with our results, similar decreases of genetic variability in oceanic islands have been reported for several reef fish species (Pinheiro et al., 2017). This suggest that, although marine species of low dispersal as O. insularis have been able to successfully colonize novel isolated habitats, such as the volcanic islands of the Atlantic, island colonization has led to a significant decrease of genetic variability in the insular range of the species, mirroring what was established in terrestrial systems (White & Searle, 2007; Boessenkool et al., 2007).

Demographic history is conditioned on oceanic currents and environmental change

Historical changes in the demographic parameters such as effective population size and gene flow are determinant for current patterns genetic diversity within species (Excoffier, Foll & Petit, 2009; Hewitt et al., 2000, Ellegren & Galtier, 2016). Although demographic changes driven by glacial cycles are well known to condition genetic diversity in terrestrial systems (Canestrelli, Sacco & Nascetti, 2012; Cheddadi et al., 2006), it is less clear how such global changes have affected marine systems that are highly dependent on shallow habitats, such asO. insularis .
Our demographic models per population detect contrasting demographic histories for coastal and oceanic lineages. For all the three lineages tested, we conclusively reject a scenario of a stable population size in favor scenarios showing one (for S-Oceanic and S-Coastal lineages) or two (for N-Coastal) instantaneous changes of effective population size. Consistent with relative values of Tajima’s D (Table S7, S8), we estimate that demographic expansions are larger for the two coastal lineages (of 11.96 and 4.17-fold for the N-Coastal and S-Coastal lineages, respectively) than for the oceanic lineage (1.38 fold for the S-Oceanic; Fig. S16). The magnitude of these expansions is consistent with studies of environmental niche modelling that suggested a strong increase of habitat suitability since the LGM along the American coast for this species (Lima et al., 2020), associated with the increase of sea level and the consequent exposure of the continental shelf and increase of sea temperature (Ludt & Rocha, 2015). In contrast, in the oceanic islands covered in this study, increases in sea level likely exposed less coastal habitat (Ávila et al., 2018), explaining the more modest expansion of population size.
Our demographic models for two populations reflect similar changes in effective population size, but reveal patterns of genetic connectivity between lineages. For the populations along the coast (N-Coastal vs S-Coastal comparison), the best model reflects a population split without gene flow, followed by a population expansion with gene flow (Fig 4A). Gene flow is strongly asymmetric, with migration from the N-Coastal to the S-Coastal lineage being 42-fold larger than in the opposite direction. This observation is coincident with the Brazil Current (BC), which moves southwards from the range of the N-Coastal lineage (Fig. 1A). This current has been associated to unidirectional gene flow in rockpool blennies (Neves et al., 2016), suggesting that it may be an important driver of diversification of marine organisms with pelagic dispersal associated to reef habitats.
For the S-Oceanic/S-Coastal comparison (Fig. 4B), the best model consists of a population split with instantaneous size change followed by constant migration. Although migration rates are asymmetric, migration from the coast to the island is only 1.5-fold larger than in the opposite direction. This mildly asymmetric gene flow is consistent with cyclonic eddies that lay north and south of the sea mountain chain connecting the coast with Trindade and Martim Vaz (Arruda et al., 2013, Mill et al., 2015), likely facilitating gene flow in both directions (Pinheiro et al., 2015). Together, these findings suggest that oceanic current not only work as major oceanographic barriers discussed above, but also mediate genetic connectivity between genetically isolated populations.

Implications for conservation biology

Our findings on the evolutionary processes shaping the diversification of O. insularis have direct implications for the conservation of this species. As the taxonomic status of the octopus species mainly targeted by American fisheries has been resolved only recently, little is known about the catch compositions of octopus fisheries in this continent. Annual catches of O. americanus Monfort, 1802 are around 15,000 t in Mexico and 2,000 t in Brazil (Jereb et al. 2014). Although O. insularis has a distribution similar to O. americanus , the reported catch of O. insularis species are much lower (around 500 t; Haimovici et al., 2014). Considering that O. insularis was described only in 2008, and since then new studies have identified a greater distribution range than originally thought, the conservation status of its stocks remains unclear, mainly due to the common misidentification between both species. With the reduction of pelagic fish stocks due to overexploitation, there is a tendency to target fishing towards cephalopods (Rosa et al. 2020), which was indeed observed by Lopes et al. (2020) on O. insularis fisheries over the last 10 years in Northeast Brazil. Thus, identifying different stocks within this species is crucial to propose management strategies that avoid overexploitation of this important fishery resource.
Our findings on population structure imply that management plans forO. insularis must consider at least six evolutionarily independent units with confined distributions (Fig. 2B). Moreover, our finding of significantly higher levels of inbreeding and lower genetic diversity in the island lineages relative to the coastal ones (Fig. 4) imply that three island lineages deserve a higher protection status relative to the three coastal ones. Currently, the insular lineages receive some protection status (Marine Protection Atlas, accessed on 24.01.2022): N-SPS is part of a Brazilian Marine Protected Areas (MPA of São Pedro/ São Paulo), S-Oceanic is a protected Brazilian military area, and N-Atlantic is partially (Marine Protection Zone of Saint Helena) or fully protected (MPA of Ascension). Apart from Rocas Atoll, which is a fully protected area, the other archipelagos are only partial no-take areas (Giglio et al., 2018), being exploited mainly by artisanal fishing in part of the territory and potentially exposing the entire stock of the local lineage of O. insularis . In other species, decreased genetic diversity of natural populations has been associated with breakdown of multiple life history traits (Mills et al., 2012; Reed & Frankham, 2003; DeWoody et al., 2021), a hypothesis that needs to be evaluated in this system by future studies. Yet, given that the adaptive potential of natural populations directly depends on the amount of genetic diversity segregating within populations, our finding of low genetic connectivity and high inbreeding in N-SPS and N-Atlantic imply that these populations deserve a higher protection status. We thus advocate that the entirety of these protected areas be classified as areas of no-take as new conservation measures of O. insularis.Restricting fishery activities to the more genetically diverse coastal lineages, educating local stakeholders on the morphological differences between sympatric octopus species, and establishing set quotas for landings of O. insularis will favor the sustainable management of this economically and ecologically important species.

Acknowledgements

This work was supported by the Brazilian National Research Council (CNPq 481492/2013- 9) and Coordination for the Improvement of Higher Education Personnel (CAPES, Projeto Ciências do Mar II-23038.004807/2014-01). They provided a grant for FDL and logistic and financial support to TSL, FDL and SM. We are thankful to the Chico Mendes Institute for Biodiversity Conservation (MMA/ICMBIO) and to the Brazilian Navy for logistics support during the field work.

Data Accessibility

Raw reads sequenced for this study were deposited in the dryad archive (doi: XXX). All the bioinformatic scripts are available in a public GitHub repository https://github.com/casparbein/OctopusInsularis.

Author Contributions

FDL, SML, and RJP conceived the project. FDL and TSL performed the specimen sampling and HL and LAR produced the genomic data. BB and RJP analysed the data and interpreted the results. BB and RJP wrote the first version of the manuscript, and FDL, TSL and SML contributed to the final version.