Global analysis of cytochrome c oxidase subunit 1 (cox1) gene variation in Dibothriocephalus nihonkaiensis (Cestoda: Diphyllobothriidae).

The cestode Dibothriocephalus nihonkaiensis (syns. Diphyllobothrium nihonkaiense and Diphyllobothrium klebanovskii), the broad fish tapeworm, is a parasitic agent of intestinal infection acquired by consumption of raw or undercooked Pacific salmon, Onchorhynchus spp. Sequencing studies conducted about a decade ago revealed the presence of two major lineages (A and B) in the broad fish tapeworm population within Asian coastal areas. However, in spite of the accumulation of sequence data on GenBank recently, no further genetic analyses of D. nihonkaiensis have been attempted. The present study assessed for the first time the global cox1 variation in D. nihonkaiensis. Novel partial cox1 sequences of 14 isolates of D. nihonkaiensis from 12 patients were generated, and a global genetic analysis was performed using the 14 novel and 79 previously published sequences for isolates from definitive and second intermediate hosts of this species was performed. A total of 48 haplotypes of three haplotype groups (Types A, B and C) were identified, and co-infections with genetically different D. nihonkaiensis were highlighted in humans and Pacific salmon.

Introduction

The broad fish tapeworm, Dibothriocephalus nihonkaiensis (Yamane, Kamo, Bylund & Wikgren, 1986) (Cestoda: Diphyllobothriidae), previously known as Diphyllobothrium nihonkaiense and Diphyllobothrium klebanovskii, is a major causative agent of intestinal parasitic diseases attributable to ingestion of second-stage larvae (plerocercoids) in raw or undercooked Pacific salmon, Onchorhynchus spp. (Waeschenbach et al., 2017; Scholz et al., 2019). Genetic diversity of D. nihonkaiensis was previously investigated in two studies conducted about a decade ago (Arizono et al., 2009a; Suzuki et al., 2010), describing two major genetic lineages (i.e. Types A and B) in phylogenetic analyses of partial cytochrome c oxidase subunit 1 (cox1) gene using about 20 isolates from humans, bears and Pacific salmon. However, despite of the accumulation of sequence data on GenBank after the two previous studies no global genetic analysis of D. nihonkaiensis has been attempted. Theoretically, humans might be infected with multiple heterogenic D. nihonkaiensis by eating Pacific salmon, but this possibility has not been confirmed because all reported cases represent single parasite infections (Yamane and Shiwaku, 2003; Arizono et al., 2009b). In addition, Pacific salmon might also be infected with multiple heterogenic D. nihonkaiensis by ingestion of the first intermediate host (copepod) harbouring procercoid larvae of this cestode, but this possibility has also not been addressed.

In the present study, we applied genetic analyses of 14 new isolates from 12 clinical cases in Japan in comparison with all isolates available in the GenBank database to assess the global genetic variability of D. nihonkaiensis. The study also confirmed co-infections with heterogenic isolates in humans and intermediate fish hosts.

Materials and methods

Isolates and DNA extraction

Proglottids of 14 specimens of D. nihonkaiensis were collected from 12 patients at the medical institutions in Osaka and Hyogo prefectures, Japan, between 2007 and 2020 (see patient data in Supplementary Table S1). Patients #3 and #12 excreted proglottids of two tapeworms each. The strobilae were naturally passed with faeces or were excreted after anthelminthic treatment, and most isolates (n = 13) were stored in 70% ethanol at 4–8 °C or −20 °C. One isolate (DnHs12-2) was fixed in 10% buffered formalin for about two months (July 8th to September 7th, 2020). Some proglottids of this isolate were transferred in distilled water overnight to remove the formalin solution.

DNA from a small piece of the proglottid was extracted and purified using DNeasy Blood and Tissue Kit or QIAamp DNA Mini Kit (Qiagen GmbH, Hilden, Germany). The isolates were identified as D. nihonkaiensis by multiplex PCR (Wicht et al., 2010b) or sequence analysis of the cytochrome c oxidase subunit 1 (cox1) (Yera et al., 2006).

