Specimens of Chaetopterus qiani Sun and Qiu, 2014, Chaetopterus sp., Phyllochaetopterus sp., Phyllochaetopterus hainanensis Wang and Li, 2017, Mesochaetopterus tingkokensis Zhang et al., 2015, Spiochaetopterus sp. A and Spiochaetopterus sp. B were collected from the intertidal and subtidal of Hong Kong and Hainan, China ( Table 1 , Figure 1 ) and preserved in 95% ethanol for DNA extraction or 4% formalin for morphological observations. Genomic DNA for each species was extracted using the CTAB method. Paired-end sequencing on the Illumina platform was performed at Novegene or Sangon to obtain 3.67 to 8.86 Gb 150 bp data ( Supplementary Table 1 ).

Assembly for all species except P. hainanensis and Spiochaetopterus sp. B was performed using CLC Genomics Workbench v7.03 (CLCbio, Arhus, Denmark). For P. hainanensis and Spiochaetopterus sp. B, SPAdes-3.15.4 (Prjibelski et al., 2020) was applied for contig assembly using multiple kmer settings (kmer = 21, 33, 55, 77, 99, 127), and NOVOPlasty v2.7.0 (Dierckxsens et al., 2017) and GetOrganelle (Jin et al., 2020) were run under default settings to build scaffolds, one of which was found to be the full mtgenome.

The chaetopterid mitochondrial genes or mtgenome or 18S/28S rRNA were detected using BLAST V2.13.0 (Camacho and Madden, 2013). The mitogenomes were annotated on the MITOS web server using genetic code 05 for invertebrates (Bernt et al., 2013c). The boundaries of protein-coding genes (PCGs) and rRNA genes were manually examined and adjusted based on alignment with published mtgenomes of this family. The GC-skew and AT-skew were defined according to Perna and Kocher (1995): AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C).

A three-gene dataset (mitochondrial cox1, nuclear 18S and 28S rRNA genes) contained 50 species, including 41 published in Moore et al. (2017) and Britayev et al. (2017), seven sequenced in this study, and two as outgroups, Eurythoe complanata (Pallas, 1766) and Owenia fusiformis Delle Chiaje, 1844, were used in our phylogenetic analyses ( Supplementary Table 2 ). The sequences of each gene were aligned using the MAFFT version 7 web server (Katoh et al., 2017), and concatenated using PhyloSuite (Zhang et al., 2020). Poorly arranged locations and very dispersive regions were removed using less stringent selection settings of Gblocks Server which include smaller final blocks, gap positions within the final blocks, and less strict flanking positions (Talavera and Castresana, 2007). The best nucleotide evolutionary model for each partition was selected based on the Akaike information criterion (AIC) (Darriba et al., 2012) of the PartitionFinder2 module in PhyloSuite. Phylogenetic analysis was performed using Maximum Likelihood (ML) and Bayesian inference (BI) were respectively executed using the IQ-TREE module in PhyloSuite with 10,000 Ultrafast Bootstrap (UFBoot) replicates with the SH-aLRT test (Minh et al., 2013) and the MrBayes module in PhyloSuite with 1000 Sampling Freq replicates.

Phylogenetic analyses were also conducted using all mitochondrial genes (13 PCGs, 2 rRNAs, and 22 tRNAs) for all mtgenomes, except in Spiochaetopterus sp. A which contained only 3 PCGs and 11 tRNAs. Eurythoe complanata (Pallas, 1766) (accession number KT726962.1) and Owenia fusiformis (Delle Chiaje, 1844) (accession number NC_028712.1) were used as the outgroups. The saturation of genetic sequences was assessed using DAMBE (Xia, 2018). The methods for sequence alignment, concatenations, removal of poorly aligned locations, model selection and phylogenetic tree construction were the same as those mentioned above.

