Parabreviscolex niepini were collected from the schizothoracine fish Schizopygopsis younghusbandi Regan, 1905 in the Yarlung Tsangpo River at Linzhi (29°39'N, 94°21'E), Tibet, China, and the specimens were fixed in 100% ethanol and stored at Freshwater Fisheries Research Center, Chinese Academy of Fishery Sciences. Total genomic DNA was extracted using a TIANamp Micro DNA Kit (Tiangen Biotech, Beijing, China), according to the manufacturer’s instructions. DNA was stored at -20 °C for further molecular analyses.

The whole mitogenome was amplified with primers designed based on closely related tapeworms (Suppl. material 1). PCR reactions were performed in a 20 µL reaction mixture, containing 7.4 µL dd H₂O, 10 µL 2×PCR buffer (Mg²⁺, dNTP plus, Takara, China), 0.6 µL of each primer (10 µM), 0.4 µL ExTaq polymerase (Takara, China), and 1 µL DNA template (200 ng/µL). Amplification was conducted as follows: initial denaturation at 95 °C for 2 min, followed by 40 cycles at 95 °C for 10 sec, 46 °C – 53 °C for 30 sec (annealing temperature depending on the primers used, see Suppl. material 1), and 68 °C for 90 sec, and final extension at 68 °C for 10 min. PCR products were sequenced bidirectionally at Sangon Biotech (Shanghai, China) using the primer walking strategy.

The amplified fragments were quality-proofed, and BLASTN (Altschul et al. 1990) to confirm the fragments were the actual target sequence. The complete mitochondrial genomic sequence of P. niepini was assembled manually in a stepwise manner using the DNAstar v7.1 program (Burland 2000). To determine the gene boundaries, it was aligned against the reference mitogenomic sequences of Atractolytocestus huronensis Anthony, 1958 (KY486754) using the program MAFFT 7.149 (Katoh and Standley 2013) integrated with Geneious (Kearse et al. 2012). The mitogenome was annotated and characterized mainly following previous descriptions (Zhang et al. 2017a, b; Zou et al. 2017; Li et al. 2018). Protein-coding genes (PCGs) were found by searching for ORFs (employing genetic code 9, an echinoderm mitochondrial genome) and checking the nucleotide alignments against the reference genome in Geneious. All tRNAs were identified, and confirmed with ARWEN (Laslett and Canback 2008) and MITOS (Bernt et al. 2013) web servers. Similarly, rrnL and rrnS were initially found using MITOS and their boundaries were determined by the alignments with the reference genome in Geneious. The NCBI submission file and tables with statistics for mitogenomes were generated using a GUI-based program, MitoTool (Zhang 2016b). Tandem Repeats Finder (Benson 1999) was employed to find tandem repeats in the non-coding regions. A nucleotide composition table was then used to make the broken line graph of A+T content in ggplot2 (Hadley 2009). Codon usage and relative synonymous codon usage (RSCU) for twelve protein-encoding genes (PCGs) of P. niepini was computed and sorted using MitoTool, and finally imported to ggplot2 to draw the RSCU figure. Ggplot2 was used to draw scatter diagrams for the principal component analysis (PCGs) and nucleotide skews. Input files for the PCGs of codon usage pattern, as well as analyses of amino acid usage pattern and nucleotide skews, were generated by MitoTool. PASW 18.0 (Allen and Bennett 2010) was used to conduct a principal component analysis and generate data for the scatter diagram.

Phylogenetic analyses were undertaken using nucleotide sequences of all 36 genes of the newly sequenced mitogenome of P. niepini and 36 selected cestodes mitogenomes available in the GenBank (Suppl. material 2). The mitogenomic sequences of Khawia sinensis (NC_034800/KR676560) and Caryophyllaeus brachycollis Janiszewska, 1953 (NC_035430/KT028770) from the common carp, sequenced and deposited in GenBank by the same researchers (Feng et al. 2017), were reassigned herein as Khawia sp. 1 and Khawia sp. 2, respectively, because the species identifications were questionable. Caryophyllaeus brachycollis mainly infests the cyprinid Barbus and Abramis in European countries, while its occurrence in China is rare (Barčák et al. 2014). We considered that the researchers have misidentified the two common tapeworms Khawia sinensis and Khawia japonensis (Yamaguiti, 1934) from the common carp.

