# A Comparative Analyses of the Complete Mitochondrial Genomes of Fungal Endosymbionts inSogatella furcifera, White-Backed Planthoppers

Authors: Nak Jung Choi; Hong Xi; Jongsun Park
Journal: International Journal of Genomics
 (2021)

---

## Abstract

Sogatella furcifera Horvath, commonly known as the white-backed planthoppers (WBPH), is an important pest in East Asian rice fields. Fungal endosymbiosis is widespread among planthoppers in the infraorder Fulgoromorpha and suborder Auchenorrhyncha. We successfully obtained complete mitogenome of five WBPH fungal endosymbionts, belonging to the Ophiocordycipitaceae family, from next-generation sequencing (NGS) reads obtained from S. furcifera samples. These five mitogenomes range in length from 55,390 bp to 55,406 bp, which is shorter than the mitogenome of the fungal endosymbiont found in Ricania speculum, black planthoppers. Twenty-eight protein-coding genes (PCGs), 12 tRNAs, and 2 rRNAs were found in the mitogenomes. Two single-nucleotide polymorphisms, two insertions, and three deletions were identified among the five mitogenomes, which were fewer in number than those of four species of Ophiocordycipitaceae, Ophiocordyceps sinensis, Hirsutella thompsonii, Hirsutella rhossiliensis, and Tolypocladium inflatum. Noticeably short lengths (up to 18 bp) of simple sequence repeats were identified in the five WBPH fungal endosymbiont mitogenomes. Phylogenetic analysis based on conserved PCGs across 25 Ophiocordycipitaceae mitogenomes revealed that the five mitogenomes were clustered with that of R. speculum, forming an independent clade. In addition to providing the full mitogenome sequences, obtaining complete mitogenomes of WBPH endosymbionts can provide insights into their phylogenetic positions without needing to isolate the mtDNA from the host. This advantage is of value to future studies involving fungal endosymbiont mitogenomes.


---

## Body

## 1. Introduction

Sogatella furcifera Horvath commonly known as the white-backed planthopper (WBPH) is a planthopper belonging to the infraorder Fulgoromorpha [1] and suborder Auchenorrhyncha [2]. It has migrated to temperate climates from subtropical regions and become a major pest in rice fields across East Asia [3–6]. In particular, migration from China to Japan via Korean peninsula has highlighted the extent of its spread across the region [7]. Sogatella furcifera has already been registered in the National Species List of Korea [8] indicating that this species has been frequently found within the country. It damages rice plants by feeding directly on them, producing a characteristic symptom, hopper burn [9]. Because of the importance of WBPH as a threat to agriculture, the mitochondrial genome (mitogenome) as well as whole genome sequences of S. furcifera has been sequenced successfully [10, 11]. The fundamental background of WBPH genomic research is, therefore, well established. For example, the complete genome sequence of the Cardinium bacterial endosymbiont of S. furcifera was also completed from the same raw reads generated by the whole genome project [12]. Another bacterial endosymbiont of WBPH, Wolbachia, which alters host reproductions by parthenogenesis, feminization, male-killing, and induction of cytoplasmic incompatibility in arthropods [13], also causes the cytoplasmic incompatibility in WBPH together with Cardinium endosymbiont [14].Besides these bacterial endosymbionts, fungal endosymbiont has been identified using PCR method in planthopper,Ricania japonica [15]. This yeast-like endosymbiont uses the enzyme uricase to recycle uric acid secreted by the host species, assisting in metabolic processes [15]. In addition, yeast-like symbionts have been identified in Nilaparvata lugens, a brown planthopper [16, 17] which also support the host’s uric acid metabolism [18]. However, there was no sequence information of this endosymbiont until the complete fungal mitogenome was obtained from the raw reads of Ricania speculum, a black planthopper [19]. This mitogenome was identified as an Ophiocordycipitaceae species by comparing already known several complete mitogenomes in this family [19]. This result suggests that next-generation sequencing technology that provides a large number of short reads can be used to provide evidence for the existence of endosymbiont species using DNA extracted from insect species. These results draw comparison to previous studies that have successfully identified a multiple number of complete organelle or bacterial genomes from one NGS library [12, 19–37].Here, we reported the first complete mitogenomes of fungal WBPH endosymbiont from five WBPH samples isolated in Korea and China. The five mitogenomes display 55,390 to 55,406 bp in length, shorter than that ofR. speculum [19]. The numbers of intraspecific variations among the five mitogenomes are fewer in number than those of the four Ophiocordycipitaceae species. Phylogenetic analysis based on conserved PCGs across Ophiocordycipitaceae mitogenomes displays that the five mitogenomes were clustered with that of R. speculum, forming an independent clade. Once additional planthopper fungal endosymbiont mitogenomes become available, their phylogenetic relationships as well as evolutionary histories based on their complete mitogenomes will become clearer.


## 2. Materials and Methods

### 2.1. DNA Preparation and Genome Sequencing of Four WBPH Samples

All four WBPH samples were captured at two places in Korea (Table1). One individual of WBPH was frozen with liquid nitrogen using 1.5 ml microtube and then was ground using a plastic pestle. The Quick-DNA Miniprep Plus Kit (Zymo Research, USA) was used for extracting DNA. Genome sequencing was performed using NovaSeq6000 at Macrogen Inc., Korea, from the extracted DNA of four WBPH samples with constructing a 350 bp pair-end library.Table 1
List of five WBPH fungal endosymbiont mitogenomes.
No.NameSample locationNCBI accessionLength (bp)GC ratio (%)No. of PCGsNo. of tRNAsNo. of rRNAsReference1KR321-2, Daesong-ri, Geumnam-myeon, Hadong-gun, Gyeongsangnam-do, Republic of KoreaMW11513155,39330.728162This study2KR.1D1411-6, Wolga-ri, Gunnae-myeon, Jindo-gun, Jeollanam-do, Republic of KoreaMW37371055,40630.728162This study3KR.5DMW37371155,39030.728162This study4KR.11DMW37686255,39330.728162This study5WGSUniversity of Science and Technology of China, Anhui province, ChinaBK05918655,39330.728162SRR3954848 [11]


### 2.2. Assembly and Annotation of the Five Fungal WBPH Endosymbiont Mitogenomes

De novo assembly, with confirmation, was accomplished with Velvet v1.2.10 [38] after filtering raw reads using Trimmomatic v0.33 [39]. After obtaining mitogenome contig sequences with the condition that sequence coverage is more than 60x, gaps were filled with GapCloser v1.12 [40], and all bases from the assembled sequences were confirmed by checking each base in the alignment (tview mode in SAMtools v1.9 [41]) against the assembled mitogenome generated with BWA v0.7.17 [42]. The circular form of mitogenomes was confirmed by the pair-end reads connecting both sides of mitogenomes. All these bioinformatic analyses were conducted under the environment of the Genome Information System (GeIS; http://geis.infoboss.co.kr/) like the previous studies of mitogenomes [19, 21–24, 26, 28, 30, 32, 33, 36, 43–91].Geneious Prime® 2020.2.4 (Biomatters Ltd, Auckland, New Zealand) was used for mitogenome annotation with referring to the mitogenome ofR. speculum fungal endosymbiont (NC_049089) [19] by transferring annotations while correcting exceptional cases, including missing start or stop codons. Also, the “FindORF” function in Geneious Prime® 2020.2.4 together with BLAST v2.2.24 [92] was also utilized to find novel PCGs including LAGLIDADG endonucleases. tRNAs were predicted and confirmed using tRNAScan-SE v2.0.6 [93].


### 2.3. Identification of Sequence Variations on the Complete Mitogenomes of WBPH Fungal Endosymbionts

Single-nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs) were identified using the “Find variations/SNP” function implemented in Geneious Prime® 2020.2.4 (Biomatters Ltd, Auckland, New Zealand) based on multiple sequence alignment of the five mitogenomes of WBPH fungal endosymbionts conducted by MAFFT v7.450 [94]. Each identified variation was manually checked to understand which mitogenome has them.


### 2.4. Identification of Simple Sequence Repeats (SSRs)

