Using genomics and morphometrics to monitor data-poor and commercially exploited octopod populations

Over 150 species of benthic octopods have been described within the ‘catch-all’ Octopus genus (Family: Octopodidae) and yet, many Octopus species harvested by fisheries remain unidentified to species-level due to a lack of distinguishing traits. Within species, there is also limited information on how populations differ genetically and the level of connectivity between populations. Therefore, we sampled octopods from commercial fisheries in southeast Australia, in order to identify the species, examine the phylogeographic relationships among species and the level of population genetic structuring within species, as well as to look for any adaptive genetic variation. The mitochondrial gene, cytochrome oxidase subunit III (COIII), was sequenced in 346 octopods along with single nucleotide polymorphisms using double digest restriction site-associated DNA sequencing (ddRADseq). Morphometric traits were also measured in mature specimens. The southern keeled octopus (‘Octopus’berrima) and pale octopus (‘Octopus’pallidus) were identified using COIII data. For ‘Octopus’berrima, we found that some populations whilst being morphologically similar were genetically distinct. In contrast, ‘Octopus’pallidus populations were both morphologically and genetically distinct across the studied regions. Our results provide key information to better inform conservation and management decisions for developing octopod fisheries in southeast Australia and highlight the importance of genomics tools in the conservation management of commercially and recreationally important species.


Introduction
Advances in technology and increasingly affordable sequencing have made genomic resources easily available for non-model organisms, including species of conservation concern. This increasing accessibility opens up the field of conservation genomics for wider use in applications such as captive breeding, reintroduction programmes, and population monitoring. While these genomic resources are usually applied to species under threat, importance should also be placed on species with unclear conservation status, especially those that are commercially harvested. Without fundamental genetic knowledge of an exploited species, genetic diversity could be lost. This loss of diversity could result in a reduction of overall fitness and persistence of a commercially harvested population especially as these threats are exacerbated by the ongoing climate crisis. Given the potentially positive impact from conservation actions and interventions (Hoffmann et al. 2010;McGowan et al. 1 3 119 Page 2 of 21 2016), it is essential that conservation genomic research is also applied to species that are potentially under threat from anthropogenic pressures.
Benthic octopods of the family Octopodidae (Mollusca: Cephalopoda) have an extensive distribution worldwide and a high commercial value (Lima et al. 2017). Despite extensive revisions over the past decades, the taxonomy of this family is still in constant flux with new species being discovered much faster than they are formally described (Norman and Hochberg 2005;Reid 2016). Within Octopodidae the 'catch-all' Octopus genus is mainly targeted by octopod fisheries, which are an emerging market worldwide as global fish stocks are declining due to climate change and overharvesting (Sauer et al. 2019). In fact, catch statistics place the 2018 global capture production of octopods at 377,358 tonnes, worth US$1.78 billion in exports (FAO 2020). This worldwide increase in octopod fisheries, including in Australia, comes despite limited information on which species are harvested (Lima et al. 2017), thus posing a risk of overharvesting genetically distinct individuals from any one species group.
A diverse range of octopod species are harvested in Australia (Martino et al. 2021), however, octopod taxonomy is poorly understood for most commercial species, and species can be difficult to identify from visual methods alone. As such, the species composition of commercial octopus catches is often unknown or uncertain. More research is thus required to verify species using genetic approaches, as well as to determine their phylogeographic relationships. This unresolved taxonomy has thus resulted in most species being placed within the Octopus genus despite this genus being polyphyletic and requiring major revision (Guzik et al. 2005;Reid 2016). Many species currently placed within this genus are similar in structural morphology, making their identification and resolution of evolutionary relationships even more challenging (Norman and Hochberg 2005). Combining both molecular methods and morphometrics has thus allowed unresolved taxa as well as the identification of cryptic species to be tackled Guerra et al. 2010;Amor et al. 2014). Mitochondrial DNA (mtDNA), such as cytochrome oxidase III (COIII), has commonly been used as a molecular marker in marine organisms, as it is useful for inferring evolutionary relationships and distinguishing closely related species (Fadhlaoui-Zid et al. 2012). In octopods, COIII has been used to delineate species and genera , as well as the phylogenetic structure of species (including those from the Octopus genus) (Guzik et al. 2005;Allcock et al. 2008), and the genetic structure of Octopus vulgaris (Fadhlaoui-Zid et al. 2012). COIII is thus used in the current study to determine which species are being commercially harvested in southeast Australia, and to infer their evolutionary history.
Within each octopod species, there is also limited information on how populations differ genetically and at what level gene flow occurs. Despite spanning vast geographical distances, the seascape is typically characterized by a lack of conspicuous physical barriers to gene flow. For instance, oceanic currents are known to contribute to population connectivity in larval organisms (Palumbi 1994;Tepolt et al. 2009;Timm et al. 2020). However, these barriers can also exist in non-physical forms, such as seasonally (Palumbi 1994;Rueda et al. 2013), environmental gradients (Palumbi 1994;Nanninga et al. 2014) and dispersal abilities (Higgins et al. 2013). The reproductive strategies of octopods may have a direct influence on population connectivity. Holobenthic octopods produce large, well-developed hatchlings (without planktonic stages) that immediately adopt a benthic lifestyle, in contrast with merobenthic species that produce small planktonic hatchlings that drift with the currents and settle later in life. The planktonic phase of these free-swimming paralarvae of merobenthic species are not well-known but are temperature-dependent and estimated to last for weeks or months. For example, Octopus rubescens paralarvae remain in the plankton for 1-2 months while captive Octopus vulgaris Type IV paralarvae remain planktonic for about 40 days (Jereb et al. 2016). Holobenthic species typically have lower rates of dispersal and thus gene flow, resulting in higher levels of population structuring compared with merobenthic species (Higgins et al. 2013;Morse et al. 2017). Since accurate descriptions of population structuring have direct implications for conservation and fisheries management, our study applied next-generation sequencing technology, specifically double digest restriction site-associated DNA sequencing (ddRADseq), which can target thousands of single nucleotide polymorphisms (SNPs) even for non-model species like octopods (Morse et al. 2017;Amor et al. 2019;Baden et al. 2023). Compared with traditional markers like microsatellites and mitochondrial DNA, these genome-wide SNPs provide higher-resolution data that enables the detection of finer-scale population patterns, thereby allowing us to correctly define populations at the genetic level to prevent local extinction and conserve population-level genetic diversity (Benestan et al. 2015;Andrews et al. 2016;Miller et al. 2019;Nowland et al. 2019). To the best of our knowledge, SNPs have been used in only three octopod studies to date-assessing the population genetics of the Southern blue-ringed octopus (Hapalochlaena maculosa) using DArTseq (Morse et al. 2017) and of the Brazil reef octopus (Octopus insularis) using ddRADseq (Bein et al. 2022), phylogenetics of Hapalochlaena using DArTseq (Whitelaw et al. 2023), and resolving the common octopus (Octopus vulgaris) species complex using ddRADseq (Amor et al. 2019).
This study focuses on evolutionary relationships at the population level. Our main objectives are to: (1) genetically identify octopod species harvested in southeast Australia and to determine their phylogeographic relationships; (2) determine key external morphometric traits and assess their ability to distinguish these species; (3) characterize the population genetic structuring of these species across their distribution range; and (4) search for any local adaptation by identifying putative loci under selection. This study presents the first genome-wide data for the conservation and fisheries management of these species across southern Australia.

