Genetic Diversity in Echinococcus multilocularis From the Plateau Vole and Plateau Pika in Jiuzhi County, Qinghai Province, China

The Qinghai-Tibet Plateau is a highly endemic area of alveolar echinococcosis where a series of intermediate hosts, especially voles and pikas, are infected with Echinococcus multilocularis. The metacestodes of E. multilocularis are fluid-filled, asexually proliferating cysts, and they are mainly found in the host's liver in the form of tumor-like growths. In this study, we investigated the genetic variations of E. multilocularis in four mitochondrial (mt) genes, namely, NADH dehydrogenase subunit 5 (nad5), adenosine triphosphate subunit 6 (atp6), cytochrome c oxidase subunit 1 (cox1), and NADH dehydrogenase subunit 1 (nad1). The complete nad5, atp6, cox1, and nad1 genes were amplified separately from each hydatid cyst isolate using polymerase chain reaction (PCR) and then sequenced. Phylogenetic trees were then generated based on the combined mt genes using MrBayes 3.1.2 and PAUP version 4.0b10. The results showed that thirty of 102 voles and two of 49 pikas were infected with E. multilocularis. The genetic variation distances among all E. multilocularis samples were 0.1–0.4%, 0.2–0.4%, 0.1–0.6%, and 0.1–0.4% for nad5, atp6, nad1, and cox1, respectively. Compared to previous studies of the genetic diversity of E. multilocularis based on the cox1 gene, the genetic distances within the same group were 1.3–1.7% (Mongolia strain), 0.6–0.8% (North American strain), 0.3–0.6% (European strain), and 0.1–0.4% (Asian strain). Based on concatenated sequences of the nad5, atp6, cox1, and nad1 genes all haplotypes were divided into two clusters. In conclusion, the genetic diversity of E. multilocularis based on mt genes on a small local area is at low level but between different regions with long distance and different ecological environment each other, the genetic diversity is at relatively high level; genetic variation is higher in the nad1 gene than that in the other three mt genes. The results on a local scale provide basic information for further study of the molecular epidemiology, genetic differences and control of E. multilocularis in Qinghai Province, China.


INTRODUCTION
Echinococcus multilocularis is a small cestode that cause the parasitic zoonosis alveolar echinococcosis (AE), which was one of the 17 neglected tropical diseases (NTDs) prioritized by the World Health Organization (WHO) in 2012 (Agudelo Higuita et al., 2016). AE was discovered in the Nineteenth century; it has a considerable socioeconomic impact and until now has been sporadically found in humans (Nakao et al., 2007;Spahn et al., 2016). E. multilocularis is mainly distributed in holarctic regions, including North America (also found in southwestern Ontario), Europe, and Asia (Massolo et al., 2014;Oksanen et al., 2016;Deplazes et al., 2017;Trotz-Williams et al., 2017). Asia, especially China, has many highly endemic areas, such as Qinghai, Gansu, and Sichuan. E. multilocularis infects two kinds of hosts: the typical intermediate host (IH), a wide spectrum of mammalian species including small herbivorous, rodents (predominantly) and pikas, and the typical definitive hosts (DH), which are canids and mammalian species, including foxes, wolves, dogs, and cats (Vuitton et al., 2003;Deplazes et al., 2011;Hegglin and Deplazes, 2013;Conraths and Deplazes, 2015;Raoul et al., 2015;Knapp et al., 2016;Eckert and Thompson, 2017). Humans play the role of an aberrant IH for E. multilocularis and can be infected when ingesting eggs released into vegetables and food by adult worms. Then, the parasite larvae travel to internal organs, mainly the liver (Torgerson et al., 2010;Conraths et al., 2017).
Mitochondrial DNA (mtDNA) has an important role in the taxonomy of Echinococcus species. According to mtDNAs, E. granulosus species include 10 genotypes (G1-G10) and genetic variation in mtDNAs has also been found within E. multilocularis (Bowles et al., 1992;Nakao et al., 2007;McManus, 2013). From early studies, E. multilocularis has been divided into two geographical genotypes, M1 (Europe) and M2 (China, Alaska and North America), based on four nucleotide substitutions in the nad1 gene (Bowles and McManus, 1993;Okamoto et al., 1995). Nakao and his co-workers found 17 regional haplotypes based on three mt protein-coding genes: cytochrome c oxidase subunit 1 (cox1), cytochrome b (cytb), and NADH dehydrogenase subunit 2 (nad2) (Nakao et al., 2009). All of those studies show a relatively low mt nucleotide diversity within E. multilocularis. To date, more genes have been used for investigating genetic variation in E. multilocularis. Although the nuclear genome and the microsatellite targets are also used for molecular markers (Knapp et al., 2008;Umhang et al., 2014;Laurimaa et al., 2015;Karamon et al., 2017), mt genes are more commonly used in the analysis of genetic variation in E. multilocularis, because they have a great copy number and evolve at a high rate (Brown et al., 1979;Ciesielski et al., 2016).
In our study, we chose mtDNAs as molecular markers to assess the diversity of E. multilocularis in plateau voles (Neodon/ Microtus fuscus) and plateau pikas (Ochotona curzoniae). In addition to the common mt cox1 and nad1 genes, NADH dehydrogenase subunit 5 (nad5) and adenosine triphosphate subunit 6 (atp6) were also included for genetic variation analysis, which provide basic information on molecular epidemiology, genetic variations or differences in E. multilocularis in China.

