Morphological diversity and a ribosomal phylogeny of Rhabdopleura (Hemichordata: Graptolithina) from the Western Pacific (Singapore and New Zealand), with implications for a re-evaluation of rhabdopleurid species diversity

The recent discovery of Rhabdopleura in Singapore and the chance collecting of fresh material from northern New Zealand (Three Kings Shelf) provided an opportunity to sequence the specimens with an aim to determine their species identity. Phylogeny reconstructions of two new Rhabdopleura taxa based on ribosomal and mitochondrial genes suggest a different identity from known samples, including putative Rhabdopleura annulata, first described from the Three Kings Shelf but sequenced from the Great Barrier Reef. Pairwise distances between rhabdopleurids for the 16S rRNA locus were several magnitudes larger than that of 18S rRNA, and might potentially be a suitable barcoding gene once sufficient samples of conspecifics are collected to determine the barcoding gaps. Type material of R. annulata was re-examined, as well as Rhabdopleura material from eight other New Zealand localities from north of subtropical Raoul Island (Kermadec Ridge) at ~29° S to the subantarctic Campbell Plateau at ~49° S. Six morphological characters, four of them new, were applied to all samples. The findings from morphology suggest (1) the holotype and cotype of R. annulata might not be conspecific; (2) there are 3–4 variants (species?) on the Three Kings Shelf; and (3) there are 2–3 additional variants (species?) elsewhere in the New Zealand region.


Introduction
The pterobranch hemichordate genus Rhabdopleura Allman, 1869, an extant graptolite taxon (Mitchell et al. 2013), is widely distributed in space and time. It is found at the present day in all major oceans, ranging across Arctic, tropical and sub-Antarctic latitudes and from the intertidal to c. 900 m (Tassia et al. 2016;Beli et al. 2018). Notwithstanding this distribution, only five named species are presently accepted out of the eight species described (Swalla and van der Land 2022). Fossil species nominally range from the Middle Cambrian to the Eocene (Chapman et al. 1995;Sennikov 2016) (but Mierzejewski (1986), Urbanek and Dilly (2000) and Mierzejewski and Kulick (2001) have interpreted the earliest true representatives of Rhabdopleura as Middle Jurassic in age). Sennikov (2016) recognizes eight fossil species. The Communicated by K. Kocot paucity of fossil species is unsurprising, given that Rhabdopleura colonies not only lack a hard skeleton, but, unlike planktonic graptolites, are found on benthic hard substrata and typically not amenable to fossilization in fine sediments. On the other hand, the low number of described living species invites explanation, especially since their geographic and environmental distribution is substantial.
The low number of described species may be attributed to the lack of study of the range of morphological variation in nominal species and limited molecular sequencing to determine genetic diversity. Ironically, even though there are so few living species, some recent workers have noted the difficulty of determining the identity of newly discovered material in remote locations based on morphological criteria. Such was the case when Dilly and Ryland (1985) encountered an intertidal Rhabdopleura species in Fiji, commenting: "The claims for any of them to be separate species are somewhat tenuous, and depend, as always, on the estimation of what are differentiating characters." The relative (or apparent) paucity of such characters led to early synonymy. For example, the type species, Rhabdopleura normani Allman, 1869, apparently widely distributed along western European coasts (Burdon-Jones 1954), was interpreted by Schepotieff (1909) to include four other eastern Atlantic species, namely Rhabdopleura mirabilis Sars, 1872, Rhabdopleura compacta Hincks, 1880, Rhabdopleura grimaldii Jullien, 1890 and Rhabdopleura manubialis Jullien and Calvet, 1903. In the same work, Schepotieff (1909) described a new species, Rhabdopleura striata, from Sri Lanka. Stebbing (1970a) subsequently clarified the separate status of Rhabdopleura compacta, removing it from synonymy with R. normani. The latter name has since been accorded to specimens from Bermuda (Dilly 1985), southern Argentina (López Gappa 1987) and, with less certainty, to the material from Fiji (Dilly and Ryland 1985). Authors agree that R. mirabilis is reliably synonymized but the status of R. grimaldii and R. manubialis, both from the Azores, is "uncertain/nomen dubium" (Swalla and van der Land 2022), underscoring the need for re-examination.
A similar example is that of Rhabdopleura annulata Norman, 1921. It was formally described from northern New Zealand, but was at the same time attributed by its author to unidentified Siboga Expedition material from Indonesia described by Harmer (1905). It was subsequently reported from southern Australia (Johnston 1937) (but Shepherd (1997) re-interpreted it as R. normani), Darwin, Northern Territory (Gordon 2009) and the Great Barrier Reef, Queensland (Ramírez-Guerrero et al. 2020). Only one other species has been described since Norman (1921), i.e. Rhabdopleura recondita Beli, Cameron & Piraino in Beli et al. 2018, a Mediterranean species hitherto overlooked because of its highly cryptic habit of growing in dead skeletal parts of several bryozoan species. Relatively few specimens of Rhabdopleura have been sequenced, including some that have not been taxonomically determined to species. The results suggest the possibility of cryptic species Ramírez-Guerrero et al. 2020).
Recently, in the course of targeted collecting of bryozoans by five of the present authors for a global molecular phylogeny (Orr et al. 2022), the genus Rhabdopleura was newly discovered in Singapore . The sole specimen was not able to be determined to species on the basis of morphology and appeared to be unrelated to Rhabdopleura previously collected from Sulawesi and West Papua by the Siboga Expedition. The same targeted bryozoan collecting also yielded a sample of Rhabdopleura from the Three Kings Shelf, northern New Zealand, geographically close to the type locality of R. annulata. We sequenced both the Singaporean and Three Kings Shelf samples to determine their phylogenetic relationships.
At the same time, owing to the evident challenge of identifying Rhabdopleura on the basis of morphological criteria, the opportunity was taken to examine a collection of Rhabdopleura from several New Zealand localities, ranging from the Kermadec Ridge at~28°S to the subantarctic Campbell Plateau at~49°S. These were not amenable to molecular sequencing in the same manner as the Singapore and Three Kings Shelf samples, having either been first preserved in seawater-formalin, followed by isopropanol, or discovered on dried unpreserved dead molluscan shell (principally Pecten and glycymeridid bivalves). A goal of the comparative morphological study was to ascertain not only if the sequenced sample from the Three Kings Shelf matched R. annulata, the type material of which was newly re-examined, but also if there exist overlooked or undervalued morphological criteria that could be useful in highlighting within-species morphological variation or discriminating among putative other species.
The purpose of this paper, therefore, is as follows: (1) to reconstruct a phylogeny including the two new Rhabdopleura samples using genome skimming to recover mitochondrial and ribosomal genes, and (2) to discuss the most useful morphological criteria applicable to species discrimination in Rhabdopleura.