Amplification and sequencing

Sequences of a partial fragment of the cox1 gene were used for genetic analysis because of the numerous available sequences from the definitive and fish intermediate host species, and because of its usefulness for genetic differentiation of D. nihonkaiensis (Arizono et al., 2009a; Suzuki et al., 2010; Autier et al., 2019). A portion of the cox1 gene (approximately 710 bp) was amplified using the primer set: forward (5′-TTG ATC GTA AAT TTG GTT C-3′); reverse (5′-AAA GAA CCT ATT GAA CAA AG-3′) (Arizono et al., 2009a). PCR amplification was performed in a volume of 25 μl using TaKaRa EX Taq Hot Start Version containing 10× PCR buffer, 20 mM MgCl2, 2.5 mM of each dNTP, 5 units/μl of Takara Ex Taq HS DNA polymerase (TaKaRa Shuzo Co. Ltd., Shiga, Japan), 0.5 μM of each primer, and 2.5 μl of DNA sample. After denaturation at 94 °C for 5 min, amplification was carried out by 35 cycles consisting of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and elongation at 72 °C for 1 min, followed by a final extension step at 72 °C for 7 min. Reactions were performed in a thermocycler (GeneAmp PCR System 9700 or 2720; Applied Biosystems, USA). Aliquots of the PCR products were separated by electrophoresis on a 3% agarose gel and were visualized under UV light after staining with ethidium bromide. Then the PCR products were purified using either the QIAquick Gel Extraction Kit or the QIAquick PCR Purification Kit (Qiagen Inc., Germany). DNA sequencing was performed using a BigDye Terminator v3.1 Cycle Sequencing kit with the primer sets used in the PCR on an automated sequencer (ABI3130; Applied Biosystems, USA). Sequence chromatograms from each strand were inspected using Sequencher DNA Sequence Analysis Software Version 4.1 (Gene Codes Corp., USA).

Phylogenetic and haplotype network analyses

A total of 79 cox1 sequences of D. nihonkaiensis (syns. D. nihonkaiense and D. klebanovskii; see Arizono et al., 2009a and Waeschenbach et al., 2017) available on GenBank were included in the present analyses (Table 1). Other available sequences of D. nihonkaiensis deposited on GenBank were excluded from the analyses because of their short length and the presence of ambiguous nucleotide positions or missing information for the hosts. Alignment and phylogenetic analysis were conducted using MEGA7 software (Kumar et al., 2016). The newly generated cox1 sequences were aligned with the 79 D. nihonkaiensis sequences (Table 1) registered to date in GenBank. Six sequences of related species (Dibothriocephalus latus, Dibothriocephalus dendriticus, Dibothriocephalus ditremus, Dibothriocephalus ursi, Diphyllobothrium balaenopterae and Spirometra mansoni) were used as the outgroup and for comparison of cox1 sequence divergence among different species. The final alignment was trimmed to match the shortest sequence and comprised 666 nucleotide positions with no gaps. The best-fitting model for nucleotide substitution (Hasegawa-Kishino-Yano model) was estimated with MEGA7 and applied in the maximum-likelihood analysis. Models with the lowest BIC scores (Bayesian Information Criterion) are considered to describe the substitution pattern the best. A discrete gamma distribution was used to model evolutionary rate differences among sites. The rate variation model allowed for some sites to be evolutionarily invariable. The phylogram reliability was tested with the bootstrap method using 1,000 replications. The following diversity indices for D. nihonkaiensis populations were calculated using DnaSP (ver. 5.10.01): number of haplotypes (h); haplotype diversity (Hd); nucleotide diversity (π); and Tajimaʼs D. The haplotype network of D. nihonkaiensis populations was constructed based on the median-joining method implemented in the Network software ver. 10.2.0.0 (Fluxus Technology Ltd., www.fluxus-engineering.com) using the data set prepared with DnaSP.

Table 1

Data for the isolates of Dibothriocephalus nihonkaiensis analyzed in this study

