Systematics and Phylogenetic Relationships of New Zealand Benthic Octopuses (Cephalopoda: Octopodoidea)

The systematics of the New Zealand octopods have only been reviewed twice in the last 100 years. In these revisions many species have been provisionally classified in the genus Octopus. Recent genetic studies have synonymized some New Zealand species with octopuses from other regions. The present study investigates the systematics and phylogeny of octopuses from New Zealand using eighty eight specimens, three mitochondrial genes (16S rRNA, cytochrome c oxidase subunit I, and cytochrome c oxidase subunit III) and one nuclear gene (Rhodopsin). Forty-four new octopod DNA sequences (belonging to 13 species) were included, adding to the 83 existing sequences from GenBank. All sequences were used to generate phylogenetic trees based on Maximum Likelihood (ML) and Bayesian inference (BI), with a data set composed by 97 species, including octopod sister groups and Vampyroteuthis infernalis as an outgroup. Gene tree and species delimitation analyses revealed a distinct genetic difference between two sympatric Graneledone subspecies, which we propose as valid species. Muusoctopus tangaroa is a sister species of M. thielei from Kerguelen; while Enteroctopus zealandicus forms a clade with E. megalocyathus from South America and E. dofleini from the North Pacific. Similarly, Octopus campbelli, O. huttoni, and O. mernoo form a monophyletic group with Robsonella fontaniana from South America, Scaeurgus unicirrhus from the Atlantic and O. pallidus from Australia. Pinnoctopus cordiformis is close to Grimpella thaumastocheir and several species of Octopus sensu lato as in previous phylogenetic studies. This study suggests that octopuses from New Zealand have different phylogenetic and biogeographic origins and represent independent radiations into this region.


INTRODUCTION
The incirrates include benthic and pelagic octopods in two superfamilies: Argonautoidea Cantraine, 1841 (pelagic octopods), and Octopodoidea (d'Orbigny, 1840) (benthic and pelagic octopods). Benthic octopuses are a group of over 200 species inhabiting all oceans of the world, from tropical to polar regions, and from the intertidal to at least 3,000 m depth (Norman, 2000;Nesis, 2003;Hoving et al., 2014;Jereb et al., 2014). Despite this diversity, the current systematic relationships within the group are still poorly understood given their variable morphology and lack of suitable characters for morphological analysis (Strugnell et al., 2014). In recent years, several changes in octopod taxonomy have been proposed, including a new phylogenetic classification that positions incirrate octopuses in the superfamily Octopodoidea, which is composed of six families: Octopodidae, Megaleledonidae, Enteroctopodidae, Amphitretidae, Eledonidae, and Bathypolypodidae (Strugnell et al., 2014).
Considering that nearly 200 species are currently incorporated within this categorization, most of the octopod phylogenies published to date have included only few species (<30) (Carlini et al., 2001;Guzik et al., 2005;Strugnell et al., 2005;Strugnell et al., 2014), with just a few studies considering more than fifty species (see Lindgren et al., 2012;Ibáñez et al., 2014Ibáñez et al., , 2018. Incorporating much more species into octopod phylogenies seems problematic, as most species are recognized only from type material that has been fixed in formaldehyde and consequently lacks color and characters seen only in living specimens, and for which DNA sequences are not available. This has hindered the analysis of phylogenetic relationships as missing taxa can significantly influence tree topology (Graybeal, 1998;Poe and Swofford, 1999;Rosenberg and Kumar, 2003); therefore, including missing and poorly studied species in new DNA sequence analyses is important to provide a more complete and updated understanding of the phylogenetic relationships among benthic octopuses.
In the specific case of octopod fauna from New Zealand, this was initially reported by Dell (1952) in a monograph describing 14 species of benthic and pelagic octopuses. O'Shea (1999) revised New Zealand octopod fauna, placing 39 octopod species in six families and 14 genera, with two new genera and 16 new species. Many new species were assigned to an unplaced genus provisionally called 'Octopus' (Norman and Hochberg, 2005;Jereb et al., 2014); however, neither Dell nor O'Shea had the benefit of obtaining genetic information for complementing their morphologic approach. The most recent review of New Zealand biodiversity includes 41 octopod taxa (Spencer et al., 2009), although recent genetic studies have synonymized two of the New Zealand species (Octopus gibbsi O'Shea, 1999;Amor et al., 2014 andO. jollyorum Reid andWilson, 2015;Gleadall, 2016).
Biogeographically, New Zealand's marine fauna comprises both subtropical and tropical species (Shears et al., 2008). The relationships among benthic octopuses from New Zealand, Australia, and South America has been hypothesized from evidence based on morphology and distribution (O'Shea, 1999). Previous biogeographic studies based on ancestral distribution inferred from phylogenies did not include the New Zealand octopod fauna (i.e., Strugnell et al., 2008Strugnell et al., , 2011Ibáñez et al., 2016), suggesting that the inclusion of those species could dramatically change not only the phylogenetic hypothesis but also our knowledge of the biogeographic events that would explain the origin of the octopus in this region.
The aims of this study were to: (i) determine the evolutionary relationships within the New Zealand octopuses and their genetic relationships in the context of octopus phylogeny and (ii) clarify the current status of some species identities. Furthermore, establishing taxonomic clarity on a regional subset of species is prerequisite to larger scale revisions of the broader group. For this purpose, we constructed a molecular phylogeny of 88 species of benthic octopuses, in addition to nine outgroups to estimate their phylogenetic relationships.

Sampling
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. 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.

DNA Extraction, PCR Amplification, and Sequencing
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 R 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 2 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). Proteincoding 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
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 hillclimbing 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 burnin, 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 1 .
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 2 using both the COI and COIII

Phylogenetic Analysis
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.

Molecular Phylogeny
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 (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).

DISCUSSION
This study evidenced that benthic octopuses from New Zealand have different phylogenetic and biogeographic origins. Our review based on museum collections and phylogenetic analyses indicated that the New Zealand octopod fauna is composed of 16 species distributed in eight genera. This is a low diversity compared to tropical regions, but higher than cold and temperate ecosystems (Rosa et al., 2019

New Zealand Octopus Systematics
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., 2015Briceño et al., , 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  (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 subspecies 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.

Octopus Phylogeny
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., 2014Ibáñez et al., , 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.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

ETHICS STATEMENT
The animal study was reviewed and approved by the Universidad Andres Bello.

AUTHOR CONTRIBUTIONS
CI and MP-G conceived the idea, designed the study, analyzed the data, and led the writing of the manuscript. MF, SC, and PR collaborated in writing and provided editorial advice. All authors have read and commented on the manuscript.

FUNDING
This study was supported by the NIWA Internal Project CDTT 1607 "New Zealand octopus phylogeny" and the Victoria University funding. This work was partially funded by FONDECYT research grants 1181153, 11170617 and 11181320, awarded to CI, SC, and MP-G, respectively. The additional support from the INACH research grant RG 50-18 awarded to MP-G is also appreciated.