A total of 88 octopuses were obtained and examined from stored collections (Table 1 and Figure 1). Of these, 20 specimens were captured by bottom trawl during fisheries research voyages aboard the National Institute of Water and Atmospheric Research, Ltd. (NIWA) vessel R/V Tangaroa. A further set of 38 specimens were collected by New Zealand Ministry for Primary Industries scientific observers program from New Zealand fishing vessels. Finally, 30 octopuses were collected by NIWA staff during the annual Bluff oyster survey in South Island (Figure 1). Specimens were captured at depths ranging from 38 to 1208 m.

All specimens are deposited at NIWA Invertebrate Collection, Wellington, New Zealand (Table 1) and are available for examination. Additionally, we reviewed some type specimens from NIWA to confirm their identification and taxonomic status.

Pinnoctopus cordiformis (Quoy and Gaimard, 1832): Neotype NIWA 43044 (H-668), ML 120 mm, 41°09.14′ S, 173°15.07′ E, 23–24 m, 21/03/1997.

Octopus mernoo (O’Shea, 1999): Holotype NIWA 7555 (H-666), ML 33.5 mm, 43°51.19′ S, 178°58.81′ E, 480 m, 13/09/1989.

Graneledone taniwha (O’Shea, 1999): Holotype NIWA 662, ML 121 mm, 44°41.90′ S, 177°23.71′ W, 1135–1157 m, 17/10/1995.

Muusoctopus tangaroa (O’Shea, 1999): Holotype NIWA 7546 (H-660), ML 97 mm, 44°06.99′ S, 178°26.01′ E, 936–999 m, 11/10/1995.

Muusoctopus clyderoperi (O’Shea, 1999): Holotype NIWA 7556 (H-667), ML 80 mm, 39°58.55′ S, 178°14.80′ E, 900 m, -/04/1994.

Muusoctopus tegginmathae (O’Shea, 1999): Holotype NIWA 7545 (H-659), ML 64.5 mm, 39°57′ S, 178°16′ E, 1020–1250 m, 25/09/1995.

Most specimens were frozen at sea prior to being shipped to NIWA where they were defrosted. In the laboratory, mantle tissue samples were taken and stored in 99% ethanol until required for the molecular analysis. Foveaux Strait specimens were shipped in plastic bags on ice, inside an insulated box to NIWA, Wellington, and processed immediately. Small whole animals were preserved in 99% ethanol. Tissue subsamples were taken from larger animals before they were fixed in a buffered 5% formaldehyde solution, then transferred to 80% ethanol, for anatomical and morphological analyses.

Total DNA was extracted from 66 specimens of the total 88 examined using a high-salt extraction protocol (Aljanabi and Martinez, 1997), the phenol/chloroform method (Sambrook et al., 1989), or DNeasy^® purification kits (mouse tail protocol, Qiagen GmbH, Germany). Polymerase Chain Reaction (PCR) amplifications were carried out in 25 μL volumes with 5 units of Platinum^TM Taq DNA polymerase (Invitrogen) with 20 mM Tris HCl (pH 8.4), 50 mM KCl, 2.5 mM dNTPs, 3 mM MgCl₂ and 0.5 μM each of primers of the mtDNA genes Cytochrome Oxidase I (COI), Cytochrome Oxidase III (COIII) and 16S rRNA (Simon et al., 1991). Primers for COI were modified from Folmer et al. (1994) to match octopus DNA sequences in GenBank (Forward: TYTCAACAAATCATAAAGAYATTG G, Reverse: TATACTTCTGGRTGACCAAARAATCA). Primers for COIII were also modified from the literature (Forward: CAATGATGACGWGAYATTATTCG; Guzik et al., 2005 and Reverse: TCTACAAAATGTCAYTATCA; Simon et al., 1994). After an initial denaturation (2 min at 94°C), the reaction mixtures were subjected to 30–40 cycles of 94°C (30 s), [40–50°C (30 s) for COI; 45–65°C (30 s) for 16S; 40–45°C (30 s) for COIII], and 72°C (60 s) followed by a final extension at 72°C (10 min) using a thermal cycler. PCR products were purified using ExoSAP-IT and the DNA sequences were determined using a 3730 ABI Genetic Analyzer at Macrogen, Inc. (Seoul, South Korea). The resultant DNA sequences were aligned by Muscle using default parameters (Edgar, 2004) implemented in MEGA ver. X software (Kumar et al., 2018). Sequences generated in this study are available from GenBank (Table 1). Protein-coding sequences (COI and COIII) were translated to amino acids using the invertebrate mitochondrial genetic code to check for errors or gaps in MEGA.

