Introduction
The genus Bursaphelenchus Fuchs, Reference Fuchs1937 (Nematoda: Aphelenchoididae, Parasitaphelenchinae) comprises over 140 described species (Kanzaki and Giblin-Davis Reference Kanzaki and Giblin-Davis2018; Ryss and Subbotin Reference Ryss and Subbotin2017; Zhao et al. Reference Zhao, Gu, Greenwood, Rogan, Ho and Taylor2024). Most species exhibit a mycophagous lifestyle, colonizing recently dead wood and utilizing various insect groups as phoretic hosts; only a minority are obligate plant parasites. Among these, two economically significant phytopathogenic species stand out: B. xylophilus (Steiner & Buhrer, 1934) Nickle, 1970, the causal agent of pine wilt disease responsible for devastating epidemics across Asia and Europe (Kikuchi et al. Reference Kikuchi, Cotton, Dalzell, Hasegawa, Kanzaki, McVeigh, Takanashi, Tsai, Assefa, Cock, Otto, Hunt, Reid, Sanchez-Flores, Tsuchihara, Yokoi, Larsson, Miwa, Maule, Sahashi, Jones and Berriman2011), and B. cocophilus (Cobb, 1919) Baujard, 1989, which causes red ring disease in multiple palm species throughout Central and South America (Subbotin and Del Prado-Vera Reference Subbotin and Del Prado-Vera2024). Both species are transmitted by insect vectors: B. xylophilus by cerambycid beetles and B. cocophilus principally by palm weevils (Ryss et al. Reference Ryss, Kazartsev and Kerchev2024).
Additionally, several other Bursaphelenchus species demonstrate pathogenicity toward conifers, including B. mucronatus Mamiya & Enda, 1979, and B. sexdentati Rühm, 1960 (Hunt 1993; Kanzaki and Giblin-Davis Reference Kanzaki and Giblin-Davis2018). Due to their severe economic impact, Bursaphelenchus species are classified as high-priority quarantine pests subject to strict international phytosanitary measures (Akbulut et al. Reference Akbulut, Yüksel, Serin and Edem2015; Gu et al. Reference Gu, Braasch, Burgermeister and Zhang2006; Kanzaki and Futai Reference Kanzaki and Futai2006). Notably, multiple species have been intercepted in imported wood packaging during phytosanitary inspections at Chinese ports. For instance, B. chengi Li, Trinh, Waeyenberge & Moens, Reference Li, Trinh, Waeyenberge and Moens2008, was intercepted in wood packaging from Taiwan at Nanjing Port (Li et al. Reference Li, Trinh, Waeyenberge and Moens2008), while B. parantoniae Maria, Fang, He, Gu & Li, Reference Maria, Fang, He, Gu and Li2015, was isolated from Pinus sp. packaging material imported from Belgium (Maria et al. Reference Maria, Fang, He, Gu and Li2015). Similarly, B. sinensis Palmisano, Ambrogioni, Tomiczek & Brandstetter, 2004, was first described in wood packaging exported from China to Austria (Palmisano et al. Reference Palmisano, Ambrogioni, Tomiczek and Brandstetter2004) and later documented in dead Pinus trees in Japan and China (Gu et al. Reference Gu, Zheng and Chen2008; Kanzaki and Futai Reference Kanzaki and Futai2007).
The description of numerous new Bursaphelenchus species, particularly cryptic and closely related taxa in recent years, highlights the necessity of robust molecular barcoding with multiple genetic markers to resolve their phylogenetic relationships and improve species delineation.
The phylogeny and classification of nematodes have long been subjects of debate among taxonomists (Coomans Reference Coomans2000). As maternally inherited, extranuclear genomes, mitochondrial genomes (mitogenomes) offer valuable insights into species identification and evolutionary genetics (Zardoya Reference Zardoya2020). Consequently, mitogenomic data have emerged as a powerful tool for resolving phylogenetic uncertainties that remain poorly supported by traditional multi-gene analyses (Cameron Reference Cameron2014; Du et al. Reference Du, Guo, Gao, Wang, Qing and Li2022).
To date, over 200 complete nematode mitogenomes have been deposited in GenBank. However, the majority represent animal-parasitic species, with plant-parasitic nematodes being significantly underrepresented (Gendron et al. Reference Gendron, Qing, Sevigny, Li, Liu, Blaxter, Powers, Thomas and Porazinska2024; Kern et al. Reference Kern, Kim and Park2020). Despite the ecological and economic importance of the genus Bursaphelenchus and its extensive species diversity, only three mitogenomes (B. mucronatus, B. okinawaensis Kanzaki, Maehara, Aikawa & Togashi, Reference Kanzaki, Maehara, Aikawa and Togashi2008, and B. xylophilus) were available (Kikuchi et al. Reference Kikuchi, Cotton, Dalzell, Hasegawa, Kanzaki, McVeigh, Takanashi, Tsai, Assefa, Cock, Otto, Hunt, Reid, Sanchez-Flores, Tsuchihara, Yokoi, Larsson, Miwa, Maule, Sahashi, Jones and Berriman2011; Sultana et al. Reference Sultana, Han and Park2013; Sun et al. Reference Sun, Shinya, Dayi, Yoshida, Sternberg and Kikuchi2020).
In this study, we sequenced and assembled the complete mitogenomes of three Bursaphelenchus species, including B. chengi, B. parantoniae, and B. sinensis, using a combination of Illumina and Oxford Nanopore sequencing platforms. We comprehensively annotated their mitogenomic architectures, including gene content and structural features. Furthermore, these new mitogenomic data were employed to elucidate the phylogenetic position of Bursaphelenchus within the Nematoda.
Materials and methods
Culture of nematodes
The B. chengi specimens were isolated from solid wood packaging materials originating from Taiwan (Li et al. Reference Li, Trinh, Waeyenberge and Moens2008), while B. parantoniae specimens were isolated from packaging wood imported from Belgium (Maria et al. Reference Maria, Fang, He, Gu and Li2015). The specimens of B. sinensis were collected from dying tree of Pinus massoniana in Ningbo, Zhejiang Province, China. All nematode specimens were extracted using a modified Baermann funnel technique and subsequently maintained in axenic cultures on potato dextrose agar (PDA) medium with Botrytis cinerea as the nutrient source.
DNA extraction and high-throughput sequencing
Genomic DNA was extracted from approximately 20,000 live nematodes per species using the Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech, Shanghai, China) following the manufacturer’s protocol. DNA concentration was quantified using a Qubit® fluorometer (Invitrogen, Carlsbad, CA, USA). For B. parantoniae, sequencing libraries were prepared with the TruSeq® DNA Sample Prep Kit (Illumina Inc., San Diego, CA, USA) and sequenced on an Illumina NovaSeq 6000 platform (150 bp paired-end reads). For B. chengi and B. sinensis, Nanopore libraries were constructed using the SQK-LSK114 Ligation Sequencing Kit (Oxford Nanopore Technologies, Oxford, UK) and sequenced on a PromethION P48 platform (Oxford Nanopore Technologies).
Genome sequence assembly
Quality filtering of raw reads was performed using FASTP (Chen et al. Reference Chen, Zhou, Chen and Gu2018) with the following criteria: (1) removal of reads less than 100 bp in length, (2) elimination of reads with Q-scores less than 20, and (3) exclusion of reads containing more than 3 uncalled bases (N). For B. parantoniae, initial de novo assembly using MitoZ (Meng et al. Reference Meng, Li, Yang and Liu2019) failed to generate complete circular mitogenomes. We therefore implemented a reference-guided approach: (1) Illumina reads were mapped to published reference nematode mitogenomes using NextGenMap (Sedlazeck et al. Reference Sedlazeck, Rescheneder and von Haeseler2013) with default parameters; (2) mapped reads were extracted with SAMtools (Li et al. Reference Li, Handsaker, Wysoker, Fennell, Ruan, Homer, Marth, Abecasis and Durbin2009); and (3) final assembly was achieved in NOVOPlasty (Dierckxsens et al. Reference Dierckxsens, Mardulyn and Smits2017) using the cox1 sequence from the initial MitoZ assembly as a ‘seed’ for the seed-extend algorithm. For B. chengi and B. sinensis, hybrid assembly was conducted using Unicycler (Wick et al. Reference Wick, Judd, Gorrie and Holt2017). Potential mitogenome contigs (10–30 kb) were identified via SeqKit (Shen et al. Reference Shen, Le, Li and Hu2016) and verified through findmitoscaf module in MitoZ (Meng et al. Reference Meng, Li, Yang and Liu2019). The complete circular mitogenome sequences were deposited to the GenBank.
Gene annotation and phylogenetic analyses
Following mitogenome assembly, protein-coding genes (PCGs), transfer and ribosomal RNAs, and non-coding regions were annotated using the MITOS2 web server (Bernt et al. Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013) with the invertebrate mitochondrial genetic code (genetic code 5). The missing genes were predicted based on sequence similarity to homologous genes from closely related species, including B. mucronatus, B. okinawaensis, and B. xylophilus. The AT-rich region was analyzed using Tandem Repeats Finder (Benson Reference Benson1999), characterising (1) region length, (2) repeat unit size, and (3) copy number within this region.
All 12 PCGs from each mitogenome were selected for phylogenetic reconstruction. Each PCG was individually aligned using TranslatorX (Abascal et al. Reference Abascal, Zardoya and Telford2010) under the invertebrate mitochondrial genetic code to ensure proper amino acid translation and codon alignment. The resulting alignments were then concatenated into a supermatrix, and optimal partitioning schemes with corresponding substitution models were determined using ModelFinder (Kalyaanamoorthy et al. Reference Kalyaanamoorthy, Bui, Wong, von Haeseler and Jermiin2017) implemented in IQ-TREE. We reconstructed the phylogeny through maximum likelihood (ML) analysis with a dataset comprising the complete mitogenomes of 35 nematode species (including our three newly sequenced Bursaphelenchus species). Thaumamermis cosgrovei, a phylogenetically distant isopod-parasitic nematode, was designated as the outgroup. The ML analysis was performed on the CIPRES Science Gateway (Miller et al. Reference Miller, Pfeiffer and Schwartz2010) using RAxML-HPC2 (Stamatakis Reference Stamatakis2014) with 1,000 bootstrap replicates to assess nodal support.
Results
Morphology characterisation
Three Bursaphelenchus species were observed morphologically with typical aphelenchid esophagus, tail terminal bursa, and two incisures in the lateral field. However, B. chengi (Figure 1A–E) is characterised by males having robust and strongly curved spicules with lamina distally angular, rostrum digitate, condylus rounded and cucullus broad and blunt, bursa truncate with posterior margin indented or rounded with a fine axial point, tail arcuate with a pointed terminus; females having a small vulval flap and conical tail gradually tapering to an almost straight or slightly recurved, pointed or rounded terminus. B. parantoniae (Figure 1F–J) is characterised by males having stout and arcuate spicules with lamina smoothly curved, condylus slightly elongate, truncated, or sometimes rounded with slight posterior recurvature, rostrum triangular with bluntly rounded tip, cucullus flattened, and bursa flap oval or irregularly shaped, possessing a rounded projection in the middle; and females having a small vulval flap and conical tail with a hyaline region, and terminus bluntly rounded or sometimes possessing a ventral depression distally. B. sinensis (Figure 1K–P) is characterised by males having a squared condylus ventrally bent, a rostrum not clearly differentiated, a ventral limb appearing anteriorly thickened and posteriorly weak and membranous, a distal end of the lamina truncate with a small thickening not expanded as a true cucullus, a tail arcuate with an obtuse or almost pointed talon-like terminus, and a small oval bursa with a small central appendage; females having a small vulval flap, a conical tail that is almost straight or slightly recurved ventrally, and a terminus that is almost pointed or rounded. Additionally, feeding apparatus dimorphism (Kanzaki et al. Reference Kanzaki, Ekino and Giblin-Davis2019) was observed – that is, the mycophagous form possessed a thin stylet with distinct basal knobs (Figure 1K), whereas the predatory form exhibited a flattened cephalic region and a longer stylet with a wider lumen, lacking basal knobs (Figure 1L).
Morphological characteristics of three Bursaphelenchus species. (A–E): Bursaphelenchus chengi; (F–J): B. parantoniae; (K–P): B. sinensis. A, F: Anterior body region; K, L: Head region; B, G, M: Male tail region; C, H, N: Male tail region in ventral view showing terminal bursa; D, I, O: Female vulval region; E, J, P: Female tail. (Scale bar = 10 μm).
Mitogenome characterisation
The complete mitogenomes of three Bursaphelenchus species – B. chengi (17,670 bp; GenBank accession number PV423453), B. parantoniae (15,021 bp; PV366404), and B. sinensis (18,386 bp; PV492208) – were successfully sequenced and assembled. All three mitogenomes comprise circular DNA molecules containing the typical nematode mitochondrial gene complement: 12 protein-coding genes (PCGs) including atp6, cob, cox1-3, nad1-6, and nad4L, 2 ribosomal RNA genes (rrnL and rrnS), and 22 transfer RNA (tRNA) genes. While the overall gene content was conserved, we observed notable specific rearrangement types in tRNA gene arrangements. Specifically, the trnV gene was located between trnP and nad6 in both B. chengi and B. parantoniae but positioned between trnE and rrnS in B. sinensis (Figure 2). This variation in tRNA gene order suggests potential evolutionary divergence within the genus.
Mitogenome structures of three Bursaphelenchus species. The protein-coding genes and RNA genes are indicated as colored blocks, with all genes encoded in the anticlockwise direction. Non-coding regions are shown as uncolored segments. Variations in trnV positions are highlighted with red boxes.
Furthermore, the mitochondrial DNA (mtDNA) of three Bursaphelenchus species exhibited a strong AT-rich bias. The mtDNA nucleotide composition of B. chengi comprised 32.5% A, 46.3% T (total A+T = 78.8%), 7.4% C, and 13.8% G. Similarly, the mtDNA nucleotide composition of B. parantoniae was 32.7% A, 50.2% T (total A+T = 82.9%), 5.9% C, and 11.2% G. In contrast, the mtDNA nucleotide composition of B. sinensis was 36.1% A, 52.3% T (total A+T = 88.4%), 3.9% C, and 7.7% G. Notably, we identified a large non-coding region between nad4 and cox1 that showed extreme AT-richness, significantly exceeding the genomic average. This region displayed substantial length variation among B. chengi (4,641 bp), B. parantoniae (1,962 bp), and B. sinensis (5,251 bp). In the non-coding region of B. chengi, the AT-rich region spanned 1,966 bp (genomic position: 13,588–15,553 bp), containing a tandem repeat structure with a 168 bp monomer unit repeated 11.7 times, with a nucleotide composition of 47.8% A and 48.7% T (total A+T = 96.5%). In B. parantoniae, the AT-rich region spanned 658 bp (position: 14,027–14,684 bp), containing a tandem repeat structure with a 187 bp monomer unit repeated 3.5 times, with a nucleotide composition of 49.9% A and 49.8% T (total A+T = 99.7%). In B. sinensis, the AT-rich region spanned 3,760 bp (position: 644–4,403 bp), containing a tandem repeat structure with a 271 bp monomer unit repeated 12.4 times, with nucleotide composition of 46.6% A and 49.4% T (total A+T = 96.0%).
The 12 protein-coding genes (PCGs) showed highly conserved lengths and genomic organization across all three species. In B. chengi, the PCGs spanned 10,156 bp in total length, accounting for 78.2% of the coding region, and the gene lengths ranged from 234 bp (nad4L) to 1,563 bp (cox1). In B. parantoniae, the PCGs covered 10,186 bp in total length, occupying 78.2% of the coding region, and the gene size range was consistent with B. chengi. In B. sinensis, the PCGs were slightly longer (10,206 bp total length, 79.0% of the coding region), and the gene lengths ranged from 240 bp (nad4L) to 1,569 bp (nad5 and cox1).
Most genes predominantly used ATT as the initiation codon, with ATA and TTG occurring less frequently. Most of the PCGs preferred using TAA as the termination codon and showed a marked preference for T-rich codons, while a small number of PCGs utilized TAG and single T as their alternative stop codons.
The rRNA genes displayed conserved locations. The large subunit ribosomal RNA (rrnL) was positioned between cox2 and nad3, while the small subunit ribosomal RNA (rrnS) was located between nad4l and nad1 (Table 1).
Mitochondrial genome organization in Bursaphelenchus chengi, B. parantoniae, and B. sinensis
All rRNA regions exhibited strong AT bias. For B. chengi, the A + T content of rrnS and rrnL were 76.0% (34.5% A, 41.5% T) and 79.0% (34.4% A, 44.6% T), respectively. For B. parantoniae, the A + T content of rrnS and rrnL were 82.7% (36.6% A, 46.1% T) and 86.4% (39.6% A, 46.8% T), respectively. Similarly, for B. sinensis, the A + T contents of rrnS and rrnL were 86.4% (37.7% A and 48.7% T) and 88.8% (38.1% A and 50.7% T), respectively.
Phylogeny and gene arrangement
The maximum likelihood phylogenetic analysis of B. chengi, B. parantoniae, and B. sinensis was conducted using translated amino acid sequences from their 12 protein-coding genes (PCGs). For comparative analysis, corresponding sequences from 32 additional nematode species (28 Chromadoreans and 4 Enopleans) obtained from GenBank were included (Figure 3). The phylogenetic reconstruction revealed that the newly sequenced species (B. chengi, B. parantoniae, and B. sinensis) formed a well-supported clade (BS = 100) with B. xylophilus and B. mucronatus. This clade showed a sister relationship (BS = 100) with a branch comprising Aphelenchoides species. Both these sister branches were robustly positioned within Aphelenchoididae family clade.
Maximum likelihood phylogenetic tree inferred from amino acid sequences of 12 concatenated protein-coding genes. Newly sequenced Bursaphelenchus species are highlighted in bold.
Gene arrangement analysis revealed a high degree of conservation in mitochondrial gene order across Bursaphelenchus species, with the exception of B. sinensis (Figure 4). A similar arrangement pattern was also observed in certain Aphelenchoides species – specifically, A. medicagus and A. besseyi exhibited identical gene orders. Notably, this conservation may extend beyond Aphelenchoididae to the suborder Rhabditina. For instance, representatives of the families Neodiplogasteridae (Pristionchus pacificus) and Rhabditidae (Caenorhabditis afar, C. elegans, and Cruznema tripartitum) exhibited identical protein-coding gene orders to Aphelenchoididae species despite some tRNA translocations.
Comparative linear mitochondrial gene arrangement of Bursaphelenchus species and related species in order Rhabditida. All protein-coding and RNA genes are shown in their relative positions, uniformly encoded in the same transcriptional direction. Non-coding regions have been omitted for clarity.
Discussion
The era of nematode mitogenomics began in 1990 with the landmark sequencing of C. elegans (Okimoto et al. Reference Okimoto, Macfarlane, Clary and Wolstenholme1992). Since then, advances in high-throughput sequencing technologies have enabled rapid expansion of nematode mitogenome data. Mitochondrial genomes have become invaluable molecular resources, offering multiple advantages for taxonomic identification, population genetics, and evolutionary studies due to their (1) maternal inheritance, (2) elevated mutation rates, (3) absence of introns, (4) limited recombination, and (5) conserved genomic structure (Boore Reference Boore1999; Cameron Reference Cameron2014). Despite these advances, current mitogenome representation remains disproportionately low compared to rRNA sequences, with only approximately 200 nematode mitogenomes available to date (Kern et al. Reference Kern, Kim and Park2020). Our study significantly expands this limited dataset by characterising three new Bursaphelenchus mitogenomes (B. chengi, B. parantoniae, and B. sinensis), thereby doubling the available mitogenomic resources for this ecologically important genus from three to six complete sequences. These novel mitogenomes provide essential genomic foundations for reconstructing the evolutionary history and phylogenetic relationships within the suborder Tylenchina.
Contemporary understanding of nematode phylogenetics relies on three primary molecular data types: (1) ribosomal gene sequences (Blaxter et al. Reference Blaxter, De Ley, Garey, Liu, Scheldeman, Vierstraete, Vanfleteren, Mackey, Dorris, Frisse, Vida and Thomas1998; De Ley and Blaxter Reference De Ley and Blaxter2004; Holterman et al. Reference Holterman, van der Wurff, van den Elsen, van Megen, Bongers, Holovachov, Bakker and Helder2006; van Megen et al. Reference van Megen, van den Elsen, Holterman, Karssen, Mooyman and Helder2009), (2) genomic/transcriptomic datasets (Ahmed et al. Reference Ahmed, Roberts, Adediran, Smythe, Kocot and Holovachov2022; Qing et al. Reference Qing, Zhang, Sun, Ahmed, Lo, Bert, Holovachov and Li2025; Smythe et al. Reference Smythe, Holovachov and Kocot2019), and (3) mitochondrial gene markers (Kern et al. Reference Kern, Kim and Park2020; Kim et al. Reference Kim, Lee, Kil and Park2020). The remarkable diversity in mitochondrial genome composition, structure, and content among nematodes suggests complex evolutionary dynamics involving multilevel selection pressures and accelerated mitogenomic evolution (Gendron et al. Reference Gendron, Qing, Sevigny, Li, Liu, Blaxter, Powers, Thomas and Porazinska2024). Mitochondrial gene order has proven particularly valuable for phylogenetic inference (Boore and Brown Reference Boore and Brown1998; Fritzsch et al. Reference Fritzsch, Schlegel and Stadler2006). Our analyses revealed a highly conserved gene arrangement pattern supporting the close evolutionary relationship between Bursaphelenchus and Aphelenchoides species, consistent with previous findings from both rRNA and mitogenome-based phylogenies (Holterman et al. Reference Holterman, van der Wurff, van den Elsen, van Megen, Bongers, Holovachov, Bakker and Helder2006; Kim et al. Reference Kim, Kim, Nadler and Park2016; Park et al. Reference Park, Sultana, Lee, Kang, Kim, Min, Eom and Nadler2011). Interestingly, we observed nearly identical gene orders across three families (Aphelenchoididae, Rhabditidae, and Neodiplogasteridae; Figure 4), suggesting strong evolutionary constraints in these lineages. In stark contrast, Enoplea exhibited substantial variation in mitochondrial gene order, with frequent rearrangements even among congeners. This pronounced variability limits the utility of gene order as a phylogenetic marker within Enoplea (Hyman et al. Reference Hyman, Lewis, Tang and Wu2011; Tang and Hyman Reference Tang and Hyman2007), highlighting fundamental differences in mitochondrial genome evolution between nematode classes.
Acknowledgements
This work was supported by the Key Research & Development Project of Ningxia Autonomous Region, China (2023BCF01022) and the Fundamental Research Funds for the Central Universities (KYPT2025001).