Specimen collection
Octopods were collected between August 2019 and November 2020 from research fishing operations within the South Australia (SA) commercial octopod fishery or purchased directly from commercial fishers in Victoria and Tasmania (see Table S1 for sampling information). These fisheries use unbaited shelter pots attached to demersal longlines (Leporati et al. 2009;Emery et al. 2016;Martino et al. 2021). Sampling locations comprised five localities in SA and one each in Victoria and Tasmania (Fig. 1). Upon capture and euthanasia, the specimens were frozen at − 20 °C. They were then thawed in batches for morphological measurements and dissection for DNA analyses. For the latter, arm tips were dissected and stored in 96% ethanol at − 80 °C without any fixing.

Morphometric trait measurements and analyses
We categorised the reproductive stages following de Lima et al. (2014) and Kivengea et al. (2014). We also measured eight traits from males and females: whole weight (WWt), intraocular distance (IOD), dorsal mantle length (ML), mantle width (MWd), length of third arm to beak (LA3B), width of the largest sucker (LSk), width of the stoutest arm (WdStout) and length of the deepest web sector (LDWebS) based on Norman and Sweeney (1997) and Amor et al. (2014). We measured three additional traits from males: length of hectocotylus (LHec), length of ligula (Lig), and length of calamus (Cal). Measurements were made according to taxonomic descriptions in Reid (2016). Specimen photos are provided in Fig S2. All analyses were conducted in R version 4.0.2 (R Core Team 2020). Statistical significance was indicated by a p value less than 0.05. We used Principal Component Analysis (PCA) to investigate if mature (reproductive stage three) octopods could be identified to species-level and population-level by using MWd, LA3B, LSk, WdStout, and LDWebS. LHec, Lig and Cal were also included for male octopods. We standardized these traits to remove the effect of size by dividing all traits by the ML. These data were also normalised. The PCA analyses were conducted using the prcomp function and ggplot2 package (Wickham 2016).
Based on the clustering from the PCA results, we conducted a decision tree analysis for a subset of the data using RapidMiner Studio 9.8.0 with the Classification and Regression Trees (CART) method. To assess if the individuals can be morphologically identified to species and determine the key distinguishing traits of each species, we split the data into two datasets (a pool of all females and a pool of all males), also using eight and eleven morphological traits respectively. To determine if each species can be morphologically assigned to a population and determine the key distinguishing traits of each population, we split the data into four datasets ('Octopus' berrima females, 'Octopus' pallidus females, 'Octopus' berrima males and 'Octopus' pallidus males). Octopus is a 'catch-all' genus for octopods (Amor et al. 2017), and the genus of these two species will hereafter be referred to in quotation marks throughout the article as we note that 'Octopus' berrima, 'Octopus' pallidus and 'Octopus' australis are phylogenetically distant to the Octopus vulgaris group (Guzik et al. 2005;Acosta-Jofré et al. 2012) and should be assigned to a new genus in future taxonomic research. For each analysis, we randomly split the source dataset at a ratio of 70:30 into test and training datasets respectively (Kotu and Deshpande 2014). The decision tree model was then nested within the 'Optimize Parameters' operator to select the best model parameters for criterion (e.g. information gain, gain ratio), confidence, minimal leaf size, minimal size for split, and maximal depth. Decision trees were thereby generated with the best distinguishing attribute situated on the top of each tree, or as the root node.
In a separate analysis, we also compared the wet weights across localities of each species by generating box plots. The Kruskal-Wallis test was conducted using the kruskal. test and pairwise.wilcox.test functions from base R when comparing more than two localities. The Mann-Whitney U test was conducted using the wilcox.test function from base R when comparing two localities.

DNA extraction, PCR and sequencing
All laboratory work was conducted at the South Australian Regional Facility for Molecular Ecology and Evolution (SARFMEE). We used a high-salt protocol (Aljanabi and Martinez 1997) to extract DNA from a small amount of muscle tissue from the arm tip of each octopod specimen.

Double-digest ddRADseq library preparation
We re-extracted DNA from 134 of the tissues, a random subset of the above samples, using Qiagen DNeasy Blood and Tissue kit, following manufacturer's protocol. 25 μg of RNAse A was then added to each sample for 30 min following overnight digestion. We then quantified the extracted DNA using Qubit (Invitrogen) and Glomax (Promega) fluorometers.
We prepared the double-digest restriction site-associated DNA sequencing (ddRADseq) reactions in two plates of 96, following the protocol of Poland et al. (2012) with some modifications. We randomly selected a subset of samples to be replicated within each plate. The first plate consisted of 88 samples (wells A1 to H11) and 8 replicates (wells A12 to H12). The second plate consisted of 46 samples (wells A1 to F6) and 10 replicates (wells G6 to H7). The remaining 40 wells contained samples from other projects and were excluded from this study. 200 ng of DNA was digested at 37 °C for 2 h using 8 U of PstI (six-base recognition site, CTG CAG ) and HpaII (four-base recognition site, CCGG) in 20 μL of 1 × CutSmart Buffer (New England BioLabs [NEB]).
We ligated barcoded adapters, unique to each plate of 96 samples, to the DNA in 40 μL consisting of 20 μL of digested DNA, 200 U of T4 ligase, 0.1 ρmol of forward (rare) and 15 ρmol of reverse (common) adapters, and 1 × T4 Buffer. Upon incubation for 2 h at room temperature, the mixture was then heated at 65 °C for 20 min to inactivate the enzymes. Within each plate, we pooled the ligation products into batches of 24 samples, which were then purified using the QIAquick PCR purification kit (Qiagen) and eluted in 120 μL of EB buffer (Qiagen) per batch.
We added full-length Illumina adapters via polymerase chain reaction (PCR), and to avoid creating a bias, we split each batch of 24 samples into eight replicates, with 10 μL of purified library in each replicate tube (leaving 40 μL per batch as a back-up). Each PCR tube contained 30 μL volumes containing 10 μL of purified library, 1 × Hot Start Taq Master Mix (NEB), and 0.66 μM each of the forward and reverse primers. The PCR conditions were: 95 °C for 30 s, 16 cycles of 95 °C for 30 s, 65 °C for 20 s, and 68 °C for 30 s, followed by 68 °C for 5 min, and 25 °C for 1 min. The eight PCR replicates per batch were re-pooled and purified using the QIAquick PCR purification kit as above, eluting in 30 μL of EB buffer (Qiagen). We then pooled all four libraries per batch into one, totalling two libraries for all samples. We used a two-step double-SPRI manufacturer's protocol to select for fragments between 100 and 300 bp using AMPure XP bead (Beckman Coulter). Both libraries were quantified using Tapestation 2200 (Agilent), and were sequenced in two runs of 2 × 150 bp (paired-end) with a 25% spike-in of PhiX on the Illumina HiSeq at Genewiz in China.

