High Mitochondrial Genome Diversity and Intricate Population Structure of Bursaphelenchus xylophilus in Kyushu, Japan

The nematode samples used in this study were collected from 12 different forests in Kyushu, Japan where PWD was reported from 2012 to 2014 (Fig. 1). Pine chips were obtained from five points in the trunk, using a drill with a diameter of 16 mm. Nematodes were extracted from the chips using the Baermann funnel method (Iwahori and Futai, 1993), cultured on Botrytis cinerea Pers. (1794) grown on barley culture medium, and incubated at 25°C for 10 d to 20 d. After successful rearing, nematodes were removed from the culture medium, and stored in 5-ml tubes containing distilled water, at 4°C until use. Mitochondrial DNA (mtDNA) sequencing was performed for each individual nematode. The number of nematode individuals in each population for which mtDNA sequences were successfully obtained and used in further analysis are listed in Table 1. The number of trees in each population that became the source of the sequenced nematodes are also listed.

Figure 1

Figure 2

A primer pair for long PCR was designed in Oligo7 (Molecular Biology Insights), based on the complete mitochondrial genome sequence of Bursaphelenchus xylophilus (14,778 bp; NCBI Accession No GQ332424). The forward (5′-TCCTCCATTAAGAACTTTAGGGCAT-3′) and reverse (5′-TACAGTCAAAGCAATAGGACGAGA-3′) primers designed produced an amplicon about 8 kb long that covered more than 50% of the mitochondrial genome.

Each nematode was used for direct amplification of mtDNA. Under a binocular stereomicroscope, each individual was separated from the nematode suspension using a needle and divided into two halves using a scalpel. Both halves were then transferred to a 25-µl PCR mixture [1× Gflex PCR Buffer (Takara), 0.43 µM each primer, 0.625 units Tks Gflex DNA Polymerase (Takara)]. The PCR reaction was carried out in the SureCycler 8800 (Agilent Technologies) device and comprised an initial denaturation at 94°C for 1 min, followed by 30 cycles at 98°C for 10 s and 68°C for 8 min. The resulting amplicon was electrophoresed on 1% agarose gel with ethidium bromide, and the ∼8-kb-long fragment obtained was excised from the gel with the MagExtractor®-PCR & Gel clean up (Toyobo), following the manufacturer’s instructions. Using the excised fragment as DNA template, a sequencing library was prepared using the TruSeq Nano DNA Library Prep Kit (Illumina) and mtDNA sequencing was carried out in the Illumina MiSeq (Illumina) platform, following the manufacturer’s protocols.