Ethics Statement
All animals were handled in strict accordance with good animal practice according to the Animal Ethics Procedures and Guidelines of the People's Republic of China, and the study was approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (No. LVRIAEC2012-007). In addition, all mice were handled in strict accordance with the animal protection laws of the People's Republic of China (A Draft of an Animal Protection Law in China released on September 18, 2009).

Biosecurity Statement
The biohazards, biological select agents, toxins, restricted materials or reagents involved in the research have been carried out according to the "Measures for the prevention and control of environmental pollution of hazardous chemicals" issued by the state environmental protection administration (No. 27_2005-03-30) and the "Regulations on the management of medical waste" promulgated by the state council of the People's Republic of China (2003-06-16).

Sample Collection
The specimens of E. multilocularis were collected from Jiuzhi County, Golog Autonomous Prefecture, Qinghai Province (33.32 • N, 100.53 • E) at an average elevation of 4,100 m above sea level. All samples were collected at random in the field amid a "Rodent control program" that was carried out by the local Center for Animal Disease Control and Prevention. Cystic lesions were collected from two intermediate host species: plateau voles and pikas. They were isolated and placed in 70% (v/v) ethanol before being sent to the laboratory at Lanzhou Veterinary Research Institute (LRVI). The lesions were used for DNA extraction, PCR amplification and sequencing to confirm the parasite species and analyze genetic variations. Additionally, one a Based on complete concatenated nad5, atp6, nad1 and cox1 nucleotide sequences.
E. multilocularis metacestode sample from the Xinjiang Uygur Autonomous Region was preserved by our laboratory.

DNA Extraction
The cystic lesions were rinsed repeatedly using phosphate buffered saline (PBS) to remove the ethanol, and then, the lesions were disrupted using TissueLysers before DNA extraction. With the use of the TIANamp Genomic DNA Kit (TianGen Biotech, Haplotype Gene GenBank accession number Beijing, China), the total DNA was extracted according to the manufacturer's protocols.

