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.