The quality of the obtained sequences was evaluated. Only the sequences with quality over Q20 were extracted and were used for assembly into the mtDNA sequence (ca. 8 kb) on GS Reference Mapper Software (Roche). The corresponding region of the complete mitochondrial genome sequence of Bursaphelenchus xylophilus (GQ332424) was used as reference. Contigs with an average depth of 451 for the 285 individuals were used for sequence alignment and for the identification of single nucleotide polymorphisms (SNPs) and insertions and deletions (indels) in Sequencher (Gene Codes). The exact coding region of each gene was determined based on the annotation listed for GQ332424 and on basic local alignment search tool (BLAST) analysis conducted on the National Center for Biotechnology Information (NCBI) web platform (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

The evaluation of genetic polymorphisms, including the number of polymorphic sites, nucleotide substitutions, and transitions and transversions without indels, was carried out in DnaSP (Librado and Rozas, 2009). This software was also used to determine the number of haplotypes (Nh) and nucleotide and haplotype diversities. Haplotype diversity (Hd: Nei and Tajima, 1981) and Nh were computed to evaluate mitogenomic diversity in each population. MEGA (Kumar et al., 2016) was used to translate protein-coding genes based on the invertebrate mitochondrial codon usage table, and synonymous and non-synonymous substitutions in protein-coding regions were identified. The number of variations, including nucleotide substitutions in genes and non-coding regions, and variations in amino acid sequences, were also determined. The significant differences in SNP densities between genes or regions were examined by Student’s t-test in Excel (Gosset, 1908).

Haplotype distribution and diversity (Hd) in each population were investigated independently. Haplotype sequence data were imported into MEGA to build a phylogenetic tree of haplotypes using the maximum likelihood (ML) method based on the Kimura 2-parameter model (Kimura, 1980), with 1,000 bootstrap replications. Phylogenetic relationships among the 12 populations in Kyushu were analyzed using haplotype frequencies in populations. Haplotype composition and frequency of haplotypes in each population were used to calculate genetic differentiation (Gst), gene flow (Nm), and genetic distances (Ds) among populations. A phenogram was also drawn for the 12 populations based on the neighbor-joining (NJ) method (Saitou and Nei, 1987).

A fragment of about 8,060 bp, corresponding to 54.5% of the PWN mitochondrial genome (ca. 14.8 kb), was sequenced. In the PWN mitochondrial genome, there is an extremely AT-rich and extensive non-coding region of about 1.7 kb (13,129-14,778 of GQ332424), while all mitochondrial genes are placed in the remaining 13.1 kb (Sultana et al., 2013). The sequence decoded in the present study corresponds to 61.4% of the latter fragment.

A total of 158 SNPs were detected in mtDNA sequences, corresponding to one SNP per 51 bp, on average (Table 2). In total, 22 genes (8 protein-coding genes, 12 transfer RNA (tRNA) genes, and 2 rRNA genes) and 6 non-coding regions were confirmed in the 8,060 bp sequence. The number of SNPs in each gene/region is shown in Table 2. Nine tRNA genes (trnD, trnG, trnH, trnA, trnP, trnV, trnW, trnE, and trnY) showed no SNPs. Thus, the 158 SNPs corresponded to 139 SNPs from protein-coding genes, 3 SNPs from the remaining 3 tRNA genes, 13 SNPs from 2 rRNA genes, and 3 SNPs from non-coding regions.

A total of 88% (139/158) of the SNPs were found in the protein-coding genes, leading to an average occurrence of one SNP/41 bp. Most of these SNPs (40 SNPs, 1 SNP/39 bp) were found in nad5, but cox1 showed the highest average occurrence of SNPs (39 SNPs, 1 SNP/30 bp). Genes coding for NADH dehydrogenase subunits 1 and 6 (nad1 and nad6) presented high number and average occurrence of SNPs (22 SNPs, 1 SNP/40 bp and 13 SNPs, 1 SNP/33 bp, respectively). Genes coding for NADH dehydrogenase subunit 3 (nad3) and ATP synthase subunit 6 (atp6) showed a relatively low average occurrence of SNPs, as only six (1 SNP/57 bp) and seven (1 SNP/60 bp) SNPs were found in each sequence, respectively. Cytochrome c oxidase subunit 2 (cox2) was the most conservative among the eight protein-coding genes (seven SNPs, 1 SNP/99 bp). Because only three SNPs were detected in 3 of the 12 tRNAs for which nucleotide sequences were decoded (i.e., trnC, trnM, and trnS), the average occurrence of SNPs in tRNA genes was extremely low (1 SNP/223 bp). According to the predicted secondary structure of mitochondrial tRNA genes, 20 of the 22 tRNA genes identified in PWN have a unique structure; a loop of variable size (TV-replacement loop; 6-12 bp) is thus displayed instead of the TΨC arm loop (Sultana et al., 2013). In the present study, the SNPs of trnC and trnS occurred in the TV-replacement loop and TΨC arm loop, respectively; the SNP of trnM occurred close to its 5′ extremity. The two rRNA genes (rrnS and large subunit rRNA gene (rrnL)) found here contained six and seven SNPs, respectively, thereby showing a lower average occurrence of SNPs than the protein-coding genes (1 SNP/116 bp for rrnS and 1 SNP/134 bp for rrnL). The total length of the six non-coding regions identified here was 50 bp (1-20 bp in each region). Generally, non-coding regions are the most variable in the genome (Van der Veer and De Vries, 2004), which was also revealed for the partial mitochondrial genome of PWN. For example, the region comprising cox1-trnC was only 15 bp but contained two SNPs, thereby presenting an extremely high occurrence rate of SNPs (1 SNP/7.5 bp). In the trnD-trnG region, one SNP was detected in three base pairs. The differences in SNP densities between protein-coding genes and tRNA genes or rRNA genes were significant (P ⩽0.01) according to the t-test.

The 158 SNPs corresponded to 122 transitions (Ts) and 36 transversions (Tv), leading to a Ts/Tv ratio of 3.4 (Table 2). In total, 114 of the 122 Ts were found in the protein-coding genes. These genes presented Ts/Tv ratios of 2.5-8.8 (average 4.6), which indicated that Ts were more frequent that Tv. The cox1 gene contained most Ts and showed the highest Ts/Tv ratio (8.8). Nad5, which harbored the largest number of SNPs, held more Tv than cox1 and therefore showed a lower Ts/Tv ratio (3.0). Only one Tv was identified among the five SNPs observed in NADH dehydrogenase subunit 4L gene (nad4L). On the other hand, no Tv appeared in the cox2. The Ts/Tv ratio of the two rRNA genes was 1.2, which means that the number of Ts and Tv was almost identical. Only Tv were detected in all of the tRNA genes, while in the non-coding region there were two Tv in cox1-trnC and one Ts in trnD-trnG (Table 2).

Amino acid variation was examined for all SNPs detected in the protein-coding genes. The numbers of synonymous and non-synonymous substitutions (Syn and Non-syn, respectively) in amino acids were 118 and 20, respectively (Table 2). Within the 118 Syn, 104 were Ts and 14 were Tv. A total of 97% (115) of the Syn were due to mutations at the third position of the codon, while the remaining three Syn were caused by mutations at the first codon position (C/T substitution in nad5, C/T substitution in nad6, and C/T substitution in nad1). In all, 10 (50%) of the 20 Non-syn were caused by mutations at the second codon position, 7 were due to mutations at the first codon position, and 2 were due to mutations at the third codon position. One Non-syn recognized in nad5 corresponded to simultaneous changes in the second and the third positions of the codon. The overall Syn/Non-syn ratio was 5:9. Thus, Syn occurred about six times more frequently than Non-syn, which means that most variations in the partial mitochondrial genome of PWN could be considered neutral (Kimura, 1983). However, the Syn/Non-syn ratio in nad3 (2.0) and nad5 (3.3) were lower than in other genes, suggesting that relatively high selective effects might have occurred in these two genes.

Because the mitochondrial genome of PWN is extremely AT-rich (Sultana et al., 2013), it is difficult to accurately estimate the length variation of long homopolymer (poly-A and poly-T) regions (Linnertz et al., 2012). Therefore, a PCR was employed in the present study to allow the occurrence of artificial length variation due to possible PCR slippage in simple sequence repeats (Hauge and Litt, 1993; Murray et al., 1993). All indels (length variation) recognized in homopolymer regions were excluded due to the low quality of the identification of nucleobases (i.e., Q score) of their last parts in many individuals. Only two indels were detected outside homopolymer regions and no indels were found in the eight protein-coding genes (Table 2). Only 1 of the 12 tRNA genes (trnY) displayed one-base deletion in its anticodon stem, which was found in 1 of the 285 nematodes investigated. In the cox1-trnC inter-genic spacer region, one four-base deletion was observed (Table 2).

Haplotype analysis was performed for the 285 PWNs sampled, and 30 haplotypes (hts-01–30 in Table 3) were detected based on 160 polymorphic sites that consisted of 158 SNPs and two indels. Polymorphic sites of haplotypes and the position of each SNP in relation to the reference sequence (GQ332424) are displayed in Table S1.

The frequency of each haplotype detected in each of the 12 investigated populations is shown in Table 3. The most frequent haplotype in Kyushu was ht-21 (34.4%, 98 nematodes), although ht-01 (15.1%, 43 nematodes) and ht-13 (14.7%, 42 nematodes) also showed relatively high frequency. On the other hand, 14 haplotypes (ht-05, hts-08–09, hts-11–12, hts-16–19, ht-23, hts-25–26, and hts-29–30) were identified as extremely rare, as they were detected in only one nematode. In total, 20 of 30 haplotypes were population-specific, while the remaining 10 haplotypes were distributed in multiple populations. The most frequent haplotype (ht-21) was observed in 7 of the 12 populations, including all 6 populations in the southern region of Kyushu (Shintomi, Miyazaki, North Nichinan, South Nichinan, Sendai, and Ibusuki) and in the central western region (Amakusa). Haplotypes 10, 13, and 14 were found in the northern region (Karatsu, Itoshima, Yukuhashi, and Chikujo) and central western region (Amakusa), and four haplotypes (hts-01–04) were detected in the northern and southern Kyushu regions. Therefore, haplotype composition notably differed among regions.

Haplotype diversity in the 12 populations varied from 0.30 to 0.83, and its average was 0.55 (Table 3); in eight populations it was above 0.5. Particularly, populations from Amakusa, Karatsu, and Ibusuki showed Hd above 0.8, which is considered high haplotype polymorphism, while populations from Sendai, Miyazaki, and Itoshima showed Hd below 0.4 and Nh below 4. The haplotype diversity for the entire Kyushu region was 0.83.

The haplotypes found in three populations from northeastern Kyushu (Itoshima, Yukuhashi, and Chikujo) were similar to each other, and common haplotypes were also detected in four populations from southeastern Kyushu (Shintomi, Miyazaki, North Nichinan, and South Nichinan).

A phylogenetic tree for the 30 haplotypes was constructed based on the ML method, using the 158 SNPs identified in the present study (Fig. 2). The closely related Bursaphelenchus mucronatus (Accession GU177865) was employed as outgroup to determine the root of the phylogenetic tree.

The dendrogram clustered the 30 haplotypes into two Clades (Clade I and Clade II). A major phylogenic difference was observed between the two clades; there was a minimum difference of 58 SNP sites (between ht-26 of Clade I and ht-25 of Clade II) and a maximum difference of 105 sites (between ht-27 of Clade I and hts-02/03/21 of Clade II). In contrast, no significant difference emerged among haplotypes within each clade. Among the 24 haplotypes within Clade I, 91 sites differed between hts-10 and 30, while only four sites differed between hts-25 and 28, among the six haplotypes within Clade II. Clade I was further divided into two subclades (Subclades I-1 and I-2) and three single haplotypes (hts-17, 26, and 10). The sequences of these three haplotypes differed from that of the other 21 haplotypes contained in Clade I. Subclades I-1 and I-2 consisted of 11 and 10 haplotypes, respectively. However, the maximum disagreement within Clade I consisted of only 16 SNP sites between ht-30 (Subclade I-1) and ht-12 (Subclade I-2), and no significant difference was identified between the two subclades. Most of the haplotypes within Subclade I-1 were distributed in the northern Kyushu region, haplotypes within Clade II were mainly distributed in southern Kyushu region, and haplotypes in Subclade I-2 were detected in both northern and southern Kyushu regions.

Phylogenetic relationships among the mitogenomic sequences of six isolates (Ibaraki, Japan: AP017463; Korea: GQ332424 and NC023208; and Portugal: JQ429761, JQ514067, and JQ514068) were also analyzed. No sequence polymorphisms were identified for isolate within the same area in Korea and Portugal and, therefore, a phylogenetic tree was built after aligning three mitogenomic sequences (Ibaraki: AP017463, Korea: GQ332424, and Portugal: JQ429761) and the 30 haplotype sequences from Kyushu (hts-01–30). As a result, the Ibaraki (East Japan) sequence was assigned to Subclade I-1, whereas sequences from Korea and Portugal were allocated to Clade II. Sequences from Ibaraki and Portugal were highly similar to ht-13 and hts-21/23, respectively. The Korean sequence was not closely related to seven of the haplotypes within Clade II (hts-02, 03, 21, 23, 25, 28, and JQ429761). Comparing the Korean sequence with the other 32 sequences (haplotype and mitochondrial sequences) revealed 35 SNPs and 8 indels in 2 rRNA genes, and no significant differences for the 8 protein-coding genes and 12 tRNA genes.

Genetic distances among populations were 0.14 to 0.92, based on haplotype frequency in each population (Table 4). The maximum distance was obtained between Miyazaki (southeastern region) and Karatsu (northwestern region), and no common haplotypes were detected. On the other hand, the Ds among three populations in northeastern Kyushu (Itoshima, Yukuhashi, and Chikujo) and the Ds among four populations in southeastern Kyushu (Shintomi, Miyazaki, North Nichinan, and South Nichinan) were 0.14 to 0.37 and 0.14 to 0.35, respectively. Geographical distances among populations within each group were very small.

The NJ tree constructed for the 12 populations based on Ds (Fig. 3) revealed that the three populations in northeastern Kyushu and the four populations in southeastern Kyushu formed cohesive clades. Three populations, located in the northwestern region (Matsuura and Karatsu) and southernmost region (Ibusuki), formed a loose clade. The other two populations (Amakusa in the central western region and Sendai in the southwestern region) were clade-independent.

Figure 3

The Gst and Nm among populations were 0.33 and 1.01, respectively.

Table S1