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

PCR amplification and DNA sequencing

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.

Sequence data analyses

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).

Mitochondrial DNA polymorphism analysis

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).

Phylogenetic analysis

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).

Results

Analysis of indels

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 distribution and diversity

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.

Phylogenic analysis of haplotype data

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 relationships among populations

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.

Discussion

In the present study, a method for performing long PCR using single nematodes was directly developed for individually sequencing each mitogenome. Sequences of about 8 kb, corresponding to 55% of the whole mitochondrial genome of PWN were obtained for 285 nematode individuals using this method. A total of 160 polymorphic sites were detected in these sequences, leading to a frequency of polymorphism as high as one polymorphic site per 51 bp. However, a low SNP frequency was detected in tRNA genes, although this is consistent with tRNAs gene sequences being more conservative than that of other genes (Van der Veer and De Vries, 2004). The mitogenomic polymorphisms of the PWN have been mainly used in interspecies comparative studies (Pereira et al., 2013; Sultana et al., 2013). Their application to intraspecific studies is limited, and a lack of polymorphisms’ evaluation is evident.

Mitogenome sequences (three complete and three incomplete) of six PWN isolates, one from Japan, two from Korea, and three from Portugal, have been reported so far. No variations were detected between the two Korean sequences (14,788 bp) or among the three Portuguese sequences (about 12 kb). However, 171 polymorphic sites containing 150 SNPs were confirmed among the sequences from Ibaraki, Korea, and Portugal, whereas sequences from Kyushu alone showed 158 SNPs.

Sequences of four PWN mitochondrial genes (cox1, cytb, nad5, and rrnS) have been revised to date. In cox1, 52 SNPs were found among 36 isolates from Japan, Korea, China, Portugal, Canada, USA, and Mexico. In total, 18 SNPs were identified in cytb sequences of 20 isolates obtained from Japan, Korea, China, Portugal, and USA, and 11 SNPs were observed in nad5 sequences from 17 isolates derived from these five countries. In total, 13 polymorphic sites containing two SNPs and 11 indels were observed in rrnS sequences from six isolates. Except for the cox1, which has been sufficiently investigated on relatively large sample sizes, sequence diversity in mitochondrial genes is generally not high. However, the present study revealed 39 SNPs in cox1, 40 SNPs in nad5, and 6 SNPs in rrnS in sequences from the Kyushu region alone, which can be considered a very limited area compared to the areas examined in the studies mentioned above.

Past genetic studies on PWN employed exclusively DNA extracted from one isolate. As an isolate consists of many nematodes, several genotypes are included and sequencing based on isolates masks many low-frequency variations. Therefore, it is presumed that only high-frequency variations could be treated as sequences of the isolate, which seriously underestimates variation. Sequencing based on the method employed here allowed reading sequences of individual mitogenomic sequences, thereby detecting rare alleles/haplotypes. A considerably large sample size (285 individuals) might also contribute to correctly estimate diversity. Unlike the worldwide PWN collections performed in previous studies, the present study targeted the Kyushu region only, and used a localized collection of PWN within an extremely small area. Nevertheless, it was evident that extremely high sequence diversity was preserved in the mitochondrial genome of PWN. Sequencing based on individuals instead of isolates is therefore a powerful approach for population genetic studies.

In the phylogenetic tree produced for haplotypes, 12 (ht-01, ht-04–14) out of the 14 haplotypes distributed in northern Kyushu belonged to Clade I. On the contrary, there was no evident trend in the haplotypes from southern Kyushu (Table 3 and Fig. 2). Sasebo (Nagasaki prefecture) was the first invasion point of PWN in Asia (Yano, 1913) from which PWN gradually spread to Shikoku and Honshu (Futai, 2008). It was then disseminated from Japan to Korea and China (Mamiya and Shoji, 2009). Furthermore, PWN populations appearing in Portugal in 1999 were assumed to have originated from East Asia (Metge and Burgermeister, 2006). Based on the method presented here, the propagation route of PWN could be further examined using mitochondrial haplotypes from invaded countries (Japan, Korea, China, and Portugal) and native countries (U.S.A and Canada).

Based on the past records of pine wilting and timber import and movement, Kishi (1988) estimated three possible routes (Sasebo route, Nichinan route, and Aira route) through which PWN invaded and expanded into Kyushu. In the Sasebo route, PWN invaded from Sasebo, northwestern Kyushu in 1905, and extensively dispersed into the northern region of Kyushu. Because ht-13 was the major haplotype in the three populations from northeastern Kyushu (Itoshima, Yukuhashi, and Chikujo), this haplotype might characterize populations dispersed through the Sasebo route. Haplotype-13 was also relatively frequent in Amakusa in central western Kyushu, it was likely held with high frequency in the PWNs of the Sasebo route and then transmitted to northern Kyushu and Amakusa regions. According to Kishi (1988), pinewoods withered by PWD in Nagasaki prefecture were transferred to a pulp mill in Yatsushiro (Kumamoto prefecture), extending PWD to the central western region of Kyushu. A PWN group that irrupted from Nichinan in Miyazaki prefecture (Nichinan route) in 1939 expanded north and south along the east coast of Kyushu (Kishi, 1988). The PWN collections from four populations (Shintomi, Miyazaki, North Nichinan, and South Nichinan) in the coastal area of Miyazaki prefecture showed similar haplotype compositions, and ht-21, which presented an extremely high frequency, might be considered as the specific haplotype of the Nichinan route. The PWNs expanding through the Aira route (third route) invaded from Aira (Kagoshima prefecture) in 1942 and dispersed to the western region of Kagoshima prefecture (Sendai and Ibusuki) (Kishi, 1988). Although four haplotypes appeared in the Sendai population, most nematodes (84%) were ht-21. Because this haplotype might be considered specific to the Nichinan route and haplotype composition in Sendai was somewhat similar to that of populations dispersing through the Nichinan route, the Sendai population might have been recently influenced by the Nichinan route. However, the three low-frequency haplotypes (hts-24–26) that were not observed in other populations might persist as a vestige of the Aira route. The Ibusuki population, which was also regarded as dispersed through the Aira route, shared only one haplotype (ht-21) with Sendai and its frequency was quite low. The Ibusuki population, however, shared three haplotypes (hts-02–04) with Matsuura, northwestern Kyushu, which is geographically distant (Table 3 and Fig. 3), suggesting that artificial long-range migration might have occurred. Thus, in the region occupied by Aira route PWNs, profound changes in haplotype composition might be rapidly progressing by natural migration or artificial transfer from other nematode strains (Fig. 1).