To assess the phylogenetic placement of Chaetopteridae in Annelida, we used all protein-coding amino acid datasets of 240 species in Struck et al. (2023) and the seven here sequenced, with Terebratulina retusa (Linnaeus, 1758) as outgroup ( Supplementary Table 2 ). The methods for model selection and ML phylogenetic tree construction were the same as the above, while BI was executed with two chains for 5,000,000 generations using Mpi-MrBayes v3.2 (Ronquist et al., 2012). In ML, a bootstrap >90 was considered strong clade support, 70–90 as moderate, and < 70 as weak support (Krenz et al., 2005). For BI, posterior probabilities > 0.95 were considered strong support (Jacobsen et al., 2010).

CREx2 was used to assess the mitogenomic rearrangement (Bernt et al., 2007), and the TreeRex (Bernt and Middendorf, 2011) allowing deducing the ancestral gene order of the inner nodes for Chaetopteridae and predicting the evolution process of present species based on the given phylogenetic tree. The parameters were set following the software recommendations: -s (strong consistency method), -w (weak consistency method), and -W (parsimonious weak consistency method).

The pairwise genetic distance analyses of amino acid sequences of each mitochondrial PCG were conducted using MEGA X (Kumar et al., 2018). Base-substitution rates for the 13 PCGs of all chaetopterid mtgenomes were calculated followed by Sun et al. (2021). In short, the sequences for each gene were aligned using default parameters of the Muscle module within MEGA X. The non-synonymous to synonymous rate ratio (Ka/Ks) was then calculated using the YN method implemented in KaKs_Calculator 2.0, with Ka/Ks indicating the strength of the selective pressure as > 1 positive selection, = 1 neutral evolution and < 1, purified selection (Yang and Nielsen, 2000). Spiochaetopterus sp. A was excluded from this analysis due to its incomplete mtgenome.

Sequence assembly resulted in six complete or nearly complete mtgenomes ranging from 15,665 bp to 21,822 bp ( Table 1 ) except for Spiochaetopterus sp. A which resulted in a mitochondrial contig with 3,388 bp with only three potential PCGs (cox3, nad2, nad3) and 11 tRNAs ( Figure 2 ; Supplementary Table 3 ). The start gene protein-coding codon for the species of Chaetopteridae was ATG, ATA and ATT, while the stop codon was either TAA, TAG, or defective TA or T ( Supplementary Table 4 ).

All chaetopterid mtgenomes exhibit a high base bias and are AT-rich (58.2–70.17%), with only P. hainanensis (59.45%) and Phyllochaetopterus sp. (58.12%) having lower values but also close to 60% A+T ( Table 2 ; Supplementary Table 5 ). The 13 PCGs (63.65%) and the third codon (63.45%) were also AT-rich. The full mtgenome and the PCGs showed negative AT-skew (-0.213–0.018, -0.296–0.06) and negative GC-skew (-0.318–0.172, -0.465–0.121) in all chaetopterid species except for the full mtgenome of P. hainanensis whose AT-skew (0.018) was greater than GC-skew (-0.314) ( Supplementary Table 5 ). The tRNAs and rRNAs GC-skew are also negative, while their AT-skew were positive in all species except Phyllochaetopterus sp. ( Supplementary Table 5 ).

Three conserved gene blocks are present in all chaetopterids except Spiochaetopterus sp. B: rrnS-trnV-rrnL-trnL1-nad3-trnS1-nad2, trnI-nad6-trnG-atp6-trnQ-trnW, and cob-trnS2-nad4L-nad4-trnH-nad5, with those in Chaetopterus and Mesochaetopterus being identical except that the location of trnY and trnF ( Figure 2 ). >100 bp non-coding regions between different genes occur in different mtgenome locations but show differences between genera ( Figure 2 ; Supplementary Table 6 ): only 0–1 in Chaetopterus, 2–4 in Mesochaetopterus and Phyllochaetopterus, two in Spiochaetopterus sp. B and five in the incomplete mtgenome of Spiochaetopterus sp. A. In addition, the non-coding regions interrupted several PCGs: cox1 in most Mesochaetopterus and Chaetopterus spp., Phyllochaetopterus sp.; nad6 and nad2 in M. tingkokensis and M. japonicus.

