The specimens of P. wilsoni were sampled from Dalian, Liaoning Province, China (38.8732°N, 121.6767°E) and identified by Deyuan Yang. All specimens were fixed directly in 95% ethanol. Two specimens with a fully-everted pharynx were used for morphological and molecular studies. They were deposited at Xiamen University (XMU) under voucher numbers 23007-1 and 23007-2, respectively. Specimens were identified based on the key in Glasby and Hsieh (2006). Methods used in the morphological study were followed by Yang et al. (2022). In summary, whole worms were photographed with an Olympus E-M1 Mark II camera with a 60 mm macro lens, detailed structures with a Zeiss Discovery V20 stereomicroscope, and a Nikon 80i compound microscope. Image stacks were obtained using Helicon Focus 7 (https://www.heliconsoft.com/heliconsoft-products/helicon-focus/) and post-processed using Adobe Photoshop. The terminology of nereidid followed Wilson et al. (2023) generally. To better describe areas II, III, and VI of paragnaths, we introduced the notations ‘L:R’ and ‘L:M:R’. ‘L:R’ is for areas II and VI, and ‘L:M:R’ to describe the paragnath patch in area III, where ‘L’ represents the number of paragnaths on the left, ‘M’ in the middle, and ‘R’ on the right. For example, for area II, ‘2:3’ signifies 2 paragnaths on the left and 3 on the right; for area III, ‘3:5,2R:7’ signifies 3 paragnaths on the left, a central patch with 5 paragnaths in 2 rows, and 7 paragnaths on the right.

Before DNA extraction, each individual was cleaned with 95% ethanol. To avoid gut contamination, two to seven parapodia were clipped from the specimens. Whole genomic DNA was extracted using the TIANamp Genomic DNA Kit (TIANGEN, Beijing, China). Genome skimming was conducted on the Illumina NovaSeq X Plus with a PE150 strategy. DNA extraction and sequencing were performed at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). Our initial sequencing effort was to obtain 5 Gb of raw data for each sample. Voucher number 23007-2 was increased to 25 Gb to explore more universal single-copy ortholog genes.

Raw paired-end reads were removed from sequence adapters, and low-quality regions were trimmed using Fastp v.0.23.4 (Chen et al. 2018). The quality of all cleaned reads was checked with FastQC v.0.12.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and subsequently, FastQC results were summarized using MultiQC v.1.8 (Ewels et al. 2016) (Suppl. material 1).

To ensure reliable assembly for the present data, various assemblers were used: 1) SPAdes v.3.15.5 (Bankevich et al. 2012), Ray v.2.3.1 (Boisvert et al. 2010), Megahit v.1.2.9 (Li et al. 2015), and IDBA-UD v.1.1.3 (Peng et al. 2012); 2) GetOrganelle v.1.7.6.1 (Jin et al. 2020); 3) MITObim v.1.9 (Hahn et al. 2013) and NovoPlasty v.4.3.1 (Dierckxsens et al. 2017), with 16S (GenBank accession numbers: LC482188) as the seed; 4) our new assembly pipeline, FastMitoAssembler (original version at https://github.com/suqingdong/FastMitoAssembler). The pipeline implemented in FastMitoAssembler involved the following steps: First, MEANGS v.1.2.1 (Song et al. 2022) was used to obtain and assemble mitochondrial genes, which then served as a seed input for NovoPlasty. Second, the results of NovoPlasty served as a seed input for GetOrganelle. Third, MitoZ v3.6 (Meng et al. 2019) was used for annotation. This workflow was managed using Snakemake (Köster and Rahmann 2012), which can handle large-scale samples with minimal input–only raw data.

For de novo assembly results (SPAdes, Megahit, Ray, and IDBA-UD), we used MitoFinder v.1.4 (Allio et al. 2020) to extract mitochondrial contigs (mt contigs) and annotations, using the Platynereis dumerilii mitogenome (GenBank: NC 000931) as a reference. Unfortunately, no mt contigs were found from the Ray assembler. If the circularization of a mt contig was not automatically completed by MitoFinder, the circules.py script (https://github.com/chrishah/MITObim/tree/master/misc_scripts) was used to assess the circularization of the contig. After obtaining mitogenomes from various assemblers, Quast v.5.2.0 (Gurevich et al. 2013) was employed for assembly quality evaluation (Suppl. material 2). A BLASTn query against GenBank in the NCBI (NCBI BLAST Analysis, https://BLAST.ncbi.nlm.nih.gov/BLAST.cgi) was also performed to check for contamination in the mitogenomes (Suppl. material 3: table S1). The GC content and read (per-base) coverage depth of mitogenome sequences were calculated using the visualize subcommand in MitoZ v.3.6 with default parameters.

The mitogenome was annotated using MITOS2 (Donath et al. 2019), MitoFinder v.1.4.1, and MitoZ v.3.6. All annotation files (in .gb or .gbf formats) were reordered with COX1 as the start gene using PhyloSuite v.1.2.3 (Zhang et al. 2020). The reordered GenBank format (gb) file was converted to a fasta file and uploaded to GeSeq (https://chlorobox.mpimp-golm.mpg.de/geseq.html) to combine the results of ARWEN v.1.2.3 and tRNAscan-SE v.2.0.7 into a single gb file. This approach allowed us to compare tRNA annotations generated from three common methods (ARWEN, tRNAscan-SE, and MiTFi). All annotations were then loaded into Geneious Prime v.2022.2.2 for manual curation. First, we utilized the MAFFT plugin in Geneious to align these gb files. Then, we checked the coding region annotations with the assistance of the ORFs function based on the invertebrate mitochondrial genetic code table 5. Second, tRNAs were evaluated using the following criteria: tRNAscan-SE was given precedence over ARWEN and MiTFi, while MiTFi had precedence over ARWEN. Thus, if tRNAscan-SE identified a tRNA, its result was chosen. In cases where both MiTFi and ARWEN identified a tRNA, the result from MiTFi was selected. Subsequently, we imported other gb files of Perinereis into Geneious to compare each gene and manually edit them.

The boundaries of each gene were determined by the following rules: 1) Protein-coding genes can overlap with each other but cannot overlap with tRNA genes. In Nereididae, only the ND4 and ND4L coding genes are allowed to overlap. 2) The boundaries of mitochondrial rRNA genes (12S and 16S) were defined by flanking genes. 3) The large non-coding region (or putative control region) was determined by the boundaries of neighboring genes and a low GC region.