Two trematode species, Dicrocoelium dendriticum (Rudolphi, 1819) (NC_025280) and Dicrocoelium chinensis Tang & Tang, 1978 (NC_025279), were used as outgroups. The nucleotide sequences for all 12 PCGs, two rRNAs and 22 tRNAs were extracted from GenBank files. The PCGs were translated into amino acid sequences (employing genetic code 9) using MitoTool, and aligned in batches with MAFFT integrated into another GUI-based program BioSuite (Zhang 2016a) using codon-alignment mode. RNAs were aligned with structural alignment mode using the Q-INS-i algorithm incorporated into MAFFT-with-extensions software. BioSuite was then used to concatenate these alignments and remove ambiguously aligned fragments from the concatenated alignments by another plug-in program, Gblocks 0.91b (Talavera and Castresana 2007). Phylogenetic analyses were conducted using maximum likelihood (ML) and Bayesian inference (BI) methods. Selection of the most appropriate evolutionary model for the dataset was carried out using ModelFinder (Kalyaanamoorthy et al. 2017). Based on the Bayesian information criterion, GTR+I+G was chosen as the optimal model for both analyses. ML analysis was performed in RaxML GUI (Silvestro and Michalak 2011) using a ML + rapid bootstrap (BP) algorithm with 1000 replicates. BI analysis was performed in MrBayes 3.2.6 (Ronquist et al. 2012) with default settings, and 6×10⁶ metropolis-coupled MCMC generations.

To determine lineage-specific positively selected sites in individual mitochondrial PCGs, a branch-site model incorporated by CodeML within PAML package (Yang 2007) was used. The resultant ML and/or BI tree (unrooted tree with outgroups removed) was employed for the analysis. The alternative model, MA fixes ω at 1 for each branch except for the specified branch leading to P. niepini (foreground branch), wherein ω is presumed to be greater than 1. The first null model MAnull fixes ω at 1 for every branch in the tree, whereas the second null model M1a fixes ω at 1 for every branch except for the foreground branch, where ω is assumed to be in the range 0 to 1. The null model and alternative model were compared via a likelihood ratio test (LRT), and positive selection was confirmed when P<0.05. Comparing MA to MAnull can estimate positive selection, while comparing MA to M1a can identify instances of relaxation of selective constraints as well as positive selection (Láruson 2017). The posterior probabilities value (≥ 95%) of Bayes Empirical Bayes (BEB) method was used to identify for positively selected sites (Yang 2005).

The closed-circular mitochondrial genome of Parabreviscolex niepini is 15,034 bp in size (GenBank accession number: MG674140). The mitogenome is composed of 12 protein-encoding genes (PCGs), 22 tRNA genes, two rRNA genes, two non-coding regions, and it lacks the atp8 gene (NCR) (Fig. 1). As is common in flatworms, all genes are transcribed from the same strand (Le et al. 2002). Eight overlapping regions and 16 intergenic regions were found in the genome (Table 1). In accordance with other caryophyllidean species, the A+T content of the whole genome (59.6%) and its elements are lower than in the segmented cestodes (Fig. 2d). The mitogenome of P. niepini exhibits G-skew and T-skew, which is also the case in other cestodes (Fig. 2a). However, the unsegmented cestodes appear to exhibit less mutation bias than segmented cestodes (lower GC-skew and higher AT-skew values, Fig. 2a).

Figure 1. 

Circular representation of the mitochondrial genome of Parabreviscolex niepini. Different colors were used to indicated protein-coding genes (12) (red), tRNAs (22) (yellow), rRNAs (2) (green), and non-coding regions (grey). Tapeworm was stained with iron acid carmine.

Figure 2. 

a The comparison of nucleotide skewness of the full genomes for the mitogenome of Parabreviscolex niepini and other cestodes b, c Principal component (PC) analysis of the codon usage and amino acid usage in the PCGs of P. niepini and other cestodes. The first PC (PC1) and the second PC (PC2) of the codon usage and amino acid usage accounted for 96.7% and 98.08% of the variability, respectively. d G+T content of complete genomes and their individual elements. The six caryophyllideans are represented by triangles in a-c. Abbreviations: AH: Atractolytocestus huronensis; BO: Breviscolex orientalis; Ksp2: Khawia sp. 2; KSK: Khawia sinensis; Ksp1: Khawia sp. 1; PN: Parabreviscolex niepini.

Coalesced PCGs were 9999 bp in size, the lowest A+T content (56.6%) in all selected eucestodes (Suppl. material 2), which is also reflected in individual PCGs from 54% (nad3) to 63.2% (nad4L) (Suppl. material 3). ATG is the most commonly used initial codon for eight PCGs; exceptions are nad1, nad3, nad6, and nad5, which use GTG. Among the terminal codons, nine out of 12 are TAG, while nad5 uses TAA, cox3, and cox1 uses abbreviated stop codons (T--) (Table 1).