Species delimitation was evaluated by using the Bayesian Poisson tree processes (bPTP) analyses (Zhang et al., 2013). Previously, Maximum Likelihood and Bayesian phylogenetic analyses were performed, including two preliminary steps on the aligned DNA sequences. First, Xia’s test for saturation of the phylogenetic signal of each gene was performed using Dambe ver. 6.0 (Xia, 2017). Second, the best substitution model for each gene was estimated with jModelTest (Posada, 2008) using the Bayesian Information Criterion (BIC).

The phylogenetic relationships of the benthic octopuses were examined using a Maximum Likelihood (ML) reconstruction via the IQ-TREE online server (Trifinopoulos et al., 2016) with hill-climbing NNI tree search strategy (Nguyen et al., 2015). The ModelFinder option (Kalyaanamoorthy et al., 2017) was used under a partition scheme including codon position for coding genes (COI, and COIII). Statistical support was estimated using 5,000 ultrafast bootstrap replicates (Minh et al., 2013). The trees were rooted using the cirrates Opisthoteuthis mero O’Shea (1999) and O. chathamensis O’Shea (1999) as outgroups, as Cirrata is well-established as the sister group of Incirrata (Voight, 1997; Young et al., 1998; Lindgren et al., 2012).

Phylogenetic reconstruction was inferred from a partitioned matrix (16S, COI, COIII) with a different substitution model for each gene. Bayesian analyses were conducted using MrBayes ver. 3.2 (Ronquist et al., 2012) with four chains, each with 10 million generations, sampled every 1,000 generations. Bayesian analyses were performed several times to compare the likelihood values of each run using Tracer ver. 1.5 (Rambaut and Drummond, 2009). The first 1,000 trees of each run were discarded as burn-in, and a consensus of the remaining trees was calculated. FigTree ver. 1.4 was used to edit the trees (Rambaut, 2009). In both phylogenetic analyses (ML and BI), we used specimens for which all mitochondrial genes were available (44 specimens, Table 1). The consensus tree was finally used as input for the species delimitation analysis with the Bayesian Poisson Tree Process method (bPTP; Zhang et al., 2013) as implemented in the web server¹.

Additional species boundaries were delimited using the Automatic Barcode Gap Discovery method (ABGD; Puillandre et al., 2012). The ABGD method recursively searches for major changes in the slope of ranked pairwise genetic distances between groups of individual sequences. Through this, ABGD proposes a distance superior to maximal intraspecific sequence divergences, as determined using a coalescent model. These distances potentially correspond to the frontiers between intra and interspecific distances, or the so-called barcode gap. ABGD analyses were performed online² using both the COI and COIII data sets independently, excluding outgroups. The p-distance with a minimum gap width of 1.5 were selected. The remaining parameters were set as default (Pmin = 0.001, Pmax = 0.100, Steps = 10, Number of bins = 20).

In a second phylogenetic analysis, mitochondrial DNA sequences obtained during the present study were combined with those of other species available at GenBank (16, COI, COIII, and Rhodopsin) to explore the phylogenetic position of the New Zealand taxa. These included Rhodopsin (RHO) sequences in a matrix with 97 species (Table 2), including 88 species of octopuses from the superfamily Octopodoidea, in addition to outgroups from two pelagic octopuses of the superfamily Argonautoidea (Argonauta argo and Argonauta nodosus), six cirrates (Opisthoteuthis mero, O. chathamensis, O. depressa, O. massyae, Cirroctopus glacialis, and Stauroteuthis gilchristi), and the vampire squid Vampyroteuthis infernalis. These analyses were performed in the same way as the previous analysis in IQ-TREE and MrBayes.

Species delimitation using bPTP agreed in defining the clusters of our dataset, with both analyses identifying 12 entities of New Zealand octopuses with posterior probabilities of conspecificity ranging from 0.90 to 1.0 (Figure 2). In this phylogeny, we found four clades within the Incirrata: Clade 1, compound by A. nodosus (Lightfoot, 1786); and Clade 2 compound by four species of deep-sea dwelling Graneledone and Thaumeledone (family Megaleledonidae), all of them recognized by a single sucker row. Within this clade, the reciprocal monophyly between G. taniwha taniwha O’Shea, 1999 and G. taniwha kubodera (O’Shea, 1999), in addition to the species delimitation analysis (bPTP and ABGD), suggest they would correspond to separate species (Figure 2). In Clade 3, species delimitation analyses evidenced two species, Enteroctopus zealandicus (Benham, 1944) and Muusoctopus tangaroa (O’Shea, 1999); whereas Clade 4 was compound by species of the family Octopodidae: P. cordiformis (Quoy and Gaimard, 1832), Octopus mernoo (O’Shea, 1999), Octopus huttoni (Benham, 1943) and Octopus campbelli (Smith, 1902) (Figure 2). Species delimitation analyses with ABGD for both coding genes (COI and COIII) identified 12 genetic groups of 11 morphological species based on p-distances on the gene tree (Figure 2).