Morphological study
The provenance of all specimens examined for this study is given in Table 1. Figure 1 plots localities for all the New Zealand samples, most of which are registered in the NIWA Invertebrate Collection of the National Institute of Water and Atmospheric Research, Wellington. Most specimens were examined by scanning electron microscopy (SEM). Selected type and other material were examined uncoated by SEM at the Natural History Museum, London (NHMUK), using a JEOL IT500 SEM, supplemented by light microscopy of slide material. Voucher specimens at the University of Oslo (NHMO) were imaged uncoated using a Hitachi TM4040PLUS Tabletop SEM. That prepared for SEM was air-dried from ethanol unless the specimen was already dry.
Materials examined at NIWA were imaged using a Hitachi TM3000 Tabletop SEM and coated using gold-palladium in a sputter coater for 30 s to avoid charging.
Historically, erect and creeping tubes of Rhabdopleura have mostly been imaged and measured using transmitted light microscopy. We have found the clarity and magnification of SEM images to be highly optimal for providing metric data, measured using Fiji (ImageJ) software (Schindelin et al. 2012). Terminology for descriptions and measurements (Fig. 2) is based on that provided by Maletz et al. (2014). A variety of metric measures were applied to all specimens examined by us to determine which morphological characters might be more informative. We settled on four discriminating characters that provide unambiguous and easy-to-measure standard metrics using SEM (Table 2). They are as follows: 1) External diameter of erect tubes between fusellar collars (Fig. 2c). Since tube diameter may vary from its base to its aperture, and sometimes just proximal to the aperture, we determined that measurements are best made in the distal half proximal to the aperture (where possiblesometimes a majority of tubes in a colony may be distally broken depending on method of collection, by hand if by diving or by dredge in deeper water, the latter being more damaging to the fragile colonies). 2) Number of fusellar collars in 500 μm of erect tubes. As with character 1, we determined that these should be counted in the distal half of tubes proximal to the aperture where there seems to be more consistency in fusellus height (which can be shorter and compressed in the proximal quarter of erect tubes in some samples). This is a new character. Its value is in smoothing across the variability in fusellus height over a standard distance. 3) Fusellus height (Fig. 2c). This is the distance between one fusellus collar and the next. As with characters 1 and 2, measurements are best made in the distal half of tubes below the aperture. In some colonies, individual fusellus height may be consistent, but mostly it can vary, so it is useful to make two measurements for each fusellus, one on both sides of the face that is presented in the SEM image. 4) Point-to-point angles of divergence in zigzag sequences of fusellar sutures on creeping tubes (Fig. 2e). Inasmuch as some of the zigzag sutures may be curved, while others are straight, straight lines were drawn between points of intersection, regardless of curvature, to obtain measurable angles that can be compared across samples. (Curvature of zigzag suture components means that initial angles are small at points of intersection and then widen. Some measurements of initial angles were made but eventually deemed to be less informative, for comparative purposes, than point-to-point angles.) This is a new character.
Two other characters, both descriptive as used here, are included in Table 2. These are as follows:  Table 2 Measurements of key characters of selected specimens examined in this study. All specimens are from New Zealand except for that from Singapore. The first two data columns pertain to sequenced specimens. Data columns 3 and 4 pertain to type specimens of Rhabdopleura annulata Norman, 1921 in the Natural History Museum, London. The remaining columns are arranged in latitudinal order from north of subtropical Raoul Island, Kermadec Ridge, to the subantarctic Campbell Plateau. For all columns, data pertain to one selected erect tube (or two short tubes if their distal ends broken, as in the two sequenced specimens) and a representative creeping tube with a clean unbroken sequence of zigzag fusellar sutures between two erect tubes. For each relevant metric (characters 1, 3, 4) n refers to the number of measurements Character  10-12 (11) 11 19-22 (21) 20, 21 9-13 (11)13 10-20 (16) 17 12-14 (14) 14) 3. Distance between fusellar collars of erect tube (distal part): range and mean (μm) *Sequenced ⁋Voucher material; only a few short broken erect tubes present, hence distal parts of tubes not available for measurements; data for character 2 extrapolated for numbers of fusellar-collars in 500 μm **Newly re-examined for measurements by SEM †To determine a modal value for a sometimes limited dataset (data columns 1-4), actual degrees were scaled up or down to the nearest 5°interval 5) The width of the fusellar-suture zigzag sequence along uninterrupted straight sections of the creeping tube ( Fig.  3). There are three main patterns. The zigzag sequence can cross the entire width of the creeping tube (e.g. Fig.  3a), occupy a narrow zone along the midline of the creeping tube (e.g. Fig. 3g-j), or be variably intermediate between these two extremes. There can be some variation within a single colony (e.g. Fig. 3g, h). This is a new character. 6) Curvature of the zigzag components of the fusellar suture of the creeping tubes. They can be all or mostly curved (e.g. Fig. 3b), mostly straight (e.g. Fig. 3j), or variable within a zigzag sequence. This is a new character.