Host	Locality	Isolate	GenBank ID	Haplotype	Genetic lineage	Reference
Human (n = 54)	Japan: Kyoto	Dn1	AB288371	H_4	A	Arizono et al. (2009a)
		Dn2	AB288372	H_2	A
		Dn3	AB288373	H_4	A
		Dn4	AB374999	H_37	A
		Dn5	AB375000	H_4	A
		Dn6	AB375001	H_38	A
		Dn7	AB375002	H_39	A
		Dn8	AB375003	H_40	B
		Dn9	AB375004	H_41	B
	Japan: Tokyo	Dnh1	AB521692	H_32	A	Suzuki et al. (2010)
		Dnh2	AB521693	H_32	A
		Dnh3	AB521694	H_4	A
		Dnh4	AB521695	H_33	A
	Japan: Nara	2009-40	AB573407	H_4	A	Nishiohuku et al. (unpublished)
		2009-18	AB573405	H_19	A
		2009-19	AB573406	H_22	A
		2009-41	AB573408	H_23	A
		2010-6	AB573409	H_3	A
	Japan	2013-16	AB821272	H_4	A	Yamasaki and Ashida (unpublished)
	Japan: Tochigi	2017-024	LC312466	H_7	B	Yamasaki and Shimada (unpublished)
	China: Shanghai	CHN-005	AB684623	H_4	A	Chen et al. (2014)
		CHN-002	AB684621	H_15	A
		CHN-003	AB684622	H_16	A
	South Korea	PCH	DQ768189	H_45	A	Yera et al. (2008)
		CSS	DQ768190	H_46	A
		KCH	DQ768188	H_44	A
	China: Heilongjiang	Sample 2	LC070677	H_9	B	Cai et al. (2017)
	South Korea	Not named 1	EF420138	H_42	A	Kim et al. (2007)
	Japan: Hokkaido	Not named 2	AB268585	H_18	A	Nakao et al. (2007)
	Japan: Tochigi	Not named 3	AB508838	H_4	A	Yanagida et al. (2010)
	Japan: Tokyo	Not named 4	AB610797	H_4	A	Yamasaki and Nakamura (unpublished)
	Japan	Not named 5	AB781787	H_17	A	Ishida and Yamasaki (unpublished)
	Japan: Nagasaki	Not named 6	AB544064	H_7	B	Ishida and Yamasaki (unpublished)
	Japan: Saitama	Not named 7	AB597273	H_1	A	Ikeda et al. (2012)
		Not named 8	AB597274	H_24	A
	Switzerland	Not named 9	AM412559	H_43	A	Wicht et al. (2007)
		Not named 10	AM412560	H_7	B
	Japan: Hamamatsu	Not named 11	AB636314	H_18	A	Ohta et al. (2011)
	Japan: Asahikawa	Not named 12	AB364645	H_7	B	Yamasaki and Kuramochi (2009)
	Japan: Tokyo	Not named 13	AB015755	H_6	A	Miyadera et al. (2001)
	Japan: Osaka	DnHs1	LC589648	H_4	A	This study
		DnHs2	LC589649	H_1	A
		DnHs3-1a	LC589650	H_2	A
		DnHs3-2	LC589651	H_3	A
		DnHs4	LC589652	H_4	A
		DnHs5	LC589653	H_4	A
		DnHs6	LC589654	H_5	A
		DnHs7	LC589655	H_4	A
		DnHs8	LC589656	H_47	A
		DnHs9	LC589657	H_21	A
		DnHs10	LC589658	H_21	A
		DnHs11	LC589659	H_48	A
	Japan: Hyogo	DnHs12–1	LC589660	H_4	A
		DnHs12–2	LC589661	H_29	B
Brown bear (n = 3)	Russia: Kamchatka	RB1	AB375660	H_7	B	Arizono et al. (2009a)
		RB2	AB375661	H_34	B
		RB3	AB375662	H_35	A