The consensus tree from the Bayesian analysis of the combined sequences (GenBank and new sequences) had high posterior probability values (>0.95) for most nodes, and a similar topology compared to the ML tree (Figure 2). The ML tree from IQ-TREE had a better topology, resolving the polytomies present in the Bayesian tree with high bootstrap support (>70, Figure 3, Mendeley Dataset: DOI: 10.17632/5vkm46hm49.2). For this reason, we present the ML tree with posterior probability values from the Bayesian analysis. Within this tree, the genus Octopus is polyphyletic, probably related to the fact that many Octopus species are poorly described (c.f. Norman and Hochberg, 2005). New Zealand benthic octopuses are placed in clades that correspond to three distinct families. Clade 1, Enteroctopodidae, is compound by the species Enteroctopus dofleini (Wülker, 1910) from Alaska, Enteroctopus megalocyathus from Chile and E. zealandicus from New Zealand (Figure 3). In the same clade, Muusoctopus thielei (Robson, 1932) from Kerguelen is closely related to M. oregonensis (Voss and Pearcy, 1990) from the North Pacific and M. tangaroa from New Zealand. Within clade 2, Megaleledonidae, Thaumeledone zeiss O’Shea, 1999 was included in a clade comprising T. gunteri (Robson, 1930), T. rotunda (Hoyle, 1885) and T. peninsulae (Allcock, Collins, Piatkowski and Vecchione, 2004) from Antarctica (Figure 3). In the same clade, Graneledone taniwha taniwha and G. taniwha kubodera are sister taxa, as G. challengeri (Berry, 1916) and G. antarctica (Voss, 1976) (Figure 3). Clade 3 comprised the family Octopodidae (Figure 3), where P. cordiformis was a sister species of Grimpella thaumastocheir (Robson, 1928), both species nearly related to the Octopus s.l. clade composed by Indo-Pacific species. The phylogenetic position of P. cordiformis and G. thaumastocheir was similar to previous phylogenetic analysis, with the absence of ink sac in G. thaumastocheir being an apparent adaptation to deep sea (Guzik et al., 2005; Strugnell et al., 2014). Octopus campbelli, O. huttoni and O. mernoo formed a monophyletic group with Robsonella fontaniana (d’Orbigny, 1834) from Chile, Scaeurgus unicirrhus (Delle Chiaje in d’Orbigny, 1841) from the Atlantic, and O. pallidus (Hoyle, 1885), from Australia. Clade 4 was composed of O. fitchi (Berry, 1953), Paroctopus digueti (Perrier and Rochebrune, 1894) from Northeastern Pacific and O. tehuelchus (d’Orbigny, 1834), from the Southwest Atlantic (Figure 3). Clade 5 included Ameloctopus litoralis (Norman, 1992) from Australia, Cistopus (Gray, 1849) and Octopus s. l. from the West Pacific. Clade 6 included Abdopus aculeatus (d’Orbigny, 1834), O. cyanea (Gray, 1849), and O. laqueus (Kaneko and Kubodera, 2005) from the West Pacific (Figure 3). Clade 7 contained the genus Amphioctopus (Fischer, 1882) and Hapalochlaena (Robson, 1929), from the West Pacific (Figure 3). Finally, clade 8 was composed by octopodid species including O. oliveri (Berry, 1914) with members of the Octopus sensu stricto clade including O. tetricus (Gould, 1852), from Australia/New Zealand and others from America (Figure 3).