DNA extraction, sequencing and genome skimming
Extraction of DNA and sequencing generally followed Orr et al. (2022). Briefly, samples were crushed with a pestle in a mixture of both the lysis buffer and proteinase K from DNeasy Blood and Tissue Kit (Qiagen). Subsequently, genomic DNA was extracted from digested samples following the manufacturer's protocol using the aforementioned kit.
Sequencing was performed at the Norwegian Sequencing Centre (Oslo, Norway) on an Illumina HiSeq 4000 lane (150×150 bp). Adapter removal and quality trimming on raw reads were performed using cutadapt v1.4.2 (Martin 2011). Assembly of trimmed reads was conducted using SPAdes v3.13.0 (Bankevich et al. 2012) under default settings. Identification of 18S rRNA genes was conducted with BLASTn by searching assembled sequences against Rhabdopleura sequences from GenBank (Table S1; see Ramírez-Guerrero et al. 2020). A single best hit was extracted for each sample sorted by highest bitscore, followed by evalue and then percentage identity. Overlapping regions with that of the subject were then extracted from the contig assembled to obtain the region of interest for downstream phylogenetic reconstruction. Mitochondrial genomes were assembled for each sample using NOVOPlasty v4.3.1 (k-mer=33) (Dierckxsens et al. 2017), with the complete mitochondrial genome of congeneric R. compacta (FN908482; Perseke et al. 2011) as reference. Annotation of the mitochondrial genome was conducted first using MITOS2 (Donath et al. 2019), followed by NCBI's ORFfinder (Wheeler et al. 2007) under translation table 24 (Rhabdopleuridae mitochondrial code) to refine the annotation for protein-coding genes (PCGs). Finally, we manually corrected PCG annotations with amino acid alignments of other hemichordate mitochondrial genomes.

