The Morphology, Mitogenome, Phylogenetic Position, and Symbiotic Bacteria of a New Species of Sclerolinum (Annelida: Siboglinidae) in the South China Sea

Sclerolinum annulatus n. sp. (Annelida: Siboglinidae) is described based on specimens collected from soft sediment of the Haima cold seep in the South China Sea. Morphologically, S. annulatus n. sp. is distinct in having a tube with transverse rings and a forepart (i.e., anterior region) containing one arched row of elongated plaques on both sides of the dorsal furrow. Genome skimming, assembly, and annotation produced a nearly complete mitogenome of S. annulatus n. sp. with 15,553 bp nucleotides that encodes 13 protein-coding genes, two rRNA, and 22 tRNA. Phylogenetic analyses based on the mitochondrial cytochrome c oxidase I (cox1) gene and a concatenated dataset comprising the mitochondrial cox1 and 16S rRNA genes along with the nuclear 18S rRNA gene both strongly support the placement of S. annulatus n. sp. in the genus Sclerolinum Southward, 1961. Based on cox1, S. annulatus n. sp. is most closely related to an undescribed siboglinid from off Kushiro in Japan (“Pogonophora” sp. Kushiro-SK-2003). Transmission electron microscopy, microbial 16S rRNA amplicon sequencing, phylogenetic reconstruction, and stable isotope analyses together indicate that S. annulatus n. sp. hosts a single phylotype of sulfur-oxidizing endosymbionts.

During a remotely operated vehicle (ROV) cruise in 2019, we discovered dense aggregations of Sclerolinum tubeworms in the Haima cold seep (Figure 2). The objectives of this study were, therefore, to describe the morphological characteristics of this Sclerolinum species, assemble its mitogenome, determine its phylogenetic relationships with other siboglinids, and characterize the composition and trophic mode of its symbiotic bacteria.

Sample Collection and Preservation
Samples were collected in May 2019 using the ROV Haima onboard the research vessel (R/V) Haiyang 6 of Guangzhou Marine Geological Survey (Guangzhou, China) from the Haima cold seep in the South China Sea (represented by a red asterisk in Figure 1; 16 • 54.04 N, 110 • 28.45 E, 1,433 m water depth). All the retrieved samples are incomplete with the anterior and/or trunk regions only. The samples were rinsed several times with chilled seawater to remove the attached sediment. Some individuals (within their tubes) were preserved in 100% ethanol or 10% formalin diluted with Milli-Q water for morphological investigation. Other individuals (also within their tubes) were frozen at −80 • C for molecular and stable isotope analyses. Type specimens have been deposited in Tropical Marine Biodiversity Collections of the South China Sea, Chinese Academy of Sciences (Guangzhou, China) as vouchers (Catalog numbers: TMBC030854-TMBC030861).