The most complete information to date on octopod fauna from New Zealand waters correspond to O’Shea’s (1999) monograph; however, considering the limitation imposed by the lack of any genetic analysis, identifying octopuses solely using O’Shea’s descriptions can be challenging. Norman et al. (2014) were critical on this revision and questioned the validity of some species (e.g., Thaumeledone, Pinnoctopus). In fact, Norman and Hochberg (2005) placed several species from New Zealand in different genera than those proposed by O’Shea (1999) (e.g., P. cordiformis and P. kermadecensis). O’Shea (1999) resurrected P. cordiformis designating a neotype in the absence of type material, but other authors argued that morphologically, this species is a senior synonym of Macroctopus maorum (Hutton, 1880) (Robson, 1929; Norman and Hochberg, 2005). Our phylogenetic results added to the review of the neotype (NIWA 43044, H-668) agreed with O’Shea (1999) and suggest placing the currently recognized species Octopus cordiformis in the genus Pinnoctopus. Similarly, previous studies also placed P. kermadecensis within the genus Callistoctopus (Norman and Hochberg, 2005; Reid and Wilson, 2015); however, molecular information would first be required in order to confirm the valid status of the species. In this context, O’Shea (1999) suggested both Macroctopus and Callistoctopus are junior synonyms of Pinnoctopus, and Voss (1981) suggested that no valid basis exist for recognizing any distinction between Callistoctopus and Octopus. Therefore, and based on this information, our study suggests that P. cordiformis and P. kermadecensis would remain as the correct names. Indeed, the name P. cordiformis has been consistently used in several studies recently carried out in New Zealand and Australia (see Carrasco, 2014; Orbach and Kirchner, 2014; Briceño et al., 2015, 2016).

In the present study, most New Zealand octopodids (P. cordiformis, O. campbelli, O. huttoni and O. mernoo) were found within Clade 3 (see Figure 3), sharing a common ancestor with other Pacific species from different genera (Grimpella, Octopus, Robsonella, and Scaeurgus). In fact, Robson (1929) observed morphological similarities between Robsonella and Scaeurgus in the terminal organ’s shape and the presence of a ligula with robust cheeks. The specimens examined from our Clade 3 (New Zealand and Chile) shared a similar hectocotylus shape, suggesting this clade require a new classification. Our genetic analyses suggest that some shallow-water, small-bodied octopuses from New Zealand (Octopus campbelli, O. huttoni, and O. mernoo) currently placed in Octopus sensu lato would belong to the Robsonella clade. In fact, Adam (1938) recognized only two species (R. fontaniana and R. campbelli), while Pickford (1955) recognized three species (R. fontaniana, R. campbelli, and R. huttoni). Clearly all these species share a recent common ancestor and have similar morphologies. Based on this information, we suggest identifying the New Zealand species O. campbelli, O. huttoni and O. mernoo as species within the genus Robsonella, as they shared several morphological features (radula, skin, hectocotylus) following the diagnosis of the genus (Ibáñez et al., 2008). Other authors have also used this identification (Sweeney and Roper, 1998; Sweeney, 2017).

Our DNA sequence data suggests that there was a close relationship (99% similarity in 16S and COIII, and 100% similarity in COI) between E. zealandicus (yellow octopus) from New Zealand and E. megalocyathus (red octopus) from Chile, relationship that deserve an improved revision as Hudelot (2000) also suggested that E. magnificus (Villanueva et al., 1992) and E. zealandicus may be conspecific. Both species (E. zealandicus and E. megalocyathus) are similar but have slight differences in morphometric, meristic measurements and coloration, which require further investigation (M.C. Pardo-Gandarillas et al. unpublished data). Norman and Hochberg (2005) and Jereb et al. (2014) classified Octopus oliveri in Octopus sensu lato provisionally, but our phylogenetic analysis suggests that this species is a close relative to the Octopus sensu stricto group to retain it in this genus. Another species not included in our genetic study, a recently described species from the Kermadec Islands (Octopus jollyorum; Reid and Wilson, 2015) has been suggested as a junior synonym of O. sinensis (d’Orbigny, 1841) (Amor et al., 2017), although Gleadall (2016) identified them as different species with clear morphological differences.

The finding of reciprocal monophyly and species delimitation analyses suggested that the subspecies G. taniwha taniwha and G. taniwha kubodera are different species. Both sub-species have similar morphological features, including suckers, gills, and wart counts (Ibáñez et al., 2012). O’Shea (1999) suggested that future descriptions of Graneledone species should provide more morphological detail (cartilaginous cluster distribution and composition, and arm sucker counts), and proposed that the only way to differentiate these two species is by the number of cartilaginous processes per cluster (i.e., G. taniwha taniwha 1–37 and G. taniwha kubodera 4–13). Similarly, specimens of each taxon were also examined here, suggesting they had different wart head counts (i.e., G. taniwha taniwha 10–14, G. taniwha kubodera 5–8) and confirming the distinction observed within the molecular phylogeny.