Phylogenetic tree and haplotype network construction (mitochondrial genetic dataset)
We trimmed the COIII sequences for quality in CodonCode Aligner 9.0.1 and imported the consensus sequences into the Basic Local Alignment Search Tool (BLAST) query search for species identification based on the lowest E-value and highest percent identity. The best multiple sequence alignment was obtained on GUIDANCE2 (Sela et al. 2015) using the MAFFT algorithm via the CIPRES Science Gateway (Miller et al. 2010). We then uploaded this alignment onto FaBox 1.5 DNA to Haplotype Collapser and Converter to obtain haplotype sequences (Villesen 2007). These were then used to generate a phylogenetic tree in IQ-TREE Web Server (Trifinopoulos et al. 2016) using the ultrafast bootstrap (UFBoot) approximation (Hoang et al. 2018) with 1000 replicates. The sH-aLRT test (Guindon et al. 2010) was also performed during the tree generation, such that a clade is reliable when both SH-aLRT > = 80% and UFboot > = 95%. We then edited the tree on the Interactive Tree of Life v6 (iTOL) (Letunic and Bork 2021).
Each species from the MAFFT alignment was subset into separate alignments, examined for missing data, trimmed to a length common across all sequences, and samples with more than 5% missing data were removed (i.e. three 'Octopus' berrima samples). We then constructed a haplotype network using the minimum spanning inference method in PopART (Bandelt et al. 1999;Leigh and Bryant 2015). We obtained the summary statistics on within-sample genetic diversity [number of segregating sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), average no. of nucleotide differences (K)] using DNAsp v 6.12 (Rozas et al. 2017). Among-sample genetic divergence, or fixation index (F ST ), was calculated with 1000 permutations using ARLEQUIN v 3.5 (Excoffier and Lischer 2010). We also examined population demographic histories using neutrality tests, Tajima's D (Tajima 1989) and Fu's F S (Fu 1997) with 10,000 simulations in ARLEQUIN v 3.5. Large negative and significant D, F S and F* indicate demographic expansion.