Pink salmon (n = 1)	USA: Alaska	US361b	KY000483	H_8	A	Kuchta et al. (2017)
Chum salmon (n = 19)	Russia: Okhotsk	Ok1	AB375672	H_4	A	Arizono et al. (2009a)
	Russia: Kamchatka	Ok2	AB375673	H_36	A
	Japan: Aomori	Dnk1	AB521674	H_4	A	Suzuki et al. (2010)
	Japan: Hokkaido	Dnk2	AB521675	H_26	A
	Japan: Hokkaido	Dnk3	AB521676	H_27	A
	Japan: Hokkaido	Dnk4	AB521677	H_28	C
	Japan: Hokkaido	Dnk5	AB521678	H_29	B
	Japan: Miyagi	Dnk6	AB521679	H_4	A
	Japan: Miyagi	Dnk7	AB521680	H_30	A
	Japan: Iwate	Dnk8	AB521681	H_30	A
	Japan: Iwate	Dnk9	AB521682	H_30	A
	Japan: Hokkaido	Dnk10	AB521683	H_4	A
	Japan: Hokkaido	Dnk11	AB521684	H_4	A
	Japan: Hokkaido	Dnk12	AB521685	H_4	A
	Japan: Hokkaido	Dn-Ok1	AB548647	H_4	A	Wicht et al. (2010b)
		Dn-Ok2	AB548648	H_25	A
		Dn-Ok3	AB548649	H_7	B
	Japan: Iwate	TD01–1	LC511596	H_20	B	Jin et al. (unpublished)
		TD01–2	LC511597	H_21	A
Cherry salmon (n = 16)	Japan: Hokkaido	Dnm1	AB521686	H_31	A	Suzuki et al. (2010)
	Japan: Iwate	Dnm2	AB521687	H_4	A
	Japan: Iwate	Dnm3	AB521688	H_26	A
	Japan: Iwate	Dnm4	AB521689	H_26	A
	Japan: Hokkaido	Dnm5	AB521690	H_4	A
	Japan: Hokkaido	Dnm6	AB521691	H_4	A
	Japan: Hokkaido	Dn-Om1	AB548650	H_4	A	Wicht et al. (2010b)
	Japan: Niigata	Om7-2b	AB924498	H_4	A	Watanabe et al. (2014)
		Om33	AB924500	H_4	A
		Om38–1	AB924505	H_13	A
		Om38–2	AB924506	H_14	B
		Om38–5	AB924503	H_11	A
		Om38–7	AB924504	H_12	A
		Om38–9	AB924499	H_4	A
		Om38–10	AB924501	H_10	A
		Om38–12	AB924502	H_4	A

Isolates and haplotypes presented in bold originate from the same host individual.