The genus Benthoctopus is not sustainable, as Gleadall et al. (2010) pointed out that Benthoctopus is a junior synonym of Bathypolypus and identified species in genus Benthoctopus as species of Muusoctopus. Therefore, there is no reason in retaining both Benthoctopus and Muusoctopus (e.g., Norman et al., 2014). In this context, Ibáñez et al. (2016) also suggested (based on morphology and genetics) that all New Zealand species of Benthoctopus should be included within the genus Muusoctopus. Future research should therefore target to compare morphometrics and genetic data to specifically determine the number of Muusoctopus species present in New Zealand waters.

The phylogeny reported here contains several groups, represented by the families Bathypolypodidae, Eledonidae, Megaleledonidae, Enteroctopodidae, and Octopodidae (Figure 3). Within Octopodidae, we recognized two monophyletic groups, one composed by Octopus sensu stricto (Clade 8) and the other by species of the genera Amphioctopus and Hapalochlaena (Clade 7). Two other groups were paraphyletic, being composed by species of Octopus sensu lato (Clades 3, 4, and 5) in addition to Abdopus, Callistoctopus, Cistopus, Pinnoctopus, Robsonella, and Scaeurgus. Our finding of Octopus as a polyphyletic group is consistent with previous studies (e.g., Strugnell et al., 2005; Lindgren et al., 2012; Ibáñez et al., 2014, 2018), suggesting that the genus Octopus need an urgent revision.

Recent studies carried out by Strugnell et al. (2014) evaluated the phylogenetic relationships of 23 octopods using four mitochondrial and three nuclear genes, and evidenced a similar topology compared to that obtained with our fewer genes (three mtDNA regions) and higher number of species (88 spp. in five families; see Figure 3). In this context, increasing the coverage of species and including additional characters is a well-recognized approach to improve phylogenetic analyses (Graybeal, 1998; Poe and Swofford, 1999; Rosenberg and Kumar, 2003), suggesting that our phylogeny is a solid estimation of the evolutionary relationships. Since our phylogenetic reconstruction was based only on three mitochondrial genes and one nuclear gene, it is plausible to expect that the inclusion of more markers would improve our understanding of the evolutionary relationships of the New Zealand octopuses. However, the polyphyletic nature of the genus Octopus is clearly an artifact of poorly described species being placed into the genus because of uncertainty about their true taxonomic positions (sensu Norman and Hochberg, 2005). As suggested by previous authors (Gleadall, 2004; Kaneko et al., 2011), Octopus systematics still requires an extensive revision in order to solve some of the difficulties in finding informative morphological characters. Based on this information, we suggest that several species included in Clade 1 (representing the family Enteroctopodidae) may not belong to the genus Octopus (e.g., O. californicus, O. conispadiceus, O. hongkongensis).

The presence in New Zealand of 16 species of benthic octopus from different genera and environments suggests a history of several radiations from tropical and cold-water ancestors. Most octopuses included in the present study inhabit the Indo-Pacific, a region that based on the high diversity of benthic octopuses is recognized as the potential origin of the family Octopodidae, and from where many species radiated worldwide (Rosa et al., 2019). Close biogeographic relationships of benthic octopuses from New Zealand, Australia and the Southern Ocean have been recently revealed (Rosa et al., 2019), confirming the taxonomic relationships proposed by O’Shea (1999). The close phylogenetic relationships of the New Zealand O. campbelli, O. huttoni and O. mernoo with O. pallidus from Australia and R. fontaniana from South America is probably related to dispersal events after the circumpolar current was established during the Cenozoic (Strugnell et al., 2008). The same pattern would be expected for the Octopus s.s. clade with species from America, Australia and New Zealand, and the Mediterranean, which suggests a classic Tethyan distribution. For the deep-sea species of the genera Graneledone and Thaumeledone, an Antarctic origin is probable based on the findings by Strugnell et al. (2008). In the case of Muusoctopus, dispersal events from the North Pacific to the Southern Ocean and Atlantic (Gleadall, 2013) and from the Atlantic to the Southern Ocean (Ibáñez et al., 2016) have been previously proposed. The molecular phylogenetic approach presented here has added important information to the current systematics of the New Zealand octopod fauna; nonetheless, further studies are still required considering larger sampling sizes and a mixture of both mitochondrial and nuclear molecular markers to properly clarify their biogeographic origin and diversification.