Phylogenetic reconstructions and pairwise distance estimations
Two separate phylogenetic reconstructions were executed. The first was a mitochondrial genome phylogeny for hemichordates based on 13 protein-coding genes (PCGs), building upon the phylogeny in Li et al. (2019). The second was a rRNA phylogeny of Rhabdopleura based on Ramírez-Guerrero et al. (2020) (Table S2).
For the hemichordate mitochondrial genome phylogeny, PCG sequences were downloaded from both GenBank and Li et al. (2019). Sequences were translated using ExPASy (Gasteiger et al. 2003) and then aligned using MAFFT v7.427 under L-INS-i setting (Katoh and Standley 2013). Amino acid alignments were then back-translated using pal2nal v14 (Suyama et al. 2006) and trimmed using trimAl v1.4 with the heuristic algorithm activated (Capella-Gutiérrez et al. 2009). Evolutionary models for each gene were determined using ModelTest-NG (Darriba et al. 2020) under default settings. To ensure that the mitochondrial genome was properly assembled, a maximum likelihood phylogeny with models specified according to ModelTest-NG (Darriba et al. 2020) was reconstructed by RAxML-NG v0.8.1 (Kozlov et al. 2019) for each mitochondrial PCG with 10 random and 10 parsimony starting trees and 100 bootstrap pseudoreplicates. We expected that both samples sequenced here would be recovered as a monophyletic group with the other congeners.
After inspecting all gene trees and ensuring the mitochondrial genomes were accurately recovered, aligned DNA sequences from 13 PCGs from all taxa were concatenated into a single matrix. A maximum likelihood phylogeny for hemichordates was reconstructed using RAxML-NG v0.8.1 (Kozlov et al. 2019), with the matrix partitioned by gene and models specified as determined above. A total of 50 random and 50 parsimony trees were specified as starting trees, with 1000 bootstrap pseudoreplicates to estimate node support. Furthermore, Bayesian inference (BI) was performed using MrBayes v3.2.6 (Ronquist et al. 2012) on the concatenated data matrix partitioned by gene, and each partition model was estimated using ModelTest-NG (-T mrbayes) (Darriba et al. 2020). Two distinct runs over four Markov chains Monte Carlo, each with 12 million generation and one tree logged every 100 generations, were performed. Convergence was determined using Tracer v1.7 (Rambaut et al. 2018), which showed that the first 20,001 trees could be discarded as burnin (average standard deviation of split frequencies = 0.000892).
For the rRNA phylogeny of Rhabdopleura, 16S and 18S rRNA genes were concatenated for analysis. Sequences were first aligned using MAFFT v7.427 under L-INS-i strategy (Katoh and Standley 2013). Subsequently, trimming and DNA evolution model selection, as well as maximum  Finally, uncorrected pairwise distances between Rhabdopleura sequences were separately computed for both 16S and 18S ribosomal genes from the trimmed alignment using MEGA X (Kumar et al. 2018) (Tables S1, S2).