All isolates originate from three cherry salmons (#7, #33, #38) that have returned to the Miomote River (Watanabe et al., 2014).

Results

A partial cox1 fragment (666–711 nt) was sequenced for the 14 clinical isolates examined in this study. The phylogenetic tree inferred from the cox1 dataset is presented in Fig. 1. The 93 sequences of D. nihonkaiensis were well separated from the other species (Fig. 1A), the congeners D. ursi (AB605762), D. dendriticus (AB374223), D. latus (AB504899) and D. ditremus (AB979518), and D. balaenopterae (AB822370) and S. mansoni (AB374543) with genetic distances (Kimura-2-parameter model) of 0.049–0.062, 0.065–0.080, 0.072–0.080, 0.112–0.124, 0.155–0.163 and 0.189–0.199, respectively (Supplementary Table S2), and the values between D. nihonkaiensis and the other Dibothriocephalus species (D. dendriticus, D. latus and D. ditremus) correspond to the values previously reported by Wicht et al. (2010a). The sequences of D. nihonkaiensis formed three clades (Types A, B and C; Fig. 1B). The genetic distances within Type A and Type B were 0.000–0.014 and 0.000–0.006, respectively, and the values between Type A and Type B, Type A and Type C and Type B and Type C were 0.012–0.026, 0.018–0.026 and 0.017–0.021, respectively. These values correspond well to the intraspecific divergence values in D. nihonkaiensis based on partial cox1 sequences (Arizono et al., 2009a).

Fig. 1

A Phylogenetic relationships of the 93 Dibothriocephalus nihonkaiensis isolates from humans, bears, and Pacific salmon with congeners (Dibothriocephalus ursi, D. dendriticus, D. latus and D. ditremus) and related species (Diphyllobothrium balaenopterae and Spirometra mansoni), based on partial cox1 sequences (666 bp). B Expanded phylogenetic tree of the Dibothriocephalus nihonkaiensis isolates from the present (14 in red) and previous (40 in black) clinical cases, bears (3 in blue), pink salmon (1 in purple), chum salmon (19 in green) and cherry salmon (16 in pink). The two isolates from cases #3 and #12 are indicated with red dashed and solid arrows, respectively. Similarly in both clinical isolates, the two isolates (TD01-1, 2) from a chum salmon and the seven isolates (Om38-1, Om38-2, Om38-5, Om38-7, Om38-9, Om38-10 and Om38-12) from a cherry salmon are indicated by green and pink arrows, respectively.

A total of 48 haplotypes were identified among the 93 partial cox1 sequences in the dataset studied (Fig. 2, Table 1). These haplotypes fell into three haplotype groups corresponding to Types A, B and C identified in the phylogenetic analysis (Fig. 2). Type A and Type B both had a central haplotype with a high frequency of appearance. Type A haplogroup was composed of 39 haplotypes among 78 sequences and these haplotypes formed a diffuse network. Type B haplogroup was composed of 8 haplotypes among 14 sequences, and Type C haplogroup was represented by a single haplotype with at least 11 nucleotide differences from Type B populations. A total of 33 haplotypes were identified among 54 sequences from patients. Three sequences for isolates from bears represented 3 haplotypes; 19 sequences for isolates from the chum salmon (Onchorhynchus keta) represented 11 haplotypes; 16 sequences for isolates from the cherry salmon (Onchorhynchus masou) represented 8 haplotypes; and one sequence for isolate from the pink salmon (Onchorhynchus gorbuscha) represented a different type from other haplotypes. Haplotype H_4 was dominant (29 out of 93 sequenced isolates, 31%). Among these haplotypes, several haplotypes originated from the same host (Table 1, isolates and their haplotypes in bold). Namely, 4 haplotypes originated from two patients examined in the present study (H_2 and H_3 from patient #3; and H_4 and H_29 from patient #12), 2 haplotypes (H_20 and H_21) were from one chum salmon (both sequences are published on GenBank only) and 6 haplotypes (H_4, H_10, H_11, H_12, H_13 and H_14) were from one cherry salmon; these sequences were published by Watanabe et al. (2014), but their differences and haplotypes have not been identified.

Fig. 2

Haplotype network for Dibothriocephalus nihonkaiensis constructed based on partial cox1 sequences. A total 93 sequences were analyzed falling into three types (A, B and C). The size of each circle corresponds to the frequency of the relevant haplotype. The minimum size of the node indicates one individual. Punctuations on branches indicate the number of mutated positions. The circles are colour-coded based on host species (green, chum salmon; pink, cherry salmon; purple, pink salmon; blue, brown bears; red, patients examined in the present study; black, patients examined in previous studies).

Table 2 shows the indices of genetic diversity in D. nihonkaiensis populations examined in the present study. Haplotype diversity for all populations was high (0.825–0.898), but the nucleotide diversity was low (0.00226–0.00692). Tajimaʼs D values were negative in both types (A, B) and all populations with statistically significant values for Type A and all populations, suggesting that D. nihonkaiensis population is genetically diverse and rapidly expanding.

Table 2

Population genetic analysis using partial cox1 sequences of Dibothriocephalus nihonkaiensis

Population	n	h	Hd	π	Tajimaʼs D	P-value
A	78	39	0.860	0.00339	−2.23538	0.01
B	14	8	0.825	0.00226	−1.51416	0.10
A, B, C	93	48	0.898	0.00692	−1.84475	0.05

Abbreviations: n, number of samples; h, number of haplotypes; Hd, haplotype diversity; π, nucleotide diversity.

Discussion

The global analysis of the genetic variation of cox1 gene revealed that the isolates identified as Type A and Type B in the previous studies by Arizono et al. (2009a) and Suzuki et al. (2010) exhibit similar clustering in the present analysis.

However, the genetic divergence between the isolate Dnk4 and Type A, or between the isolate Dnk4 and Type B was similar to that between Type A and Type B (Supplementary Table S2), and the isolate Dnk4 formed a separate clade from Type B with statistical support (Fig. 1B). Moreover, the isolate Dnk4 was clearly different from Type A and Type B haplogroups in the haplotype network (Fig. 2). Therefore, in the present study, the isolate Dnk4 originating from a chum salmon landed at Hokkaido was identified as a new Type C (or haplotype group C). To the best of our knowledge, no haplotype network analysis has been carried out for D. nihonkaiensis. The present analysis demonstrates the usefulness for detecting sequence variation among D. nihonkaiensis population. Two haplotypes (H_33 and H_35) were slightly more mutated than the other haplotypes in Type A (Fig. 2). Since the bootstrap value of Type A was low in the phylogenetic analysis (Fig. 1B), further analysis of the isolates closely related to both haplotypes is required to determine whether these two haplotypes belong to a new type.

There has been no report about co-infection with genetically different D. nihonkaiensis in humans and salmon, and the present study is the first confirmation for this phenomenon in the definitive and second intermediate hosts. As in the present clinical cases, genetic analysis of the tapeworms made it possible to clearly identify multiple infections. Especially, in the case of excretion of one scolex and two strobilae, it is also effective for subsequent treatment and follow-up care to distinguish whether the strobilae are derived from the same specimen.

Species of Pacific salmon, the second intermediate hosts of D. nihonkaiensis, migrate widely across the northern Pacific including the Okhotsk and the Bering seas, which is an endemic area of this cestode, until returning to their respective rivers where they were born. It remains unclear where the Pacific salmon are becoming infected with this cestode. The genetic separation of D. nihonkaiensis populations might be related to regional differences in endemic areas, although no molecular evidence about this has been shown yet. There would be at least two ecological populations if the theory about brackish-water origin of D. nihonkaiensis is correct (Muratov, 1992; Kuchta et al., 2015, 2017): one of the Asian coastal areas and one of the North American coastal areas. In the present analysis, one isolate (Table 1, US361b) from a wild pink salmon captured at south-central Alaska (Kuchta et al., 2017) was included in Type A. Moreover, two isolates (Table 1, Not named 9 and 10) from two Swiss clinical cases suspected to be caused by eating wild Pacific salmon imported from Canada or North America (Wicht et al., 2007) were also included in Type A, suggesting that D. nihonkaiensis population does not exhibit regional genetic differences based on the cox1 sequence analysis.

Further comparative haplotype analysis of the isolates from definitive and second intermediate hosts collected in Pacific and Atlantic areas may be necessary to clarify the ecology of D. nihonkaiensis.

Conclusion

The present updated assessment of the genetic diversity of D. nihonkaiensis populations by phylogenetic and haplotype network analyses targeting the partial cox1 gene showed the presence of three genetic lineages in this cestode species. The present study also provides the first published evidence of co-infection with genetically different D. nihonkaiensis isolates in the same definitive or second intermediate host individual.

Funding

This study was supported by a Grant-in-Aid for Scientific Research (C) (20K05902) from the Japan Society for the Promotion of Sciences ().

CRediT author statement

Niichiro Abe: Conceptualisation, Methodology, Validation, Formal Analysis, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualisation, Supervision, Investigation. Takashi Baba: Methodology, Formal Analysis, Data Curation, Writing - Review & Editing. Yoshitaka Nakamura: Methodology, Investigation, Data Curation, Writing - Review & Editing. Shintaro Murakami: Methodology, Investigation, Data Curation, Writing - Review & Editing. All authors read and approved the final manuscript.

Data availability

The sequences generated in the present study were deposited in DNA Data Bank of Japan (DDBJ) under the accession numbers LC589648-LC589661.

Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.