Simple sequence repeats (SSRs) were identified on the mitogenome sequence using the pipeline of the SSR database (SSRDB;http://ssr.pe.kr/; Park et al., in preparation). Based on the conventional definition of an SSR on an organelle genome, monoSSR (1 bp) to hexaSSR (6 bp), the total length of SSRs on the mitogenome exceeds 10 bp. Owing to the different criteria of SSRs on organelle genomes [95–101], we adopted the criteria used in various organelle genome analyses [21, 44, 102–104], as follows: the monoSSR (unit sequence length of 1 bp) to hexaSSR (6 bp) are used as normal SSRs, and heptaSSR (7 bp) to decaSSR (10 bp) are defined as extended SSRs. Among the normal SSRs, pentaSSRs and hexaSSRs for which the number of unit sequences is 2 are classified as potential SSRs.


### 2.5. Construction of Phylogenetic Trees

Five conserved PCGs, includingATP8, CO2, NAD3, NAD4, and NAD4L, from 26 fungal mitogenomes including the five mitogenomes assembled in this study and one outgroup species, Fusarium graminearum, were aligned independently using MAFFT v7.450 [94] and concatenated using the Perl script, one of the component of GenomeArchive® (http://www.genomearchive.info) [105]. The model test was conducted with jModelTest v2.1.5 [106]. The neighbor-joining (NJ) and maximum-likelihood (ML) trees were reconstructed in MEGA X [107]. In the ML analysis, a heuristic search was used with nearest-neighbor interchange (NNI) branch swapping, general time reversible (GTR) model, and uniform rates among sites. All other options used the default settings. Bootstrap analyses with 1000 pseudoreplicates were conducted with the same options. The posterior probability of each node was estimated by Bayesian inference (BI) using the MrBayes v3.2.7 [108] plug-in implemented in Geneious Prime® 2020.2.4. The HKY85 model with gamma rates was used as a molecular model. A Markov chain Monte Carlo (MCMC) algorithm was employed for 1,100,000 generations, sampling trees every 200 generations, with four chains running simultaneously. Trees from the first 100,000 generations were discarded as burn-in.


## 2.1. DNA Preparation and Genome Sequencing of Four WBPH Samples

All four WBPH samples were captured at two places in Korea (Table1). One individual of WBPH was frozen with liquid nitrogen using 1.5 ml microtube and then was ground using a plastic pestle. The Quick-DNA Miniprep Plus Kit (Zymo Research, USA) was used for extracting DNA. Genome sequencing was performed using NovaSeq6000 at Macrogen Inc., Korea, from the extracted DNA of four WBPH samples with constructing a 350 bp pair-end library.Table 1
List of five WBPH fungal endosymbiont mitogenomes.
No.NameSample locationNCBI accessionLength (bp)GC ratio (%)No. of PCGsNo. of tRNAsNo. of rRNAsReference1KR321-2, Daesong-ri, Geumnam-myeon, Hadong-gun, Gyeongsangnam-do, Republic of KoreaMW11513155,39330.728162This study2KR.1D1411-6, Wolga-ri, Gunnae-myeon, Jindo-gun, Jeollanam-do, Republic of KoreaMW37371055,40630.728162This study3KR.5DMW37371155,39030.728162This study4KR.11DMW37686255,39330.728162This study5WGSUniversity of Science and Technology of China, Anhui province, ChinaBK05918655,39330.728162SRR3954848 [11]


## 2.2. Assembly and Annotation of the Five Fungal WBPH Endosymbiont Mitogenomes

De novo assembly, with confirmation, was accomplished with Velvet v1.2.10 [38] after filtering raw reads using Trimmomatic v0.33 [39]. After obtaining mitogenome contig sequences with the condition that sequence coverage is more than 60x, gaps were filled with GapCloser v1.12 [40], and all bases from the assembled sequences were confirmed by checking each base in the alignment (tview mode in SAMtools v1.9 [41]) against the assembled mitogenome generated with BWA v0.7.17 [42]. The circular form of mitogenomes was confirmed by the pair-end reads connecting both sides of mitogenomes. All these bioinformatic analyses were conducted under the environment of the Genome Information System (GeIS; http://geis.infoboss.co.kr/) like the previous studies of mitogenomes [19, 21–24, 26, 28, 30, 32, 33, 36, 43–91].Geneious Prime® 2020.2.4 (Biomatters Ltd, Auckland, New Zealand) was used for mitogenome annotation with referring to the mitogenome ofR. speculum fungal endosymbiont (NC_049089) [19] by transferring annotations while correcting exceptional cases, including missing start or stop codons. Also, the “FindORF” function in Geneious Prime® 2020.2.4 together with BLAST v2.2.24 [92] was also utilized to find novel PCGs including LAGLIDADG endonucleases. tRNAs were predicted and confirmed using tRNAScan-SE v2.0.6 [93].


## 2.3. Identification of Sequence Variations on the Complete Mitogenomes of WBPH Fungal Endosymbionts

Single-nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs) were identified using the “Find variations/SNP” function implemented in Geneious Prime® 2020.2.4 (Biomatters Ltd, Auckland, New Zealand) based on multiple sequence alignment of the five mitogenomes of WBPH fungal endosymbionts conducted by MAFFT v7.450 [94]. Each identified variation was manually checked to understand which mitogenome has them.


## 2.4. Identification of Simple Sequence Repeats (SSRs)

Simple sequence repeats (SSRs) were identified on the mitogenome sequence using the pipeline of the SSR database (SSRDB;http://ssr.pe.kr/; Park et al., in preparation). Based on the conventional definition of an SSR on an organelle genome, monoSSR (1 bp) to hexaSSR (6 bp), the total length of SSRs on the mitogenome exceeds 10 bp. Owing to the different criteria of SSRs on organelle genomes [95–101], we adopted the criteria used in various organelle genome analyses [21, 44, 102–104], as follows: the monoSSR (unit sequence length of 1 bp) to hexaSSR (6 bp) are used as normal SSRs, and heptaSSR (7 bp) to decaSSR (10 bp) are defined as extended SSRs. Among the normal SSRs, pentaSSRs and hexaSSRs for which the number of unit sequences is 2 are classified as potential SSRs.


## 2.5. Construction of Phylogenetic Trees

Five conserved PCGs, includingATP8, CO2, NAD3, NAD4, and NAD4L, from 26 fungal mitogenomes including the five mitogenomes assembled in this study and one outgroup species, Fusarium graminearum, were aligned independently using MAFFT v7.450 [94] and concatenated using the Perl script, one of the component of GenomeArchive® (http://www.genomearchive.info) [105]. The model test was conducted with jModelTest v2.1.5 [106]. The neighbor-joining (NJ) and maximum-likelihood (ML) trees were reconstructed in MEGA X [107]. In the ML analysis, a heuristic search was used with nearest-neighbor interchange (NNI) branch swapping, general time reversible (GTR) model, and uniform rates among sites. All other options used the default settings. Bootstrap analyses with 1000 pseudoreplicates were conducted with the same options. The posterior probability of each node was estimated by Bayesian inference (BI) using the MrBayes v3.2.7 [108] plug-in implemented in Geneious Prime® 2020.2.4. The HKY85 model with gamma rates was used as a molecular model. A Markov chain Monte Carlo (MCMC) algorithm was employed for 1,100,000 generations, sampling trees every 200 generations, with four chains running simultaneously. Trees from the first 100,000 generations were discarded as burn-in.


## 3. Results and Discussions

### 3.1. Complete Mitogenome of Fungal WBPH Endosymbionts

We successfully assembled fungal endosymbiont mitogenomes from four WBPH samples isolated in Korea and China and one public dataset of NGS raw reads (Table1). This is the first WBPH fungal endosymbiont mitogenome identified. Their lengths ranged from 55,390 bp to 55,406 bp (Table 1), which is shorter than that of R. speculum (66,785 bp) [19]. In these mitogenomes, there were 28 protein-coding genes (PCGs), 12 tRNAs, and 2 rRNAs (Table 2). Some of the PCGs found were LAGLIDADG endonucleases, which are usually found in intronic regions of various fungal mitogenomes, contributing to the expansion of their length [19, 55, 108–112]. In comparison to the previously sequenced mitogenome of the fungal endosymbiont of R. speculum, there were slightly fewer PCGs and tRNAs found in the WBPH endosymbiont mitogenomes. There were three fewer PCGs for three reasons: the smaller number of LAGLIDADG endonucleases, the absence of one endonuclease and a GIY-YIG endonuclease, and the presence of two additional PCGs—a hypothetical protein and a LAGLIDADG/HNH endonuclease. This particular configuration of PCGs is usually identified in other fungal mitogenomes; for example, two mitogenomes of Fusarium oxysporum (GenBank accessions are MN259514 and MN259515) display two completely different PCGs in each mitogenome [54, 56]. There are also five fewer tRNAs because of the different configurations: tRNA-Asp, tRNA-Cys, tRNA-Ile, and two tRNA-Ser (also found in the mitogenome of the fungal symbiont of R. speculum [19]). This difference in configuration of tRNAs between two different fungal symbionts suggests that tRNA configuration may not be critical because essential tRNAs absent in the fungal mitogenome can be supported from the nuclear genome [113].Table 2
List of genes in WBPH fungal endosymbiont mitogenomes.
No.NameTypeStart positionEnd positionLength (bp)StrandNo. of exons1Cytochrome bCDS1608182975Reverse52LAGLIDADG/HNH endonucleaseCDS74297776348Reverse13NADH dehydrogenase subunit 5CDS843210,4051974Reverse14NADH dehydrogenase subunit 4 LCDS10,40510,680276Reverse15Cytochrome c oxidase subunit IICDS10,89613,952894Reverse66ATP synthase F0 subunit cCDS14,05214,237186Reverse17NADH dehydrogenase subunit 3CDS14,31514,968447Reverse28LAGLIDADG endonucleaseCDS14,68314,880198Reverse19NADH dehydrogenase subunit 2CDS14,96918,1931431Reverse210Hypothetical proteinCDS15,62616,390765Reverse111tRNA-MettRNA18,21018,28374Reverse112tRNA-HistRNA18,32518,39874Reverse113tRNA-LeutRNA18,48518,56985Reverse114tRNA-LystRNA18,57018,64273Reverse115tRNA-PhetRNA18,64318,71573Reverse116tRNA-LeutRNA19,67319,75583Reverse117tRNA-MettRNA19,75819,83073Reverse118tRNA-GlutRNA19,90419,97673Reverse119Large subunit rRNArRNA20,06727,3374793Reverse320Ribosomal protein S3CDS20,67221,9161245Reverse121tRNA-ProtRNA27,40027,47172Reverse122NADH dehydrogenase subunit 6CDS27,67728,441765Reverse123tRNA-GlytRNA28,54428,61471Reverse124Cytochrome c oxidase subunit IIICDS29,22432,387795Reverse325Intron-encoded endonuclease aI1CDS29,93230,138207Reverse126Hypothetical proteinCDS30,39530,664270Reverse127Alpha-beta-hydrolaseCDS30,42330,686264Forward128tRNA-AsntRNA32,43732,50872Reverse129Small subunit rRNArRNA34,25435,6151362Reverse130ATP synthase F0 subunit aCDS35,91838,001762Reverse231ATP synthase F0 subunit 8CDS38,06138,237177Reverse132NADH dehydrogenase subunit 4CDS38,62040,0351416Reverse133NADH dehydrogenase subunit 1CDS40,11942,5821107Reverse434tRNA-ArgtRNA42,72642,79974Reverse135Cytochrome c oxidase subunit ICDS42,80955,2931314Reverse736LAGLIDADG endonucleaseCDS43,66044,8801221Reverse137LAGLIDADG endonucleaseCDS45,12546,054930Reverse138LAGLIDADG endonucleaseCDS46,85748,2481392Reverse139LAGLIDADG endonucleaseCDS48,62049,483864Reverse140LAGLIDADG endonucleaseCDS49,52550,8141290Reverse141LAGLIDADG endonucleaseCDS50,86551,734870Reverse142Hypothetical proteinCDS54,29054,850561Forward1Several PCGs in the fungal mitogenomes have been invaded by introns multiple times. For example,COX1 contains three introns, and COB has five introns in the Hirsutella thompsonii mitogenome [114]. This phenomenon contributes to increased fungal mitogenome: Aspergillus pseudoglaucus and Aspergillus egyptiacus are longer than the other Aspergillus mitogenomes because of the presence of many introns on major PCGs [55, 115]. The fungal mitogenomes examined in this study also present many introns on PCGs including COB, COX1, NAD1, ATP8, COX3, COX2, and NAD2 (Figure 1), which is a major reason for the expansion of fungal mitogenomes together with endonucleases.Figure 1
Complete mitogenome of WBPH fungal endosymbionts. Blue circle indicates fungal endosymbiont mitogenome, yellow arrows are protein-coding genes, purple arrows are tRNAs, and red arrows mean rRNAs. Direction of arrows indicates direction of transcription. Each gene name was displayed with lines directing to the corresponded arrows. Numbers displayed outside of blue circle mean base pair.The gene order of WBPH andR. speculum fungal symbiont mitogenomes was the same when PCGs except endonucleases and rRNAs are considered. However, intron structures of COX1, COX2, NAD2, NAD3, NAD5, and ATP synthase F0 subunit present different configurations between the two mitogenomes (Figure 2). The intron structures of NAD5 and NAD2 present reduce of a reduction in the number of exons via removal of intron regions in the WBPH fungal endosymbiont mitogenome (Figures 2(a) and 2(d)), whereas those of COX2, NAD3, and the ATP synthase F0 subunit display insertions of one intron into the WBPH fungal endosymbiont mitogenome (Figures 2(b), 2(c), and 2(e)). This indicates that the reduction in the total length of the WBPH fungal symbiont mitogenome is not primarily caused by reducing the number of exons, unlike in Aspergillus mitogenomes [55, 116]. In addition, COX1, which contains the largest number of exons in these mitogenomes, lost the sixth and seventh exons of the R. speculum fungal endosymbiont mitogenome in the mitogenome of WBPH endosymbiont (Figure 2(f)). However, the total length of COX1, including the introns of WBPH fungal endosymbionts, is longer than that of R. speculum fungal endosymbionts by 1 kb (Figure 2(f)), reflecting complex events that occurred during the evolution of both mitogenomes. Additional studies are required to identify the correct exons of the COX1 gene of this fungal endosymbiont. For example, alignment of RNA-Seq raw reads against this mitogenome could provide expressed regions in this mitogenome.Figure 2
Schematic diagram of exon/intron structure of protein-coding genes displaying different configurations between fungal endosymbiont mitogenomes of WBPH andR. speculum. Orange-colored diagrams indicate components of WBPH fungal endosymbionts, and blue-colored diagrams are those of R. speculum fungal endosymbionts. Exon/intron structures of (a) NAD5, (b) COX2, (c) NAD3, (d) NAD2, (e) ATP synthase F0 subunit, and (f) COX1 are displayed.
(a)(b)(c)(d)(e)(f)Once more fungal symbiont mitogenomes are available, patterns of presence and absence of tRNAs, additional endonucleases, and intron structures of PCGs in endosymbiont mitogenomes will elucidate a detailed evolutionary history of these genes.


### 3.2. Identification of Intraspecific Variations on Fungal WBPH Endosymbiont Mitogenomes

We identified two SNPs, three insertions, and two deletions via multiple sequence alignments of the five fungal mitogenomes (Table3). One of two SNPs was identified in KR.5D WBPH and changed leucine (L) to glutamine (Q) in the ATP synthase F0 subunit (Table 3). One 10 bp insertion in the intergenic space was found in KR.1D WBPH, while the remaining two insertions and all three deletions were 1 to 3 bp in length (Table 3).Table 3
List of intraspecific variations identified on the five WBPH fungal endosymbiont mitogenomes.
No.TypeCoordination of multiple sequence alignmentsStrainsBase changesPosition1Insertion4209-4210KR.1D, KR.5D- to CCIntergenic2Insertion27,476-27,486KR.1D- to TGGGCCCCCCIntergenic3SNP27,487KR.1DA to CIntergenic4Deletion32,727-32,728KR.5DCC to -Intergenic5SNP37,574KR.5DA to TL to Q in ATP synthase F0 subunit6Insertion38,727KR.1D- to GIntergenic7Deletion38,728-38,730KR.5DGGG to -IntergenicThe proportions of these intraspecific SNPs, insertions, and deletions in these fungal mitogenomes were 0.0036%, 0.020%, and 0.012%, respectively. The proportion of insertions and deletions was higher than that of SNPs. Interestingly, there is geographical variation in the fungal symbiont mitogenomes. The mitogenome of WBPH endosymbionts used in the whole genome sequencing (WGS) and the KR.11D isolate were identical to that of KR, while the other three WBPHs captured in other locations in Korea displayed intraspecific variations. The sample used in the WGS originated from the University of Science and Technology of China (Anhui province, China), indicating that KR 11D and KR WBPH samples obtained in Korea have migrated from the similar region to the WGS sample. However, further analyses of their complete mitogenomes or whole genomes will be needed to provide more supportive data for identifying their origins.There is a relatively small number of intraspecific SNPs and INDELs identified from these fungal mitogenomes in comparison to those of other fungal mitogenomes, for example, 16 to 17 SNPs (0.055% to 0.0582%) and 22 to 27 INDELs (0.075% to 0.092%) onAspergillus flavus [52, 53] and 62 SNPs (0.15%) and 181 INDELs (0.43%) on Fusarium oxysporum f.sp. lactucae [56]. They are also fewer than those identified in insect mitogenomes [10, 22, 23, 43, 45–51].Based on 25 available complete fungal mitogenomes in Ophiocordycipitaceae, four species,Ophiocordyceps sinensis, Hirsutella thompsonii, Hirsutella rhossiliensis, and Tolypocladium inflatum, contain more than one complete fungal mitogenome (Table 4). We investigated intraspecific variations in the mitogenomes of these four species (Table 5). There are significantly more INDELs than SNPs identified in the four fungal species, a trend identical to that observed in the four mitogenomes of fungal endosymbiont WBPH with the exception of their absolute amounts. Moreover, there were at least three times more SNPs and INDELs in these fungal mitogenomes than that in the fungal symbiont of WBPHs. This phenomenon can be explained by two major factors: first, the geographical distribution or genetic background of WBPH samples is relatively limited in comparison to those of the four fungal species, and second, the surroundings of fungal endosymbionts are less dynamic than those of normal fungal species, causing low selection pressure from the environment. This second factor is supported by two studies: first, the bacterial genome of aphid endosymbiont Buchnera aphidicola (Aphis gossypii) displays a low level of intraspecific variation in comparison to those of host mitogenome (Bae et al., under revision), and second, the whole genome of endosymbiont of Pediculus humanus capitis also shows low-level intraspecific variations in comparison to those of their whole genomes [117].Table 4
List of available complete fungal mitogenomes in Ophiocordycipitaceae.
No.SpeciesNCBI accessionLength (bp)GC ratio (%)No. of PCGsNo. of tRNAsNo. of rRNAsReference1Ophiocordyceps sinensisNC_034659157,53930.288272Unpub.2Ophiocordyceps sinensisMH400233157,58430.276272Unpub.3Ophiocordyceps sinensisKP835313157,51030.26∗232[123]4Hirsutella thompsoniiNC_04016562,50929.830272[114]5Hirsutella thompsoniiMH36729665,33230.332272[114]6Hirsutella thompsoniiMH36729560,36230.029272[114]7Hirsutella rhossiliensisMG97907162,94928.333262[129]8Hirsutella rhossiliensisNC_03016462,48328.224262Unpub.9Hirsutella vermicolaNC_03661053,79325.327252[130]10Hirsutella minnesotensisNC_02766052,24528.430252[131]11Tolypocladium sp.MN58326546,46626.115262[132]12Tolypocladium guangdongenseMT47126746,10226.130272[133]13Tolypocladium ophioglossoidesNC_03138435,15927.519252[134]14Tolypocladium cylindrosporumNC_04683934,69827,024262[135]15Tolypocladium inflatumNC_03638225,32827.815252[136]16Tolypocladium inflatumKY92488025,23827.815252[136]17Tolypocladium inflatumKY92488125,32827.815252[136]18Tolypocladium inflatumKY92488225,32827.815252[136]19Tolypocladium inflatumKY92488324,79327.815252[136]20Ophiocordycipitaceae sp.NC_04908966,78530.631172[19]∗Mitogenome annotation of this genome (KP835313) seems not to be complete because several major genes, such as COX1, NAD1, NAD5, and COB, which have many introns on fungal mitogenomes that were not annotated as CDS.Table 5
Number of intraspecific variations identified from four fungal species belonging to Ophiocordycipitaceae.
No.SpeciesNo. of mitogenomesAligned length (bp)No. of SNPsSNP coverage (%)No. of INDELsINDEL coverage (%)1Ophiocordyceps sinensis3157,606160.0101440.0912Hirsutella thompsonii366,6352810.4264899.743Hirsutella rhossiliensis264,85870.0120083.104Tolypocladium inflatum525,338300.123751.48


### 3.3. Identification and Comparative Analysis of Simple Sequence Repeats on the Five WBPH Fungal Endosymbiont Mitogenomes

Simple sequence repeats (SSRs) identified from organellar genomes have been utilized as molecular markers in various species such as plant species [99, 118–122], suggesting that SSRs on fungal endosymbiont mitogenomes can be used as molecular markers to identify the geographical origins of WBPH. In total, 23 normal and 6 extended SSRs were identified from fungal endosymbiont mitogenomes (Figure 3(b)), with the exception of the fungal endosymbiont mitogenome of WBPH KR.1D which displays 24 normal and 6 extended SSRs (Table 6). The fungal endosymbiont mitogenome of WBPH KR.1D has one more monoSSR (Table 6) with a unit sequence of C and length of 15 bp caused by one insertion (Table 3). In addition, 140 potential SSRs were also identified in the five mitogenomes (Table 6). SSRs identified in the mitogenome were distributed evenly (Figure 3(a)), suggesting that there was no hot spot of SSRs in these fungal mitogenomes.Figure 3
Numbers and distribution of SSRs on the mitogenome of fungal endosymbiont of WBPH KR. (a) Black circle indicates fungal endosymbiont mitogenome, yellow arrows are protein-coding genes, purple arrows are tRNAs, and red arrows mean rRNAs. SSRs, extended SSRs, and potential SSRs are displayed with yellow, green, and light green colors, respectively. (b) The number of SSRs along with SSR types is displayed. Blue color indicates SSRs and extended SSRs, and orange color means potential SSRs.Table 6
List of SSRs identified from the five fungal mitogenomes of WBPH endosymbionts.
SSR typeKRKR.1DKR.5DKR.11DWGSMonoSSR89888DiSSR88888TriSSR22222TetraSSR44444PentaSSR00000HexaSSR11111HeptaSSR55555OctaSSR11111NonaSSR00000DecaSSR00000Subtotal2930292929PentaPotentialSSR101101101101101HexaPotentialSSR3939393939Subtotal140140140140140The length of the identified SSRs is relatively short (a maximum length of 18 bp; Figure4(a)) in comparison to those of other fungal species in the same family: Ophiocordyceps sinensis (up to 24 bp) [123], as well as fungal species in the other families, such as Pestalotiopsis fici (up to 45 bp) [124]. Moreover, the maximum length of SSRs identified from the mitogenome of R. speculum (NC_049089) [19] was 18 bp, suggesting that this short SSR length can be linked to the evolution of endosymbiont mitogenomes.Figure 4
Comparative analysis of simple sequence repeats on mitogenome of WBPH fungal endosymbiont. (a) Number of SSRs, extended SSRs, and potential SSRs along with its length.x-axis indicates SSR length (bp), and y-axis means the number of SSRs. Blue-colored bars indicate SSRs and extended SSRs, and orange bars mean potential SSRs. (b) Number of SSRs, extended SSRs, and potential SSRs based on positions, genic, intronic ORF, intergenic, and intronic. (c) Number of genic SSRs along with genes (blue bars) as well as length of genes (orange line). x-axis is genes containing SSRs, and y-axis indicates the number of SSRs (left) and gene length (bp; right).
(a)(b)(c)Out of 191 normal SSRs, extended SSRs, and potential SSRs, 84 (43.98%) are located in the genic region (genic and intronic ORF categories in Figure4(b); Table 7). The intronic ORF position indicates the location of the PCGs placed at the introns of other PCGs, most of which are LAGLIDADG endonucleases (Table 2). Nearly half of the SSRs are in PCGs, which are conserved in comparison to intron and intergenic regions, indicating that these SSRs can be utilized for distinguishing species level or even higher rank. In the intergenic region, there were 61 SSRs (31.94%), and in comparison, only 24 SSRs (12.57%) in the intergenic region (Figure 4(b); Table 7). These SSRs are located in relatively nonconserved regions in comparison to PCG regions, suggesting that these SSRs can be used to distinguish intraspecific levels, such as population or geographical origins. Once more endosymbiont mitogenomes are available in the near future, these SSRs can be evaluated for their use in identification of species and their geographical origin as well as evolutionary history of their mitogenomes.Table 7
List of SSRs identified on fungal endosymbiont mitogenome of WBPH KR.
No.NameSSR typeTypeStartEndUnit sequenceRepeat numberGenes1M0000001Normal SSRMonoSSR42094219C11(Intron)Cob2M0000002Normal SSRMonoSSR55485558A11Cob3M0000003Normal SSRMonoSSR2231422323T10Large subunit ribosomal RNA4M0000004Normal SSRMonoSSR2248222492T11Large subunit ribosomal RNA5M0000005Normal SSRMonoSSR2844228451T106M0000006Normal SSRMonoSSR3271532728C147M0000007Normal SSRMonoSSR3831538327G138M0000008Normal SSRMonoSSR4003540044T10NAD49D0000001Normal SSRDiSSR57885797TA5(Intron)Cob10D0000002Normal SSRDiSSR1370813721AT7(Intron)COX211D0000003Normal SSRDiSSR2006820077AT5Large subunit ribosomal RNA12D0000004Normal SSRDiSSR2064620657AT6(Intron)large subunit ribosomal RNA13D0000005Normal SSRDiSSR2726627275TA5Large subunit ribosomal RNA14D0000006Normal SSRDiSSR3799538004TA5ATP synthase F0 subunit a15D0000007Normal SSRDiSSR4043040441AT6NAD116D0000008Normal SSRDiSSR5535255361TA517T0000001Normal SSRTriSSR3582935840ATT418T0000002Normal SSRTriSSR5278052791ATA4(Intron)COX119Te0000001Normal SSRTetraSSR2472124732ATTT3Large subunit ribosomal RNA20Te0000002Normal SSRTetraSSR4259842609TTTA321Te0000003Normal SSRTetraSSR4847748488ATAA3(Intron)COX122Te0000004Normal SSRTetraSSR5270952720AATA3(Intron)COX123P0000001Potential SSRPentaSSR586595TTGT2(Intron)Cob24P0000002Potential SSRPentaSSR19231932TAATA2(Intron)Cob25P0000003Potential SSRPentaSSR35233532TAAAA2(Intron)Cob26P0000004Potential SSRPentaSSR41164125TTGTC2(Intron)Cob27P0000005Potential SSRPentaSSR50685077ATAAT2(Intron)Cob28P0000006Potential SSRPentaSSR65406549TAATG2(Intron)Cob29P0000007Potential SSRPentaSSR66466655ATTTT2(Intron)Cob30P0000008Potential SSRPentaSSR67146723TTTT2(Intron)Cob31P0000009Potential SSRPentaSSR76367645AGCAA2LAGLIDADG/HNH endonuclease, (Intron)Cob32P0000010Potential SSRPentaSSR91099118AAGTT2NAD533P0000011Potential SSRPentaSSR92649273ATAA2NAD534P0000012Potential SSRPentaSSR1030510314AGACA2NAD535P0000013Potential SSRPentaSSR1085710866ATTCA236P0000014Potential SSRPentaSSR1136911378AGATA2COX237P0000015Potential SSRPentaSSR1251112520TTATA2(Intron)COX238P0000016Potential SSRPentaSSR1260012609TAAGA2(Intron)COX239P0000017Potential SSRPentaSSR1271012719AAGCG2(Intron)COX240P0000018Potential SSRPentaSSR1279612805TTAAC2(Intron)COX241P0000019Potential SSRPentaSSR1330313312TAATA2COX242P0000020Potential SSRPentaSSR1494614955AAAAG2NAD343P0000021Potential SSRPentaSSR1674716756TCGAG2(Intron)NAD244P0000022Potential SSRPentaSSR1738317392TCATT2NAD245P0000023Potential SSRPentaSSR1742717436AATAA2NAD246P0000024Potential SSRPentaSSR1751417523AAATG2NAD247P0000025Potential SSRPentaSSR1768517694AATAA2NAD248P0000026Potential SSRPentaSSR1808318092AATAC2NAD249P0000027Potential SSRPentaSSR1810918118TATT2NAD250P0000028Potential SSRPentaSSR1812118130ATAGA2NAD251P0000029Potential SSRPentaSSR1840018409TTATG252P0000030Potential SSRPentaSSR1881918828GATA253P0000031Potential SSRPentaSSR1907819087ATTTT254P0000032Potential SSRPentaSSR1919019199TTGTA255P0000033Potential SSRPentaSSR1930219311ATAAT256P0000034Potential SSRPentaSSR1961319622AACT257P0000035Potential SSRPentaSSR1962919638TATT258P0000036Potential SSRPentaSSR1990419913TAGAC2tRNA-Glu59P0000037Potential SSRPentaSSR2038420393TTATT2Large subunit ribosomal RNA60P0000038Potential SSRPentaSSR2096820977TTATT2RPS3, (Intron)large subunit ribosomal RNA61P0000039Potential SSRPentaSSR2114621155TGTAT2RPS3, (Intron)large subunit ribosomal RNA62P0000040Potential SSRPentaSSR2117321182TATTA2RPS3, (Intron)large subunit ribosomal RNA63P0000041Potential SSRPentaSSR2173021739TTATT2RPS3, (Intron)large subunit ribosomal RNA64P0000042Potential SSRPentaSSR2203822047TTTTA2(Intron)large subunit ribosomal RNA65P0000043Potential SSRPentaSSR2212122130TTATT2(Intron)large subunit ribosomal RNA66P0000044Potential SSRPentaSSR2260622615TAATA2Large subunit ribosomal RNA67P0000045Potential SSRPentaSSR2351223521AAGAC2Large subunit ribosomal RNA68P0000046Potential SSRPentaSSR2396223971TTTTC2Large subunit ribosomal RNA69P0000047Potential SSRPentaSSR2445324462AATTA2Large subunit ribosomal RNA70P0000048Potential SSRPentaSSR2450124510ATTTA2Large subunit ribosomal RNA71P0000049Potential SSRPentaSSR2521125220TTTAC2Large subunit ribosomal RNA72P0000050Potential SSRPentaSSR2535725366TTTT2Large subunit ribosomal RNA73P0000051Potential SSRPentaSSR2618126190CATTT2(Intron)large subunit ribosomal RNA74P0000052Potential SSRPentaSSR2722427233ATTTC2Large subunit ribosomal RNA75P0000053Potential SSRPentaSSR2770027709TTAAG2NAD676P0000054Potential SSRPentaSSR2784427853ATAAT2NAD677P0000055Potential SSRPentaSSR2801328022TAAAA2NAD678P0000056Potential SSRPentaSSR2911929128TCCCC279P0000057Potential SSRPentaSSR2927329282CAGTA2COX380P0000058Potential SSRPentaSSR2997329982TGAT2Intron-encoded endonuclease aI1, (Intron)COX381P0000059Potential SSRPentaSSR3086030869AGTG2(Intron)COX382P0000060Potential SSRPentaSSR3261832627TCCCC283P0000061Potential SSRPentaSSR3339733406TAAAT284P0000062Potential SSRPentaSSR3341633425ATGGT285P0000063Potential SSRPentaSSR3391633925AGAGA286P0000064Potential SSRPentaSSR3484234851AATT2Small subunit ribosomal RNA87P0000065Potential SSRPentaSSR3668036689TTAAA2(Intron)ATP synthase F0 subunit a88P0000066Potential SSRPentaSSR3699337002TTAAA2(Intron)ATP synthase F0 subunit a89P0000067Potential SSRPentaSSR3702137030ATTTT2(Intron)ATP synthase F0 subunit a90P0000068Potential SSRPentaSSR3707037079AAGGA2(Intron)ATP synthase F0 subunit a91P0000069Potential SSRPentaSSR3773637745ATTTG2ATP synthase F0 subunit a92P0000070Potential SSRPentaSSR3824638255TATTT293P0000071Potential SSRPentaSSR3882038829ACAAT2NAD494P0000072Potential SSRPentaSSR3894038949ATAAA2NAD495P0000073Potential SSRPentaSSR4020940218TTCAG2NAD196P0000074Potential SSRPentaSSR4049840507AATAC2NAD197P0000075Potential SSRPentaSSR4072440733GTTA2(Intron)NAD198P0000076Potential SSRPentaSSR4111541124AATGG2(Intron)NAD199P0000077Potential SSRPentaSSR4146341472AATAT2NAD1100P0000078Potential SSRPentaSSR4186741876TACAA2(Intron)NAD1101P0000079Potential SSRPentaSSR4197341982ATATT2NAD1102P0000080Potential SSRPentaSSR4241542424TAGTT2(Intron)NAD1103P0000081Potential SSRPentaSSR4316343172TACAC2(Intron)COX1104P0000082Potential SSRPentaSSR4380843817TATTT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1105P0000083Potential SSRPentaSSR4400744016AATTT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1106P0000084Potential SSRPentaSSR4407944088ATAT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1107P0000085Potential SSRPentaSSR4416044169TTATA2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1108P0000086Potential SSRPentaSSR4435944368TAATT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1109P0000087Potential SSRPentaSSR4771747726TGTTT2LAGLIDADG endonuclease (QPC56054.1), (Intron)COX1110P0000088Potential SSRPentaSSR4841148420ATATA2(Intron)COX1111P0000089Potential SSRPentaSSR4896548974TATAT2LAGLIDADG endonuclease (QPC56060.1), (Intron)COX1112P0000090Potential SSRPentaSSR4985249861TATTT2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1113P0000091Potential SSRPentaSSR5003850047ATAAA2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1114P0000092Potential SSRPentaSSR5057250581AAATA2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1115P0000093Potential SSRPentaSSR5068650695CATAG2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1116P0000094Potential SSRPentaSSR5090150910TATTT2LAGLIDADG endonuclease (QPC56059.1), (Intron)COX1117P0000095Potential SSRPentaSSR5210552114ATAG2(Intron)COX1118P0000096Potential SSRPentaSSR5216552174TATTT2(Intron)COX1119P0000097Potential SSRPentaSSR5315053159TTTAC2(Intron)COX1120P0000098Potential SSRPentaSSR5321453223ATAT2(Intron)COX1121P0000099Potential SSRPentaSSR5326153270TTATA2(Intron)COX1122P0000100Potential SSRPentaSSR5352553534ATATT2(Intron)COX1123P0000101Potential SSRPentaSSR5508555094ATAT2COX1124H0000001Potential SSRHexaSSR14091420ATTTAG2(Intron)Cob125H0000002Potential SSRHexaSSR15441555GAATTA2(Intron)Cob126H0000003Potential SSRHexaSSR18191830TTAATC2(Intron)Cob127H0000004Potential SSRHexaSSR23532364ATTTT2(Intron)Cob128H0000005Normal SSRHexaSSR25482565AAATAT3Cob129H0000006Potential SSRHexaSSR29963007TTTTTA2(Intron)Cob130H0000008Potential SSRHexaSSR59355946TTTATT2(Intron)Cob131H0000009Potential SSRHexaSSR65126523TAAATC2(Intron)Cob132H0000011Potential SSRHexaSSR75067517GATTA2LAGLIDADG/HNH endonuclease, (Intron)Cob133H0000012Potential SSRHexaSSR89658976AACTA2NAD5134H0000013Potential SSRHexaSSR99729983ATCCC2NAD5135H0000014Potential SSRHexaSSR1234812359TAAAT2(Intron)COX2136H0000015Potential SSRHexaSSR1247512486AAAGT2(Intron)COX2137H0000016Potential SSRHexaSSR1347913490ATTTA2(Intron)COX2138H0000017Potential SSRHexaSSR1794917960GTTAAT2NAD2139H0000018Potential SSRHexaSSR1797517986TAAAAA2NAD2140H0000019Potential SSRHexaSSR1935319364TAATAC2141H0000020Potential SSRHexaSSR2111021121TTTTAA2RPS3, (Intron)large subunit ribosomal RNA142H0000023Potential SSRHexaSSR2240322414TATGCC2Large subunit ribosomal RNA143H0000024Potential SSRHexaSSR2382423835TCCGCA2Large subunit ribosomal RNA144H0000025Potential SSRHexaSSR2458524596GAACT2Large subunit ribosomal RNA145H0000026Potential SSRHexaSSR2651026521AAATA2(Intron)large subunit ribosomal RNA146H0000027Potential SSRHexaSSR2704027051TATTTT2Large subunit ribosomal RNA147H0000028Potential SSRHexaSSR2766927680TTTAT2NAD6148H0000029Potential SSRHexaSSR2825328264TATTAA2NAD6149H0000030Potential SSRHexaSSR3100831019TCTGA2(Intron)COX3150H0000031Potential SSRHexaSSR3419634207TAGTT2151H0000032Potential SSRHexaSSR3680236813GTGTA2(Intron)ATP synthase F0 subunit a152H0000034Potential SSRHexaSSR3788537896AGATAA2ATP synthase F0 subunit a153H0000035Potential SSRHexaSSR4128741298ATTTAA2NAD1154H0000036Potential SSRHexaSSR4442544436TCCATC2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1155H0000037Potential SSRHexaSSR4577945790TCCATC2LAGLIDADG endonuclease (QPC56058.1), (Intron)COX1156H0000038Potential SSRHexaSSR4617546186TATTTA2(Intron)COX1157H0000039Potential SSRHexaSSR4634546356TTATT2(Intron)COX1158H0000040Potential SSRHexaSSR4660946620TTAATA2(Intron)COX1159H0000041Potential SSRHexaSSR4735847369ATAAAC2LAGLIDADG endonuclease (QPC56054.1), (Intron)COX1160H0000042Potential SSRHexaSSR5088950900TTTTAA2LAGLIDADG endonuclease (QPC56059.1), (Intron)COX1161H0000043Potential SSRHexaSSR5348353494CTTAT2(Intron)COX1162H0000044Potential SSRHexaSSR5410554116TTACCC2(Intron)COX1163H0000045Potential SSRHexaSSR5536455375TTCT2164cHp0000001Extended SSRHeptaSSR898911AATTATA2(Intron)Cob165cHp0000002Extended SSRHeptaSSR1397913992AATAATA2166cHp0000003Extended SSRHeptaSSR1590915922GGTATTT2Hypothetical protein, (Intron)NAD2167cHp0000005Extended SSRHeptaSSR3424234255TTATAA2Small subunit ribosomal RNA168cHp0000006Extended SSRHeptaSSR4493044943ATTATT2(Intron)COX1169O0000001Extended SSROctaSSR4066240677TTCATAT2(Intron)NAD1In the genic region, 84 SSRs were distributed in 24 different genes consisting of 21 PCGs, 2 rRNAs, and 1 tRNA (Figure4(c); Table 7). The large subunit RNA contained the most SSRs and the genes COX1, COX3, NAD3, two LAGLIDADG endonucleases, intron-encoded nuclease aI1, hypothetical protein, and tRNA-Glu contained the fewest (Figure 4(c); Table 7). Considering the length of these genes, some, including large submit RNA, NAD2, LAGLIDADG endonuclease (QPC56057.1), NAD1, NAD6, ATP synthase F0 subunit a, and LAGLIDADG/HNH endonuclease, displayed a relatively large number of SSRs (Figure 4(c); Table 7). Meanwhile, the remaining genes have a relatively low number of SSRs. This inequality of SSR distribution in PCGs can be another useful characteristic for developing efficient molecular markers. In addition, SSRs in PCGs are known to affect the functions of those PCGs especially for adaptation to environmental factors in fungi [125–127], suggesting that these SSRs can also affect the functions of mitochondrial PCGs.


### 3.4. Phylogenetic Analysis of 25 Fungal Mitogenomes of Ophiocordycipitaceae

We constructed bootstrapped maximum-likelihood (ML) and Bayesian inference (BI) phylogenetic trees using 26 fungal mitogenomes consisting of 5 mitogenomes used in this study, 25 mitogenomes in the Ophiocordycipitaceae family, and 1 outgroup species (Fusarium graminearum) [128]. Due to the incomplete annotation of the Ophiocordyceps sinensis fungal mitogenome (KP835313), five PCGs, NAD5, COB, COX1, NAD1, and NAD4, containing introns are not correctly annotated. Only five conserved PCGs, ATP8, COX2, NAD2, NAD3, and NAD4L, were selected and aligned individually. Subsequently, this alignment was concatenated to construct three phylogenetic trees.Five fungal endosymbiont mitogenomes of WBPH were well clustered with another fungal symbiont mitogenome ofR. speculum (NC_049089) [19] with high supportive values (Figure 5). This indicates taxonomic similarity between the R. speculum endosymbiont and the five WBPH endosymbionts, suggesting that other fungal endosymbionts may also be independently clustered with other fungal species in the sample family, Ophiocordycipitaceae. In terms of evolution, it can be explained by the two hypotheses: (i) independent evolution once this endosymbiont entered the host insect species or (ii) independent taxonomic groups of Ophiocordycipitaceae entering into the host insect species multiple times during evolution. To determine which hypothesis is more likely, we would need more endosymbiont mitogenomes from various host insect species of infraorder Fulgoromorpha and suborder Auchenorrhyncha as well as mitogenomes from neighboring noninsect endosymbiont fungal species.Figure 5
Phylogenetic trees of 26 fungal complete mitogenomes in Ophiocordycipitaceae. Neighbor-joining (bootstrap repeat is 10,000) and maximum-Likelihood (bootstrap repeat is 1000) phylogenetic trees as well as Bayesian inference tree (1,100,000 generations) of 26 fungal mitogenomes of Ophiocordycipitaceae. Phylogenetic tree was displayed based on the maximum-likelihood tree. The numbers above branches indicate bootstrap support values of maximum-likelihood and neighbor-joining phylogenetic trees and posterior possibility value of Bayesian inference tree, respectively. Scientific names inside the parenthesis indicate those of host species.Four fungal species used to investigate intraspecific variations in mitogenomes,Hirsutella thompsonii, Hirsutella rhossiliensis, Ophiocordyceps sinensis, and Tolypocladium inflatum, also display rigid clades covering all mitogenomes of each species with high supportive values (Figure 5). Three mitogenomes of Ophiocordyceps sinensis were clustered with the longest branch length among the four species, of which Hirsutella thompsonii had the second longest (Figure 5). These branch lengths were not proportional to the ratio of SNPs and INDELs (Table 4). The topology of the Tolypocladium genus in the trees was not congruent between the ML and BI trees with low bootstrap values (Figure 5), indicating that additional conserved gene sequences are required to resolve this clade properly.


## 3.1. Complete Mitogenome of Fungal WBPH Endosymbionts

We successfully assembled fungal endosymbiont mitogenomes from four WBPH samples isolated in Korea and China and one public dataset of NGS raw reads (Table1). This is the first WBPH fungal endosymbiont mitogenome identified. Their lengths ranged from 55,390 bp to 55,406 bp (Table 1), which is shorter than that of R. speculum (66,785 bp) [19]. In these mitogenomes, there were 28 protein-coding genes (PCGs), 12 tRNAs, and 2 rRNAs (Table 2). Some of the PCGs found were LAGLIDADG endonucleases, which are usually found in intronic regions of various fungal mitogenomes, contributing to the expansion of their length [19, 55, 108–112]. In comparison to the previously sequenced mitogenome of the fungal endosymbiont of R. speculum, there were slightly fewer PCGs and tRNAs found in the WBPH endosymbiont mitogenomes. There were three fewer PCGs for three reasons: the smaller number of LAGLIDADG endonucleases, the absence of one endonuclease and a GIY-YIG endonuclease, and the presence of two additional PCGs—a hypothetical protein and a LAGLIDADG/HNH endonuclease. This particular configuration of PCGs is usually identified in other fungal mitogenomes; for example, two mitogenomes of Fusarium oxysporum (GenBank accessions are MN259514 and MN259515) display two completely different PCGs in each mitogenome [54, 56]. There are also five fewer tRNAs because of the different configurations: tRNA-Asp, tRNA-Cys, tRNA-Ile, and two tRNA-Ser (also found in the mitogenome of the fungal symbiont of R. speculum [19]). This difference in configuration of tRNAs between two different fungal symbionts suggests that tRNA configuration may not be critical because essential tRNAs absent in the fungal mitogenome can be supported from the nuclear genome [113].Table 2
List of genes in WBPH fungal endosymbiont mitogenomes.
No.NameTypeStart positionEnd positionLength (bp)StrandNo. of exons1Cytochrome bCDS1608182975Reverse52LAGLIDADG/HNH endonucleaseCDS74297776348Reverse13NADH dehydrogenase subunit 5CDS843210,4051974Reverse14NADH dehydrogenase subunit 4 LCDS10,40510,680276Reverse15Cytochrome c oxidase subunit IICDS10,89613,952894Reverse66ATP synthase F0 subunit cCDS14,05214,237186Reverse17NADH dehydrogenase subunit 3CDS14,31514,968447Reverse28LAGLIDADG endonucleaseCDS14,68314,880198Reverse19NADH dehydrogenase subunit 2CDS14,96918,1931431Reverse210Hypothetical proteinCDS15,62616,390765Reverse111tRNA-MettRNA18,21018,28374Reverse112tRNA-HistRNA18,32518,39874Reverse113tRNA-LeutRNA18,48518,56985Reverse114tRNA-LystRNA18,57018,64273Reverse115tRNA-PhetRNA18,64318,71573Reverse116tRNA-LeutRNA19,67319,75583Reverse117tRNA-MettRNA19,75819,83073Reverse118tRNA-GlutRNA19,90419,97673Reverse119Large subunit rRNArRNA20,06727,3374793Reverse320Ribosomal protein S3CDS20,67221,9161245Reverse121tRNA-ProtRNA27,40027,47172Reverse122NADH dehydrogenase subunit 6CDS27,67728,441765Reverse123tRNA-GlytRNA28,54428,61471Reverse124Cytochrome c oxidase subunit IIICDS29,22432,387795Reverse325Intron-encoded endonuclease aI1CDS29,93230,138207Reverse126Hypothetical proteinCDS30,39530,664270Reverse127Alpha-beta-hydrolaseCDS30,42330,686264Forward128tRNA-AsntRNA32,43732,50872Reverse129Small subunit rRNArRNA34,25435,6151362Reverse130ATP synthase F0 subunit aCDS35,91838,001762Reverse231ATP synthase F0 subunit 8CDS38,06138,237177Reverse132NADH dehydrogenase subunit 4CDS38,62040,0351416Reverse133NADH dehydrogenase subunit 1CDS40,11942,5821107Reverse434tRNA-ArgtRNA42,72642,79974Reverse135Cytochrome c oxidase subunit ICDS42,80955,2931314Reverse736LAGLIDADG endonucleaseCDS43,66044,8801221Reverse137LAGLIDADG endonucleaseCDS45,12546,054930Reverse138LAGLIDADG endonucleaseCDS46,85748,2481392Reverse139LAGLIDADG endonucleaseCDS48,62049,483864Reverse140LAGLIDADG endonucleaseCDS49,52550,8141290Reverse141LAGLIDADG endonucleaseCDS50,86551,734870Reverse142Hypothetical proteinCDS54,29054,850561Forward1Several PCGs in the fungal mitogenomes have been invaded by introns multiple times. For example,COX1 contains three introns, and COB has five introns in the Hirsutella thompsonii mitogenome [114]. This phenomenon contributes to increased fungal mitogenome: Aspergillus pseudoglaucus and Aspergillus egyptiacus are longer than the other Aspergillus mitogenomes because of the presence of many introns on major PCGs [55, 115]. The fungal mitogenomes examined in this study also present many introns on PCGs including COB, COX1, NAD1, ATP8, COX3, COX2, and NAD2 (Figure 1), which is a major reason for the expansion of fungal mitogenomes together with endonucleases.Figure 1
Complete mitogenome of WBPH fungal endosymbionts. Blue circle indicates fungal endosymbiont mitogenome, yellow arrows are protein-coding genes, purple arrows are tRNAs, and red arrows mean rRNAs. Direction of arrows indicates direction of transcription. Each gene name was displayed with lines directing to the corresponded arrows. Numbers displayed outside of blue circle mean base pair.The gene order of WBPH andR. speculum fungal symbiont mitogenomes was the same when PCGs except endonucleases and rRNAs are considered. However, intron structures of COX1, COX2, NAD2, NAD3, NAD5, and ATP synthase F0 subunit present different configurations between the two mitogenomes (Figure 2). The intron structures of NAD5 and NAD2 present reduce of a reduction in the number of exons via removal of intron regions in the WBPH fungal endosymbiont mitogenome (Figures 2(a) and 2(d)), whereas those of COX2, NAD3, and the ATP synthase F0 subunit display insertions of one intron into the WBPH fungal endosymbiont mitogenome (Figures 2(b), 2(c), and 2(e)). This indicates that the reduction in the total length of the WBPH fungal symbiont mitogenome is not primarily caused by reducing the number of exons, unlike in Aspergillus mitogenomes [55, 116]. In addition, COX1, which contains the largest number of exons in these mitogenomes, lost the sixth and seventh exons of the R. speculum fungal endosymbiont mitogenome in the mitogenome of WBPH endosymbiont (Figure 2(f)). However, the total length of COX1, including the introns of WBPH fungal endosymbionts, is longer than that of R. speculum fungal endosymbionts by 1 kb (Figure 2(f)), reflecting complex events that occurred during the evolution of both mitogenomes. Additional studies are required to identify the correct exons of the COX1 gene of this fungal endosymbiont. For example, alignment of RNA-Seq raw reads against this mitogenome could provide expressed regions in this mitogenome.Figure 2
Schematic diagram of exon/intron structure of protein-coding genes displaying different configurations between fungal endosymbiont mitogenomes of WBPH andR. speculum. Orange-colored diagrams indicate components of WBPH fungal endosymbionts, and blue-colored diagrams are those of R. speculum fungal endosymbionts. Exon/intron structures of (a) NAD5, (b) COX2, (c) NAD3, (d) NAD2, (e) ATP synthase F0 subunit, and (f) COX1 are displayed.
(a)(b)(c)(d)(e)(f)Once more fungal symbiont mitogenomes are available, patterns of presence and absence of tRNAs, additional endonucleases, and intron structures of PCGs in endosymbiont mitogenomes will elucidate a detailed evolutionary history of these genes.


## 3.2. Identification of Intraspecific Variations on Fungal WBPH Endosymbiont Mitogenomes

We identified two SNPs, three insertions, and two deletions via multiple sequence alignments of the five fungal mitogenomes (Table3). One of two SNPs was identified in KR.5D WBPH and changed leucine (L) to glutamine (Q) in the ATP synthase F0 subunit (Table 3). One 10 bp insertion in the intergenic space was found in KR.1D WBPH, while the remaining two insertions and all three deletions were 1 to 3 bp in length (Table 3).Table 3
List of intraspecific variations identified on the five WBPH fungal endosymbiont mitogenomes.
No.TypeCoordination of multiple sequence alignmentsStrainsBase changesPosition1Insertion4209-4210KR.1D, KR.5D- to CCIntergenic2Insertion27,476-27,486KR.1D- to TGGGCCCCCCIntergenic3SNP27,487KR.1DA to CIntergenic4Deletion32,727-32,728KR.5DCC to -Intergenic5SNP37,574KR.5DA to TL to Q in ATP synthase F0 subunit6Insertion38,727KR.1D- to GIntergenic7Deletion38,728-38,730KR.5DGGG to -IntergenicThe proportions of these intraspecific SNPs, insertions, and deletions in these fungal mitogenomes were 0.0036%, 0.020%, and 0.012%, respectively. The proportion of insertions and deletions was higher than that of SNPs. Interestingly, there is geographical variation in the fungal symbiont mitogenomes. The mitogenome of WBPH endosymbionts used in the whole genome sequencing (WGS) and the KR.11D isolate were identical to that of KR, while the other three WBPHs captured in other locations in Korea displayed intraspecific variations. The sample used in the WGS originated from the University of Science and Technology of China (Anhui province, China), indicating that KR 11D and KR WBPH samples obtained in Korea have migrated from the similar region to the WGS sample. However, further analyses of their complete mitogenomes or whole genomes will be needed to provide more supportive data for identifying their origins.There is a relatively small number of intraspecific SNPs and INDELs identified from these fungal mitogenomes in comparison to those of other fungal mitogenomes, for example, 16 to 17 SNPs (0.055% to 0.0582%) and 22 to 27 INDELs (0.075% to 0.092%) onAspergillus flavus [52, 53] and 62 SNPs (0.15%) and 181 INDELs (0.43%) on Fusarium oxysporum f.sp. lactucae [56]. They are also fewer than those identified in insect mitogenomes [10, 22, 23, 43, 45–51].Based on 25 available complete fungal mitogenomes in Ophiocordycipitaceae, four species,Ophiocordyceps sinensis, Hirsutella thompsonii, Hirsutella rhossiliensis, and Tolypocladium inflatum, contain more than one complete fungal mitogenome (Table 4). We investigated intraspecific variations in the mitogenomes of these four species (Table 5). There are significantly more INDELs than SNPs identified in the four fungal species, a trend identical to that observed in the four mitogenomes of fungal endosymbiont WBPH with the exception of their absolute amounts. Moreover, there were at least three times more SNPs and INDELs in these fungal mitogenomes than that in the fungal symbiont of WBPHs. This phenomenon can be explained by two major factors: first, the geographical distribution or genetic background of WBPH samples is relatively limited in comparison to those of the four fungal species, and second, the surroundings of fungal endosymbionts are less dynamic than those of normal fungal species, causing low selection pressure from the environment. This second factor is supported by two studies: first, the bacterial genome of aphid endosymbiont Buchnera aphidicola (Aphis gossypii) displays a low level of intraspecific variation in comparison to those of host mitogenome (Bae et al., under revision), and second, the whole genome of endosymbiont of Pediculus humanus capitis also shows low-level intraspecific variations in comparison to those of their whole genomes [117].Table 4
List of available complete fungal mitogenomes in Ophiocordycipitaceae.
No.SpeciesNCBI accessionLength (bp)GC ratio (%)No. of PCGsNo. of tRNAsNo. of rRNAsReference1Ophiocordyceps sinensisNC_034659157,53930.288272Unpub.2Ophiocordyceps sinensisMH400233157,58430.276272Unpub.3Ophiocordyceps sinensisKP835313157,51030.26∗232[123]4Hirsutella thompsoniiNC_04016562,50929.830272[114]5Hirsutella thompsoniiMH36729665,33230.332272[114]6Hirsutella thompsoniiMH36729560,36230.029272[114]7Hirsutella rhossiliensisMG97907162,94928.333262[129]8Hirsutella rhossiliensisNC_03016462,48328.224262Unpub.9Hirsutella vermicolaNC_03661053,79325.327252[130]10Hirsutella minnesotensisNC_02766052,24528.430252[131]11Tolypocladium sp.MN58326546,46626.115262[132]12Tolypocladium guangdongenseMT47126746,10226.130272[133]13Tolypocladium ophioglossoidesNC_03138435,15927.519252[134]14Tolypocladium cylindrosporumNC_04683934,69827,024262[135]15Tolypocladium inflatumNC_03638225,32827.815252[136]16Tolypocladium inflatumKY92488025,23827.815252[136]17Tolypocladium inflatumKY92488125,32827.815252[136]18Tolypocladium inflatumKY92488225,32827.815252[136]19Tolypocladium inflatumKY92488324,79327.815252[136]20Ophiocordycipitaceae sp.NC_04908966,78530.631172[19]∗Mitogenome annotation of this genome (KP835313) seems not to be complete because several major genes, such as COX1, NAD1, NAD5, and COB, which have many introns on fungal mitogenomes that were not annotated as CDS.Table 5
Number of intraspecific variations identified from four fungal species belonging to Ophiocordycipitaceae.
No.SpeciesNo. of mitogenomesAligned length (bp)No. of SNPsSNP coverage (%)No. of INDELsINDEL coverage (%)1Ophiocordyceps sinensis3157,606160.0101440.0912Hirsutella thompsonii366,6352810.4264899.743Hirsutella rhossiliensis264,85870.0120083.104Tolypocladium inflatum525,338300.123751.48


## 3.3. Identification and Comparative Analysis of Simple Sequence Repeats on the Five WBPH Fungal Endosymbiont Mitogenomes

Simple sequence repeats (SSRs) identified from organellar genomes have been utilized as molecular markers in various species such as plant species [99, 118–122], suggesting that SSRs on fungal endosymbiont mitogenomes can be used as molecular markers to identify the geographical origins of WBPH. In total, 23 normal and 6 extended SSRs were identified from fungal endosymbiont mitogenomes (Figure 3(b)), with the exception of the fungal endosymbiont mitogenome of WBPH KR.1D which displays 24 normal and 6 extended SSRs (Table 6). The fungal endosymbiont mitogenome of WBPH KR.1D has one more monoSSR (Table 6) with a unit sequence of C and length of 15 bp caused by one insertion (Table 3). In addition, 140 potential SSRs were also identified in the five mitogenomes (Table 6). SSRs identified in the mitogenome were distributed evenly (Figure 3(a)), suggesting that there was no hot spot of SSRs in these fungal mitogenomes.Figure 3
Numbers and distribution of SSRs on the mitogenome of fungal endosymbiont of WBPH KR. (a) Black circle indicates fungal endosymbiont mitogenome, yellow arrows are protein-coding genes, purple arrows are tRNAs, and red arrows mean rRNAs. SSRs, extended SSRs, and potential SSRs are displayed with yellow, green, and light green colors, respectively. (b) The number of SSRs along with SSR types is displayed. Blue color indicates SSRs and extended SSRs, and orange color means potential SSRs.Table 6
List of SSRs identified from the five fungal mitogenomes of WBPH endosymbionts.
SSR typeKRKR.1DKR.5DKR.11DWGSMonoSSR89888DiSSR88888TriSSR22222TetraSSR44444PentaSSR00000HexaSSR11111HeptaSSR55555OctaSSR11111NonaSSR00000DecaSSR00000Subtotal2930292929PentaPotentialSSR101101101101101HexaPotentialSSR3939393939Subtotal140140140140140The length of the identified SSRs is relatively short (a maximum length of 18 bp; Figure4(a)) in comparison to those of other fungal species in the same family: Ophiocordyceps sinensis (up to 24 bp) [123], as well as fungal species in the other families, such as Pestalotiopsis fici (up to 45 bp) [124]. Moreover, the maximum length of SSRs identified from the mitogenome of R. speculum (NC_049089) [19] was 18 bp, suggesting that this short SSR length can be linked to the evolution of endosymbiont mitogenomes.Figure 4
Comparative analysis of simple sequence repeats on mitogenome of WBPH fungal endosymbiont. (a) Number of SSRs, extended SSRs, and potential SSRs along with its length.x-axis indicates SSR length (bp), and y-axis means the number of SSRs. Blue-colored bars indicate SSRs and extended SSRs, and orange bars mean potential SSRs. (b) Number of SSRs, extended SSRs, and potential SSRs based on positions, genic, intronic ORF, intergenic, and intronic. (c) Number of genic SSRs along with genes (blue bars) as well as length of genes (orange line). x-axis is genes containing SSRs, and y-axis indicates the number of SSRs (left) and gene length (bp; right).
(a)(b)(c)Out of 191 normal SSRs, extended SSRs, and potential SSRs, 84 (43.98%) are located in the genic region (genic and intronic ORF categories in Figure4(b); Table 7). The intronic ORF position indicates the location of the PCGs placed at the introns of other PCGs, most of which are LAGLIDADG endonucleases (Table 2). Nearly half of the SSRs are in PCGs, which are conserved in comparison to intron and intergenic regions, indicating that these SSRs can be utilized for distinguishing species level or even higher rank. In the intergenic region, there were 61 SSRs (31.94%), and in comparison, only 24 SSRs (12.57%) in the intergenic region (Figure 4(b); Table 7). These SSRs are located in relatively nonconserved regions in comparison to PCG regions, suggesting that these SSRs can be used to distinguish intraspecific levels, such as population or geographical origins. Once more endosymbiont mitogenomes are available in the near future, these SSRs can be evaluated for their use in identification of species and their geographical origin as well as evolutionary history of their mitogenomes.Table 7
List of SSRs identified on fungal endosymbiont mitogenome of WBPH KR.
No.NameSSR typeTypeStartEndUnit sequenceRepeat numberGenes1M0000001Normal SSRMonoSSR42094219C11(Intron)Cob2M0000002Normal SSRMonoSSR55485558A11Cob3M0000003Normal SSRMonoSSR2231422323T10Large subunit ribosomal RNA4M0000004Normal SSRMonoSSR2248222492T11Large subunit ribosomal RNA5M0000005Normal SSRMonoSSR2844228451T106M0000006Normal SSRMonoSSR3271532728C147M0000007Normal SSRMonoSSR3831538327G138M0000008Normal SSRMonoSSR4003540044T10NAD49D0000001Normal SSRDiSSR57885797TA5(Intron)Cob10D0000002Normal SSRDiSSR1370813721AT7(Intron)COX211D0000003Normal SSRDiSSR2006820077AT5Large subunit ribosomal RNA12D0000004Normal SSRDiSSR2064620657AT6(Intron)large subunit ribosomal RNA13D0000005Normal SSRDiSSR2726627275TA5Large subunit ribosomal RNA14D0000006Normal SSRDiSSR3799538004TA5ATP synthase F0 subunit a15D0000007Normal SSRDiSSR4043040441AT6NAD116D0000008Normal SSRDiSSR5535255361TA517T0000001Normal SSRTriSSR3582935840ATT418T0000002Normal SSRTriSSR5278052791ATA4(Intron)COX119Te0000001Normal SSRTetraSSR2472124732ATTT3Large subunit ribosomal RNA20Te0000002Normal SSRTetraSSR4259842609TTTA321Te0000003Normal SSRTetraSSR4847748488ATAA3(Intron)COX122Te0000004Normal SSRTetraSSR5270952720AATA3(Intron)COX123P0000001Potential SSRPentaSSR586595TTGT2(Intron)Cob24P0000002Potential SSRPentaSSR19231932TAATA2(Intron)Cob25P0000003Potential SSRPentaSSR35233532TAAAA2(Intron)Cob26P0000004Potential SSRPentaSSR41164125TTGTC2(Intron)Cob27P0000005Potential SSRPentaSSR50685077ATAAT2(Intron)Cob28P0000006Potential SSRPentaSSR65406549TAATG2(Intron)Cob29P0000007Potential SSRPentaSSR66466655ATTTT2(Intron)Cob30P0000008Potential SSRPentaSSR67146723TTTT2(Intron)Cob31P0000009Potential SSRPentaSSR76367645AGCAA2LAGLIDADG/HNH endonuclease, (Intron)Cob32P0000010Potential SSRPentaSSR91099118AAGTT2NAD533P0000011Potential SSRPentaSSR92649273ATAA2NAD534P0000012Potential SSRPentaSSR1030510314AGACA2NAD535P0000013Potential SSRPentaSSR1085710866ATTCA236P0000014Potential SSRPentaSSR1136911378AGATA2COX237P0000015Potential SSRPentaSSR1251112520TTATA2(Intron)COX238P0000016Potential SSRPentaSSR1260012609TAAGA2(Intron)COX239P0000017Potential SSRPentaSSR1271012719AAGCG2(Intron)COX240P0000018Potential SSRPentaSSR1279612805TTAAC2(Intron)COX241P0000019Potential SSRPentaSSR1330313312TAATA2COX242P0000020Potential SSRPentaSSR1494614955AAAAG2NAD343P0000021Potential SSRPentaSSR1674716756TCGAG2(Intron)NAD244P0000022Potential SSRPentaSSR1738317392TCATT2NAD245P0000023Potential SSRPentaSSR1742717436AATAA2NAD246P0000024Potential SSRPentaSSR1751417523AAATG2NAD247P0000025Potential SSRPentaSSR1768517694AATAA2NAD248P0000026Potential SSRPentaSSR1808318092AATAC2NAD249P0000027Potential SSRPentaSSR1810918118TATT2NAD250P0000028Potential SSRPentaSSR1812118130ATAGA2NAD251P0000029Potential SSRPentaSSR1840018409TTATG252P0000030Potential SSRPentaSSR1881918828GATA253P0000031Potential SSRPentaSSR1907819087ATTTT254P0000032Potential SSRPentaSSR1919019199TTGTA255P0000033Potential SSRPentaSSR1930219311ATAAT256P0000034Potential SSRPentaSSR1961319622AACT257P0000035Potential SSRPentaSSR1962919638TATT258P0000036Potential SSRPentaSSR1990419913TAGAC2tRNA-Glu59P0000037Potential SSRPentaSSR2038420393TTATT2Large subunit ribosomal RNA60P0000038Potential SSRPentaSSR2096820977TTATT2RPS3, (Intron)large subunit ribosomal RNA61P0000039Potential SSRPentaSSR2114621155TGTAT2RPS3, (Intron)large subunit ribosomal RNA62P0000040Potential SSRPentaSSR2117321182TATTA2RPS3, (Intron)large subunit ribosomal RNA63P0000041Potential SSRPentaSSR2173021739TTATT2RPS3, (Intron)large subunit ribosomal RNA64P0000042Potential SSRPentaSSR2203822047TTTTA2(Intron)large subunit ribosomal RNA65P0000043Potential SSRPentaSSR2212122130TTATT2(Intron)large subunit ribosomal RNA66P0000044Potential SSRPentaSSR2260622615TAATA2Large subunit ribosomal RNA67P0000045Potential SSRPentaSSR2351223521AAGAC2Large subunit ribosomal RNA68P0000046Potential SSRPentaSSR2396223971TTTTC2Large subunit ribosomal RNA69P0000047Potential SSRPentaSSR2445324462AATTA2Large subunit ribosomal RNA70P0000048Potential SSRPentaSSR2450124510ATTTA2Large subunit ribosomal RNA71P0000049Potential SSRPentaSSR2521125220TTTAC2Large subunit ribosomal RNA72P0000050Potential SSRPentaSSR2535725366TTTT2Large subunit ribosomal RNA73P0000051Potential SSRPentaSSR2618126190CATTT2(Intron)large subunit ribosomal RNA74P0000052Potential SSRPentaSSR2722427233ATTTC2Large subunit ribosomal RNA75P0000053Potential SSRPentaSSR2770027709TTAAG2NAD676P0000054Potential SSRPentaSSR2784427853ATAAT2NAD677P0000055Potential SSRPentaSSR2801328022TAAAA2NAD678P0000056Potential SSRPentaSSR2911929128TCCCC279P0000057Potential SSRPentaSSR2927329282CAGTA2COX380P0000058Potential SSRPentaSSR2997329982TGAT2Intron-encoded endonuclease aI1, (Intron)COX381P0000059Potential SSRPentaSSR3086030869AGTG2(Intron)COX382P0000060Potential SSRPentaSSR3261832627TCCCC283P0000061Potential SSRPentaSSR3339733406TAAAT284P0000062Potential SSRPentaSSR3341633425ATGGT285P0000063Potential SSRPentaSSR3391633925AGAGA286P0000064Potential SSRPentaSSR3484234851AATT2Small subunit ribosomal RNA87P0000065Potential SSRPentaSSR3668036689TTAAA2(Intron)ATP synthase F0 subunit a88P0000066Potential SSRPentaSSR3699337002TTAAA2(Intron)ATP synthase F0 subunit a89P0000067Potential SSRPentaSSR3702137030ATTTT2(Intron)ATP synthase F0 subunit a90P0000068Potential SSRPentaSSR3707037079AAGGA2(Intron)ATP synthase F0 subunit a91P0000069Potential SSRPentaSSR3773637745ATTTG2ATP synthase F0 subunit a92P0000070Potential SSRPentaSSR3824638255TATTT293P0000071Potential SSRPentaSSR3882038829ACAAT2NAD494P0000072Potential SSRPentaSSR3894038949ATAAA2NAD495P0000073Potential SSRPentaSSR4020940218TTCAG2NAD196P0000074Potential SSRPentaSSR4049840507AATAC2NAD197P0000075Potential SSRPentaSSR4072440733GTTA2(Intron)NAD198P0000076Potential SSRPentaSSR4111541124AATGG2(Intron)NAD199P0000077Potential SSRPentaSSR4146341472AATAT2NAD1100P0000078Potential SSRPentaSSR4186741876TACAA2(Intron)NAD1101P0000079Potential SSRPentaSSR4197341982ATATT2NAD1102P0000080Potential SSRPentaSSR4241542424TAGTT2(Intron)NAD1103P0000081Potential SSRPentaSSR4316343172TACAC2(Intron)COX1104P0000082Potential SSRPentaSSR4380843817TATTT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1105P0000083Potential SSRPentaSSR4400744016AATTT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1106P0000084Potential SSRPentaSSR4407944088ATAT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1107P0000085Potential SSRPentaSSR4416044169TTATA2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1108P0000086Potential SSRPentaSSR4435944368TAATT2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1109P0000087Potential SSRPentaSSR4771747726TGTTT2LAGLIDADG endonuclease (QPC56054.1), (Intron)COX1110P0000088Potential SSRPentaSSR4841148420ATATA2(Intron)COX1111P0000089Potential SSRPentaSSR4896548974TATAT2LAGLIDADG endonuclease (QPC56060.1), (Intron)COX1112P0000090Potential SSRPentaSSR4985249861TATTT2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1113P0000091Potential SSRPentaSSR5003850047ATAAA2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1114P0000092Potential SSRPentaSSR5057250581AAATA2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1115P0000093Potential SSRPentaSSR5068650695CATAG2LAGLIDADG endonuclease (QPC56055.1), (Intron)COX1116P0000094Potential SSRPentaSSR5090150910TATTT2LAGLIDADG endonuclease (QPC56059.1), (Intron)COX1117P0000095Potential SSRPentaSSR5210552114ATAG2(Intron)COX1118P0000096Potential SSRPentaSSR5216552174TATTT2(Intron)COX1119P0000097Potential SSRPentaSSR5315053159TTTAC2(Intron)COX1120P0000098Potential SSRPentaSSR5321453223ATAT2(Intron)COX1121P0000099Potential SSRPentaSSR5326153270TTATA2(Intron)COX1122P0000100Potential SSRPentaSSR5352553534ATATT2(Intron)COX1123P0000101Potential SSRPentaSSR5508555094ATAT2COX1124H0000001Potential SSRHexaSSR14091420ATTTAG2(Intron)Cob125H0000002Potential SSRHexaSSR15441555GAATTA2(Intron)Cob126H0000003Potential SSRHexaSSR18191830TTAATC2(Intron)Cob127H0000004Potential SSRHexaSSR23532364ATTTT2(Intron)Cob128H0000005Normal SSRHexaSSR25482565AAATAT3Cob129H0000006Potential SSRHexaSSR29963007TTTTTA2(Intron)Cob130H0000008Potential SSRHexaSSR59355946TTTATT2(Intron)Cob131H0000009Potential SSRHexaSSR65126523TAAATC2(Intron)Cob132H0000011Potential SSRHexaSSR75067517GATTA2LAGLIDADG/HNH endonuclease, (Intron)Cob133H0000012Potential SSRHexaSSR89658976AACTA2NAD5134H0000013Potential SSRHexaSSR99729983ATCCC2NAD5135H0000014Potential SSRHexaSSR1234812359TAAAT2(Intron)COX2136H0000015Potential SSRHexaSSR1247512486AAAGT2(Intron)COX2137H0000016Potential SSRHexaSSR1347913490ATTTA2(Intron)COX2138H0000017Potential SSRHexaSSR1794917960GTTAAT2NAD2139H0000018Potential SSRHexaSSR1797517986TAAAAA2NAD2140H0000019Potential SSRHexaSSR1935319364TAATAC2141H0000020Potential SSRHexaSSR2111021121TTTTAA2RPS3, (Intron)large subunit ribosomal RNA142H0000023Potential SSRHexaSSR2240322414TATGCC2Large subunit ribosomal RNA143H0000024Potential SSRHexaSSR2382423835TCCGCA2Large subunit ribosomal RNA144H0000025Potential SSRHexaSSR2458524596GAACT2Large subunit ribosomal RNA145H0000026Potential SSRHexaSSR2651026521AAATA2(Intron)large subunit ribosomal RNA146H0000027Potential SSRHexaSSR2704027051TATTTT2Large subunit ribosomal RNA147H0000028Potential SSRHexaSSR2766927680TTTAT2NAD6148H0000029Potential SSRHexaSSR2825328264TATTAA2NAD6149H0000030Potential SSRHexaSSR3100831019TCTGA2(Intron)COX3150H0000031Potential SSRHexaSSR3419634207TAGTT2151H0000032Potential SSRHexaSSR3680236813GTGTA2(Intron)ATP synthase F0 subunit a152H0000034Potential SSRHexaSSR3788537896AGATAA2ATP synthase F0 subunit a153H0000035Potential SSRHexaSSR4128741298ATTTAA2NAD1154H0000036Potential SSRHexaSSR4442544436TCCATC2LAGLIDADG endonuclease (QPC56057.1), (Intron)COX1155H0000037Potential SSRHexaSSR4577945790TCCATC2LAGLIDADG endonuclease (QPC56058.1), (Intron)COX1156H0000038Potential SSRHexaSSR4617546186TATTTA2(Intron)COX1157H0000039Potential SSRHexaSSR4634546356TTATT2(Intron)COX1158H0000040Potential SSRHexaSSR4660946620TTAATA2(Intron)COX1159H0000041Potential SSRHexaSSR4735847369ATAAAC2LAGLIDADG endonuclease (QPC56054.1), (Intron)COX1160H0000042Potential SSRHexaSSR5088950900TTTTAA2LAGLIDADG endonuclease (QPC56059.1), (Intron)COX1161H0000043Potential SSRHexaSSR5348353494CTTAT2(Intron)COX1162H0000044Potential SSRHexaSSR5410554116TTACCC2(Intron)COX1163H0000045Potential SSRHexaSSR5536455375TTCT2164cHp0000001Extended SSRHeptaSSR898911AATTATA2(Intron)Cob165cHp0000002Extended SSRHeptaSSR1397913992AATAATA2166cHp0000003Extended SSRHeptaSSR1590915922GGTATTT2Hypothetical protein, (Intron)NAD2167cHp0000005Extended SSRHeptaSSR3424234255TTATAA2Small subunit ribosomal RNA168cHp0000006Extended SSRHeptaSSR4493044943ATTATT2(Intron)COX1169O0000001Extended SSROctaSSR4066240677TTCATAT2(Intron)NAD1In the genic region, 84 SSRs were distributed in 24 different genes consisting of 21 PCGs, 2 rRNAs, and 1 tRNA (Figure4(c); Table 7). The large subunit RNA contained the most SSRs and the genes COX1, COX3, NAD3, two LAGLIDADG endonucleases, intron-encoded nuclease aI1, hypothetical protein, and tRNA-Glu contained the fewest (Figure 4(c); Table 7). Considering the length of these genes, some, including large submit RNA, NAD2, LAGLIDADG endonuclease (QPC56057.1), NAD1, NAD6, ATP synthase F0 subunit a, and LAGLIDADG/HNH endonuclease, displayed a relatively large number of SSRs (Figure 4(c); Table 7). Meanwhile, the remaining genes have a relatively low number of SSRs. This inequality of SSR distribution in PCGs can be another useful characteristic for developing efficient molecular markers. In addition, SSRs in PCGs are known to affect the functions of those PCGs especially for adaptation to environmental factors in fungi [125–127], suggesting that these SSRs can also affect the functions of mitochondrial PCGs.


## 3.4. Phylogenetic Analysis of 25 Fungal Mitogenomes of Ophiocordycipitaceae

We constructed bootstrapped maximum-likelihood (ML) and Bayesian inference (BI) phylogenetic trees using 26 fungal mitogenomes consisting of 5 mitogenomes used in this study, 25 mitogenomes in the Ophiocordycipitaceae family, and 1 outgroup species (Fusarium graminearum) [128]. Due to the incomplete annotation of the Ophiocordyceps sinensis fungal mitogenome (KP835313), five PCGs, NAD5, COB, COX1, NAD1, and NAD4, containing introns are not correctly annotated. Only five conserved PCGs, ATP8, COX2, NAD2, NAD3, and NAD4L, were selected and aligned individually. Subsequently, this alignment was concatenated to construct three phylogenetic trees.Five fungal endosymbiont mitogenomes of WBPH were well clustered with another fungal symbiont mitogenome ofR. speculum (NC_049089) [19] with high supportive values (Figure 5). This indicates taxonomic similarity between the R. speculum endosymbiont and the five WBPH endosymbionts, suggesting that other fungal endosymbionts may also be independently clustered with other fungal species in the sample family, Ophiocordycipitaceae. In terms of evolution, it can be explained by the two hypotheses: (i) independent evolution once this endosymbiont entered the host insect species or (ii) independent taxonomic groups of Ophiocordycipitaceae entering into the host insect species multiple times during evolution. To determine which hypothesis is more likely, we would need more endosymbiont mitogenomes from various host insect species of infraorder Fulgoromorpha and suborder Auchenorrhyncha as well as mitogenomes from neighboring noninsect endosymbiont fungal species.Figure 5
Phylogenetic trees of 26 fungal complete mitogenomes in Ophiocordycipitaceae. Neighbor-joining (bootstrap repeat is 10,000) and maximum-Likelihood (bootstrap repeat is 1000) phylogenetic trees as well as Bayesian inference tree (1,100,000 generations) of 26 fungal mitogenomes of Ophiocordycipitaceae. Phylogenetic tree was displayed based on the maximum-likelihood tree. The numbers above branches indicate bootstrap support values of maximum-likelihood and neighbor-joining phylogenetic trees and posterior possibility value of Bayesian inference tree, respectively. Scientific names inside the parenthesis indicate those of host species.Four fungal species used to investigate intraspecific variations in mitogenomes,Hirsutella thompsonii, Hirsutella rhossiliensis, Ophiocordyceps sinensis, and Tolypocladium inflatum, also display rigid clades covering all mitogenomes of each species with high supportive values (Figure 5). Three mitogenomes of Ophiocordyceps sinensis were clustered with the longest branch length among the four species, of which Hirsutella thompsonii had the second longest (Figure 5). These branch lengths were not proportional to the ratio of SNPs and INDELs (Table 4). The topology of the Tolypocladium genus in the trees was not congruent between the ML and BI trees with low bootstrap values (Figure 5), indicating that additional conserved gene sequences are required to resolve this clade properly.


## 4. Conclusions

We successfully elucidated the five complete mitogenomes of the fungal endosymbiont of WBPH from various sources of NGS raw reads obtained from WBPH samples. These five complete mitogenomes show common and their own characteristics in comparison to the previously elucidated complete mitogenome of theR. japonica fungal endosymbiont [19]. There were fewer intraspecific variations in the five WBPH endosymbiont mitogenomes in comparison to those identified from the four Ophiocordycipitaceae fungal species, Ophiocordyceps sinensis, Hirsutella thompsonii, Hirsutella rhossiliensis, and Tolypocladium inflatum. This can be explained by the narrow geographical distribution and/or genetic background and the low selection pressures of endosymbionts. We identified 191 SSRs were from each WBPH fungal symbiont complete mitogenomes, except for the WBPH_KR.1D mitogenome, which presented an additional SSR. These SSRs are relatively short in length (a maximum length of 18 bp) compared to those of other fungal mitogenomes. Nearly half of the SSRs are in the genic region, suggesting that these SSRs may be more conserved and they may affect the functionality of PCGs. Based on the phylogenetic trees of 5 conserved PCGs of 26 fungal mitogenomes, including one outgroup species, WBPH fungal endosymbiont mitogenomes were clustered with that of R. speculum with high supportive values. This suggests that these insect-hosted fungal endosymbionts have been evolved independently from the other fungal species in the Ophiocordycipitaceae family. Owing to the advantages of NGS raw reads, which can detect sequences from unknown or unexpected organisms [12, 19–37], we successfully identified the complete mitogenomes of WBPH fungal endosymbionts within the NGS raw reads, suggesting that we can understand their phylogenetic positions of fungal symbiont with high resolution without the need to isolate the symbiont from the host. Furthermore, our study shows that NGS raw reads of insects generated in the future can be used to pinpoint further fungal endosymbionts that have previously been difficult to identify. This method could provide novel insights into their phylogenetic positions as well as interactions with their host species.



---

*Source: 6652508-2021-06-09.xml*