The CGView Comparison Tool (CCT) (https://github.com/paulstothard/cgview_comparison_tool) (Grant et al. 2012) was employed to compare all accessible mitogenomes of Nereididae in NCBI (as of 15 January 2024) to evaluate the quality of annotations further.

The high-copy nuclear markers (18S to 28S genes and histone genes) were assembled using GetOrganelle v.1.7.6.1. Considering the lack of effective software or pipelines to annotate these genes, we manually annotated nuclear rRNA genes in Geneious using primer pairs to define gene boundaries: 18SA and 18SB (Medlin et al. 1988) for 18S, F63.2 and R3264.2 (Struck et al. 2006) for 28S, and ITS18SFPOLY and ITS28SRPOLY (Nygren et al. 2009) for the ITS1-5.8S-ITS2 region. After completing one sequence annotation, we utilized the “Transfer annotations” function within Geneious Prime v.2022.2.2 to annotate other sequences. The histone genes were annotated using the same function, with GenBank accession number X58895 (Urechis caupo) as a reference at an 85% similarity threshold. If the number of gene annotations was incomplete in one sequence, we imported the corresponding cleaned sequencing reads into Geneious and reassembled the sequence employing the “map to reference” function in Geneious Prime v.2022.2.2 (mapper: Geneious, fine-tuning: iterate up to 25 times) and then re-annotated it. Once annotations were completed, each gene was extracted for NCBI BLAST analysis. The read (per-base) coverage depth of nuclear genes was calculated using a custom script, which used BWA (Li 2013) and SAMtools (Li et al. 2009).

BUSCO v.5.6 was employed to generate universal single-copy orthologs (USCOs) (Simão et al. 2015) using only the assembly results (contigs) from SPAdes with the following parameters: -l metazoan_odb10, -m genome, --augustus, -c 40. A total of 618 single-copy genes were searched.

The phylogenetic analysis utilized two types of datasets: complete or near-complete mitogenomes of Nereididae and DNA barcoding sequences (COX1, 16S, and ITS) of Perinereis. All publicly available Nereididae mitogenomes from GenBank (as of 15 January 2024), along with two outgroup species (Craseoschema thyasiricola NC 060815 and Leocrates chinensis NC 066969, belonging to Chrysopetalidae and Hesionidae, respectively), were included in the phylogenetic analysis. The outgroup taxa were selected based on previous hypotheses that Chrysopetalidae and Hesionidae are sister groups to Nereididae (Glasby 1993; Pleijel and Dahlgren 1998; Tilic et al. 2022). To ensure consistent annotation criteria and to correct potential errors in previous annotations of mitogenomes, we re-annotated all available Nereididae mitogenomes using the methods described above. Ultimately, Platynereis massiliensis (NC 051996) and Alitta succinea (NC 051993) from the Nereididae were excluded from this study (Suppl. material 3: table S2).

To identify potential taxonomically mislabeled sequences in the mitogenomes dataset, we extracted the COX1 and 16S genes and conducted NCBI BLAST analysis. The results were downloaded as CSV files. For each sequence, the top 10 matching sequences were sorted based on identity using WEKits v.1.0 (https://github.com/GP-sir/wekits/releases). Next, we retrieved the source information (subjectAcc) for these matched sequences from NCBI. The species identification of each sequence was then manually verified in Microsoft Excel. Specifically, we focused on matched sequences with identity over 97% and checked whether they were supported by reliable taxonomic references (Suppl. material 3: table S3, S4).

All manually curated mitogenomes were imported into PhyloSuite v.1.2.3. The extracted protein-coding genes (PCGs) and two rRNAs (2R) from these sequences were aligned using MAFFT in normal mode. MACSE v.2 (Ranwez et al. 2018) was used to improve the multiple sequence alignment of 13 PCGs in the “refinement strategy.” After completing the alignment tasks, each alignment was manually checked in Geneious. Ambiguously aligned regions were removed using Gblocks with the default settings or trimAl v.1.2 (Capella-Gutiérrez et al. 2009) with the automated1 setting (Suppl. material 3: table S5). Five mitogenomic datasets were used in the phylogenetic analyses: (i) 13PCGs: PCGs with all three codon positions; (ii) 13PCGs12: PCGs without 3^rd codon positions; (iii) 13PCGs_2R: PCGs dataset plus two rRNA genes; (iv) 13PCG12_2R: PCG12 dataset plus two rRNA genes; and (v) 13PCGsAA: amino acid sequences translated from 13PCGs. The heterogeneity of each dataset was assessed using AliGROOVE v.1.08 (https://github.com/PatrickKueck/AliGROOVE) with the default settings.

ModelFinder v.2.2.0 was used to select the best substitution models for each partition of the maximum likelihood (ML) and Bayesian inference (BI) analyses (Suppl. material 3: table S6). ML analysis was carried out in IQ-TREE v.2.2.2 under the edge-linked partition model with 20,000 ultrafast bootstraps. BI analysis was performed using MrBayes v.3.2.7a (Ronquist et al. 2012), with two parallel runs for at least 2,000,000 generations, ensuring the average standard deviation of split frequencies was below 0.01 (Ronquist et al. 2012).

Considering the mitogenomes of Perinereis are limited, shorter DNA sequences (DNA barcodes) were also used to explore the phylogenetic position of Perinereis wilsoni. All available Perinereis genes from NCBI (https://www.ncbi.nlm.nih.gov/nuccore/?term=Perinereis) were downloaded. The COX1, 16S, 18S, 28S, H3, and ITS genes were further extracted because these molecular markers are widely employed in Annelida. To organize the data more efficiently, we used a custom script (Suppl. material 4) to categorize these data based on voucher numbers (treating isolates, clones, strains, and haplotypes as synonymous with voucher numbers). Ultimately, the COX1, 16S, and ITS genes were selected for phylogenetic analysis due to their higher representation among taxa, and new sequences from this study were added to explore the phylogenetic position of Perinereis wilsoni (Suppl. material 3: table S7). The trees constructed with the COX1, 16S, and ITS genes used the same methods as the mitogenome data.

Multiple alignment summary statistics were generated using the alignment summary function in BioKIT (Steenwyk et al. 2022), including the number of taxa (sequences), alignment lengths, number of constant sites, number of parsimony informative sites and variable sites, and the frequency of all character states.

Pairwise tree structure comparison was conducted using the all.equal.phylo function in the ape v.5.7.1 package (Paradis and Schliep 2019), and topological differences were assessed using TreeSpace (Jombart et al. 2017), both implemented in R v.4.3.1 (R Core Team 2023). Finally, iTOL v.6 (Letunic and Bork 2021) was used to visualize the trees.

Strand asymmetries were calculated using the following formulae (Perna and Kocher 1995): AT-skew = (A - T) / (A + T); GC-skew = (G - C) / (G + C). Codon usage and relative synonymous codon usage (RSCU) of 13 PCGs were computed in PhyloSuite and visualized using the ‘ggplot2’ package (Wickham 2016) in R v.4.1.3 (R Core Team 2023). DnaSP v.6.0 (Rozas et al. 2017) was used to calculate the non-synonymous (Ka)/synonymous (Ks) substitution rates among the 13 PCGs of Nereididae and the nucleotide diversity (Pi) with a sliding window of 100 bp and a step size of 25 bp. A Ka/Ks ratio < 1 indicates purifying selection, while a ratio > 1 indicates positive selection.

To determine which molecular marker is reliable for accurate species delimitation, all labeled P. wilsoni and P. mictodonta COX1, 16S, and ITS genes in NCBI were used. We used the following criteria: 1) whether the presence of a barcode gap: the minimum interspecific genetic distance is greater than the maximum intraspecific genetic distance; 2) whether each species is recovered as monophyletic; and 3) whether there is a small overlap between intra- and interspecific distances and a large barcode gap, as recommended by dos Santos Vieira et al. (2020). Genetic distances (p-distance and Kimura 2-parameter) over sequence pairs were calculated in MEGA X (Kumar et al. 2018). The heat maps were generated in TBtools v.2.080 (Chen et al. 2020).

The mitogenomes (23007-1 and 23007-2) generated from various assemblers yielded consistent results, except for the Ray assembler. They are 15,817 bp with an average coverage depth of 315×, and read mapping is 0.09% (Suppl. material 2). The genomes contain 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and one putative control region measuring 1160 bp. All genes are distributed on the heavy (H-) strand, similar to other Nereidinae species (Fig. 2A, Suppl. material 3: table S8). The nucleotide identity of 13 PCGs between P. wilsoni and the other 31 Nereididae species showed that the other species of Perinereis and Neanthes acuminata (OQ729916) have a high nucleotide identity with P. wilsoni (Fig. 2D).

Figure 2. 

The gene map of Perinereis wilsoni (23007-1), 23007-2, is the same as 23007-1 but only shows one; A. Mitogenome; B. Nuclear rRNA cluster; C. And histone genes; D. CCT (CGView Comparison Tool) map and sequence identity compare the mitogenome between P. wilsoni and the other nereidids. A–C. The green lines depict the distribution of coverage depth; D. Starting from the outermost ring, the feature rings depict: 1. COG (Clusters of Orthologous Groups of Proteins) functional categories for forward strand coding sequences; 2. Forward strand sequence features; 3. The remaining rings show regions of sequence similarity detected by BLAST comparisons between CDS translations from the reference genome and 31 comparison genomes. BLAST identities are organized from high to low, with higher values closer to the outer ring.

The NCBI BLAST results of the mitogenome revealed that the sequence identity with the published 16S sequences of P. wilsoni (LC482171–LC482183, Tosuji et al. 2019) ranged from 98.34% to 99.52%. However, the NCBI BLAST results for the full-length 16S indicated over 98% sequence identity with both P. mictodonta (e.g., LC482161–LC482168, Tosuji et al. 2019) and P. wilsoni (e.g., LC482171–LC482183, Tosuji et al. 2019). These results show that the 16S marker is not suitable for distinguishing these two species.

The nuclear rRNA contigs for the two specimens (23007-1 and 23007-2) yielded 10,983 bp and 10,820 bp with average coverages of 218.3× and 306×, respectively. Both results contained the full 18S, ITS1, 5.8S, ITS2, and 28S regions, with gene lengths of 1,849 bp, 358 bp, 160 bp, 305 bp, and 3,840 bp, respectively (Fig. 2B). The nuclear histone genes for both specimens were 6,399 bp, with coverage depths of 1,547.5× and 2,620×, respectively, incorporating complete sequences of H3 (411 bp), H2A (375 bp), H2B (372 bp), and H4 (312 bp) genes (Fig. 2C).

The NCBI BLAST analysis of nuclear rRNA and histone genes showed that the 18S sequence had over 99% sequence identity with eight species belonging to five genera, indicating that 18S may not be suitable for species identification: 91.43% for 28S with Alitta virens (OW028578); 98.48%-98.54% for ITS with P. wilsoni (AF332158–AF332162, Chen et al. 2002); 96%-98% for H3 genes with Perinereis sp. A ZW-2022 (OL546356), P. aibuhitensis (MW622055), and P. suluana (JX443591); approximately 92% for H2A and H2B with A. virens (OW028584), P. dumerilii (X53330); and 89.4% for H4 (Suppl. material 3: table S1).

The mitochondrial genes in P. wilsoni exhibit a high A + T content of 64.6% (34.3% T and 30.3% A) and lower levels of C and G at 21.6% and 13.8%, respectively (Suppl. material 3: table S8). AT-skew and GC-skew are negative for the mitogenome (-0.062 and -0.221) and PCGs except for COX2 (-0.084 and -0.250), and AT-skew is positive for tRNAs and rRNAs (0.041 and 0.043) (Suppl. material 3: table S8). The mitogenome is compact, with a total of three gene overlaps ranging in length from 2 to 7 bp. There are 12 gene spacers (259 bp in total) from 1 bp to 78 bp. All 13 PCGs start with the ATG codon and stop with T, TAA, or TAG (Suppl. material 3: table S9). A total of 22 tRNA genes with lengths ranging from 53 to 68 bp were identified in the mitogenome of P. wilsoni.

The most frequently utilized amino acids in the mitogenome of P. wilsoni are Leu (14.83%), Ile (9.94%), Ser (8.60%), and Ala (7.84%). The least common amino acids are Cys (1.01%), Arg (1.78%), Asp (1.78%), and Gln (1.94%) (Fig. 3, Suppl. material 3: table S10). Relative synonymous codon usage (RSCU) values for the 13 PCGs showed that UCU (Ser) and UUA (Leu) are the two most frequent codons, whereas UGG (Trp) and CCG (Pro) have the lowest frequencies (Fig. 3).

Figure 3. 

The relative synonymous codon usage (RSCU) in the mitogenome of Perinereis wilsoni. The numbers on the bars represent amino acid composition.

Summary statistics for multiple alignments of various datasets are available in Suppl. material 3: table S11. The heterogeneity analysis of five mitogenome datasets revealed that datasets excluding the third codon position (13PCGs12, 13PCGs12_2R) exhibited lower heterogeneity compared to those including the third codon position (13PCGs, 13PCGs_2R). The lowest heterogeneity was observed in the 13PCGsAA dataset (see Suppl. material 5).

A total of ten ML and BI trees were inferred from five mitogenome datasets of Nereididae (see Suppl. material 6). These trees exhibit eight different topologies based on pairwise tree comparisons, indicating the significant influence of the chosen mitogenomic dataset and tree-building method on the inferred phylogenetic relationships. In general, ML and BI trees inferred from the same dataset are generally consistent, except for some obvious differences observed in the 13PCGs_2R and 13PCGs12_2R datasets (see Suppl. material 6). In the 13PCGs12 dataset, we could not re-root the ML and BI trees using the outgroup, which was excluded from TreeSpace analysis. The remaining eight trees were clustered into four groups, representing four types of trees (see Fig. 4B–F). The most common topology is depicted in Fig. 4A, D.

Figure 4. 

An analysis of phylogenies from five mitogenome datasets. A. Maximum likelihood (ML) and Bayesian inference (BI) tree of Nereididae based on the dataset 13PCGs. The GenBank accession numbers used are listed after the species names. The scale bar (0.5) corresponds to the estimated number of substitutions per site. Numbers at nodes are statistical support values for ML bootstrap support. Asterisks denote 100% bootstrap support. “-” indicates no support value. Color-coded clades are four subfamilies within Nereididae. The gene order is shown to the right. “?” indicates the deletion of a gene; B. A two-dimensional MDS plot of eight trees (excluding 13PCGs12_ML and 13PCGs12_BI), colored by different clusters. C–F. For each cluster identified in B., a representative tree was selected. Note: Neanthes glandicincta (NC 035893) is an incorrectly identified taxon, which should be the genus Dendronereis. Neanthes acuminata (OQ729916) and Perinereis fayedensis (OQ729919) should be Perinereis suezensis and Perinereis damietta, respectively. Perinereis aibuhitensis (NC 023943) should be Perinereis linea (NC 063944).

Phylogenetic analyses based on five datasets of the available mitochondrial genomes (32 taxa in Nereididae, including 8 taxa in Perinereis) provide support for the monophyly of the genus Perinereis. Perinereis was either sister to the genus Nereis and the species Cheilonereis cyclurus (MF538532) (Fig. 4C, E, F) or sister to the genera Platynereis, Hediste, Alitta, Nectoneanthes, and the species Laeonereis culveri (KU992689) (Fig. 4D). In all trees, P. wilsoni was sister to P. vancaurica and P. nuntia with high nodal support values (BS ≥ 94%, PP = 1) (Fig. 4 and Suppl. material 6). The P. nuntia species group was found not to be monophyletic. Specifically, P. linea (NC 063944), P. aibuhitensis (NC 023943), and P. vancaurica (ON611802) have two or three paragnaths on Area VI of the pharynx, which belong to Perinereis Group 2 (Hutchings et al. 1991); all other Perinereis species (except Perinereis cultrifera) in the trees belong to the P. nuntia group. Yet, the species of Group 2 were nested within the P. nuntia group.

Phylogenetic trees (ML and BI) based on COX1 genes recovered that P. wilsoni (23007-1, 23007-2) have a closer relationship with P. mictodonta (KC800632, KC800630, KC800628), with nodal support values (BS = 83%, PP = 0.94). All taxa labeled P. wilsoni and P. mictodonta did not each form a monophyletic group, respectively; instead, they were divided into five distinct clades. This suggests the potential of cryptic species present within P. wilsoni and P. mictodonta or that the specimens were sampled from geographically distant localities with some degree of isolation. These five clades were also supported by genetic distances, which have a distinct barcode gap from each other (Fig. 5A, D, and Suppl. material 7). The phylogenetic trees (ML and BI) based on 16S genes suggested that P. wilsoni and P. mictodonta were each not monophyletic (Fig. 5B, E, and Suppl. material 8). Genetic distances also showed that the sequences of P. wilsoni and P. mictodonta had no barcode gap (Fig. 5B). The phylogenetic trees (ML and BI) based on ITS genes supported P. wilsoni as a sister to P. mictodonta (BS = 96%, PP = 1). Although P. wilsoni and P. mictodonta formed a monophyletic cluster with each other, the genetic distance analyses showed that these two species have no distinct barcode gap (Fig. 5C, F, and Suppl. material 9). Based on the dos Santos Vieira et al. (2020) methodology, optimal markers were selected based on an overlap criterion of less than 20%. However, the overlap between P. wilsoni and P. mictodonta exceeded 35%, indicating that COX1, 16S, and ITS genes are not optimal molecular markers (Fig. 5A–C).

Figure 5. 

The heatmap of COX1 A, 16S B., and ITS C. p-distance for Perinereis wilsoni and P. mictodonta, with the barcode gap and distance overlap of them on the left. Bayesian inference (BI) phylogenetic trees for the two species, based on sequences of COX1 D, 16S E, and ITS F, are excerpted from Suppl. materials 7–9. Sequences from P. wilsoni are depicted in black, those from P. mictodonta in green, with sequences from this study highlighted in bold.

The positions of Paraleonnates uschakovi (NC 032361) and Laeonereis culveri (KU992689) are unstable, jumping across different phylogenies (Fig. 4, Suppl. material 6). Laeonereis culveri (KU992689) is closer to P. uschakovi (NC 032361) in four trees (13PCGs_2R_ML, 13PCGs12_ML, 13PCGs12_2R_ML, and 13PCGsAA_BI in Suppl. material 6), all with high nodal support. In contrast, in six other trees, L. culveri (KU992689) is more closely related to the genera Hediste, Alitta, and Nectoneanthes, also with high nodal support, except in the 13PCGs12_ML and 13PCGs12_BI trees (Fig. 4, Suppl. material 6). Paraleonnates uschakovi is closer to the subfamilies Gymnonereidinae, Dendronereinae, and Nereidinae (8 of 10 trees), with low nodal support in ML trees but with high nodal support in BI trees. In 13PCGs12_ML and 13PCGs12_BI trees, P. uschakovi was clustered with the outgroup, all with low nodal support (Suppl. material 6).

There are two primary types of mitochondrial gene order in the known mitochondrial genomes of Nereididae, except for L. culveri (KU992689). The first type of gene order is observed in the subfamilies Gymnonereidinae, Dendronereinae, and P. uschakovi. The second type is found in the subfamily Nereidinae, except for L. culveri (see Fig. 4).

The nucleotide diversity (Pi) analysis was conducted using concatenated alignments of 13 PCGs and two rRNAs of 32 Nereididae species. The sequence variation ratio exhibits variable nucleotide diversity between the Nereididae, with Pi values for the 100 bp windows ranging from 0.061 to 0.563 (Fig. 6A). COX1 (Pi = 0.211), COX3 (0.235), CYTB (0.238), and COX2 (0.239) exhibit a comparatively low sequence variability, whereas ATP8 (0.405), ND2 (0.386), and ND6 (0.378) have a comparatively high sequence variability. This is corroborated by the non-synonymous/synonymous (Ka/Ks) ratio analysis, which shows that COX1 (Ka/Ks = 0.036), COX3 (0.038), CYTB (0.053), and COX2 (0.065) are evolving comparatively slowly, whereas ATP8 (0.270), ND2 (0.196), and ND6 (0.173) are evolving comparatively fast (Fig. 6B). All genes are under purifying selection.

Figure 6. 

Nucleotide diversity analysis: A. Of 13 PCGs + two rRNAs and Ka/Ks rates; B. Of 13 PCGs based on 32 Nereididae species. The Pi values for the 13 PCGs + two rRNAs are shown on the graph. The red line represents the value of nucleotide diversity (Pi) (window size = 100 bp, step size = 20 bp). The pink, purple, and green columns represent the values of Ka, Ks, and Ka/Ks, respectively.

Before discussing this topic, it is crucial to specify the species concept employed in this study. Despite approximately 30 species concepts having been proposed (Hong 2020), the scientific community has not reached a consensus on defining a species. Our study adopts the Gen-morph species concept, initially proposed by Deyuan Hong (Hong 2020), which requires at least two each of independent morphological characteristics and genetic markers for species definition. According to this concept, morphological characters include quantitative morphological features that have been shown to be statistically different between species.

Perinereis wilsoni and P. mictodonta were initially established as separate species based on statistically validated morphometric differences and the results of the ITS genes (Glasby and Hsieh 2006). Subsequent studies by Park and Kim (2007) and Tosuji et al. (2019) demonstrated that COX1, 16S, and ITS genes can distinguish these two species, respectively. However, the COX1 sequences used in Park and Kim (2007) have not been made publicly available. Moreover, their study relied solely on morphological characteristics for species identification and did not validate these identifications with ITS gene analysis, casting doubt on the reliability of their findings. Our re-analysis of sequences from (Tosuji et al. 2019) showed that the 16S was insufficient to distinguish the two species. To date, only ITS genes have proven effective in distinguishing P. wilsoni and P. mictodonta (Chen et al. 2002; Glasby and Hsieh 2006; Tosuji et al. 2019). Considering the Gen-morph species concept, P. wilsoni may not be a “good” species until a second genetic marker (and a more optimal one than ITS; see below) can be found. Therefore, comprehensive specimen sampling from various geographical locations is still required to assess its taxonomic status thoroughly (Deyuan Yang et al. in preparation).

Partial mitochondrial genes, like COX1 and 16S, and nuclear genes, like 18S, 28S, ITS, and H3, are widely used in Nereididae for species discovery and phylogenetic studies (Glasby et al. 2013; Elgetany et al. 2022; Teixeira et al. 2022a; Teixeira et al. 2022b; Teixeira et al. 2024). Among these genes, the 18S and 28S genes are typically used in higher-taxon phylogenetic analyses. A comprehensive assessment of the effectiveness of these molecular markers in differentiating species, both closely and distantly related across various genera and families, remains limited (Halanych and Janosik 2006).

Given that species within the P. nuntia complex are morphologically similar and not easily distinguishable based on morphology alone, molecular-based identification can offer a faster and more reliable method for species identification when reliable molecular references are available. Our analyses were unable to confirm the discriminative capability of the COX1 gene between P. wilsoni and P. mictodonta, as the sequences sourced from public databases were not linked to morphological and ITS gene data, thus limiting our further study. Additionally, the 16S gene was proven to be insufficient to distinguish these two species. Although ITS genes have been found effective in distinguishing cryptic species (Chen et al. 2002; Pleijel et al. 2009; Nygren and Pleijel 2011; Nygren 2014), they may not be the optimal molecular marker in distinguishing P. wilsoni and P. mictodonta due to the lack of a distinctive barcoding gap between these two species. Furthermore, there are only a few available sequences of Polychaeta ITS genes in the public database, possibly due to difficulties in amplifying these sequences by PCR. We found that the primer sets ITS18SFPOLY and ITS28SRPOLY (Nygren et al. 2009) may be more effective than the one provided by Chen et al. (2002) for amplifying ITS genes in Perinereis.

With advances in sequencing technologies, high-throughput sequencing (HTS) is more efficient than PCR-based Sanger sequencing in obtaining molecular markers. Recently, some new molecular markers have been proposed, such as nearly universal single-copy nuclear protein-coding genes (Eberle et al. 2020) and organelle genomes (Margaryan et al. 2021). In this paper, we found the genome-skimming approach makes it easier to explore the mitogenome, complete nuclear ribosomal DNA (18S-ITS1-5.8S-ITS2-28S), and complete histone genes. However, it is more difficult to explore universal single-copy nuclear protein-coding genes, even with sequencing data increased to 25 Gb.

Currently, a robust phylogenetic backbone of the genus Perinereis, based on phylogenomic methods and extensive taxon sampling encompassing major species in all five informal grouping schemes proposed by Hutchings et al. (1991), is still lacking. Traditionally, the genus Perinereis has been considered polyphyletic (Bakken and Wilson 2005; Glasby et al. 2013; Alves et al. 2020; Elgetany et al. 2022). However, this view is largely based on morphological evidence or a limited number of loci and/or may include assembly errors in their datasets. In this study, although our curated mitogenome datasets support the monophyly of the genus Perinereis, the dataset includes only 8 species, whereas the genus currently comprises 103 species (Prajapat et al. 2024). Furthermore, morphological diversity within the genus has been unevenly sampled. Hutchings et al. (1991) proposed an informal grouping scheme with five groups. Our mitogenome dataset includes Group 1A, represented by 1 taxon (P. cultrifera); Group 2A, represented by 4 taxa (P. linea, P. aibuhitensis (should be P. linea), P. camiguinoides, and P. vancaurica); and Group 3A, represented by 4 taxa (P. fayedensis, N. acuminata (=P. suezensis), P. wilsoni, and P. nuntia). However, Groups 1B and 3B are not represented. To confidently establish the monophyly of the genus, future analyses should include multiple representatives from each of these informal morphological groupings.

Based on the available mitogenome datasets, Perinereis wilsoni is a sister group to P. vancaurica and P. nuntia, with high nodal support in all phylogenetic trees. In contrast, single-gene phylogenetic trees suggest that P. wilsoni is more closely related to P. mictodonta, with low nodal support. Although including more taxa, single-gene trees do not provide sufficient resolution in phylogeny. In summary, the phylogenetic relationships of P. wilsoni remain poorly understood due to the limited number of Perinereis species or other Nereididae for which genomic data are available.

A deep discussion of the phylogeny of Nereididae is beyond the scope of our current work. Here, we provide a brief discussion. The positions of P. uschakovi and L. culveri (KU992689) were observed to be unstable in the trees, inferred from different mitogenome datasets in this study.

In detail, L. culveri always nested within the subfamily Nereidinae, which was also found in previous studies using different datasets, such as COX1, 16S, and 18S (Wang et al. 2021; Alves et al. 2023; Villalobos-Guerrero et al. 2024), as well as mitogenome datasets (Alves et al. 2020). L. culveri is also a sister group to the Gymnonereidinae, Dendronereinae, and Nereidinae subfamilies, along with P. uschakovi (this study; Villalobos-Guerrero et al. 2024). P. uschakovi was always found as a sister group to other subfamilies (Alves et al. 2020; Wang et al. 2021; Alves et al. 2023; Villalobos-Guerrero et al. 2024; this study). However, it occasionally grouped with the outgroup Craseoschema thyasiricola (NC 060815, see Suppl. material 6: 13PCGs12_ML tree and 13PCGs12_BI tree) or formed a single clade with the Gymnonereidinae and Dendronereinae subfamilies (see Suppl. material 6: 13PCGsAA_ML tree), which appears to be supported by mitochondrial genome gene order. Alves et al. (2020) also found this situation when removing outgroups. Inspired by a morphological phylogenetic tree proposed by Wu et al. (1981), which is based on the characteristics of pharyngeal armature, we propose a hypothesis that Paraleonnates and Laeonereis (KU992689) may represent early branching members of Nereidinae. The lack of sufficient data on such taxa may be causing their uncertain placements within the phylogenetic framework.

In this study, we uncovered potential errors in assembly and annotation within GenBank. However, the absence of corresponding Sequence Read Archive (SRA) data in public databases hampers accurate confirmation of these errors. Consequently, we filtered these sequences based solely on our expertise. Therefore, I (Deyuan Yang) advocate for the uploading of original data (raw data) to public databases, such as NCBI. Even if some authors are unwilling to upload their data, various assembly methods should be employed to ensure the accuracy of the assemblies.

Additionally, we emphasize that carefully curated datasets (verifying the taxonomic identification of the species used in the phylogenetic study), especially those from public databases, are crucial before conducting phylogenetic studies, as taxonomic misidentifications can lead to incorrect conclusions. For example, Alves et al. (2020) concluded that Gymnonereidinae was non-monophyletic, a finding that was impacted by the inclusion of an incorrectly identified taxon, Neanthes glandicincta (NC 035893), which should have been classified under the genus Dendronereis Peters, 1854 (Zhen et al. 2022; this study). Additionally, in this study, if we had not questioned the name Neanthes acuminata (OQ729916), we would have concluded that Perinereis is non-monophyletic. However, N. acuminata should be re-identified as Perinereis suezensis (Elgetany et al. 2022).