PCR Amplification and Sequencing of MT Genes
Four pairs of primers were designed for nad5, atp6, nad1, and cox1 genes (

Sequence Alignment and Analysis
The open reading frames of all raw nucleotide sequences of the nad5, atp6, nad1, and cox1genes of each isolate were assembled, edited and concatenated into a total sequence using the software package SeqMan (DNAStar, Inc., Madison, WI, www.dnastar.com/t-megalign.aspx) and Clustal W2 (online software, http://www.ebi.ac.uk/Tools/msa/clustalw2/, serviced by EBI, the European Bioformatics Institute). The total Indicates the nucleotide is same as the E. multilocularis reference mtDNA sequence (AB018440).
sequences (one of the sequences having 100% similarity with the others was chosen) were aligned using BioEdit v7.2.3 and ClustalX 1.83 (Thompson et al., 1997;Hall, 1999). The Megalign software (DNAStar, Inc., Madison, WI) was used to calculate the genetic divergence between different haplotypes. In phylogenetic reconstruction, E. shiquicus (GenBank accession no. AB159136) was used as an out group, because it was considered the sister species of E. multilocularis (Nakao et al., 2007). Bayesian analysis was performed to combine four datasets by using MrBayes 3.1.2 with a general time reversible (GTR) model of DNA substitution and a gamma distribution rate variation across sites (Ronquist and Huelsenbeck, 2003). In this model, a Metropolis-coupled Markov chain Monte Carlo analysis was run for 1 million generations, and trees were sampled every 100 generations. The first one-fourth of 10,000 trees were treated as burn-in. The command was executed until the average standard deviation of the split frequencies was lower than 0.01 (Nakao et al., 2009;Zhao et al., 2013). Phylogenetic reconstruction was completed using PAUP (phylogenetic analysis using parsimony) version 4.0b10 (Cummings, 2014) by the neighbor-joining (NJ) and maximum parsimony (MP) methods. The NJ tree was performed from HKY85 distances with inverse-squared weighting (power = 2). The MP tree was executed using a heuristic search with TBR branch swapping options and 1,000 random sequence additions. The clade of trees was estimated with 1,000 replicates of bootstrap (BT) analysis. The descriptive tree statistics, tree length (TL), consistency index (CI), retention index (RI), rescaled consistency index (RC), and homoplasy index (HI) for Maximum Parsimonious Tree were engendered (Zhao et al., 2013).

E. multilocularis cysts
A total of 102 voles (N. fuscus) and 49 pikas (O. curzoniae) were captured in Jiuzhi County. Of these animals, 32 (30 voles, infection rate: 29.41%; 2 pikas, infection rate: 4.08%) had cystic lesions in their viscera. From these 32 individuals, 39 cysts were isolated. Twenty-seven cysts were found only in the liver, and the others were found not only in the liver but also in the lungs, the intestinal wall and the muscles ( Table 2).

Characterization of Haplotypes and Phylogenetic Trees Analysis
The combined sequences of nad5 (1,575 bp), atp6 (516 bp), nad1 (894 bp) and cox1 (1,608 bp) genes were distributed into 15 haplotypes including one from Xinjiang Uygur Autonomous Frontiers in Microbiology | www.frontiersin.org Region and one from Yushu Tibetan Autonomous Prefecture (Table 2), which were designated as JZ01 to JZ13 (haplotypes from Jiuzhi), XJ (haplotype from Xinjiang), and YS (haplotype from Yushu), respectively. When coming from a single gene, the number of haplotypes decreasedto 7 in nad5, 2 in atp6, 4 in nad1, and 11 in cox1 (Table 3). Phylogenetic analyse showed that the genetic diversity of individual gene changes was 0.1-0.4% (nad5), 0.2-0.4% (atp6), 0.1-0.6% (nad1), and 0.1-0.4% (cox1). Within the four genes, cox1 had the most mutation sites (with 9 sites), followed by 7 sites in nad5, 2 sites in atp6 and 2 sites in nad1 (Table 4). When we compared the sequences with a reference sequence for partitioned nad5, atp6, nad1, and cox1 genes, there were only two haplotypes in the four genes. The genetic distance for JZ01 wa 99.7, 99.8, 99.9, and 99.8%, and for JZ03, it was 99.7, 99.6, 100, and 99.9%. The sequence analysis with the reference sequence in this study further predicted that all of the samples or specimens isolated from plateau voles and pikas belonged to E. multilocularis.
The phylogram based on the maximum parsimony method precisely revealed that the 14 haplotypes were divided into two large haplogroups: the JZ03, JZ05-JZ10, and JZ12 haplotypes had a close relationship with each other (classified as C1 haplogroup) and the rest of the haplotypes gathered together (classified as C2 haplogroup) (Figure 1). The haplotype XJ from Xinjiang Province is at the base of the phylogenetic tree. In this phylogenetic tree, the bootstrap value of three nodes was almost the same and was low. The maximum genetic diversity between the XJ haplotype and C1 haplogroup reached 0.2%, and the divergence between haplotype XJ and C2 was 0.2%. When compared haplogroup C1 with C2, the rate of pairwise divergence was the same. Haplotype cladograms of E. multilocularis were carried by Bayesian and NJ methods using the combined nucleotide sequences of the cox1, nad1, nad5, and atp6 genes. The two methods both showed a similar result with two clades (Figure 2). Furthermore, the network of those haplotypes also divided into the two same groups (Figure 3). Compared to previous studies of the genetic diversity of E. multilocularis based on the cox1 gene, we found that the sequences could be divided into 4 groups; Mongolia strain, North American strain, European strain and Asian strain. The highest genetic divergence is 1.3-1.7% (samples collected from Inner Mongolia, Mongolia and Russia), followed by 0.6-0.8% (samples collected from the USA: Alaska and Indiana), 0.3-0.6% (samples collected from Austria, Estonia, France, Poland, Slovakia and Russia) and 0.1-0.4% (samples collected from China: Sichuan, Kazakhstan, Japan: Hokkaido). The genetic distance within the same group is 0.1-0.3% (Mongolia strain), 0.4% (North American strain), 0.1-0.3% (European strain), and 0.1-0.4% (Asian strain) (Nakao et al., 2009;Ito et al., 2010;Konyaev et al., 2012;Laurimaa et al., 2015;Karamon et al., 2017). Similarly, according to the phylogenetic tree based on the cox1 gene constructed by MEGA6, all haplotypes clustered together with the haplotypes downloaded from NCBI (samples collected from China: Sichuan, Kazakhstan, Japan: Hokkaido), belonging to the Asian clade (not shown).

DISCUSSION
The genus Echinococcus has a unique reproductive cycle, including a sexual reproduction in the adult stage and asexual proliferation in the larval phase. Because of this reproductive model, Echinococcus has genetic monomorphism in local populations (Lymbery and Thompson, 1996;Haag et al., 1999;Nakao et al., 2007). In early studies, E. multilocularis was recognized as a subspecies of E. granulosus, but Rausch and Nelson thought that E. multilocularis had unusual morphological FIGURE 2 | Haplotype cladograms of E. multilocularis inferred by the neighbor-joining (NJ) and Bayesian methods using the concatenated nucleotide sequences of the nad5, atp6, nad1, and cox1 genes. and biological peculiarities and should be a separate genus. Then, phylogenetic analyses based on the mt genes also indicated that E. multilocularis was different compared with E. granulosus (Lukashenko and Zorikhina, 1961;Rausch and Nelson, 1963;Tappe et al., 2010). Now the genus Echinococcus includes at least 10 species: E. granulosus, E. equinus, E. ortleppi, E. Canadensis, E. intermedius, E. felidis, E. multilocularis, E. vogeli, E. oligarthra, and E. shiquicus (Thompson, 2008(Thompson, , 2017Nakao et al., 2013;Thompson and Jenkins, 2014).
Until now, a majority of published data on the gene diversity and phylogenetic analysis of the genus Echinococcus have been based on short mt genes, especially the cox1 and nad1 genes. From early studies, nine genotypes were classified in E. granulosus s. l. based on the diversity of the mt genes cox1 and nad1, and these two genes were applied in E. multilocularis genetic variation analysis (Bowles et al., 1995;Tappe et al., 2010;Romig et al., 2017). In our study, we used the 4 fulllength sequences of mt genes to analyze the diversity of E. multilocularis on the local scale. Herein, the sequence divergence within E. multilocularis wase 0.1-0.4% for nad5, 0.2-0.4% for atp6, 0.1-0.6% for nad1, and 0.1-0.4% for cox1. The result showed that the nad1 gene had more diversity than those of the other genes, which is in contrast with previous studies showing that atp6 had a comparatively more diversity within the Echinococcus spp. and E. multilocularis (Yang et al., 2005;Nakao et al., 2007). However, our results were in keeping with a previous study based on the cox1 gene. The variation in the nad5 gene showed a consistent results with the published data of the other parasites (Zhao et al., 2014;Lou et al., 2015). Regarding amino acid, cox1 exhibited the highest frequency of substitution (0.7476%, 4 substitution sites/535 sites), followed by apt6 (0.5814%, 1/172), nad5 (0.3817%, 2/524), and nad1 (0.3367%, 1/297). According to the results, even though the analysis of the cox1 gene was more meaningful from an evolutionary standpoint, the phylogenetic tree based on only one gene was not accurate. The phylogenetic trees analysis revealed that all 15 haplotypes divided into 3 groups. The 14 samples isolated from Qinghai-Tibet plateau divided into two clades, FIGURE 3 | The network of E. multilocularis generated using mt concatenated genes collected from voles and pikas. Each link between haplotypes indicates one mutational difference and the unlabeled nodes point out that inferred steps not found in the samples. and the haplotype isolated from Xinjiang province was in a separate group. In C1, the JZ09 haplotype was isolated from the lungs of voles, the JZ12 haplotype was isolated from the intestines of voles, and the remaining haplotypes were found in the liver of voles. In C2, the JZ13 haplotype was isolated from the lungs in voles, the JZ11 haplotype was isolated from the lungs of pikas, and the others were found in the liver of voles. The phylogram indicted that the different intermediate hosts and parasitic sites had no effect on the genetic diversity of E. multilocularis on a local scale. From the results of the genetic divergence within haplotypes based on the cox1 gene, we observed that the European strain has the lowest genetic diversity, and the Asian strain and Mongolia strain have the highest diversity. This also implied that the relationship maybe not important between geographical distances and genetic distances.
In conclusion, our study, based on the mitochondrial complete coding genes nad5, atp6, nad1, and cox1, shows a relatively low genetic variation among the samples from voles and pikas on a local scale. The genetic variation is higher in nad1 than that in the other three mt genes, and the concatenated mt sequences of the cox1, nad1, nad5, and atp6 genes are useful in phylogenetic reconstruction within E. multilocularis isolates.

AUTHOR CONTRIBUTIONS
JL, LL, HY, and WJ conceived and designed the experiments. JL, LL, and YF performed the experiments and the data analyses. JL prepared the figures and wrote the manuscript, WJ and HY provided improving paragraphs and BF and XZ provided very constructive suggestions for revisions.