Phylogenetic analysis with 48 ingroups was conducted based on concatenated fragments of cox1, 18S, and 28S with a total length of 3,179 bp, including 649 bp for cox1, 1,653 by for 18S, and 877 bp for 28S ( Supplementary Table 7 ). Chaetopteridae appears as a well-supported clade with ML/BI = 100/1 ( Figure 3 ). Chaetopterus and Mesochaetopterus are both monophyletic sister groups, while Phyllochaetopterus and Spiochaetopterus form a well-supported clade (ML/BI: 100/1) with two moderately-supported sub-clades (ML/BI/83/0.95; ML/BI: 85/#) including species from two genera so that both are paraphyletic.

The phylogenetic analyses were based on a data matrix containing 16,939 characters, including 9,900 bp from 13 PCGs, 2,687 bp from two rRNA genes and 1,630 bp from 22 tRNA genes ( Figure 4 ; Supplementary Table 7 ). Five out of 13 PCGs showed significant saturation but were included because they contained phylogenetic information ( Supplementary Table 8 ). The best models selected for the phylogenetic analyses were shown in Supplementary Table 9 and they produced similar tree topologies to those based on the three-gene dataset ( Figure 4 ). All species of Chaetopteridae form a well-supported clade (ML/BI: 100/1) with two sub-clades, one with three species of Chaetopterus and two of Mesochaetopterus, both well supported (ML/BI: 100/1) consistently with the three-gene dataset results, in which Chaetopterus and Mesochaetopterus are also monophyletic. Another well-supported sub-clade consisted of species from Phyllochaetopterus and Spiochaetopterus, but this showed the three species of Phyllochaetopterus as paraphyletic with species of Spiochaetopterus and Phyllochaetopterus from the same location clustered together.

There were seven gene order patterns in the chaetopterid mtgenomes, including three in Chaetopterus and Mesochaetopterus, three in Phyllochaetopterus and one in Spiochaetopterus ( Figure 5 ). The gene order distance ranges from 1 to 7 ( Supplementary Table 10 ). Those between Chaetopterus and Mesochaetopterus are similar (1 to 2) but they are substantially larger than with Phyllochaetopterus (4 to 7).

The gene order evolutionary patterns were predicted by TreeREx based on the phylogenetic tree ( Figure 5 ; Supplementary Table 11 ). Reversions were detected in the most recent common ancestor of the M. tingkokensis-M. japonicus clade (trnE and trnF), in C. qiani (trnY and trnM), when compared with its most recent common ancestor with Chaetopterus. Two transpositions were detected in the most recent common ancestor of two Phyllochaetopterus sp. (OR637234 and KT726961.1) (trnN and trnK), and P. hainanensis and Spiochaetopterus sp. B (trnD and trnN), compared with its most recent hypothetical common ancestral gene order. One reversion (trnA-trnP and trnR-trnD)was detected in Phyllochaetopterus sp (KT726961.1) and two transpositions (trnC and trnM-rrnS) and one reversion (trnP and trnA) in Spiochaetopterus sp. B, compared with that of its most recent common ancestor.

Pairwise genetic distance analysis was performed using PCG amino acids. Cob, cox1, cox2 and cox3 exhibit low genetic divergence, with the lowest value in cox1 (0.014–0.283) and the highest in atp8 (0.255–0.804), followed by nad2 (0.363–0.712) ( Figure 6 ). Compared to Chaetopterus and Mesochaetopterus, Phyllochaetopterus has a higher genetic divergence in several PCGs, including atp8, cob, cox1 and nad2.

The Ka/Ks of the 13 chaetopterid mitochondrial PCGs are all lower than 0.5, suggesting that they have undergone strong purification selection. Nevertheless, the Ka/Ks of atp8 (0.378) and nad4 (0.423) are higher than those of the other PCGs ( Figure 7A ), with the smallest in the complex IV (i.e. cox1 = 0.018, cox2 = 0.051, cox3 = 0.035) and nad5 (0.070) ( Figure 7A ). Compared with Chaetopterus and Mesochaetopterus, Ka/Ks was higher in Phyllochaetopterus for all mitochondrial PCGs except nad4 ( Figure 7B ), indicating a general relaxation of purifying selection in the mtgenomes of this genus.

The ML/BI phylogenetic tree was constructed using the concatenated datasets of the amino acid sequences of 13 PCGs of mitogenome, with 247 annelid terminals ( Supplementary Table 2B , Supplementary Figure 1 ). Our results show the Oweniidae and Magelonidae as sister groups, forming a sister clade to Chaetopteridae + all other annelid families (ML/BI: 95/#).

The mitochondrial intergenic gaps of metazoans are usually composed of very short non-coding regions or very few overlapping bases (Zhong et al., 2008). Therefore, the length of metazoan mtgenome is generally stable, especially in vertebrates. Conversely, annelid mtgenomes vary substantially from 14kb (Aguado et al., 2016) to 25kb (Seixas et al., 2017; Sun et al., 2021), although chaetopterid mtgenomes show a smaller range from 15kb (Chaetopterus) to 21kb (Phyllochaetopterus). Longer annelid mtgenomes obey either an enormous control region (D-loop) (as in Spirobranchus and Siboglinum) or non-coding region within the cox1 gene (the group II intron) and many huge intergenic gaps (as in Glycera, Hydroides and Polynoidae). Longer mitogenomes generally contain more long non-coding regions and a non-conserved gene order, as in Hydroides spp. (Sun et al., 2021), Siboglinidae (Siboglinum fiordicum Webb, 1963, 19,502 bp; Li et al., 2015). However, there are some exceptions, such as among Glyceridae (Glycera unicornis Lamarck, 1818, 20,366 bp; Richter et al., 2015). In Phyllochaetopterus and Mesochaetopterus, large mitogenomes were either due to the group II introns among the cox1, to the non-coding region within the nad2, or numerous huge intergenic non-coding regions. Long non-coding regions were produced probably during changes in gene structure in Chaetopteridae, since intergenic non-coding spaces have been suggested to facilitate the inception of replication (Pons et al., 2014), or as remnants of gene order change (such as TDRL and transpositions) (Seixas et al., 2017), strands of the template mispairing, or imprecise termination during replication (Bernt et al., 2013b; Aguado et al., 2016).

The mtgenome gene orders of Chaetopteridae greatly differ from those of the ancestral annelid and most other annelids, although it is relatively conserved and shows three conserved gene blocks in most species ( Figure 2 ). Among the ten mtgenomes, those of Chaetopterus and Mesochaetopterus show very similar gene orders, with only one reversion in the position of trnY and trnM or trnE and trnF. The Phyllochaetopterus and Spiochaetopterus mtgenomes are also very similar, with several unstable positions of tRNA (i.e. trnR, trnD, trnP, trnN, trnK, trnC) and a fragment of trnM-rrnS. These results are consistent with the close phylogenetic relationships between Chaetopterus and Mesochaetopterus, and between Spiochaetopterus and Phyllochaetopterus, respectively.

Our study shows the Chaetopteridae as a sister family to most annelid families, being recovered as sisters to the Oweniidae/Magelonidae clade. Transcriptome and mtgenome studies have all recovered Chaetopteridae as one of the basal annelid families, although its relationships with other basal families are unstable (Rousset et al., 2007; Weigert et al., 2016; Struck et al., 2023). Our ML result (i.e. ((Oweniidae + Magelonidae) + (Chaetopteridae + other annelid families)) agrees with the result in Struck et al. (2023), while the BI analyses did not favor these relationships (Struck et al., 2023). Our Chaetopteridae clade was well supported, in agreement with the previous results based on universal biomarker genes (cox1, 18S, 28S) (Martin et al., 2008; Osborn et al., 2007; Morineaux et al., 2010; Weigert et al., 2014; Andrade et al., 2015; Zhang et al., 2015; Britayev et al., 2017; Moore et al., 2017; Martin et al., 2022; Yang et al., 2022), while differing from morphological studies based on larvae and adults, which classified Chaetopteridae as a family within Sabellida (Rouse and Fauchald, 1997; Brown et al., 1999; Rouse, 1999). (Apistobranchidae + (Psammodrilidae + Chaetopteridae)) formed the clade of Chaetopteriformia based on phylogenetic analysis using transcriptomes (Helm et al., 2018), however, only Chaetopteridae was included in the later report on mtgenome evolution in Annelida (Struck et al., 2023). The mtgenomes of Apistobranchidae and Psammodrilidae are needed to further explore the ancestral mtgene order of Chaetopteridae considering their close evolutional relationship with Chaetotperidae.

The monophyly of Chaetopterus and Mesochaetopterus was not supported by morphology and BI based on 28S (Osborn et al., 2007; Morineaux et al., 2010). In this study, Chaetopterus appears as a well-supported sister taxon to Mesochaetopterus, in agreement with Moore et al. (2017); Martin et al. (2022) and Zhang et al. (2015). The paraphyletic relationship of Phyllochaetopterus + Spiochaetopterus is consistent with the results of the three scholars above. Major revisions are needed to define more distinguishing characteristics of these two genera instead of the present/absent of notopodial cirrus of A1.

Among the PCGs, atp8 has the highest pairwise genetic distance ( Figure 6 ), followed by nad2 and nad6, suggesting faster substitution rates compared to other PCGs, as previously shown in deep-sea Polynoidae (Zhang et al., 2018), in which atp8 is the fastest-evolving mitochondrial gene (Sun et al., 2021). Additionally, the PCG genetic distance in Phyllochaetopterus was generally greater than that of the other genera, suggesting different selective pressures and faster mutation rates in its mitogenomes. In addition, the branch length of the Phyllochaetopterus + Spiochaetopterus clade was longer than that of the Chaetopterus + Mesochaetopterus clade ( Figures 3 , 4 ), indicating a faster substitution rate in the former.

Chaetopteridae exhibits distinct mtgenome gene order arrangement patterns compared to other annelids. Gene orders were initially considered as conserved among annelid mtgenomes (Jennings and Halanych, 2005), which is still true for many families (Weigert et al., 2016). However, substantial gene order arrangements have been found among Syllidae (Aguado et al., 2015, 2016), deep-sea Polynoidae (Zhang et al., 2018) and Serpulidae (Sun et al., 2021). Chaetopteridae shows different mtgenomes than the ancestral annelids (excluding Oweniidae and Magelonidae), mainly in the atp6, nad1, nad6, cox3, cox1, and several tRNAs positions ( Supplementary Figure 2 ). However, its gene order is relatively conserved, with three conserved regions in all species except Spiochaetopterus ( Figure 2 ). The uncertain evolutional relationships between Chaetopteridae and all other polychaetes, the gene order of its most recent ancestor cannot be determined. Therefore, it is not clear how the mitochondrial gene order of the common ancestor of Chaetopterus and Mesochaetopterus or Phyllochaetopterus and Spiochaetopterus evolved from that of the chaetopterid common ancestor.

The diversity of gene orders is positively related to the substitution rate (Bernt et al., 2013a, b; Luo et al., 2015). The conserved functions in the respiration of mitochondrial genes explain why they are undergoing purifying selection (Stewart et al., 2008). Our study shows the Ka/Ks for all chaetopterid protein-coding genes as being lower than 1, indicating purifying selection ( Figure 7 ), while being quite fast evolving compared to other polychaete (Zhang et al., 2018) or invertebrates (Ren et al., 2010). Phyllochaetopterus shows higher base-substitute rates for most mitochondrial PCGs than Chaetopterus or Mesochaetopterus, indicating faster mutation rate that is consistent with their longer branch length in phylogenetic trees. Different mutation rates have also been shown among Polynoidae, with the deep-sea genera Branchipolynoe and Branchinotogluma, having high mutation rates that suggest adaptation values (Zhang et al., 2018).

Overall, our results show different mitochondrial gene orders in Chaetopteridae compared to the conserved Annelida pattern, despite most species in this family sharing three singular conserved regions. Chaetopterus and Mesochaetopterus are sister groups forming well-supported monophyletic clades, which, together, are sister to a paraphyletic Phyllochaetopterus clade with Spiochaetopterus as subclade. The mitochondrial PCG genes in all examined chaetopterids have undergone purifying selection, but in most cases, those of Phyllochaetopterus show a higher base-substitution rate than those of the other genera.