SNP calling, filtering and processing (nuclear genetic dataset)
Raw sequencing reads were demultiplexed by each sample's unique barcode using GBSX v1.3 (Herten et al. 2015), with BBDuk used to chastity filter the reads, and BBDuk2 used to remove the barcodes, reverse adapter (allowing for 1 sequencing error), remove contaminants longer than 21 bp, trim regions at either end of the read with low quality (Q < 10), and trim the resulting reads to 80 bp in length. Each sample's raw read file was renamed, taking into account the plate and well number, the barcode sequence, the sample name and population, and the species, so that any sample mix-up could be quickly identified. STACKS v1.40 (Catchen et al. 2011(Catchen et al. , 2013 was used to process the ddRADseq data and call SNPs, using parameters recommended by Mastretta-Yanes et al. (2015) to maximize loci and minimize errors, while also exploring other values for these parameters to establish whether they changed the results. ddRAD loci were formed for each sample using the ustacks script, requiring a minimum stack read depth of three (m = 3) and a maximum of two nucleotide mismatches (M = 2) between stacks at a locus. Loci with more than three stacks (mls = 3) and more reads than two standard deviations above the mean were filtered out as they were more likely to map to multiple locations within the genome. Also, the deleveraging algorithm was used to attempt to resolve over-merged loci.
As we wanted to confirm the COIII-identified species, we analysed the full dataset (all 134 samples), before then analysing each ddRADseq-confirmed species separately (refer to methods below). Therefore, a catalogue of consensus loci was generated for the full dataset, as well as each species separately ('Octopus' berrima All Marker dataset = BAM, 'Octopus' pallidus All Marker dataset = PAM) using the cstacks script with the ustacks output files. Loci were considered as homologous within a species (BAM or PAM) or for the full dataset if they mismatched at four or fewer bases (n = 5). Alleles were identified in each individual by mapping back against this catalogue using the sstacks script. To generate a SNP dataset (i.e. only polymorphic markers), after initial processing and SNP calling, we filtered loci with heterozygosity > 0.8 (to remove potential paralogs), with more than 20% missing data and with minor allele frequencies of < 0.05 using the populations script (Fig. 2). The full dataset was output as genepop and plink format files for further processing to confirm the species designation. As a preliminary check on the full dataset, Plink v1.90b6.20 was used to identify individuals with more than 60% missing data (these three samples were removed from the dataset) and the similarity between the replicates of each replicated sample was confirmed using their closeness on a PCA.

Population assignment and structure (nuclear genetic dataset)
To examine the population genetic patterns from the 'Octopus' berrima and 'Octopus' pallidus species, we focused on only the neutrally evolving loci. We did this by first confirming the species designations from the full dataset ( Fig. S3), using PCA, AMOVA and fixed difference analysis. PCA was run using the scaleGen function in the Adegenet package in R (Jombart and Ahmed 2011) and then the dudi.pca function in the ade4 package in R (Dray and Dufour 2007), using the mean allele frequency to replace any missing values. AMOVA was conducted using the gl.amova function from the dartR package (Gruber et al. 2017), whereas fixed difference analysis was run using the dartR package. Then we explored each species dataset individually (BAM and PAM) To do this, the species-specific SNP datasets were rerun through cstacks, sstacks and populations STACKS scripts separately and output as genepop and plink format files, after removing individuals with more than 60% missing data at the species level (BAM or PAM). These steps resulted in an initial 'Octopus' berrima dataset of 77 individuals that comprised 6817 loci and an 'Octopus' pallidus dataset of 54 individuals that comprised 4690 loci). The R packages Adegenet (Jombart and Ahmed 2011) and ade4 (Dray and Dufour 2007) were used to import the genepop file for each dataset, convert it to a genind object and genlight objects for further analysis.
Each species dataset (BAM and PAM) was then examined to establish that the standard population genetic model assumptions were not violated within each population, i.e.
(1) obeying HWE (any loci not under HWE at the 5% level were identified using plink, and the Stacks populations script was re-run after removal), (2) random mating (i.e. any individuals with kinship values > 0.1 were identified using the R package hierfstat (Goudet 2005) and removed from the output files), and (3) only neutrally evolving loci (i.e. any putative loci under selection were identified using a combination of BayeScan and PCAdapt in the R package pcadapt v4.3.3, and the Stacks populations script re-run after removal; see methods below). We removed 2300 loci and 1066 loci not in HWE in the dataset ('Octopus' berrima and 'Octopus' pallidus, respectively), and eight 'Octopus' berrima samples (EP2103, EP2106, EP2157, MD2151, MD2156, VB1104, VB2116, and VB2154) and two 'Octopus' pallidus samples (PL2040 and PL2029) were possibly half-siblings with kinship values > 0.1 and so were removed. Then 349 loci identified as putatively under selection (in either BayeScan or PCAdapt analysis) from the BAM dataset and 387 loci identified as putatively under selection (in the PCAdapt analysis) from the PAM dataset, were removed as well. These loci under selection will hereafter be referred to as the BSM ('Octopus' berrima Markers under Selection dataset) and PSM ('Octopus' pallidus Markers under Selection dataset). Removal of loci under selection then resulted in an 'Octopus' berrima dataset of 4168 loci and 69 individuals ('Octopus' berrima Neutral Marker dataset = BNM), and an 'Octopus' pallidus dataset of 3235 loci and 52 individuals ('Octopus' pallidus Neutral Marker dataset = PNM) (Fig. 2). The R package dartR (Gruber et al. 2017) was then used to convert the genind object to a genlight object.
We explored each species' dataset individually using PCA, admixture plots using the software sparse non-negative matrix factorizations (sNMF), and Structure to identify populations. sNMF allows population structures to be inferred at much shorter runtimes than Structure without loss of accuracy and like PCA, it is a model-free approach that is robust to departures from classical population genetic assumptions (Frichot et al. 2014). Structure, on the other hand, is a Bayesian modelling approach with prior assumptions of the data (e.g. absence of genetic drift, ancestral populations are in Hardy-Weinberg and linkage equilibrium) (Frichot et al. 2014;Zurn et al. 2022). Hence, three methods (PCA, Structure and sNMF) were adopted to obtain a clearer picture of the population structures. The PCA was generated using the scaleGen function in the Adegenet package in R (Jombart and Ahmed 2011) and then the dudi.pca function in the ade4 package in R (Dray and Dufour 2007), using the mean allele frequency to replace any missing values. An sNMF admixture coefficient analysis was conducted for 1000 runs of subpopulations (K) ranging from 1 through 10 with parameters: 500 iterations, tolerance = 0.0001, α = 10 and 5% masked genotypes (Frichot et al. 2014). Structure was run at K = 1-10 for 10 repetitions each and for 100,000 Markov Chain Monte Carlo (MCMC) iterations after removing 20,000 iterations of burn-in. The groups identified from the PCA and sNMF analyses were further examined by estimating the pairwise fixed differences between the groups using the R package hierfstat. Both PCA and sNMF analyses are robust to departures from standard population genetic model assumptions. These analyses allowed us to define each species' populations and meta-populations for further examination.
Statistics for the BAM and PAM datasets as well as the BNM and PNM datasets, such as fixed differences, private alleles, allelic richness, F ST , and identity-by-descent were generated using the dartR package. Tajima's D was calculated for each population, meta-population and species using the TajimaD function in the R package r2vcftools (Danecek et al. 2011;Pope 2019), which is a wrapper that uses vcftools to calculate Tajima's D but performs simulations from the neutral model to conduct a significance test (with 10,000 Fig. 2 Schematic of filtering and processing steps for nuclear dataset (single nucleotide polymorphisms [SNPs] and autosomal) for 'Octopus' berrima and 'Octopus' pallidus, from the full dataset to subse-quent species-specific datasets containing markers under selection or neutral markers, each with different sets of criteria. Blue-grey arrows show the resulting plots and tables generated from each dataset simulations). We also calculated the effective population size for each population, meta-population and species using the gl.LDNE function in the R package dartR, which is a wrapper for NeEstimator v2, with default parameters.
As levels of heterozygosity are influenced by missing data, Plink v1.90b6.20 was used to remove loci with any missing data from the BNM and PNM datasets, which resulted in 175 loci for the BNM NoMiss dataset and 247 loci for the PNM NoMiss dataset. Plink was then used to calculate a SNP-based measure of observed and expected heterozygosity. As levels of heterozygosity are also biased by filtering decisions, population structure, and when only polymorphic markers are used (Schmidt et al. 2021), we also generated an autosomal version of each species dataset while taking into account population structure. To estimate this unbiased measure of heterozygosity (i.e. autosomal) (Schmidt et al. 2021), we reran the populations script allowing for no missing data, removing the maximum heterozygosity filtering and minimum minor allele frequency filtering, and including all monomorphic bases for each marker (not just the polymorphic SNP sites). This resulted in datasets with a reduced number of loci, designated the BAM auto and PAM auto datasets (Fig. 2). Autosomal measures of observed and expected heterozygosity and F IS were obtained from the sumstats summary file from Stacks.

Outlier detection and mapping to genome
We wanted to examine signatures of selection for both 'Octopus' berrima and 'Octopus' pallidus species. To do this, we focused on putative loci under selection that were identified by the BayeScan and PCAdapt analysis i.e. BSM and PSM datasets (Fig. 2). BayeScan was run with default parameters and prior odds of the neutral model of 100. PCAdapt was run using the Mahalanobis method and keeping the first 20 K-values with a p-value of 0.1 using the Benjamini and Hochberg (1995) adjustment method to deal with multiple tests. We generated a PCA (using the same methods as above) to visualise how loci under selection may be influencing populations differently compared to neutral loci.

Nuclear genetic dataset
For the full SNP dataset, a total of 609,146,289 short reads were obtained from the sequencing data of 80 'Octopus' berrima and 54 'Octopus' pallidus (Fig. 2). After demultiplexing and filtering,199,196,920 and 188,900,295 reads were retained, which corresponded to an average [± standard deviation (SD)] of 1,011,151 ± 456,047 and 939,802 ± 430,086 reads per individual (for 'Octopus' berrima and 'Octopus' pallidus respectively). After the initial processing and SNP calling using the STACKS scripts, three 'Octopus' berrima individuals were removed due to large amounts of missing data (> 60%). In the resulting ddRAD SNP dataset (131 samples comprised of 77 'Octopus' berrima and 54 'Octopus' pallidus), the number of markers ranged from 2555 to 6741 (mean 5973 ± 959 SD) for 'Octopus' berrima and 2753 to 4660 (mean 4220 ± 394 SD) for 'Octopus' pallidus, and mean depth coverage of 12.9× (± 7.7 SD) and 11.3× (± 5.7 SD) per individual, respectively. Raw and processed reads were deposited in Zenodo (https:// doi. org/ 10. 5281/ zenodo. 75394 57). From this dataset comprising 2791 loci across all 131 samples, the Structure analysis' optimum K (= 2) also confirmed the two species (Fig. S3a). Likewise, two well-differentiated groupings were evident, based on PCA (Fig. S3b), each containing samples we identified as 'Octopus' berrima and 'Octopus' pallidus via the mitochondrial COIII dataset, with no misplaced individuals. No individual SNP played a major role in separating these two groupings in the PCA (maximum loading is 0.015618; see loadings plot for this PCA in Fig. S5). There were 1858 fixed differences between the species (67%), with 2532 private alleles in 'Octopus' berrima and 2105 private alleles in 'Octopus' pallidus. The AMOVA indicated that 99.9% of the variation was between species and < 1% within species. After confirming the species identities for the dataset, two species-specific SNP datasets were generated separately.

Morphology
Out of 346 specimens in total, 170 were found to be mature-97 out of 143 males were mature, 73 out of 203 females were mature, and 2 out of the 203 females could not be assessed for reproductive stages. Morphometric analyses of the mature specimens also showed strong differentiation between the two species, but within-species variation was only found for 'Octopus' pallidus. Individuals from both species formed two distinct clusters in the PCA plots, more notably in males (Fig. 3). For both sexes, the morphological variables WdStout, LSk and LA3B contributed most to the principal components. Within males, LHec and Lig also contributed strongly to PC1. CART analysis showed that this distinction between species was mainly attributed to WdStout for females and LSk for males, where larger traits likely indicate 'Octopus' pallidus (Fig. S6). Within 'Octopus' berrima, individuals could not be geographically differentiated for either sex (Fig. 4 and Fig. S7). Distinct clusters or populations could be identified for 'Octopus' pallidus which was largely attributed to higher WdStout and LdWebS values in the Victoria/Tasmania populations (Fig. 4). Likewise, CART analysis showed that female 'Octopus' pallidus were largely differentiated by WdStout but this was less clear for males.
Significant differences were not found for wet weights of female 'Octopus' berrima across any localities, but were found in male 'Octopus' berrima between Venus Bay and Tumby Bay, and between Mount Dutton Bay and Tumby Bay (Fig. S8). Wet weights of female 'Octopus' pallidus significantly differed across all three localities (Port Lincoln in South Australia, Victoria, and Tasmania) (Fig. S8). Likewise, male 'Octopus' pallidus significantly differed in wet weights between Port Lincoln and Victoria.

Mitochondrial genetic dataset
Based on the COIII dataset, a total of 18 and 12 haplotypes were identified for 'Octopus' berrima (n = 223) and 'Octopus' pallidus (n = 120). The most common haplotype HB1 in 'Octopus' berrima ( Fig. 5a) was not found in the westernmost SA locality, Venus Bay. All individuals from Venus Bay exhibited 3 unique haplotypes (HB10, HB11, HB12) none of which were found in other populations. In contrast, the most common haplotype HP1 for 'Octopus' pallidus ( Fig. 5b) was found in all sampling sites and was also connected to almost all the other haplotypes by one mutational step.
The lowest haplotype and nucleotide diversities were found in Venus Bay for 'Octopus' berrima and in Port Lincoln for 'Octopus' pallidus, while the highest values were found in Eely Point and Tasmania for the respective species (Table 1). The highest F ST values in 'Octopus' berrima ( Table 2) were found in the pairwise comparisons between the westernmost SA site, Venus Bay, and all other populations (0.702-0.846) while the lowest were between the closest sites to one another, Mount Dutton Bay and Eely Point (0.0607) as well as Tumby Bay and Port Lincoln (0.0429). All F ST values were significant for 'Octopus' berrima except between Tumby Bay and Port Lincoln.
For 'Octopus' pallidus, F ST values were the highest between the two furthest sites, Port Lincoln and Tasmania (0.0969), and lowest between the closest sites, Victoria and Tasmania (− 0.00393). Only the F ST between Port Lincoln and Tasmania was significant. Within populations based on COIII, the test statistic for population demographic history, Tajima's D and Fu's F S , were negative and significant for 'Octopus' berrima in Venus Bay and Tumby Bay and for 'Octopus' pallidus in Victoria, suggesting population expansion. However, species-level Tajima's D was negative but not significant for both 'Octopus' berrima and 'Octopus' pallidus whereas Fu's F S was only negative and significant for 'Octopus' pallidus, indicating no signs of population expansion for both species as a whole. On the other hand,

Nuclear genetic dataset
All population genetic analyses were performed on each species dataset using only neutral loci (i.e. BNM: 69 'Octopus' berrima individuals for 4168 neutral loci, with average depth of coverage: 15× ± 8 SD; and PNM: 52 'Octopus' pallidus individuals for 3235 neutral loci, with average depth of coverage: 13.6× ± 6.2 SD), except where specified (i.e. autosomal loci for calculating unbiased heterozygosity levels). Within 'Octopus' berrima, populations geographically close to each other genetically resembled each other i.e. Eely Point and Mount Dutton Bay (within 11 km), and Port Lincoln and Tumby Bay (within 51 km), as reflected by the lower levels of F ST between these pairs of populations (Table 2). Therefore we consider the Eely Point and Mount Dutton Bay populations as one gene pool (i.e. a meta-population), and Port Lincoln and Tumby Bay as one gene pool (i.e. another meta-population). Similarly, the Victoria and Tasmania populations of 'Octopus' pallidus (within 370 km of each other) also genetically resembled each other with similarly low F ST values. Hence we consider the Victoria and Tasmania populations as one gene pool (i.e. a meta-population). For both species, observed heterozygosity values per population were consistently lower than expected, indicating an excess of homozygosity. This coincides with the positive values of the inbreeding coefficient F IS .
Diversity (observed heterozygosity and allelic richness) and differentiation levels were generally higher in the 'Octopus' pallidus populations compared with that in 'Octopus' berrima, although that was driven mostly by the 'Octopus' pallidus Victoria and Tasmania populations (  'Octopus' berrima did not cluster by region while male (c) and female (d) 'Octopus' pallidus were found to cluster by region and Port Lincoln/Tumby Bay, and 2039-2110 private alleles between Venus Bay and each other population). Similarly, the AMOVA showed 77% of the variation was between the 'Octopus' berrima populations and 23% within 'Octopus' berrima populations compared with 52% between, and 48% within, the 'Octopus' pallidus populations.

Population structure (nuclear genetic dataset)
PCAs were performed on each dataset after loci not under Hardy-Weinberg Equilibrium (HWE) and putatively under selection using BayeScan and PCAdapt, and related individuals were removed. These steps resulted in final datasets of 69 'Octopus' berrima individuals for 4168 neutral loci and 52 'Octopus' pallidus individuals for 3235 neutral loci (average depth of coverage: 15× ± 8 SD and 13.6× ± 6.2 SD, respectively). The PCA for 'Octopus' berrima (Fig. 6a) showed a geographically divergent pattern where the westernmost SA locality, Venus Bay, was clearly different from the other populations by the first axis, which explained 11.8% of the total variation. This was congruent with the COIII haplotype network. Pairs of the closest sites clustered together (Mount Dutton Bay and Eely Point populations as well as Port Lincoln and Tumby Bay populations). The PCA for 'Octopus' pallidus ( Fig. 6d) also showed a geographically divergent pattern where Port Lincoln, the westernmost site of 'Octopus' pallidus sampling range, was clearly different from the other two populations along the first axis, Fig. 5 A minimum spanning haplotype network for a 'Octopus' berrima and b 'Octopus' pallidus generated using Pop-ART with default settings. Size of circles indicate the number of individuals per haplotype which explained 25.5% of the total variation. The Victoria and Tasmania populations showed the greatest dispersion in the second axis, however the lack of clear differentiation between these two populations suggested gene flow could still be occurring between Victoria and Tasmania, or that they share a more recent ancestral connection that has not been completely erased over time. The same patterns were also seen in the sNMF plots (Fig. 6b, e) which showed some levels of admixture in a portion of individuals for both species. Some 'Octopus' pallidus individuals showed admixture between the Port Lincoln population and Victoria/Tasmania metapopulation. The best K was 3 for 'Octopus' berrima and 2 for 'Octopus' pallidus using sNMF. The Structure plot for 'Octopus' berrima had an optimal K = 2, with Venus Bay separate from all other populations (Fig. 6c). When Venus Bay individuals were removed and Structure analysis repeated, K = 2 was the optimal structural split into Port Lincoln/Tumby Bay metapopulation versus Eely Point/Mount Dutton Bay metapopulation (Fig. S9). The Structure plot for 'Octopus' pallidus remained congruent with the PCA and sNMF plots in having an optimal K = 2 (Fig. 6f).
When we examined isolation by distance, the 'Octopus' berrima dataset had significant correlation between genetic and geographic distance (mantel r = 0.8484; p-value = 0.0333), but the 'Octopus' pallidus dataset did not have significant correlation between genetic and geographic distance (mantel r = 0.9849; p-value = 0.1667). We note that our ability to detect correlation for 'Octopus' pallidus was limited by a lack of sampling at intervening localities.

Outlier detection and mapping to genome (nuclear genetic dataset)
BayeScan identified 123 loci putatively under selection for 'Octopus' berrima but no loci under selection for 'Octopus' pallidus. When PCAdapt was also used to identify putative loci under selection, the same 123 loci were identified for 'Octopus' berrima but also another 226 loci (PCAdapt identified a total of 349 loci in 'Octopus' berrima), while PCAdapt identified 387 putative loci under selection for 'Octopus' pallidus. When we examined the PCA of each dataset under selection, only the 123 putative loci identified by BayeScan in 'Octopus' berrima showed a different pattern to that of the neutral loci (Fig. 7). Unlike the PCA for neutral loci in 'Octopus' berrima ( Fig. 6a), the first axis separates the eastern (Port Lincoln and Tumby Bay) and western (Mount Dutton Bay and Eely Point) populations, and the second axis separates the far western population (Venus Bay).
We attempted to map the adapter-trimmed forward and reverse reads to the chromosome-level Octopus reference genome and transcriptome datasets published online (https:// www. ncbi. nlm. nih. gov/ assem bly/ GCF_ 00634 5805.1 for O. sinensis, the East Asian common octopus). However, we did not recover any significant BLAST hits when mapping to either the genome or transcriptome dataset to examine genes that might be driving the signals of selection.

Discussion
We combined morphometric and genetic approaches to identify octopod species and characterise the intraspecific populations of commercially exploited populations. Firstly, we identified two species: 'Octopus' berrima and 'Octopus' pallidus. We show that both species can be morphologically differentiated but key morphological traits were less clear for 'Octopus' berrima than for 'Octopus' pallidus. Population genetic structure was found within each species based on both neutral and adaptive loci, both of which reveal different genetic patterns for 'Octopus' berrima.

Conservation genomics to monitor population structure
The spatial and temporal structures of populations shape the way fishery resources are managed and assessed and yet, the population genetic structures of the commercially important 'Octopus' berrima and 'Octopus' pallidus are still understudied. 'Octopus' berrima is a moderate-sized (up to 35 cm total length) and nocturnal species endemic to southeast Australia, where it is normally found in sand and mud habitats feeding mainly on isopods, and other crustaceans and bivalves (Stranks and Norman 1992;Jereb et al. 2016;Hua et al. 2023). Found in the same geographical range, 'Octopus' pallidus is also a nocturnal species that grows up to a total length of 54 cm with a more robust appearance covered in regular oval patches, and with short, stout arms (Stranks 1988;Jereb et al. 2016). They are also found in sand and mud habitats, feeding on bivalves and other shellfish and crustaceans. Both 'Octopus' berrima and 'Octopus' pallidus are commercially important as they are harvested in developing and/or small scale fisheries in all three Australian states across their geographical range (South Australia, Victoria and Tasmania) (Leporati et al. 2008;Jereb et al. 2016;Conron et al. 2020;Krueck et al. 2021). However, the stock status for 'Octopus' pallidus has been either classified as 'negligible' (in South Australia) or 'undefined' (in Victoria) due to a lack of scientific data, whereas stocks in Tasmania are listed as 'Depleting' (Conron et al. 2020;Krueck et al. 2021). Stock status for 'Octopus' berrima on the other hand has yet to be formally defined for any region. According to IUCN, both species are listed as 'least concern' despite the lack of population data and monitoring while both continue to be commercially harvested. We found genetic structuring was present in the populations of both species across southeast Australia in the nuclear genetic dataset. Both 'Octopus' berrima and 'Octopus' pallidus are holobenthic species, producing well-developed hatchlings that settle immediately on the seafloor therefore some structuring would be expected. Table 2 Pairwise  Table 3 Population-level statistics based on neutral single nucleotide polymorphisms (SNP) loci for 'Octopus' berrima and 'Octopus' pallidus at the population, meta-population, and species level (meta-populations were defined based on the F ST , fixed differences, sNMF and Structure analyses) Observed and expected heterozygosity based on autosomal markers were also included (i.e. with no filtering for minor allele frequency nor maximum heterozygosity, no missing data, and taking into account population structure, and both monomorphic and polymorphic bases). Locality names are abbreviated as follows: This genetic differentiation was not detected at the mitochondrial DNA level except for the 'Octopus' berrima population at Venus Bay. Thus, our study highlights the importance of integrating different molecular methods, whereby COIII served as a rapid tool for species identification and provided phylogeographical information and ddRADseq informed population genetic structures. This study also underscores the ability of next-generation technology to reveal finer-scale population patterns. Octopus maorum, a merobenthic species found in the same region, was shown to be genetically homogenous between South Australia and southern Tasmania, as facilitated by the Leeuwin Current (Doubleday et al. 2009), whereas genetic structuring was found across the range of 'Octopus' pallidus although gene  (Higgins et al. 2013). Hence our data support the conclusion that the life history of octopods is a main factor contributing to the genetic structuring and low level of connectivity between meta-populations. Levels of gene flow in 'Octopus' pallidus were high between Victoria and Tasmania populations, indicating movement of individuals across Bass Strait, ~ 350 km. Individuals in the Victoria/Tasmania metapopulation and Port Lincoln were unlikely to be interbreeding due to their holobenthic nature but more sampling from locations across southern Australia is required to confirm this. High levels of gene flow in 'Octopus' berrima were also found between Port Lincoln and Tumby Bay (metapopulation Port Lincoln/ Tumby Bay), and between Eely Point and Mount Dutton Bay (metapopulation Eely Point/Mount Dutton Bay), suggesting movement of individuals over small spatial scales (~ 50 km) within each metapopulation that maintains gene flow. Both Structure and sNMF suggested that the optimal number of subpopulations was 2-3 for 'Octopus' berrima and 2 for 'Octopus' pallidus. Given that octopods are known to be sensitive to fluctuations in salinity (Vaisakha et al. 2021), differences in salinity could potentially explain the lack of gene flow in 'Octopus' berrima among these two metapopulations and Venus Bay, where Venus Bay has the highest salinity (mean 44 ppt), followed by Eely Point/ Mount Dutton Bay (mean 42 ppt) and lastly Port Lincoln/ Tumby Bay (mean 38 ppt) which are more exposed to an open ocean environment. In fact, Martino et al. (2021) found that 'Octopus' berrima from Venus Bay was the most abundant relative to other populations, and is the first known octopod species thriving in hypersaline conditions. While the association between genetic differentiation and environmental parameters has been observed in other marine species (Limborg et al. 2012;Sjoqvist et al. 2015;Miller et al. 2019;Nowland et al. 2019), to the best of our knowledge, this has not been studied in cephalopods. However, the PCA of loci putatively under selection (Fig. 7) suggests the Port Lincoln/Tumby Bay and Eely Point/Mount Dutton Bay metapopulations may be under different selective pressures. This differed from the purely phylogeographic pattern seen in the neutral loci where Venus Bay was the most divergent. The presence of a density front and summertime thermal front at the mouth of Spencer Gulf, where Port Lincoln and Tumby Bay are situated, limits the exchange of Gulf and oceanic waters (Nunes-Vaz 2014). This, in addition to, the geographical differences in salinity and depth may have contributed to the different clustering pattern in the PCA. Over time, the increasing divergence among conspecific populations as a result of geographical heterogeneity can drive the evolution of cryptic subspecies or speciation (Doebeli and Dieckmann 2003). Whether or not this has or is occurring in 'Octopus' berrima or 'Octopus' pallidus will require more widespread geographic sampling effort. In addition, Tajima's D values using the SNP dataset indicated population expansion for all 'Octopus' berrima populations except for Venus Bay (contraction), and population contraction for all 'Octopus' pallidus populations. However, based on the COIII dataset for 'Octopus' berrima, only Venus Bay and Tumby Bay showed signs of population expansion, whereas for 'Octopus' pallidus, population contraction was found only in Victoria. This discrepancy between methods was likely because Tajima's D analysis of a single mitochondrial loci is unable to distinguish between selective sweep or demographic changes whereas that of multiple unlinked loci in a nuclear DNA study is a true indicator of demographic changes (Winkelmann et al. 2013). While we obtained an effective population size of 511 [95% CI (325, 1158)] for the Victoria/Tasmania metapopulation of 'Octopus' pallidus, the small sample size of the other populations and metapopulations led to uncertainty in their effective population sizes. Nonetheless, our study highlights the importance of conservation genomics especially in commercially targeted species in order to prevent the loss of genetic diversity that can be caused by overfishing. With increasing interest in ocean-sourced proteins, it has become ever more essential in monitoring these populations. While fisheries are a useful source of specimen collection especially for understudied species and those found in remote locations, sampling bias may be an issue as seen in the lack of genomic information from intervening locations. Hence, this study was limited by the lack of access to intervening localities, which is important for the continuous monitoring of commercially harvested species.

Single-gene method in initial population analyses of understudied species
The COIII-based haplotype network highlighted the common ancestry that many populations of each species shared (i.e. sharing of a common haplotype between metapopulations Eely Point/Mount Dutton Bay and Port Lincoln/ Tumby Bay of 'Octopus' berrima, and among all populations of 'Octopus' pallidus). It also indicated that 'Octopus' berrima had a phylogeographic structure (with Venus Bay being the divergent population), which was supported by multiple approaches. However, 'Octopus' berrima also exhibited isolation by distance, which means that the highly divergent nature of Venus Bay may be biased by a lack of sampling from intervening populations. Nevertheless, such single-gene methods may still be useful in obtaining a first look at populations especially where funds and resources are lacking. For cephalopods that lack fundamental data such as species identity and ancestral gene flow, nuclear or mitochondrial genes combined with morphometrics may be useful in initial studies. However, ddRADseq may be more important for in-depth studies of gene flow and a continual monitoring of the impacts of commercial harvests, given the increased costs involved.

Genetic health of octopod populations and management implications
There are few SNP-based population genetic studies of cephalopods. Nevertheless, even with this limited comparison with other taxa, we found that based on nuclear SNPs, both O. berrima and O. pallidus showed similar heterozygosity and inbreeding levels to other octopods in the same region i.e. the southern blue-ringed octopus (Hapalochlaena maculosa) sampled from South Australia (Morse et al. 2017) (Table S10). This similarity is expected as these three species share similar life history traits and exhibit genetically structured populations.
Tumby Bay had never been fished before 2019 while Port Lincoln and Venus Bay have only been fished for 3-5 years, compared with Mount Dutton Bay and Eely Point that have been fished for over 30 years (Martino et al. 2021). Given that the Eely Point/Mount Dutton Bay metapopulation has the lowest autosomal heterozygosity level compared to other 'Octopus' berrima populations, this suggests particular attention be paid to this genetically unique locality in terms of fisheries management. It is also important to note that despite having slightly higher heterozygosity values, Venus Bay may also be vulnerable to loss of a stock given that it is geographically isolated from other localities or populations. The Victoria/Tasmania population of 'Octopus' pallidus also showed lower levels of heterozygosity compared to the South Australia (Port Lincoln) population. Since the former has been fished for at least two decades (Conron et al. 2020;Krueck et al. 2021) and previous genetic studies were based on microsatellites that prevents proper comparisons of how this might have changed over time (Higgins et al. 2013), it is possible that the Victoria/Tasmania population may have experienced a loss of genetic diversity over time but it is not clear whether the South Australia (Port Lincoln) population had lower or higher genetic diversity in the past. Nonetheless, it has been shown that overfishing causes a reduction in allelic diversity and heterozygosity (Hauser et al. 2002), and fishing pressure that concentrates on one population reduces overall genetic diversity (Cheng et al. 2021). The type of fishing gear employed in these fisheries also has significant impacts on their genetic health. Brooding females seek sheltered areas to lay their eggs such as under rocks and even inside human-made structures like glass jars. Fishers then make use of this behaviour to lay artificial octopod pots, resulting in primarily females at different brooding stages being harvested. As eggs are discarded in this process and thus lost from the gene pool, this in turn leads to reduced recruitment and hence genetic diversity. Currently in South Australia, there is a limit on the number of pots that fishers are allowed to set, and shelter pots must be used as opposed to trigger traps. Tasmanian and Victorian octopod fisheries are managed by pot limits and total allowable commercial catch (TACC) (Leporati et al. 2009;VFA 2022). This implementation of pot limits and TACC is in line with the recommended management strategy of holobenthic, semelparous octopods by Emery et al. (2016) and Martino et al. (2021). However, spatial closures (such as through MPAs) can result in a shift and therefore an increased fishing effort outside of the closures to compensate for lost opportunities (Emery et al. 2016). For genetically structured populations in particular, brooding females and juveniles from subpopulations may not be adequately protected, leading to localised depletion of individuals. Therefore on top of the current regulations, spatially rotating of fishing effort, to reduce the loss of genetically unique populations, could reduce impact of fishing on 'Octopus' berrima and 'Octopus' pallidus (Emery et al. 2016). Overall, our results emphasise the importance of continual assessment of their genetic health as only sampling at the advanced stages of exploitation makes it more difficult to detect significant changes in allelic diversity over time (Purcell et al. 1996).
The fact that these conclusions can be drawn using genome-wide analyses also reiterates the usefulness of SNPs in investigating and monitoring the levels and patterns of genetic differentiation among the populations of these two commercially harvested octopods. The use of pairwise F st as an a priori approach in identifying metapopulations coupled with a posteriori methods using PCA, sNMF and Structure analyses of SNPs not only provided consistent results of the genetic structures for both species, but were also advantageous in determining localities that are genetically homogenous. These analyses could in turn be useful for fisheries management whereby stock populations tend to be easier to manage at broader spatial scales (e.g. Eely Point/Mount Dutton Bay managed as one metapopulation). Knowing their genetic health and adaptive potential could then inform genomics-based management decisions to ensure the sustainability of these populations. This involves continual collaboration and development with both existing and new fisheries for not only obtaining genetic data but also to potentially expand the geographical spread of sampling sites that might reveal new genetic patterns. Furthermore, results from this study also provide an important baseline upon which future studies can be compared to investigate genetic health over time.

Comparative morphometrics shows trait differences in similar octopods
Many commercially harvested octopods are not identified at the species level, both in Australia and globally, which undermines their sustainable management (Sauer et al. 2019;Martino et al. 2021). Octopod species are difficult to identify by eye, as they change colors and shape, and species are often similar in appearance, particularly once killed (Amor et al. 2017). Here we demonstrate how multivariate morphometrics may be able to distinguish species of similar appearance. 'Octopus' pallidus appear to have shorter distances of beak to third arm, thicker widths of the stoutest arms, and larger widths of the largest suckers as compared to 'Octopus' berrima. Within males, 'Octopus' pallidus also have shorter hectocotylus and ligula than 'Octopus' berrima. Although the decision tree for male 'Octopus' pallidus could not reliably distinguish both species, it corroborated the fact that mature male 'Octopus' pallidus possess enlarged suckers on all arms, which are absent in 'Octopus' berrima. However, these may not be practical traits for fishers to distinguish each species. This method has also been used effectively in resolving the O. tetricus species complex, whereby male morphology can successfully differentiate O. tetricus and O. djinda as the latter has significantly higher sucker numbers (Amor et al. 2014;Amor and Hart 2021). The geographical populations of 'Octopus' pallidus were also able to be differentiated using the multivariate morphometric analysis as individuals in Victoria and Tasmania have deeper web sectors and thicker widths of the stoutest arms than those in South Australia. In contrast, the 'Octopus' berrima populations lack reliable morphological indicators to be distinctly differentiated. Whether this remains to be the case when compared to intraspecific populations further east of their distribution calls for further sampling in future studies.
Interestingly, female 'Octopus' pallidus were significantly heavier in Tasmania than Victoria, and individuals in Tasmania and Victoria were also significantly heavier than those in South Australia. The same patterns were found for male 'Octopus' pallidus when comparing Victoria and South Australia (Port Lincoln). This is possibly due to differences in sea temperatures whereby larger specimens tend to be found in colder high-latitude waters (i.e. Tasmania). Moreover it is known that, provided food is not a limiting factor, octopods in colder waters grow slower during the exponential phase and reach larger sizes at maturity than those in warmer waters (Forsythe and Hanlon 1988;Forsythe 2004). Future studies could verify and investigate if this wet weight trend is directly attributed to differences in sea temperatures.
In summary, this study provided the first comprehensive suite of genome-wide information on 'Octopus' berrima and 'Octopus' pallidus. We used a combination of morphometrics, mitochondrial gene COIII, and genomewide SNP discovery and genotyping via next-generation sequencing to provide information for important conservation and management decision-making. Our results provide evidence for genetic structuring in the populations of both species, which has implications for fisheries management. This is especially relevant to holobenthic octopods with highly structured populations, which are more vulnerable to localised overfishing or even commercial extinction. Although 'Octopus' berrima was morphologically similar within South Australia, they were genetically different between the geographical populations studied. In 'Octopus' pallidus, there was little genetic connectivity between South Australia and Victoria/Tasmania. We also showed that morphometrics can be a quick and affordable tool to gain useful insights into intra-and interspecific differences. We strongly recommend continual genetic monitoring of populations especially in commercially targeted species to evaluate diversity levels at new, existing, and potential fishing zones.