Codon usage, RSCU, and codon family proportion (corresponding to the amino acid usage) of P. niepini was investigated (Suppl. material S5). The four most abundant codon families (Phe, Val, Leu2, and Gly) encompass 38.41% of all codon families. Among these codon families, G+T-rich codons are favored over synonymous codons with lower G+T content in P. niepini (Suppl. material S5). This G+T preference corresponds well with the relatively high G+T content (Suppl. material S3) as well as G and T preference in the skewness analysis for PCGs (Suppl. material S2). Additionally, the principal component analyses (PCGs) suggested that the overall amino acid usage patterns of the unsegmented cestodes (except for Khawia sinensisKY486753) were apparently different from segmented cestodes (Fig. 2c). Noteworthy, in contrast to segmented cestodes, which have notably heightened A+T content at the 3^rd codon position, these unsegmented cestodes (except Khawia sinensisKY486753) exhibit lower and/or similar A+T content to other elements of the mitogenome (Fig. 2d).

The two rRNAs, rrnL, and rrnS are 953 and 707 bp in size, with 59.6% and 60.4% A+T content, respectively (Suppl. material S3). All 22 commonly found tRNAs are present in the mitochondrial genome of P. niepini, ranging from 58 bp (trnQ and trnR) to 66 bp in size (trnN, trnW, trnT and trnL1), and adding up to 1381 bp in total coalesced length (Table 1 and Suppl. material S2). All of the secondary structures (predicted by MITOS and ARWEN) exhibit the conventional cloverleaf structure, except for trnS1^(AGN) and trnR, which lack DHU arms. The unorthodox trnS1^(AGN) and trnR were also found in the Caryophyllidea (Li et al. 2017) and the Anoplocephalidae (Guo 2016). Additionally, the anti-codon of trnS1^(AGN) in the mitogenome of P. niepini is TCT, in contrast to other eucestodes, which use GCT, except for Khawia sinensis (KY486753) (Suppl. material S4).

The two major non-coding regions (NCR), 567 bp (NCR1) and 1428 bp (NCR2) in size, are located between trnL2 and cox2 and between trnG and cox3, respectively. The positions of the two NCR are consistently reported in other unsegmented tapeworms (see fig. s3 of Li et al. 2017). They have apparently higher A+T content (76.5% for NCR1 and 72.6% for NCR2) than other parts of the genome (Suppl. material S3). The NCR1 contain five tandem repeats (TRs), with two truncated TRs (repeat unit 1 and 5) and one T insertion in repeat unit 2 (Fig. 3). Two highly repetitive regions (HRR) were found in NCR2. HRR1 possess seven TRs, identical in size (40 bp). Repeat units 1–3 are identical in nucleotide composition. In comparison to the repeat units 1–3, repeat unit 4 differed in three nucleotides, while unit 5 and 6 differed in five nucleotides (Fig. 3). HRR2 possess 20 TRs, repeat units 1–19 are identical in size (57 bp) and nucleotide composition, whereas unit 20 is 46 bp long (Fig. 3).

Figure 3. 

Tandem repeats in two main non-coding regions of Parabreviscolex niepini.

The phylogenetic topology constructed using BI and ML methods show concordant branches and high statistical support. All bootstrap support values (BS) are higher than 68 and Bayesian posterior probabilities (BPP) are higher than 0.96. Parabreviscolex niepini exhibits the closest phylogenetic relationship with Breviscolex orientalis Kulakovskaya, 1962 and Atractolytocestus huronensis with robust support (Fig. 4). Moreover, the similarity of the codon usage pattern (Fig. 2b) lends further support to the phylogenetic affinity of P. niepini, B. orientalis and A. huronensis. The mitochondrial gene arrangement of P. niepini (Fig. 1) shows a consistent pattern in the Caryophyllidea, and obviously differs from the segment tapeworms (fig. S3 of Li et al. 2017).

Figure 4. 

Phylogenetic tree of five cestode orders inferred from maximum likelihood analysis with concatenated nucleotide sequence of all 36 genes (12 PCGs, 2 rRNAs, and 22 tRNAs). Bootstrap (BS)/bayesian posterior probability (BPP) support values are shown above the nodes, only BS < 100 and BPP < 1 are displayed. Scale bar represents the estimated number of substitutions per site.

The branch-site model tests based on the criteria of posterior probabilities ≥ 95% in the BEB analyses and in the likelihood ratio test (LRT) (P<0.05), found the amino acid positions V(6) and H(49) of P. niepinicytb (Suppl. material S6, Table 2) were under positive selection. Moreover, several sites in nad4, nad5, and cox3 were also identified to exhibit relaxed selective pressure (Table 2).