Genetic Differentiation and Spatial Structure of Phellinus noxius, the Causal Agent of Brown Root Rot of Woody Plants in Japan

Phellinus noxius is a pathogenic fungus that causes brown root rot disease in a variety of tree species. This fungus is distributed in tropical and sub-tropical regions of Southeast and East Asia, Oceania, Australia, Central America and Africa. In Japan, it was first discovered on Ishigaki Island in Okinawa Prefecture in 1988; since then, it has been found on several of the Ryukyu Islands. Recently, this fungus was identified from the Ogasawara (Bonin) Islands, where it has killed trees, including rare endemic tree species. For effective control or quarantine methods, it is important to clarify whether the Japanese populations of P. noxius are indigenous to the area or if they have been introduced from other areas. We developed 20 microsatellite markers from genome assembly of P. noxius and genotyped 128 isolates from 12 of the Ryukyu Islands and 3 of the Ogasawara Islands. All isolates had unique genotypes, indicating that basidiospore infection is a primary dissemination method for the formation of new disease foci. Genetic structure analyses strongly supported genetic differentiation between the Ryukyu populations and the Ogasawara populations of P. noxius. High polymorphism of microsatellite loci suggests that Japanese populations are indigenous or were introduced a very long time ago. We discuss differences in invasion patterns between the Ryukyu Islands and the Ogasawara Islands.