Mitochondrial genomes and gene order
A circularized mitochondrial genome was recovered for both individuals by NOVOPlasty v4.3.1 (Dierckxsens et al. 2017). The genome length of the New Zealand sample was 16,307 bp, and that of the Singapore sample was 16,500 bp. We recovered all 13 PCGs and both mitochondrial rRNA genes. Curiously, for both samples, tRNA-Arg remained unannotated, which was the only missing tRNA flagged by MITOS2 (Donath et al. 2019). In addition, the position for cytochrome c oxidase subunit II (COII) for both samples lay between ATPase 6 (ATP6) and NADH dehydrogenase subunit 6 (NAD6), whereas that in R. compacta was between cytochrome c oxidase subunit III and ATPase 8 (ATP8) (

Phylogeny reconstructions and pairwise distances
Rhabdopleura was recovered as a monophyletic clade in all reconstructions with maximum node support, and both samples in this study were nested within Rhabdopleura (Fig. 4). In both the concatenated rRNA and mitogenome trees, samples from Singapore and New Zealand were recovered as sister taxa. In the rRNA tree, weak support for this relationship was observed (bootstrap = 37, unresolved in Bayesian inference) (Fig. 4a). The phylogeny based on the mitochondrial PCGs performed better, with strong support for the close relationship between the two taxa (bootstrap = 73, posterior probability = 0.98), which were together sister to the remaining Rhabdopleura spp. (R. annulata and R. compacta) with maximal support (Fig. 4b).

Morphological study
Descriptions are given below for the voucher material of the sequenced species from Singapore (NHMO H2136) and the Three Kings Shelf (NHMO H2137) at the University of Oslo. Partial descriptions only are given of re-examined type material of R. annulata, based on some new measurements of selected characters, in order to supplement the original description (Norman 1921). Following that, an overview is given of the other specimens examined for this comparative study, sourced from several other localities in the New Zealand region. Comparative metrics of selected characters from this range of specimens are provided in Table 2.
A. Description of sequenced species from Singapore and Three Kings Shelf The descriptions below are based on small voucher specimens retained from fragments used for sequencing.
Rhabdopleura sp. from Pulau Tekukor, Singapore (Fig. 3c, 5a-d; Table 1) Description of NHMO H2136 (University of Oslo): Colony on eroded coral surface, comprising adherent creeping tubes from which arise, at right angles, 31 erect tubes of varying length (most broken) up to 980 μm, the longest with 32 irregular full fusellae. Diameter of erect tubes 125-178 μm between fusellar collars (mean 148 μm), with tube openings subcircular to elongate-oval. Fusellus height variable (19-44 μm, mean 31 μm), fusellar collars irregular, projecting 5-13 μm (mean 8 μm) circumferentially from the tube. Creeping tubes comprising a generally parallel-sided convex median linear portion that is flanked by a flat basal layer, of varying width and unpatterned, on either side (see Fig. 2d). Convex part of creeping tube 144-263 μm (mean 212 μm). Surface fusellar pattern of creeping tubes variable; in places fusellar sutures are irregularly annular but they are more often oblique, forming a zigzag pattern in which the points of the zigzag are either at or very near the edges of the median linear part of the tube, or nearer its middle, forming parallel oblique strips in a criss-cross pattern (see Jain et al. 2022, Fig. 6e-g). Spacing between the latter parallel strips 42-106 μm (mean 75 μm). Zigzag angles (point-to-point angle of divergence) 28-44°( mean 36°). Width of stolon within creeping tube 19-57 μm (mean 35 μm).

Remarks
The voucher specimen was re-examined to obtain additional measurements (cf. Jain et al. 2022). The identity of the Singapore species remains in doubt. It is not identical to the holotype specimen of R. annulata, described from near the Three Kings Islands, New Zealand, at 549 m. Although the point-to-point angles of divergence of the zigzag sutures on the creeping tubes (Char 4) are practically identical, the number of fusellar collars in 500 μm (Char 2) is non-overlapping. Erect tubes of NHMO H2136 also have a smaller average diameter. Ramírez-Guerrero et al. (2020) compared their material from Heron Island, Great Barrier Reef, with Norman's (1921) description of R. annulata and concluded that they had that species. Indeed, their figure 1C, which shows a wide zigzag suture on a creeping tube, is very close to that of Norman (1921, Fig. 3). Our phylogeny (Fig. 4a) shows that NHMO H2136 is not conspecific with R. annulata from Heron Island. On the other hand, it is closely related to a new specimen from The Three Kings Shelf, sequenced for this study (below). Rhabdopleura sp. from Three Kings Shelf, New Zealand (Fig. 3i, Fig. 4 Rhabdopleura phylogeny reconstructions based on a concatenated 16S+18S rRNA genes and b mitochondrial protein-coding genes. Numbers adjacent to nodes represent bootstrap support values followed by posterior probabilities of creeping tube 119-244 μm (mean 176 μm). Surface fusellar pattern of creeping tubes variable in vicinity of erect tubes and where branching of creeping tubes takes place; on long stretches of creeping tube fusellar strips are oblique and parallel, forming a longitudinal series of intersecting chevrons. Spacing between parallel strips 46-80 μm (mean 64 μm). Zigzag angles (point-to-point angle of divergence) 58-70°(mean 65°). Width of stolon within creeping tube 25-27 μm (mean 26 μm). Zooid individuals not seen within erect tubes.

Remarks
Although NHMO H2137 is from the same general area (Three Kings Shelf) as R. annulata, it appears not to be conspecific. Metric characters 1-4 (Table 2) are all non-overlapping and the zigzag suture of the former is confined to the middle part of the creeping tube. In comparing NHMO H2137 to the other specimens studied here, it appears closest to NIWA 90266 from Cavalli Seamount c. 200 km to the east. Although the zigzag pattern appears very similar, there is in fact no overlap in Char 4 (point-to-point angle of divergence) between the two morphologies, which suggest they might also not be conspecific.

Remarks
Surface debris obscured the fusellar patterning of the creeping tubes in SEM images, but this pattern was able to be seen in a light micrograph (Fig. 3a, 6b). It shows very wide curved zigzag sutures that occupy the full width of the creeping tube, with the zigzag angles extending onto the flat basal layer. What is surprising is that cotype specimen NHMUK 1921.7.28.1 (Fig. 2a, 3d, 6c), from a different station on the Three Kings Shelf, differs from the holotype in some respects (compare data columns 1-4 in Table 2). Character 2 (number of fusellar collars in 500 μm) is non-overlapping. Furthermore, the zigzag pattern on the creeping tube of the cotype is not as wide as in the holotype. While the amount of material that could be separated from the holotype and cotype of R. annulata for SEM examination is admittedly limited, the metric data and difference in fusellar zigzag patterning suggest that the two type specimens might not be conspecific. NHMO H2137 from the Three Kings Shelf may represent a third species from the same general area.
C. Other Rhabdopleura specimens from the New Zealand region (Fig. 2c, e, 3b, e-h, j, 7a-h, 8a-h; Table 1) Data pertaining to these specimens show a singular lack of conformity for all characters. While some characters clearly overlap or coincide, others do not. For example, NIWA 158518 from north of subtropical Raoul Island on the Kermadec Ridge (c. 29°S) is closest to NIWA158791 from the Otago Shelf (c. 46°S), and the two entities may be conspecific. They differ most in the average distance  Fig. 2a; note that the tips of the zigzags are set back from the tube margins. Scale bars: a-c 200 μm between fusellar collars (Char 3), but the ranges overlap. NIWA 90264 from Spirits Bay (SE part of Three Kings Shelf) differs from the material from near Raoul Island in character 2, which is non-overlapping, but it could be conspecific with NIWA 158517 from the Snares Plateau (c. 48°S); all characters overlap and the zigzag sutures are similar (cf. Fig. 3g, h (Snares) and Fig. 6b (Spirits  Bay)). NIWA 90266 from Cavalli Seamount has the highest average number of fusellar collars (28) in 500 μm; it could be related to NHMO H2137 from the Three Kings Shelf but character 4 (point-to-point angle of divergence in fusellar zigzags) is non-overlapping. Finally, NIWA 158156 from the subantarctic Campbell Plateau differs in several significant respects, having the largest erect-tube diameter of all other samples in Table 2. It matches NHMO H2137 in having larger zigzag angles (Char 4), but the latter has a much narrower tube diameter.

Remarks
Comparison of metric and non-metric data in Table 2 for all eleven samples suggests either that morphological diversity in Rhabdopleura is extremely variable or that there might be many cryptic species. Clearly, a major constraint in studies of Rhabdopleura, including our own, is the availability of specimens from a wide range of localities amenable to molecular sequencing to compare with morphological results.

Discussion
This study compares Rhabdopleura specimens from a wide range of latitudes (c. 29-49°S) and depths (57-630 m) in the New Zealand region, plus Singapore, using six morphological characters. We consider four of these (2, 4-6) to be effectively new. Character 2, pertaining to the number of fusellar collars over a fixed distance of 500 μm in the distal part of an erect tube proximal to the aperture, is supported by the statement in  that early fuselli of an erect tube can be of a lesser height than distal or apertural fuselli but a lesser height may also be seen at thecal apertures. Concerning characters 4-6, although authors have historically noted the zigzag pattern within fusellar sutures on creeping tubes, they have not taken this further by measuring zigzag angles or describing the range of variation shown by the sutures (which can also vary in thickness). The overlooked potential usefulness of the pattern of creeping-tube sutures may be attributable to paucity of material in many instances, though it would be useful to apply it and some of the other new characters to populations of The tube opening in a more magnified; c a section of the tube in a just distal to midlength; d a short tube (possibly broken), with curled annuli in the part near the substratum; e opening of another unbroken erect tube; f stolon emergent from a cellular mass, with some globular components, in a slightly swollen part of a creeping tube; g, h variations in creeping-tube width and fusellar patterning, the latter typically weakly and thinly developed throughout the whole colony. Scale bars: a 1 mm; b-f 100 μm, g, h 200 μm R. annulata along European coasts, for example. The new characters used here are suggested to be potentially useful for morphological discrimination within the genus, the results of which could be tested in parallel with molecular sequencing pending the discovery of fresh material. When we began our study, it was anticipated that examination of NIWA collection samples from a range of localities in the New Zealand region would suggest at least one widespread species. The sequence data and morphology indicate that there are at least two species on the Three Kings Shelf (assuming that Heron Island R. annulata is indeed this species-see below). We cannot conclusively say if there is one highly variable species or that there are instead several cryptic species. The latter hypothesis seems the more likely, but what is now needed is preserved new material of each of the variants for sequencing, which can be accompanied by detailed measurements of within-colony variation in the same samples. Norman (1921) had suggested that New Zealand R. annulata was conspecific with Siboga Expedition material from Indonesia that was described by Harmer (1905) but left in open nomenclature. This conclusion seemed supported by the attribution of Heron Island material to R. annulata by Ramírez-Guerrero et al. (2020). In our study, molecular data for the Singapore specimen does not match that from Heron Island. Instead, it is sister to new material from northern New Zealand (NHMO H2137) in the mitochondrial genome phylogeny (Fig. 4b), and both are closely related to R. normani (Fig. 4a) in the ribosomal phylogeny.
Our study finds that the length of mitochondrial genomes of our material is similar to R. compacta but that gene order is slightly different. Specifically, for the two samples sequenced in this study, the position of COII was between ATP6 and NAD6, whereas in R. compacta it was between COIII and ATP8 (Fig. S1). While some invertebrates exhibit extensive mitochondrial gene rearrangements even between congenerics (Chen et al. 2018;Tempestini et al. 2020), others demonstrate remarkable conservation of gene order over millions of years (Ren et al. 2010;Quek et al. 2021). Mitochondrial gene arrangements may serve as important phylogenetic characters, particularly in deep phylogenies (Dowton et al. 2002;Zhang et al. 2018). The similarity in gene arrangement between the two mitochondrial genomes suggests that the two taxa diverged more recently when compared to the other members of Rhabdopleura, which was corroborated by our mitochondrial genome phylogeny (Fig. 4b). Furthermore, interspecific pairwise distances for both ribosomal genes (16S and 18S) were lower between New Zealand and Singapore samples (Tables S2 and S3), cementing the close evolutionary relationship between the two samples. Nevertheless, intraspecific distances for 16S rRNA in R. recondita were lower than 4%, whereas the smallest interspecific distance stood at 13.5% between Rhabdopleura species (Tables S2 and S3). In contrast, 18S rRNA was much less variable, with intraspecific distance below 0.5%, whereas interspecific distances never went above 2.5%. Taken together, this suggests that 16S might be a useful barcoding gene for this clade, although more samples are required in order to determine an appropriate barcoding gap distance, particularly from conspecifics (reviewed in DeSalle and Goldstein 2019).
If, as appears the case from the morphological data, there is more than one species on the Three Kings Shelf, that would be consistent with what is already known about the biodiversity of the Three Kings Shelf/Spirits Bay area, which is considered New Zealand's marine-biodiversity hotspot (Cryer et al. 2000;Taylor and Gordon 2003). For example, Cryer et al. (2000) reported~300 species of Bryozoa and~215 species of sponges from an area of about 200 km 2 north of Spirits Bay, including many multispecies genera. Indeed, five closely related species of the bryozoan genus Steginoporella, all sequenced (Orr et al. 2022), are found in the vicinity of the Three Kings Islands. Simakov et al. (2015) estimated the origin of Rhabdopleura and the rest of the hemichordates to be at 559 Ma, placing it in deep evolutionary time. Fossils of putative Rhabdopleura and other pterobranchs have been discovered as early as Cambrian Series 3 (509 Ma) (Durman and Sennikov 1993;Maletz and Steiner 2015;Sennikov 2016). The small number of nominal species described so far could be an artefact of the difficulty in identifying them, given their unostentatious morphology (Chapman et al. 1995). Of course, the question remains-is Rhabdopleura as highly conserved genetically as it is morphologically? Based on our results, it appears that 16S rRNA could be a useful marker for identifying different species once sufficient sampling has been conducted for this enigmatic clade. If morphological variability suggests cryptic or otherwise unrecognized species, hinted at by differences in widths of erect-and creeping tube fusellar bands and the zigzag angles on the creeping tubes, then it is interesting to speculate on what might generate these differences. Inasmuch as fusellar strips are formed from secretions of the head shield of the zooids stretching from the erecttube and creeping-tube apertures (Dilly 1986(Dilly , 1993Rigby and Dilly 1993;Maletz and Cameron 2016), then head-shield size, and individual zooid movements, may be correlated with fusellar-strip width. A study of tube-building in R. normani led Dilly (1986) to postulate that the "beautiful symmetry and regularity" of fusellar patterns may be the result of "a restricted and limited programme of tube-building available to the zooid". Hence, size and behavioural differences among rhabdopleurid zooids may be species-specific. If fusellar zigzag patterns (inter alia) have taxonomic significance at the species level, then such characters will have applicability to fossil rhabdopleurids as well as living ones, as they can be preserved in thin organic films (Maletz and Steiner 2015;Sennikov 2016, pl. 1) and by bioimmuration (Taylor and Todd 2001, Fig. 3.2.5.2d).
It should be noted, however, that variations in zooid behaviour and tube morphometry might be caused by different environmental conditions (e.g. water currents, food availability, and sedimentation rates). For example, Stebbing (1970b) observed in R. compacta that fusellus height in erect tubes can reflect degeneration and regeneration of zooids that may in part be seasonally related. In this species, a number of fuselli of normal height may be succeeded distally by shorter fuselli, which then gradually increase in size towards the distal aperture. He attributed this increase in annulus height to increasing zooid size following regeneration. Degenerating zooids were commonly seen in autumn, but Stebbing (1970a) also noted that embryos can be produced throughout the year, presumably related to continuing food availability and/or or nutrient reserves in so-called dormant buds. Ecological studies and growth studies are desirable for the commonest Rhabdopleura species. For example, if R. normani is as widespread as reported, then it may be fairly eurythermal and less subject to seasonal degeneration than range-restricted R. compacta, which in any case may depend less on the physiological effects of seasonal temperature than the seasonal availability of planktonic food.
Larval traits indicate an important role in dispersal and offspring survival in Rhabdopleura. The larva is yolky, well-ciliated, non-planktonic and relatively short-lived (Stebbing 1970b(Stebbing , 1973Dilly 1973;Lester 1988;Sato et al. 2008;Strano et al. 2019) and relatively few are produced per zooid. Stebbing reported that larvae of R. compacta are ready to settle up to 2 days after release, but they could delay metamorphosis for up to 10 days (Sato et al. 2008), after which they apparently have reduced ability to metamorphose. According to Lester (1988), swimming larvae of R. normani settle within 24 h after being released from the parent tube and metamorphosis takes place over 7-10 days. Strano et al. (2019) found that metamorphosis in R. recondita took 10 days in laboratory conditions at 21°C. It is clear that Rhabdopleura larvae are lecithotrophic and have only enough nutrient reserves to sustain a period of swimming for 1-2 days followed by a period of metamorphosis of up to 10 days. Accordingly, reproductive effort per offspring is higher, and larval mortality is assumed to be lower than it would be if planktotrophic. In reviewing what was known about larval ecology in marine invertebrates and the implications for paleobiology, Jablonski and Lutz (1983) observed that while species with lecithotrophic larvae appear to be typically geographically more restricted, geologically short-ranging and exhibit high speciation and extinction rates as a consequence of their characteristically low larval dispersal capabilities, there are exceptions, as subsequent studies have confirmed. For example, whereas all ascidians are lecithotrophic, those taxa with the highest speciation rates are those with the largest range sizes, while those with the least dispersive larval development have the lowest speciation rates and smallest geographic ranges (Maliska et al. 2013), contrary to findings from echinoids and gastropods in the fossil record. A comprehensive study of sacoglossan gastropods determined that lecithotrophy evolved at least 27 times in this group and that while lecithotrophic lineages arose often, they rarely persisted long enough to diversify (Krug et al. 2015). Rhabdopleuridae differs from such clades as Ascidiacea and Sacoglossa, and indeed most fossilizable clades, for being far less morphologically diverse over geological time, which means that lecithotrophy has likely been the sole reproductive mode throughout the Phanerozoic.
Remarkably, this low-diversity clade has persisted. It remains to be demonstrated if the clade is truly of such low diversity as appears at the present day. Our study suggests more morphological diversity exists than is generally appreciated but how this relates to genetic diversity remains to be seen. Previous studies sequencing the genomes of low-diversity, ancient lineages such as tuatara, coelacanth and horseshoe crab, found slow rates of evolution (Amemiya et al. 2013;Gemmell et al. 2020;Shingate et al. 2020). Nevertheless, mitochondrial rearrangements in pterobranchs were previously suggested in Li et al. (2019) to possibly correlate with higher evolutionary rates such as in ascidian tunicates. Even between congeners in Rhabdopleura, we found a minor rearrangement of mitochondrial gene order (Fig. S1). Taken together, we suggest that this clade is likely more speciose than is currently sampled, and more extensive sampling efforts of both morphology and genomics (i.e., whole genome sequencing) may reveal previously undiscovered species of this inconspicuous but highly distinctive creature.

Declarations
Conflict of interest The authors declare no competing interests.
Ethical approval All applicable international, national and/or institutional guidelines for animal testing, animal care and use of animals were followed by the authors.
Sampling and field studies All necessary permits for sampling and observational field studies have been obtained by the authors from the competent authorities and are mentioned in the Acknowledgements, if applicable. The study is compliant with CBD and Nagoya protocols.
Data availability Numerical data are tabulated in the article. Genetic data are published in GenBank and accession numbers are available in Table S1. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.