Analysis of the genetic structure of the PWN population in Kyushu revealed high genetic differentiation among the 12 populations (Gst = 0.331). While 33.1% of the total genetic variability was due to variation among populations, 66.9% was due to variation within populations (Mouhaddab et al., 2015). The Nm among populations was nearly 1 (1.01) evidencing that PWN migration among populations had little effect on the genetic differentiation within the Kyushu region (Wright, 1951; Lowe et al., 2004). In the NJ tree of the 12 populations (Fig. 3) based on Ds (Table 4), the three populations within northeastern Kyushu (Itoshima, Yukuhashi, and Chikujo) and the four populations within southeastern Kyushu (Shintomi, Miyazaki, North Nichinan, and South Nichinan) formed cohesive clades. Average Gst values for the three populations of northeastern Kyushu and four populations of southeastern Kyushu were 5.7% and 1.9%, and their Nm values were 8.3 and 24.4, respectively. Thus, genetic differentiation among populations within these areas was much smaller than within the whole Kyushu area. Additionally, due to the large value of Nm, it could be inferred that genetic differentiation within the referred local regions occurred due to genetic drift (Wright, 1951; Lowe et al., 2004).

When the three population in northeastern Kyushu and four populations in southeastern Kyushu were grouped as metapopulations and analyzed with the remaining five populations (Matsuura, Karatsu, Amakusa, Sendai, and Ibusuki), average Gst and Nm for these seven populations were recalculated as 32.1% and 1.06, respectively. Because Nm was nearly 1, the influence of gene flow on the genetic structure of those areas did not appear to be substantial (Wright, 1951; Lowe et al., 2004). In nature, PWD was expanded because PWN was propagated by vector insects. Due to the severe decline of the damaged pine forests in the past, the fragmentation of PWN populations has progressed and this geographical separation might have decreased the natural long-distance gene flow. Matsuura and Karatsu populations in the northwestern region and the Ibusuki population in the southernmost region formed a loose clade in the NJ tree (Fig. 3). Matsuura and Karatsu are geographically close, whereas Ibusuki is geographically distant from the other two (Fig. 1). As Ds among the three populations were not too large (Table 4), a long-range migration of PWN might have occurred by human activity.

According to Nose et al. (2009), the two populations of Ibusuki and Sendai in southwestern Kyushu were similarly affected by the PWN expanded from Aira in Kagoshima prefecture, as these two populations presented similar haplotype composition and shared a common dominant haplotype that was not very frequent in other areas. However, in the present study, the haplotype composition of the two populations was quite distinct. A total of 84% of the individuals in Sendai showed the haplotype (ht-21) appearing in the four populations (Shintomi, Miyazaki, North Nichinan, and South Nichinan) of the southeastern region. In Ibusuki, which is located in southernmost area of Kyushu, the only haplotype common to Sendai was ht-21, but its frequency was quite low.

Population size of nematodes invading an area seems to drastically decrease following the decline of the pine forest, and small populations are maintained or nearly extinguish if the pine forest is extensively destroyed. Therefore, through migration from other areas, the genetic composition of these populations will drastically change. In extremely small populations, drastic changes might occur by random drift (Miko et al., 2009), which might be responsible by the dynamics observed in the distribution of PWN strains in the Kyushu region. In the present study, PWN sampling was performed in 2012-14, whereas that of Nose et al. (2009) was performed in 2006-07. Therefore, population composition changes seem to have occurred within a single decade. Because there is a distance limit for vector insects’ natural travel, human activities might be heavily involved in genetic dynamics.

Investigations on PWN genetic variations have been conventionally concentrated on the nuclear genome. In recent years, the number of studies on the mtDNA of PWN gradually increased, although its utilization mainly involved comparing Bursaphelenchus xylophilus with close relative species (Pereira et al., 2013; Sultana et al., 2013; Moreira et al., 2014). In Nematode, it has been widely confirmed that there are intraspecific mutations in mitogenomes (Thomas and Wilson, 1991), and intraspecific mutations are also recognized in some genes of pinewood nematodes (Valadas et al., 2013). In the present research, intraspecific variation surveys were quantitatively evaluated by the mitochondrial genome scale.

In previous intraspecific variation studies, a few mitochondrial genes were used, and sample (isolate) sizes were small. Furthermore, there was no analysis of individual nematodes. All these facts hampered the detection of detailed intraspecific variation (Valadas et al., 2013). In the present study, we clarified the high diversity of PWN mitochondrial genome by using a wide range of nucleotide sequence information and sufficient number of individuals (285 individuals) that enabled genetic population analysis based on haplotype analysis. Furthermore, the present study evidenced that the haplotype analysis of a single mitochondrial genome sequence is a valid approach for several aspects of PWN studies, including population phylogenetic analysis, genetic diversity evaluation, artificial migration possibility, and invasion and expansion pathways’ elucidation.