Morphological Characterization
Morphological features were examined under a Motic SMZ-171 Stereo Microscope (Motic, Xiamen, China), and photographed with a Canon 700D Camera (Canon, Tokyo, Japan) mounted on the microscope. To visualize morphological details, different parts of the tubes and the extracted soft body were carefully dissected for scanning electron microscopy (SEM) in Hong Kong Baptist University. In brief, the dissected tubes and tissues were treated with a gradient of hexamethyldisilazane-ethanol solutions (i.e., 50, 75, and 100%, each for 10 min), dried inside a fume hood, mounted on conductive carbon adhesives, sputter-coated with gold, and observed under a LEO 1530 Field Emission Scanning Electron Microscope (LEO Elektronenmikroskopie GmbH, Oberkochen, Germany).
To examine the symbiont morphology, the trophosome tissue (along with the tube) of one specimen fixed in 10% formalin was cut into small blocks of 2 mm in length, and sent to FIGURE 1 | A map showing the known distributions of seven described Sclerolinum species, the type locality of Sclerolinum annulatus n. sp. (i.e., the Haima cold seep in the SCS as indicated by a red asterisk), along with an undescribed species of Sclerolinum "Pogonophora" sp. Kushiro-SK-2003. The map was created using Ocean Data View (ODV) v.5.0 (https://odv.awi.de). Color coding represents the water depth (unit: m). Ant, Antarctic; CS, Caribbean Sea; EI, East Indies; EQ, equator; GoM, Golf of Mexico; JT, Java Trench; Kus, Kushiro; NF, Norwegian fiords; NS, Norwegian Sea; SCS, South China Sea. Notes: the geographic coordinates of Sclerolinum javanicum inhabiting the Java Trench and "Pogonophora" sp. Kushiro-SK-2003 from off Kushiro in Japan are not available, therefore, their distributions are indicative only. Electron Microscope Unit of The University of Hong Kong for transmission electron microscopy (TEM). Briefly, the tissue blocks were changed to cacodylate buffer (0.1 M sodium cacodylate-HCl buffer pH 7.4) containing 0.1 M sucrose to remove the fixative, post-fixed with 1% osmium tetroxide in cacodylate buffer for an hour at room temperature, washed in three changes of cacodylate buffer, then dehydrated with a gradient of ethanol solutions (i.e., 50, 70, and 90% all diluted in Milli-Q water, each for 5 min; 100%, three changes with each for 10 min). The samples were afterward exposed to propylene oxide (two changes with each for 5 min), infiltrated in epoxy resin/propylene oxide 1:1 mixture for one and a half hour followed by epoxy resin for an hour both at 37 • C, and then embedded in plastic capsules that filled with epoxy resin. After polymerization at 60 • C overnight, the embedded tissue blocks were cut into ultrathin sections in 100-nm thickness via a Leica Ultracut UCT Ultramicrotome (Leica, Wetzlar, Germany), mounted on the 150-mesh hexagonal copper grids, stained with 2% aqueous uranyl acetate for 20 min, washed well in running Milli-Q water, stained with Reynold's lead citrate for 15 min, and thereafter examined under a Philips CM 100 Transmission Electron Microscope (Philips, Amsterdam, Netherlands).

DNA Extraction, Sequencing, and Assembly
Genomic DNA was extracted from individuals T1-T3 (along with their tubes due to the difficulties to extract the soft body from the tube) using the CTAB method (Stewart and Via, 1993). Integrity of DNA was monitored with 1% agarose gel electrophoresis. Purity of DNA was checked using a NanoPhotometer R Spectrophotometer (IMPLEN, Calabasas, CA, United States). Concentration of DNA was measured using the Qubit R dsDNA BR Assay Kit (Invitrogen, Carlsbad, CA, United States) together with a Qubit R 2.0 Fluorometer (Invitrogen, Carlsbad, CA, United States). One DNA sequencing library with an insert size of approximately 350 bp was constructed with the NEBNext R Ultra TM DNA Library Prep Kit for Illumina R (NEB, Ipswich, MA, United States) for each individual, and was then sequenced on an Illumina R NovaSeq 6000 Sequencer (Illumina, San Diego, CA, United States) in Novogene (Beijing, China) to generate 150 bp paired-end reads. Raw reads were filtered using Trimmomatic v.0.38 (Bolger et al., 2014) to remove adapters and low-quality sequences with the following settings: ILLUMINACLIP: Adapters.fas:2:30:10, LEADING = 10, TRAILING = 10, SLIDINGWINDOW = 4:15, and MINLEN = 40.
Obtained high-quality reads of each individual were de novo assembled using CLC Genomics Workbench v.20.0.4 (QIAGEN Bioinformatics, Aarhus, Denmark) under the default settings. The putative mitochondrial contig and the nuclear 18S rRNA gene were identified by applying the BLASTn option of BLAST + v.2.10.1 (Camacho et al., 2009) to search against the mitogenome and 18S rRNA gene sequences of other siboglinids retrieved from GenBank 1 with an E-value threshold of 1e-50. The putative mitochondrial contig of each specimen obtained via CLC Genomics Workbench v.20.0.4 was then served as the bait for its mitogenome assembly using GetOrganelle v.1.7.5 (Jin et al., 2020) under the default settings. Alignment of the obtained mitogenomes of these three specimens was performed using MUSCLE (Edgar, 2004)

Mitogenome Annotation
Annotation of the mitogenome was performed using MITOS2 Web Server (Donath et al., 2019) under the default settings, except to change the genetic code to "5 invertebrate". Resultant boundaries of each protein-coding gene (PCG) and rRNA were further examined and adjusted by aligning with the PCGs and rRNAs extracted from several published mitogenomes of siboglinids, including S. brattstromi (Li et al., 2015), using MUSCLE implemented in MEGA X under the default settings. A gene map of the mitogenome was constructed using CGView Server (Grant and Stothard, 2008). Base constitutions, codon compositions, as well as relative synonymous codon usage (RSCU) patterns were computed using MEGA X. Alignment of each mitochondrial PCG of our specimen with that derived from the mitogenome of S. brattstromi (Li et al., 2015) was carried out using MUSCLE in MEGA X under the default settings. The Kimura-2-parameter (K2P) model (Kimura, 1980) in MEGA X was then adopted to calculate the pairwise genetic distance of each mitochondrial PCG between these two species.

Genetic Distance Estimation and Phylogenetic Analyses of Siboglinids Based on the Cox1 Gene
Sequences of the mitochondrial cytochrome c oxidase I (cox1) gene of 38 siboglinids representing 14 genera and two cirratuliformians (serving as the outgroup) were downloaded from GenBank for genetic distance estimation and phylogenetic analyses, with accession numbers provided in Supplementary Table 1.
To estimate the genetic distances, alignment of the cox1 gene was conducted using MUSCLE in MEGA X under the default settings. Poorly aligned positions were eliminated by Gblocks Server 2 under the default settings. The K2P model in MEGA X was then used for genetic distance calculation.
To analyze the phylogenetic relationships, alignment of the cox1 gene was also performed using MUSCLE in MEGA X under the default settings. Nevertheless, to retain as much nucleotide information as possible for phylogenetic reconstruction, poorly aligned positions were removed using Gblocks Server under the relaxed settings, which allow smaller final blocks, gap positions within the final blocks, and less strict flanking positions. Afterward, phylogenetic analyses were carried out using the maximum-likelihood (ML) approach implemented in IQ-TREE v.1.6.12 (Nguyen et al., 2015) and the Bayesian inference (BI) approach implemented in MrBayes v.3.2.7 (Ronquist et al., 2012). The best-fitting nucleotide-substitution model was determined based upon the Bayesian information criterion (BIC) via ModelFinder (Kalyaanamoorthy et al., 2017) implemented in IQ-TREE v.1.6.12. If the model selected by ModelFinder was not available in MrBayes v.3.2.7, the closest one available was used instead.
In particular, the ML analysis was performed by applying the GTR + F + I + G4 nucleotide-substitution model to run 1,000 replicates of the Shimodaira-Hasegawa-like (SH-like) approximate likelihood ratio tests (Guindon et al., 2010) together with 1,000 standard non-parametric bootstraps (NB) using IQ-TREE v.1.6.12. The BI analysis was executed by utilizing the GTR + I + G nucleotide-substitution model, running four independent Markov chains for 10 million generations, and sampling every 1,000 generations with the first 25% discarded as the burn-in using MrBayes v.3.2.7.

Phylogenetic Analyses of Siboglinids Based on a Concatenated Dataset of Three Genes
Sequences of the mitochondrial 16S rRNA and nuclear 18S rRNA genes of 26 siboglinids representing 11 genera and two cirratuliformians (serving as the outgroup) were also downloaded from GenBank for phylogenetic analyses based on the concatenated dataset of three genes (i.e., cox1 + 16S rRNA + 18S rRNA), with accession numbers summarized in Supplementary Table 2. Alignment and trimming of each gene followed the same methods detailed in the molecular analyses of cox1. Concatenation of the three gene alignments was achieved using SequenceMatrix v.1.7.8 (Vaidya et al., 2011). Methods for phylogenetic reconstruction of the concatenated three genes were identical to those of cox1, except that the bestfitting nucleotide-substitution model was assessed for each gene alignment and then applied to each gene partition. Particularly, the nucleotide-substitution model applied to the ML analysis was GTR + F + I + G4 for cox1, GTR + F + I + G4 for 16S rRNA, and TIM3e + I + G4 for 18S rRNA, and that applied to the BI analysis was GTR + I + G for cox1, GTR + I + G for 16S rRNA, and GTR + I + G for 18S rRNA.

Microbial 16S rRNA Amplicon Sequencing and Analyses
Genomic DNA was extracted from individuals T4-T6 (along with their tubes due to the difficulties to extract the small soft body from the tube) for sequencing the V3-V4 region of the microbial 16S rRNA gene. Integrity, purity, and concentration of DNA were examined using 1% agarose gel electrophoresis, a NanoPhotometer R Spectrophotometer (IMPLEN, Calabasas, CA, United States), and the Qubit R dsDNA BR Assay Kit (Invitrogen, Carlsbad, CA, United States) together with a Qubit R 2.0 Fluorometer (Invitrogen, Carlsbad, CA, United States), respectively. One DNA sequencing library was constructed with the NEBNext R Ultra TM DNA Library Prep Kit for Illumina R (NEB, Ipswich, MA, United States) with each individual labeled using a unique barcode, and then sequenced on an Illumina R NovaSeq 6000 Sequencer (Illumina, San Diego, CA, United States) in Novogene (Beijing, China) to generate 250 bp paired-end reads. Data filtering was carried out using the same method depicted in the section of DNA Extraction, Sequencing, and Assembly to remove adapters and low-quality sequences.
Thereafter, the microbial operational taxonomic unit (OTU) identification, quantification, and taxonomy prediction were performed using USEARCH v.11.0.667 (Edgar, 2010) according to the following steps: (1) high-quality read pairs were merged using the -fastq_mergepairs command under the default settings; (2) merged reads were further filtered using the -fastq_filter command, with the -fastq_maxee parameter set to 1; (3) unique reads were obtained after removing the replicate and singleton reads with the -fastx_uniques command; (4) OTU assignment was achieved based on the UPARSE algorithm implemented in the -cluster_otus command with a similarity cutoff of 97%, and chimeras detected in this step were also excluded from further analyses; (5) an OTU table was generated using theotutab command by mapping the merged reads to the obtained OTUs, with the mapping efficiency applied to quantify the relative abundance of each OTU; and (6) taxonomy of theses OTUs was predicted using the SINTAX algorithm (Edgar, 2016) implemented in the -sintax command to research against the Ribosomal Database Project (RDP) 16S rRNA gene training set v.16 with a bootstrap confidence cutoff of 0.80.
Sequences of the symbiont 16S rRNA gene reported from 28 siboglinids representing 12 genera and two Methylococcaceae (serving as the outgroup) were downloaded from GenBank for phylogenetic analyses. Additionally, the gill symbiont 16S rRNA gene reported from eight bivalve species representing eight genera were also included for phylogenetic inferences due to their genetic affinities to the symbionts of siboglinids based on the online BLASTn search 3 and/or previous studies (e.g., Eichinger et al., 2014;Reveillaud et al., 2018). Methods for sequence alignment, trimming, phylogenetic reconstruction, and K2P genetic distance estimation were identical to those applied to the molecular analyses of cox1. The nucleotide-substitution model utilized in the ML analysis was GTR + F + R3 and in the BI analysis was GTR.

Stable Isotope Analyses
Individuals T7-T9 (along with their tubes due to the difficulties to extract the small soft body from the tube) and individuals T10-T12 (along with their tubes due to the difficulties to extract the small soft body from the tube) were combined into two samples (i.e., Sample 1: containing individuals T7-T9; Sample 2: containing individuals T10-T12) for stable isotope analyses of carbon (C), nitrogen (N), and sulfur (S) to understand the trophic mode of their symbionts. For the C and N stable isotope analyses, the two samples were freeze-dried, homogenized, weighted using a Mettler Toledo R Microbalance (Mettler Toledo, Columbus, OH, United States), and then sealed into tin capsules. For the S stable isotope analysis, the two freeze-dried and homogenized samples were further screened with the 200mesh sieves, dried inside an oven at 60 • C, weighted using a Mettler Toledo R Microbalance (Mettler Toledo, Columbus, OH, United States), mixed with 0.2 mg vanadium pentoxide, and thereafter sealed into tin capsules. These processed samples were analyzed using a Sercon Integra2 Elemental Analyzer Isotope Ratio Mass Spectrometry (Sercon Instruments, Crewe, United Kingdom) in The Third Institute of Oceanography, Ministry of Natural Resources (Xiamen, China). Values were reported in permille ( ) with the δ notation, relative to the standards Vienna Pee Dee Belemnite (VPDB) for C, air N 2 for N, and the Canyon Diablo Triolite (CDT) for S. The precision for δ 13 C, δ 15 N, and δ 34 S determinations were all ± 0.2 .

Type Locality
The Haima cold seep, 1,433 m water depth, off southern Hainan Island, on the northwestern slope of the South China Sea.

Type Material
Type specimens are incomplete, with the anterior and/or trunk regions but lacking the opisthosomal region.

Diagnosis
Tube straight, often possessing pronounced, irregularly transverse rings, in anterior and middle regions. Cephalic lobe with a rounded triangular tip pointing anteriorly. Dorsal plaques elongated, three on each side of dorsal furrow, forming a caret shape pointing anteriorly. Ventral plaques oral, present in two irregular rows.

Description
Tubes brownish in color, straight, flexible, and elastic, protruding as dense aggregations from sediment surface to water column (Figure 2). Tube outer surface with clearly defined transverse rings in anterior (Figures 3A,F) and middle regions ( Figure 3B); rings sparsely spaced in anterior region (Figures 3A,F), closely spaced in middle region (Figure 3B), and absent in posterior region ( Figure 3C). Tube walls composed of multi-layered fibers (Figures 3D,E). Tube diameter 0.35-0.53 mm. Tube thickness 23-41 µm. Tube length up to 87.85 cm. Two wrinkled tentacles (Figures 3F,G, 4A), varying greatly in length from 1.4 to 12.5 mm according to different states of contraction among measured specimens. Forepart 75-190 µm in diameter, with a deep dorsal furrow extending posteriorly from the base of two tentacles ( Figure 4A). Epidermal glands scattered on forepart, more densely distributed anterior than posterior to bridle (i.e., thicken cuticular ridge) (Figures 4A,B). Cephalic lobe located on ventral side, with a rounded triangular tip pointing anteriorly ( Figure 4D). Bridle located approximately 0.7 mm posterior to cephalic lobe, formed by 15 elongated to oval plaques with diameter ranging from 32 to 52 µm (Figures 4A-E).
Dorsal plaques elongated, three on each side of dorsal furrow, forming a caret shape pointing anteriorly (Figures 4A-C).
Ventral plaques present in two irregularly rows, each with 4-5 plaques (Figures 4D,E). Transition between the relatively smooth forepart and the highly wrinkled and densely papillated trunk region clear. Oval plaques scattered on top of epidermal papillae in trunk region ( Figure 4F).

Distribution
Within the South China Sea, S. annulatus n. sp. has only been reported from the Haima cold seep-the type locality situated on the northwestern continental slope of the South China Sea.

Etymology
The species epithet "annulatus" refers to "ring" in Latin, which reflects the rings on the tube of this species.

Remarks
Sclerolinum annulatus n. sp. has a relatively large tube (0.35-0.53 mm diameter) among the seven described Sclerolinum species. The only two Sclerolinum species with comparable tube sizes are S. major (0.20-0.84 mm) (Southward, 1972) and S. contortum (0.20-0.61 mm) (Smirnov, 2000;Eichinger et al., 2013). Sclerolinum annulatus n. sp. stands out from all the seven described Sclerolinum species in that its tube is relatively straight and has clearly defined transverse rings in the anterior and middle regions (Figures 3A,B,F), which resembles the tubes of frenulates (Southward, 1972;Smirnov, 2000;Sen et al., 2020). In comparison, the tubes of most other Sclerolinum species are rather smooth (Southward, 1961). Although the tube of S. sibogae was described as deeply wrinkled (Southward, 1961) and those of some S. magdalenae as weakly wrinkled (Southward, 1972), they do not form complete rings. The tube of S. contortum is also different in typically having several bends (although it is smooth in one population) (Eichinger et al., 2013) or having only faint transverse rings (Georgieva et al., 2015). Furthermore, the bridle shape of S. annulatus n. sp. is unique among all the seven described Sclerolinum species in that it is composed of one obliquely arranged row of elongated plaques on each side of the furrow on the dorsum, and two rows of transversely arranged oval plaques on the ventrum (Figures 4A-E).

General Features of the Mitogenome
The assembled mitogenomes of S. annulatus n. sp. T1-T3 do not form a circle. They vary in length from 15,565 to 15,590 bp, with 15,553 bp being identical. Therefore, only the consensus nucleotides, with an overall base composition of 30.51% for A, 23.55% for C, 12.49% for G, and 33.45% for T, were retained as the mitogenome of S. annulatus n. sp. for downstream analyses. It contains 13 PCGs (i.e., atp6, atp8, cox1-3, cytb, nad1-6, and nad4l), two rRNA, and 22 tRNA which are typical of bilaterian mitogenomes (Boore, 1999), with potential missing nucleotides in the control region that lies between trnR and trnH ( Figure 5 and Supplementary Table 3). Gene order of the S. annulatus n. sp. mitogenome is identical to that of other siboglinids reported to date (Li et al., 2015;Sun et al., 2018;Zhou et al., 2020). All 13 mitochondrial PCGs use ATG as the start codon. Most mitochondrial PCGs use TAA as the stop codon, except that cox1 uses TA, cox2 uses TAG, and nad2, nad5, as well as nad6 use a single T. Comparing the 13 mitochondrial PCGs between S. annulatus n. sp. and S. brattstromi uncovered that their K2P genetic distances range from 15.00% (cox2) to 28.66% (atp8), with an average value of 23.75% (Figure 6A). Utilization of the start and stop codons of the 13 mitochondrial PCGs of S. annulatus n. sp. and S. brattstromi is almost consistent, except for the stop codons of nad1 (TAA for S. annulatus n. sp. and TAG for S. brattstromi) and nad6 (a single T for S. annulatus n. sp. and TAA for S. brattstromi). Lengths of 12 out of the 13 mitochondrial PCGs of these two species are identical. The only exception occurs in nad6, with 472 bp for S. annulatus n. sp., whereas 471 bp for S. brattstromi. With the exclusion of the stop codons, the 13 mitochondrial PCGs of S. annulatus n. sp. and S. brattstromi exhibit a total of 3,687 and 3,686 codons, respectively. The three most frequently used codons of these two Sclerolinum mitogenomes are consistent, including Leu1 (CUN), Ile (AUY), and Ser2 (UCN), each with a codons per thousand codons (CDsp T) value > 80 ( Figure 6B). Results of the RSCU patterns further illustrate that the predominant synonymous codons in the 13 mitochondrial PCGs of these two Sclerolinum species are also conserved, which generally have an over-usage of A and T at the third codon positions (Figure 6C).
Genetic Distances of Siboglinids Based on the Cox1 Gene Alignment using MUSCLE followed by trimming using Gblocks Server (under the default settings) results in a 514-bp alignment of the cox1 gene from 41 siboglinids. The K2P genetic distance of cox1 between the examined frenulates ranges from 18.02% to 24.54%, between the examined Osedax ranges from 11.39% to 24.79%, between the examined vestimentiferans ranges from 0.39% to 21.30%, and between the examined Sclerolinum ranges from 1.98% to 18.78% (Supplementary Table 4). Among them, S. annulatus n. sp. is most closely related to "Pogonophora" sp. Kushiro-SK-2003 collected via deep-sea trawling from an unknown habitat off Kushiro in Japan (1,500 m water depth) (Kojima et al., 1997), with a K2P genetic distance of 1.98%. By contrast, S. annulatus n. sp. is distantly related to S. brattstromi and S. contortum with a K2P genetic distance ranging from 17.16% to 18.78%, respectively. Furthermore, the K2P genetic distance between S. contortum collected from four populations from the Arctic to the Antarctic ranges from 0.19% to 1.38%.   (Figure 7). Specifically, in the Sclerolinum clade, S. annulatus n. sp. T1-T3 form a small clade that is sister to "Pogonophora" sp. Kushiro-SK-2003. These four sequences form a well-supported clade (ML-SH/ML-NB/BI: 100/100/1.00) that is sister to another well-supported clade (ML-SH/ML-NB/BI: 89.9/95/0.99) of Sclerolinum, which includes S. brattstromi from Norwegian fiords along with four S. contortum from the Arctic to the Antarctic. Alignment using MUSCLE followed by trimming using Gblocks Server (under the relaxed settings) resulted in a 1,211bp alignment of cox1, a 532-bp alignment of 16S rRNA,  and a 1,735-bp alignment 18S rRNA from 29 siboglinids and two cirratuliformians. Phylogenetic reconstruction based on this 3,478-bp concentrated dataset also unveiled four major clades of siboglinids (Figure 8). Among the Sclerolinum clade (ML-SH/ML-NB/BI: 100/100/1.00), S. annulatus n. sp. T1-T3 form a well-supported clade (ML-SH/ML-NB/BI: 100/100/1.00) that is sister to another well-supported Sclerolinum clade (ML-SH/ML-NB/BI: 98.8/100/1.00) comprising S. brattstromi from Norwegian fiords and two S. contortum from the Antarctica and the Haakon Mosby Mud Volcano. Compared to the cox1 tree (Figure 7), the support values of the four major siboglinid clades in the three-gene tree have in general increased (Figure 8). Nevertheless, phylogenetic relationships among the four siboglinid clades remain unsolved, with the support values between Osedax and the other three clades being quite low (Figures 7, 8).
Alignment using MUSCLE followed by trimming using Gblocks Server (under the relaxed settings) resulted in a 1,434-bp alignment of the microbial 16S rRNA gene from 29 siboglinids, eight bivalves, along with two Methylococcaceae outgroups. Phylogenetic analyses based on the microbial 16S rRNA gene elucidated the symbionts of Sclerolinum as a well-supported clade (ML-SH/ML-NB/BI: 98.4/100/1.00), which contains OTU-1 of S. annulatus n. sp. T4-T6 and two sulfuroxidizing symbionts hosted by S. contortum from the Gulf of Mexico and the Haakon Mosby Mud Volcano (Figure 10 and Supplementary Figure 1). Similarly, the symbionts of frenulates (ML-SH/ML-NB/BI: 99.9/100/1.00) and Osedax (ML-SH/ML-NB/BI: 100/100/1.00) each form a well-supported clade (Figure 10 and Supplementary Figure 1). However, the symbionts of vestimentiferans are paraphyletic, with 24 forming a clade that is sister to a clade of three bivalve symbionts and four being nested between the bivalve and Sclerolinum symbionts (Figure 10 and Supplementary Figure 1). Alignment using MUSCLE followed by trimming using Gblocks Server (under the default settings) resulted in a 429bp alignment of microbial 16S rRNA gene of five Sclerolinum. Among them, OTU-1 of S. annulatus n. sp. T4-T6 exhibit high genetic similarities, with a K2P genetic distance ranging from 0 to 0.47% (Supplementary Table 9). In addition, these OTU-1 are genetically closely related to the sulfur-oxidizing symbionts of S. contortum, with a K2P genetic distance ranging from 0.47% to 0.94% (Supplementary Table 9).

DISCUSSION
In the present study, we described the morphological features of Sclerolinum annulatus n. sp., obtained its nearly complete mitogenome, determined its phylogenetic position in the family Siboglinidae Caullery, 1914, and provided evidence for its endosymbionts inside the trophosome being sulfur oxidizers. Morphologically, S. annulatus n. sp. is distinguishable from its congeneric species in that its tube has clearly defined transverse rings similar to the rings on the frenulate tubes (Figure 3) and the unique arrangement of plaques in its bridle (Figure 4). Phylogenetic reconstruction inferred from two datasets (Figure 7: cox1 and Figure 8: a concatenated dataset of cox1 + 16S rRNA + 18S rRNA) and the K2P genetic distance estimation based on the cox1 gene (Supplementary Table 4) support the recognition of this new species as three individuals are nested within a clade of two other described Sclerolinum species, while they exhibit substantial interspecific K2P genetic distances with these two species. Besides, phylogenetic analyses based upon cox1 indicate that "Pogonophora" sp. Kushiro-SK-2003(Kojima et al., 1997 belongs to Sclerolinum (Figure 7). However, it has a K2P genetic distance of 1.98% from S. annulatus n. sp. (Supplementary Table 4). In a large-scale DNA barcoding study that involved amplification of cox1 from 1,876 specimens representing 333 provisional species, Carr et al. (2011) found an average intraspecific K2P genetic distance of 0.38%, and in most species, this distance is much less than 2%. Therefore, it is likely that "Pogonophora" sp. Kushiro-SK-2003 represents an undescribed species of Sclerolinum in the Northwest Pacific.
The complete mitogenomes of siboglinids are circular (Li et al., 2015), whereas that of S. annulatus n. sp. is not, indicating the latter may be incomplete. Nevertheless, given that the mitogenome of S. annulatus n. sp. contains all 37 genes typical in bilaterian mitogenomes (Boore, 1999) with a length close to that of S. brattstromi (15,553 bp for S. annulatus n. sp. vs. 15,383 bp for S. brattstromi), it should be nearly complete, with a small gap probably located in its control region. As having been reported in some siboglinid mitogenomes, the control region is difficult to be assembled using short reads alone due to the occurrence of repetitive sequences (Boore and Brown, 2000;Jennings and Halanych, 2005;Li et al., 2015). Further longread sequencing based assembly approach would help to address this issue in the future. The mitogenome of S. annulatus n. sp. possesses the same gene order as S. brattstromi, which is also identical to that of other siboglinid mitogenomes (Li et al., 2015;Sun et al., 2018;Zhou et al., 2020). In addition, the codon usage patterns along with RUSC of the S. annulatus n. sp. and S. brattstromi mitogenomes are highly conserved, despite the large K2P genetic distances between all their 13 mitochondrial PCGs (Figure 6).
Compared to the symbionts of other siboglinids, very little is known about the symbionts of Sclerolinum species, with only that of S. contortum having been characterized (Pimenov et al., 1999(Pimenov et al., , 2000Lösekann et al., 2008;Eichinger et al., 2014). Earlier studies based on radiotracer determination (Pimenov et al., 1999) and TEM (Pimenov et al., 2000) suggested that S. contortum inhabiting seep sediment of the Haakon Mosby Mud Volcano may host methane-oxidizing symbionts in its trophosome. By contrast, more recent phylogenetic inferences, stable isotope analyses, and fluorescent in situ hybridization experiments unveiled a single phylotype of sulfur-oxidizing symbionts inside the trophosome of S. contortum dwelling in seep sediment of the Haakon Mosby Mud Volcano (Lösekann et al., 2008) and the Gulf of Mexico (Eichinger et al., 2014). These different lines of evidence underline that the symbionts of S. contortum reported in Pimenov et al. (1999Pimenov et al. ( , 2000 might also be thioautotrophs with an unusually light δ 13 C signature and an exceptional morphology that are atypical in sulfur-oxidizing bacteria. Alternatively, the S. contortum examined by Pimenov et al. (1999Pimenov et al. ( , 2000 might be another siboglinid species that is morphologically similar to S. contortum whereas hosts methanotrophs.
In the present study, the presence of intracellular symbionts inside the trophosome of S. annulatus n. sp. was detected via TEM. These symbionts have no internal stacked membranes (Figure 9), indicating that they are not methane-oxidizing bacteria (Cavanaugh et al., 1987). Stable isotope analyses of the tubeworm fragments detected less negative δ 13 C values of S. annulatus n. sp. (i.e., −38.57 and −33.70 ) than that of S. contortum collected from seep sediment of the Haakon Mosby Mud Volcano (i.e., −47.9 ) (Lösekann et al., 2008), but within the range of δ 13 C values (i.e., −40 to −30 ) of other deep-sea organisms, including some vestimentiferans, that harbor chemoautotrophic sulfur oxidizers (Nelson and Fisher, 1995;Feng et al., 2018). The symbiont morphological characteristics and the δ 13 C signature thus together illustrate that the symbionts of S. annulatus n. sp. are sulfur-oxidizing bacteria, which may mainly utilize a carbon source less depleted in 13 C, such as inorganic carbon (Lösekann et al., 2008). The low δ 34 S values of S. annulatus n. sp. (i.e., −3.78 and −2.40 ) also imply that the assimilation of hydrogen sulfide by its symbionts might be primarily derived from the sulfate-dependent anaerobic oxidation of methane (Feng et al., 2015). Besides, the δ 15 N values of S. annulatus n. sp. (i.e., 2.28 and 3.16 ) are lower than the sediment δ 15 N values (i.e., 5 -6 ) in the South China Sea, indicating the assimilation of some locally occurring and isotopically light nitrate or ammonium by the symbionts of S. annulatus n. sp. (Feng et al., 2015).
Application of the microbial 16S rRNA amplicon sequencing resulted in more than 300 OTUs from S. annulatus n. sp. (Supplementary Table 5), nonetheless, only the OTU-1 of each individual is present in high abundance (Supplementary Tables 6-8), suggesting that OTU-1 represents the symbionts, or at least the most dominant symbionts, of S. annulatus n. sp. Co-occurrence of multiple phylotypes of symbionts has been well known for macrobenthos inhabiting chemosynthesis-based ecosystems, including those belonging to different siboglinid lineages (e.g., Kubota et al., 2007;Dubilier et al., 2008;Verna et al., 2010;Zimmermann et al., 2014). In comparison, previous studies of S. contortum revealed only one phylotype of sulfuroxidizing symbionts inside its trophosome (Lösekann et al., 2008;Eichinger et al., 2014). Although we were unable to separate the trophosome of S. annulatus n. sp. from its tube, given the high abundance of OTU-1 and the low abundance of all the other OTUs identified herein, we consider that S. annulatus n. sp. likely only hosts OTU-1 as its symbionts. Moreover, the close phylogenetic relationships between the symbionts of S. annulatus n. sp. (i.e., OTU-1) and S. contortum (Figure 10 and Supplementary Figure 1) support that OTU-1 belongs to sulfur-oxidizing bacteria. The small genetic distances between the symbionts of S. annulatus n. sp. (i.e., OTU-1) and S. contortum (Supplementary Table 9) further indicate that the symbionts of these two Sclerolinum species may be conspecific or have diverged only recently. Nevertheless, whole-genome sequencing of the siboglinid symbionts should be conducted in future studies to better understand their intra-and inter-specific diversity and differences in metabolic capabilities.
In the Haima cold seep, S. annulatus n. sp. forms huge aggregations consisting of hundreds to thousands of individuals. The posterior of the tube is inserted into sediment, which may allow the tubeworm to gain access to the reducing compounds available in porewater (Freytag et al., 2001). Numerous empty shells of the vesicomyid clam A. marissinica along with a few empty shells of the bathymodioline mussel G. haimaensis have been observed surrounding S. annulatus n. sp. (Figure 2). This phenomenon indicates that the S. annulatus n. sp. aggregations may have developed after the decline of the local chemosymbiotrophic bivalve populations. A similar successional pattern from mollusks to tubeworms has been proposed for cold seep communities in the Gulf of Mexico (Cordes et al., 2009). Additionally, S. annulatus n. sp. aggregations develop a complex three-dimensional habitat for various macrobenthos, such as the squat lobster Munidopsis sp. and the brittle star Ophiophthalmus sp. (Figure 2). A resembling association between Sclerolinum and other macrobenthos has been observed for S. contortum inhabiting the Storegga Slide seep area (Vanreusel et al., 2009) and the Loki's Castle vent field (Kongsrud and Rapp, 2012) in the Norwegian Sea. It has been suggested that S. contortum might be capable of releasing sulfate and hydrogen ions into sediment to facilitate sustained sulfide production as other siboglinid tubeworms (Dattagupta et al., 2006). Besides, S. contortum may also have the ability to modify its local environments by releasing iron and manganese from sediment to ambient water (Aquilina et al., 2014). These multiple lines of evidence imply Sclerolinum tubeworms as an ecologically important taxon in worldwide chemosynthesis-based ecosystems. Therefore, further in-depth studies are desired to enhance our knowledge on their species diversity, biogeographical distribution, adaptive evolution, as well as symbiotic relationships with sulfur-oxidizing bacteria in the coming future.

DATA AVAILABILITY STATEMENT
Raw Illumina sequencing data of Sclerolinum annulatus n. sp. have been deposited in the National Centre for Biotechnology Information (NCBI) under the BioProject PRJNA771719. Sequences of the mitogenome (accession numbers OL681896-OL681898) and the nuclear 18S rRNA gene (accession numbers OL681870-OL681872) of S. annulatus n. sp. T1-T3 have been deposited in the GenBank. Sequences of the microbial 16S rRNA gene of all the identified OTUs from S. annulatus n. sp. T4-T6 have been summarized in Supplementary Tables 6-8.

AUTHOR CONTRIBUTIONS
J-WQ, P-YQ, and TX conceived and designed this project. J-WQ and TX collected the samples. YS, ZW, TX, and AS performed morphological examination and description. TX conducted molecular analyses and led manuscript writing. All authors edited the manuscript and approved the submission.