Introduction
Phellinus noxius (Corner) G. Cunn. (Hymenochaetaceae) is a pathogenic fungus that causes brown root rot disease in a variety of tree species [1][2][3][4][5][6][7]. The fungus is distributed in tropical and sub-tropical regions in Southeast and East Asia, Oceania, Australia, Central America and Africa [6,[8][9][10][11][12][13][14]. Infection causes slow and reduced growth in trees, discolouration and wilting of leaves, defoliation, and dieback of branches [14,15]. Most affected trees eventually die, and in some cases, the fungus causes the rapid wilt and death of the tree within a few months of infection [8,11,12]. The host range of the fungus is very wide [1,6], showing little host specificity [5,7,16], and to date more than 200 woody plant species representing 59 families have been recorded as host plants [12]. The life cycle of P. noxius is similar to that of other important forest pathogens, such as Phellinus sulphurascens Pilát that causes laminated root rot of conifers [17] and Armillaria spp. that cause Armillaria root rot of woody plants [18]. The fungus infects host trees via root-to-root contact from adjacent infected trees, or from wood debris of dead trees where P. noxius can persist saprophytically more than ten years [19]. Basidiospores may function to establish new disease foci [1,14], but remain undocumented in P. noxius.
In Japan, brown root rot was first found in windbreaks composed of Casuarina equisetifolia L. on Ishigaki Island in Okinawa Prefecture in 1988 [20]. Since then, the disease has gained increasing attention as it has appeared on several islands of the Ryukyu Islands in both Okinawa and Kagoshima Prefectures, causing serious problems for shade, windbreak, and ornamental or landscape trees [14,[21][22][23]. Amami-Oshima Island currently represents the northernmost distribution point of the disease [23]. In 2012, this fungus was identified on the Ogasawara (Bonin) Islands, oceanic islands located approximately 1,000 km south of Tokyo, where it killed trees, including rare species endemic to the islands (Sahashi et al. personal communication).
Phellinus noxius is suspected to be indigenous to many tropical or subtropical areas throughout the world [8]; however, whether P. noxius in Japan is indigenous or has been introduced from other areas remains unclear. Phellinus noxius on the Ogasawara Islands was possibly introduced from the Ryukyu Islands, as several tree species, including Bischofia javanica Blume and Pinus luchuensis Mayr were introduced to these islands from the Ryukyus as timber or fuel trees in the early 1900s [24]. To establish effective control or quarantine methods for brown root rot, it is important to first determine whether the Japanese populations are indigenous to the area or introduced from other areas.
Simple sequence repeats (SSRs) or microsatellites are a group of DNA sequences with repeating units of 2-6 base pairs (bp) that are abundant in most genomes exhibiting high levels of polymorphism [25,26]. Hence, SSRs are useful molecular markers for analysing genetic diversity and have recently been used for a robust assessment of population structure in various plant pathogens [27][28][29][30]. In this study, we developed microsatellite markers based on the de novo sequencing assembly of a Japanese isolate and then analysed genetic diversity or genetic structure in P. noxius in Japan. Although population genetics studies using microsatellite markers have been conducted for other similar root rot pathogenic fungi distributed in cool-or warm-temperate areas, including Armillaria spp. [31][32][33] and Heterobasidion spp. [29,34], this study is the first to examine the population genetics of wood-decay and tree pathogenic fungi in a tropical or subtropical area.

Isolates
The isolates of Phellinus noxius used in this study are listed in Table 1. We collected infected root samples or basidiocarps (for isolate KPN246 only) from 12 of the Ryukyu Islands in 1990-2010 [14] and from 3 of the Ogasawara Islands in 2012-2013 (Fig 1). Phellinus noxius was isolated using the methods described in Sahashi et al. (2012) [14]. Isolates were cultured on potato dextrose agar (PDA; Nissui, Tokyo, Japan) in test tubes and were maintained at 25°C using periodical subculture at the Forestry and Forest Products Research Institute (FFPRI, Tsukuba, Japan). Unless two isolates from the same region were found genetically incompatible, one isolate per disease foci was used for subsequent analysis. Moreover, one isolate (P919-02W.1) of P. noxius from Pohnpei Island, Federated Stated of Micronesia, isolated by Y. Ota and N. Sahashi in 2013 was included in the development of the microsatellite marker to guarantee the robustness of the markers for future worldwide analyses.

DNA extraction
Fungal DNA was extracted from mycelia as described in Ota et al. (2014) [35]. Cultures were grown in 10 mL MYG medium (2% malt extract, 0.2% yeast extract, and 2% glucose) at 25°C in the dark and were harvested 7 days after inoculation. DNA was extracted from frozen mycelia using a DNeasy Plant Mini kit (Qiagen, Valencia, California) according to the manufacture's instructions after grinding mycelia into a fine, dry powder using a mortar and pestle in liquid nitrogen. Genetic Sifferentiation of Phellinus noxius in Japan protocol to produce 100 bp paired-end reads. Assemblies of P. noxious genome sequences were constructed from Illumina reads using an MaSuRCA assembler [36] with the following options: GRAPH_KMER_SIZE = auto, ovlMerSize = 30, cgwErrorRate = 0.15, utgErrorRate = 0.015, and KMER_COUNT_THRESHOLD = 1. MISA (http://pgrc.ipk-gatersleben.de/misa/) was used to identify di-to tri-nucleotide microsatellite loci from the genome assemblies of P. noxius KPN92 with at least eight repeats of di-and tri-nucleotides and the maximum number of bases between two microsatellite loci set to 100 bp. The total number of microsatellites identified was 334 (232 di-nucleotide microsatellite with 8-33 repeats and 102 tri-nucleotide microsatellite with 8-23 repeats). Specific primer pairs to amplify those microsatellite loci with four classes of product size were designed using Primer3 2.3.6 (http://primer3.sourceforge.net/) with the following options: the ranges of product size were 100-200, 200-300, 300-400, and 400-500. Within designated primer pairs, 50-60 pairs for each of the four product size classes (220 pairs in total) were selected arbitrarily and synthesised with the tail sequence on the 5' end of the forward primer: tail A (GCC TCC CTC GCG CCA) for product size classes of 100-200 and 300-400 and tail B (GCC TTG CCA GCC CGC) for product size classes of 200-300 and 400-500 [37]. These 220 primer pairs were tested for amplification on eight P. noxius isolates: KPN1, 21, 92, 121, 145, 157, 256, and P919-02W.1 for the first screening. PCR amplification was performed using a BIO-RAD iCycler (Hercules, California) with the following conditions: 5 min at 95°C followed by 35 cycles of 30 s at 95°C, 30 s at 55°C, and 45 s at 72°C, and a final extension of 10 min at 72°C. Each PCR reaction contained approximately 5 ng template DNA, 0.2 μM of each primer, and 1X Go Taq Green Master Mix (Promega, Madison, Wisconsin) in 25 μL total volume. PCR products were separated by electrophoresis on a 1% agarose gel in TAE buffer and visualised using ethidium bromide staining on a UV transilluminator.

Microsatellite marker development
For the second screening, 20 primer pairs in each of the four product size classes (80 primer pairs in total) were arbitrarily selected from the primer pairs that generated clear PCR products in all eight isolates in the first screening; these were combined into 20 multiplex PCR panels that included one primer set of each of the four classes. To test the amplification of multiplex PCR panels and the polymorphism of each microsatellite loci, 15 P. noxius isolates were used: the 8 isolates used in the first screening, plus KPN28, 80, 132, 141, 280 and 289. PCR amplifications using a Qiagen Multiplex PCR Kit were performed with a BIO-RAD iCycler following the manufacturer's recommended conditions: 15 min at 95°C followed by 35 cycles of 30 s at 94°C, 90 s at 55°C, and 60 s at 72°C, and a final extension of 30 min at 60°C. Each PCR reaction contained approximately 5 ng template DNA, 0.1 μM forward primer, 0.2 μM reverse primer, 0.2 mM of each universal primer labelled with fluorescent dye (Tail A with 6-FAM and Tail B with VIC, [37]), and 5 μL Master Mix (from the kit) in 10 μL total volume. Amplified products were loaded on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA), and genotype scoring was performed using the GeneScan 600 LIZ dye size standard (Applied Biosystems) and GeneMapper version 4.1 software (Applied Biosystems). Finally, we selected 20 primer sets in five multiplex PCR panels for further analysis of the 128 Japanese isolates of P. noxius.

Genetic diversity analysis
We defined isolates from the same island as a "population." From the genotyping data, 102 isolates were observed as diploid, whereas 26 isolates were haploid. The number of alleles at each microsatellite locus was calculated using the program GenAlEx version 6.5.0.1 [38] for all 128 isolates. For further analysis of genetic diversity or genetic structure, haploid isolates were excluded. The number of alleles, Shannon's information index (I), observed heterozygosity (Ho), expected heterozygosity (He), and Nei's unbiased expected heterozygosity (uHe) for each locus were calculated using GenAlEx. Deviation from Hardy-Weinberg equilibrium (HWE) for each locus and linkage disequilibrium between loci were tested using the program Fstat version 2.9.3.2 [39] under the infinite allele model (IAM), and multiple testing with the Holm-Bonferroni method [40] was performed. Weir and Cockerham's estimate of F IS [41] was calculated using Fstat. These analyses were conducted among populations (each island) as well as between the two groups of populations (the Ryukyu Islands and the Ogasawara Islands) inferred from the STRUCTURE analysis.

Genetic structure analysis
A Bayesian-based clustering method was applied to infer the genetic structure of Japanese P. noxius isolates using the program STRUCTURE version 2.3.4 [42]. An admixture model with correlated allele frequencies assuming no prior information of population origin was used. Twenty independent runs for K = 1 to 10 were performed at 100,000 Markov Chain Monte Carlo (MCMC) repetitions after a burn-in period of 50,000 iterations. The appropriate number of clusters (K) based on the ad hoc statics ΔK was determined using the method of Evanno et al. (2005) [43] with the program Structure Harvester [44].
Subsequently, analysis of molecular variance (AMOVA) was performed using the program Arlequin version 3.5.1.3 [45] to calculate the hierarchical distribution of genetic variation in Japanese isolates. All populations were initially combined into one hierarchical group, and then divided into two groups (the Ryukyu Islands and the Ogasawara islands) based on STRUCTURE analysis. The significance of components that showed variance was tested by performing 9,999 permutations.
Finally, the relationship between genetic structure and the isolation-by-distance (IBD) model was tested [46]. Values of F ST between populations were calculated using Fstat. Mantel's test was performed using GenAlEx with 9,999 permutations and assuming a linear relationship between pairwise values of F ST /(1-F ST ) and the natural logarithm of geographic distance (km) between all population pairs [47]. The central value between the maximum and minimum latitude and longitude of the isolates on the same island was used as the location of the population.

Characteristics of microsatellite markers
Of the 220 microsatellite primer pairs designed from the assembly of the P. noxius genome (isolate KPN92), 20 primer pairs in five multiplex panels were selected for use in the population analysis ( Table 2). The 20 microsatellite markers were distributed in 19 distinct scaffolds in the genome assembly and exhibited high polymorphisms at each locus ( Table 3). The sequences of these microsatellite loci inferred from the KPN92 genome assembly have been deposited in DDBJ (accession numbers are shown in Table 2). The number of alleles at each locus was 21.7 on average and ranged from 7 at Pn155 to 45 at Pn111 (Table 3).
In all, 4 to 22 isolates of P. noxius from each island were tested using 20 microsatellite markers, and remarkably, all 128 isolates exhibited different genotypes. A total of 102 multilocus genotypes were interpreted as being diploid from the microsatellite analysis; however, 25 of 83 isolates from the Ryukyu Islands and 1 of 45 isolates from the Ogasawara Islands were judged to be haploid, because one single allele was detected at all loci in these isolates. Only diploid isolates were used for further analysis of genetic structure and diversity. In addition, the Yoron Island population was excluded because all of its isolates were haploid.

Genetic diversity
A summary of the genetic diversity for the 14 populations is presented in Table 4. Expected heterozygosity across all populations was 0.67 (±0.01 SD), ranging from 0.48 (±0.07) on for the three populations in the Ogasawara Islands. There was a low number of isolates for some populations; therefore, for further analysis, each population was combined into two groups of islands (the Ryukyu Islands and the Ogasawara Islands) based on the STRUCTURE analysis. For the Ryukyu Islands (N = 58), the average observed and expected heterozygosity at each locus was 0.648 ± 0.043 and 0.832 ± 0.031, respectively. F IS at each locus ranged from 0.004 to 0.824, and significant deviation from HWE was detected at 9 of 20 loci after sequential Bonferroni correlation (α = 0.05). For the Ogasawara Islands (N = 44), the average observed and expected heterozygosity at each locus was 0.495 ± 0.053 and 0.623 ± 0.059, respectively. F IS at each locus ranged from -0.001 to 0.594, and significant deviation from HWE was detected at three loci after sequential Bonferroni correlation (α = 0.05). Allelic richness was 16.5 ± 1.7 in the Ryukyu Islands and 7.9 ± 1.3 in the Ogasawara Islands. No significant linkage disequilibrium was detected between each locus in any population after sequential Bonferroni correlation (α = 0.05). Na, Na di , Ho, He, and Rs refer to as the total number of alleles per locus in all isolates, the total number of alleles per locus in diploid isolates, the observed heterozygosity, the expected heterozygosity, and allelic richness respectively. F IS was calculated by Weir & Cockerham. * indicates that the HWE test is significant after the Holm-Bonferroni correction method (α = 0.05).
-indicates that Fis was not calculated because the loci was monomorphic.

Genetic structure
Evanno's method using Structure Harvester clearly indicated that ΔK at K = 2 was at a maximum and two was an appropriate number of clusters (Fig 2). The two clusters clearly exhibited structure between isolates of the Ryukyu Islands (composed of 11 islands from Amami-Oshima Island to Iriomote Island) and those of the Ogasawara Islands, containing the islands of Anijima, Chichi-jima, and Haha-jima (Fig 3). When all populations were combined in one hierarchical group, AMOVA analysis indicated that most of the genetic variation could be explained by differences in individual isolates within populations (85.52%) rather than by variation among populations (14.48%, P < 0.0001, Table 5). When the populations were partitioned into two groups (the Ryukyu Islands and the Ogasawara Islands) established from the STRUCTURE analysis, most of the genetic variance could be explained by differences in individual isolates within a population (79.18%, P < 0.0001). Differences in isolates among groups and among populations within groups explained 16.65% (P = 0.0029) and 4.17% (P < 0.0001), respectively.
A pairwise analysis of IBD among populations indicated a significant positive correlation between genetic distance and geographic distance (R 2 = 0.47193; P = 0.003, Fig 4). High F ST /(1-F ST ) values above 1200 km (= ln7.1 km) were consistent with the pairwise analysis between the Ryukyu Islands and the Ogasawara Islands.

Discussion
We developed 20 microsatellite markers for P. noxius. Multiplex PCR for these markers (four of each primer pair in one reaction) successfully genotyped Japanese isolates as well as an isolate from Pohnpei Island, Federated States of Micronesia, indicating that these markers are useful to genotype isolates from other geographic regions. Moreover, these markers showed enough polymorphism to analyse the genetic or clone composition of P. noxius in local populations. For some Ryukyus Islands, we were unable to obtain a sufficient number of isolates to analyse differences between islands, because the number of sites where the disease occurred was too low and some isolates were haploid. However, the number of isolates was sufficient for comparisons between the Ryukyu Islands and the Ogasawara Islands. Diploidy is the main ploidy for vegetative hyphae in Hymenochaetaceae, which is consistent with our microsatellite data for P. noxius. However, 26 of 128 isolates obtained from decayed woods or basidiocarps were judged to be haploid. Among the isolates from the Ryukyu Islands obtained between 1999 and 2014, 30.1% were haploid, whereas only one isolate was haploid among those isolated from the Ogasawara Islands during or after 2012 (Table 1). During the maintenance of the isolate cultures, we encountered cases in which some sub-cultures from a diploid isolate showed haploid microsatellite signals that had only one of two alleles of each microsatellite loci of diploid isolates. In Pholiota nameko, an edible basidiomycetous fungus, diploid mycelia often become haploid during storage via a mechanism known as monokaryotisation or dedikaryotisation [48,49]. Haploid mycelia of P. noxius may occur as primary mycelia that are derived from germinated basidiospores, however, basidiocarp formation of this fungus is very rare [3,12] and the primary mycelia are usually short lived [50]. Therefore, haploid P. noxius isolates in this study may have changed from diploid during periodical subculturing by monokaryotisation.
Phellinus noxius has two dissemination methods: asexual root-to-root contact from a diseased tree to a living tree and dispersal of sexually produced basidiospores [1,14]. Using somatic incompatibility tests, Hattori et al. [51] examined the clone distribution of P. noxius in windbreak trees on the Ishigaki Islands of Japan. They concluded that infection via both  basidiospores and root-to-root contact occurred in the area. Some researchers, however, have suspected that infection by basidiospores is rare because basidiocarps are seldom seen in areas of disease propagation [3,12]. We found that all of the isolates exhibited unique genotypes, strongly indicating that basidiospore infection is the main dissemination method for the formation of new disease foci. In the Ryukyu Islands, basidiocarps were rarely seen in areas where the disease was spreading and forming forest gaps [14]; however, they were occasionally seen on dead or fallen trees in natural forests, where the disease was not spreading. Meanwhile, basidiocarps are more frequently observed on the Ogasawara Islands, although the reason is unclear (Hattori personal observation). In such areas, basidiospores may function to produce new disease foci. Although many unique genotypes have been observed, it is possible that a small number of genotypes dominates within a single disease focus, because P. noxius spreads asexually within disease foci like Phellinus sulphurascens and Armillaria spp [1,15,17,18]. More extensive sampling within disease foci are needed to clarify the clone distribution pattern of P. noxius.
Whether P. noxius is indigenous in Japan or introduced from other areas is unknown. Because brown root rot was first recognised in Japan as recently as the 1980s, the possibility of P. noxius as an introduced pathogen has been expected [23]. Ann et al. [12] suggested that P. noxius was likely introduced to Taiwan on diseased roots of exotic trees, based on observations that the distribution of P. noxius in Taiwan is limited to areas of human activity and the disease has never been found in undisturbed forests. In general, introduced pathogens have lower genetic diversity than indigenous pathogens because of the founder effect of small population sizes and subsequent bottlenecks [52]. Introduced diseases that have had devastating effects include chestnut blight caused by an ascomycetous fungus Cryphonectria parasitica (Murrill) M.E. Barr [27], ash dieback by Hymenoscyphus fraxineus (T. Kowalski) Baral, Queloz & Hosoya [53], sudden oak death by an oomycete Phytophthora ramorum Werres, De Cock & Man in't Veld [54], and alder decline due to P. alni Brasier & S.A. Kirk [28]. Population genetics studies using microsatellite markers have indicated that the genetic diversities of these species in the area of introduction are lower than in native areas [27,28,53]. In terms of rootrotting basidiomycetous fungi, Heterobasidion irregulare Garbel. & Otrosina in Italy [29] and Armillaria mellea (Vahl) P. Kumm in South Africa [55] are known as introduced pathogens. Heterobasidion irregulare in Italy, which is suspected to have been introduced by the US military during World War II, exhibits fewer alleles (1-7) at each microsatellite locus than native populations in North America. In our study, Japanese P. noxius isolates exhibited a high number of alleles per loci (21.7 on average), suggesting that P. noxius is indigenous to Japan or was introduced to the country a very long time ago. Further studies using isolates collected from other geographic region are needed to confirm the conclusion. Although the occurrence of brown root rot in Japan was only recognised recently (i.e., in the 1980s on the Ryukyu Islands and the 2010s on the Ogasawara Islands), basidiocarps of P. noxius were recorded on a broadleaved tree on the Ogasawara Islands in 1916 [56]. This suggests that P. noxius was present on these islands without causing a conspicuous decline of resident trees. The causes of the recent outbreak of this disease in Japan has not yet been determined, although several environmental changes, including irregular climatic events such as typhoons and droughts, as well as human disturbances may have contributed to the outbreak.
The STRUCTURE analysis strongly indicated genetic differentiation between the Ryukyu and Ogasawara populations of P. noxius. Additionally, the AMOVA and IBD analysis also supported the conclusion. These findings suggest minimal gene flow between the two island chains over a long period of time or a different origin of the two populations. The Ryukyu Islands and Taiwan are continental islands that were once connected to the Eurasian continent; thus, P. noxius was able to spread to and from the continent similar to other flora and fauna [57]. In contrast, the Ogasawara Islands are oceanic islands formed by volcanic activity and were never connected to a continent or other larger islands such as the main Japanese islands. Therefore, the origins of all flora and fauna on them are thought to be introductions from other continents or islands followed by their unique evolution. Because the dispersal modes for plants on the Ogasawara Islands are by air, bird, and oceanic drift [57], P. noxius was probably introduced via one of these methods. In general, basidiospores are ephemeral, and the majority of basidiospores fall within a short distance of the basidiocarps [58,59]. However, long-distance dispersal (1000 km) has also been reported in some wood-inhabiting basidiomycetous fungi [60]. The basidiospores might have been introduced from the Mariana Islands, the nearest oceanic islands to the Ogasawara Islands, by a typhoon, as many typhoons form around the Marianas and move to the Ogasawaras.
Phellinus noxius could serve as a suitable model for studying the evolutionary history of fungi and forest diseases on oceanic islands, because it is distributed in three geographically different categories: continents, continental islands, and oceanic islands. Further population genetics studies using isolates collected from around the world will be useful for understanding the evolutional history of P. noxius and its worldwide routes of dispersal.