Population genetic structure of the sheath blight pathogen Rhizoctonia solani AG-1 IA from rice fields in China, Japan and the Philippines

Sheath blight, caused by the fungal pathogen Rhizoctonia solani AG-1 IA is one of the most important rice diseases worldwide. The objetives of this study was to determine the predominant reproductive system and the genetic structure of 18 rice-infecting populations of R. solani sampled from China, Japan and the Philippines, the most important rice production countries in Asia. Knowledge about the population genetic structure of the pathogen in Asia is useful in identifying sources of infection and formulating sustainable management strategies for rice sheath blight. From a total of 717 isolates, 423 unique multilocus genotypes were detected based on nine microsatellite loci. The three country populations of R. solani AG-1 IA exhibited a mixed reproductive system, which included both sexual and asexual components. A moderate degree of clonality indicated that the asexual sclerotia represent important source of inoculum. Population subdivision varied within and among countries, fitting the isolation by distance model. While no subdivision was found among populations within Japan or within the Philippines, subdivision was detected among populations within China. Historic migration indicated high influx of immigrants from Japan into Northern, Central and Eastern China populations. Southern China contributed a high number of immigrants to the populations from the Philippines.


Introduction
The basidiomycetous fungus Rhizoctonia solani Kühn [sexual stage: Thanatephorus cucumeris (Frank) Donk] anastomosis group 1 IA (AG-1 IA) is an extremely important soilborne pathogen causing the rice sheath blight disease, which has a worldwide distribution (Lee & Rush, 1983;Singh, Sunder, & Kumar, 2016). Economically, sheath blight ranks as the second most important rice disease in China, Japan and in Southeast Asia countries. Yield losses due to the rice sheath blight has been estimated as 24 -38 thousands tons per year in Japan but can reach 6 million tons of rice grains per year in China. In the Philippines, rice sheath blight can cause 5 -80% reduction in rice grain yield and the yield losses per hectare ranges from 0.27 -1.29 tons in the dry season to 0.23 -1.37 tons in the wet season (Ou, 1985;Ren, Gao, & He, 2001). So far, no rice varieties with considerable levels of resistance against R. solani AG-1 IA have been developed (Eizenga, Lee, & Rutger, 2002;Jia, Liu, Costanzo, Lee, & Dai, 2009;Zuo, Zhang, Chen, Chen, & Pan, 2010) especially because of the breeding challenges due the quantitative nature of the disease (Liu et al., 2013). Disease management practices with fungicides increases the production costs, are not fully effective and may select for resistance in the pathogen's populations, besides causing environmental impact (Gnanamanickam, 2009;Singh et al., 2016).
On a global scale, the pathogen is genetically diverse, and populations of R. solani AG-1 IA have restricted gene flow, eventually broken by man-mediated long distance dispersal via infected seeds and plant materials (Bernardes de Assis et al., 2009;González-Vera et al., 2010). Because Asia accounts for 92% of the world's total rice production and the majority of rice is produced in China, Japan and the Philippines (Das, 2017), strict quarentine measures should be implemented to restrict the movement of infected seed lots carrying the pathogen among countries, and to ensure a stable food security.
Despite the fact that the genetic structure of contemporaty populations of R. solani AG-1 IA from China has been extensively described, no information is currently available regarding the genetic structure of the pathogen from Japan and the Philippines, the two other most important rice producing countries in the region. A single field study focusing on the phenotypic diversity of R. solani AG-1 IA over a five-year period in Los Baños, Laguna Province, in the Philippines indicated that a single somatically compatible population (SCP) of R. solani AG-1 IA dominated a large upland field planted also with corn, cowpea, soybean and cotton, supporting a high degree of clonality . In contrast with the phenotypic diversity based on somatic compatibility, molecular markers indicated that around 30 haplotypes of R. solani AG-1 IA occurs in the Philippines and also in Japan (Pascual, Toda, Raymondo, & Hyakumachi, 2000). Similarly, information about the historical and current levels of gene flow among country populations of the pathogen is lacking.
Although the field occurrence of hymenia and basidiospores (the fungal sexual reproducing structures) had been reported as early as 1960's in Asia (Fukazu, Kashizaki, & Hirayama, 1960) there is a prevailing perception that the reproductive system of R. solani AG-1 IA is essencially assexual in rice-growing areas in Asia. This is because the role of basidiospores in the fungal disease cycle is not yet fully understood (Rosewich et al., 1999;Bernardes de Assis et al., 2009). However, strong evidence that basidiospores can be an important source of inoculum shaping the genetic structure of the fungal populations towards a mixed reproductive system has been reported in Iran, in the Middle East (Padasht-Dehkaei et al., 2010). In addition, basidiospores of R. solani AG-1 IA have been commonly detected in Asia, in Japanese rice paddy fields (Miyasaka & Nakajima, 2009). Therefore, the main objectives of this study were: (i) to determine the main reproductive system shaping the structure of R. solani AG-1 IA populations from Japan and the Philippines and compare with the structure described for the Chinese populations of the pathogen (Bernardes de Assis et al., 2009); (ii) to determine the current levels of gene flow and differentiation among regional populations of the pathogen; and (iii) to determine the main historical migration path associated with the current genetic struture of regional and country populations of R. solani AG-1 IA. Knowledge about the genetic structure of R. solani populations in Asia will be useful in identifying sources of infection and formulating sustainable disease management strategies for rice sheath blight.

Sampling, fungal isolation and DNA extraction
The isolates were obtained from infected plants sampled from rice-growing areas in Japan (nine populations), and in the Philippines (four populations), from 2008 to 2013. Samples of infected plants with characteristic symptoms of sheath-blight disease were collected from each province, from three to six infected foci at 10 meter apart in each distinct row for a total of 10 to 15 rows per field, totaling 30 to 90 isolates per field population. Isolations were done by transferring pieces of infected leaves to selective medium and incubating at 25 o C in the dark. Pure cultures were obtained by transfer of hyphal tips from the isolation medium to potato dextrose agar (PDA). Sclerotia from young cultures from each isolate were harvested after 10 days incubation in PDA at 25 o C and placed into 1.8 mL cryotubes (Nunc Cryoline System) containing silica gel (Fluka Chemie GmbH, Germany) for long-term storage at 4 o C. Mycelia were harvested after five-day incubation on PD broth at 25 o C and lyophilised for DNA extraction. DNA was extracted as described previously (Linde et al., 2005) and part of the 28S ribosomal DNA (rDNA) region was selectively amplified using primers specific for R. solani AG-1 IA to confirm the fungus identification: R. solani AGcommon primer (forward) 5′-CTCAAACAGGCA-TGCTC-3′ and R. solani AG 1-IA specific primer (reverse) 5′-CAGCAATAGTTGGTGGA-3′ (Matsumoto, 2002). We also included data obtained with the genotyping of isolates from five geographical populations of R. solani AG-1 IA from China, characterized by Bernardes de Assis et al. (2009). In total, our sample size included 717 isolates of R. solani AG-1 IA.

Microsatellite genotyping
Each isolate of R. solani AG-1 IA was genotyped using a specific set of nine polymorphic co-dominant microsatellite (SSR) markers  labeled with fluorescent primers. The PCR products were electrophoresed in an ABI PRISM3100 (Applied Biosystems) automated sequencer along with size standard (GeneScan-500; Applied Biosystems). The binning of the alleles was estimated using GenneMapper Software 5 (Thermo Fischer Scientific). Eight positive controls of standard R. solani DNA spanning the range of microsatellite size were included in each run to ensure consistency of allele scoring for every SSR marker across runs.

Fungal Reproductive system: Genotypic diversity, clonal fraction, Hardy-Weinberg and gametic equilibrium
Based on its multilocus microsatellite genotype (MLMG), each isolate of R. solani AG-1 IA were initially screened for clones using the program GenoDive 2.0b7 (Meirmans & Van Tienderen, 2004). The following indices of clonal diversity were calculated: (i) The number of different multilocus genotypes; (ii) clonal fraction (the proportion of fungal isolates originating from asexual reproduction), calculated as 1-(G/n), where G is the number of different genotypes and is the total number of isolates (Zhan, Pettway, & McDonald, 2003). Using the same program, we performed clone correction, in which only one individual of each MLMG was included per population. All subsequent analyses were performed using clone-corrected data.
A test for Hardy-Weinberg equilibrium (HWE) was performed for each locus with the program Arlequin 3.11 (Excoffier, Laval, & Schneider, 2005). This test is analogous to Fisher's exact test on a two-by-two contingency table but extended to a contingency table of arbitrary size (Guo & Thompson, 1992;Raymond & Rousset, 1995). P values were obtained using a Markov Chain Monte Carlo (MCMC) approach, generating an exact probability distribution not biased by rare alleles (Raymond & Rousset, 1995). Gametic equilibrium (GE) was assessed using a multilocus association test (Brown, Feldman, & Nevo, 1980). The hypothesis that genotypes at one locus are independent from genotypes at another locus was tested using an MCMC algorithm implemented in the program GenePop (Raymond & Rousset, 1995) available at http://genepop.curtin.edu.au. We also measured the multilocus index of association (IA) (Maynard Smith, Smith, O'rourke, & Spratt, 1993) using the program Multilocus v1.3 (Agapow & Burt, 2001). An IA significantly different from zero means disequilibrium. Inbreeding, a possible cause of deviation from HWE and GE, was quantified based on values of fixation index (FIS) calculated using Arlequin 3.11 (Excoffier et al., 2005) with 1,023 permutations of alleles between individuals.

Contemporary gene flow and historical migration
Distribution of genetic variation was assessed by pairwise comparisons of populations using analysis of molecular variance (AMOVA) with the program ARLEQUIN 3.11 (Excoffier et al., 2005). The sum of squared size differences (RST) between two haplotypes was used as the distance measure assuming stepwise mutation model (SMM) (Balloux & Lugon-Moulin, 2002). Significance was assessed using 1,000 permutations. We also tested whether population subdivision followed the isolation-by-distance model proposed by Mantel (1967). Mantel's test was carried out using GenoDive 2.0b7 (Meirmans & Van Tienderen, 2004) assuming a linear relationship between the pairwise genetic differentiation measure RST /(1-RST) and the natural logarithm of geographic distances (kilometers) between all population pairs (Rousset, 1997). The significance of the relationship was accessed with 1,000 permutations.
To determine whether any individuals in a sample were immigrants with respect to their reference geographical population, we used an assignment test performed with a Bayesian statistical model implemented by STRUCTURE 2.2 (Pritchard, Stephens, & Donnelly, 2000). This program calculates the membership coefficients (Q) of every MLMG to each of the populations. The MLMGs were assigned a priori to their reference population. By including this prior information, genotypes in the sample are assigned probabilistically to their reference population, or jointly to two or more populations if their genotypes are consistent with admixture. We performed 10 runs of MCMC simulations, with an initial burn-in of 10,000 followed by 100,000 iterations. Parameters were set using nine clusters (K = 9), admixture model, and RST values previously calculated with Arlequin 3.2 (Excoffier et al., 2005). n For migration analyses, populations were grouped according to the levels of population differentiation inferred from the population structure analyses. Therefore, because no differentiation was detected among regional populations within Japan and within the Philippines, they were grouped into two country populations. Populations from China were grouped into five regional popu lations: Northern China, Southern China, Central China, Western China and Eastern China. Historical migration among these seven geographical populations was assessed by using maximum likelihood test based on MCMC (Beerli & Felsenstein, 1999). A total of seven different migration models between populations were tested: China as recipient of immigrants from Japan (a) or from the Philippines (b); d) Japan as recipient of immigrants from China (c) or from the Philippines (d); the Philippines as recipient of immigrants from China (e) or from Japan (f); g) immigration within China, among the five Chinese geographical populations. Estimates of historical migration were obtained using five runs, and the run with the highest likelihood was chosen to represent each of the migration model. The likelihood values of the eleven migration models were compared based on the log of the Bayes Factor (LBF). LBF was calculated as 2 [ln(Prob(Data | ModelX)) -ln (Prob(Data | best of the eleven models))]. Higher LBF values reflect the best fit of the migration model to the data (Beerli, 2006;Beerli & Palczewski, 2010) (Additional file 1). For all migration analyses the data type chosen was microsatellite data with Brownian motion and assuming SMM. Each of the five runs had ten short initial chains, one long final chain, a static heating scheme with five temperatures (1.0, 100, 1,000, 10,000, and 100,000), and swapping interval of 1. The initial chains were performed with 500-recorded steps, a sampling increment of 100, with 2,500 trees recorded per short sample. The long chain was carried out with 8,334 recorded steps, a sampling increment of 500, six concurrent chains (replicates) and 500 discarded trees per chain (burn-in). The final number of sampled parameter values was 25,002,000 iterations. The values and confidence intervals for the migration rate (M), and the effective population size (θ = 2Neμ for diploids, where Ne = effective population size and μ = mutation rate inferred for each locus) were calculated using a percentile approach (Additional file 2). Migration analyses were implemented in MIGRATE-n v. 3.6.11 (Beerli & Palczewski, 2010) at the CIPRES Science Gateway v. 3.3 (https://www.phylo.org).

Fungal reproductive system: Genotypic diversity, clonal fraction, HWE and GE tests
A total of 423 different miltilocus genotypes were identified among the 717 isolates analysed from the 18 regional rice-infecting populations of R. solani AG-1 IA sampled from China, Japan and the Philippines ( Table 1). The average clonal fractions was not significantly different among countries, varying from 0.38 to 0.51. The majority of the microsatellite loci used for genotyping field populations of the pathogen were in Hardy-Weinberg equilibrium proportions and gametic equilibrium in all the three rice-infecting populations from Asia. This is evidence of random mating and a sexual recombining population structure. However, the proportions of loci under HWE varied among the three country populations (Table 1). The highest proportions of loci under HWE (0.97) was detected in the Japanese populations, which was significantly different from the Chinese populations (0.87). Therefore, most of the loci were in HWE in the rice-infecting populations from Japan and China, while in the Philippines we found HWE in 6.1 out of 9 loci. Significant inbreeding (FIS) was found in two of the five populations from China (Eastern and Western), a single population in Japan (Northeastern), and in all populations from the Philippines. These populations departed from the random mating predictions towards assortative mating. Within China, gametic equilibrium (GE) was detected in a single Chinese population (Northern) (IA = 0.12, p = 0.208), with the proportion of locus pairs at disequilibrium at 10.7%. The remaining four Chinese populations were in gametic disequilibrium, with significant IA values ranging from 0.75 in the Eastern population to 1.12 in the Western (Table 1). Within Japan, the Northeastern and the Western populations were in GE. Significant disequilibrium, with IA varying from 0.42 to 0.66 and percentage of locus pairs at disequilibrium ranging from 3.6 to 14.3% were detected in the other seven populations from Japan. Within the Philippines, all four regional populations were in gametic disequibilibrium (Table 1). In general, there was no significant difference among country populations considering the average IA and the proportion of locus pairs in disequilibrium (Table 1).
Acta Scientiarum. Agronomy, v. 42, e42457, 2020  province; Nueva Ecija (NE) and Tarlac (TA) provinces; Pangasinan (PA). b One or two monomorphic loci in the population were not included in the HWE test. c FIS is the population inbreeding coefficient. d IA is the index of multilocus gametic disequilibrium. e Pairs of loci at significant disequilibrium according to the Fisher exact test (probability test), implemented by GenePop 3.4 (Raymond & Rousset, 1995), after Bonferroni correction for multiple comparisons (Bonferroni, 1935). f Country means followed by the same capital letters are not significantly different by t-test at p < 0.05.
The reproductive system across the 18 populations ranged from completely random mating to varying degrees of inbreeding and gametic disequilibrium. Our overall interpretation was that both asexual and sexual reproductive systems were present in R. solani AG-1 IA populations. Asexual reproduction, recognized as the persistence of clonal lineages within the clone corrected fraction, varied across geographical regions and countries. The evidence for recombination (most of the loci in HWE and gametic disequilibrium) was found in all three country populations, indicating that the sexual cycle likely occurs on rice in China, Japan and the Philippines and thus, basidiospores may play an important role in the epidemiology of the sheath blight, despite the fact that they are considered not easily detectable in nature. This finding is consistent with earlier studies of rice infecting populations from India (Linde et al., 2005), USA (Rosewich et al., 1999) and Iran (Padasht-Dehkaei et al., 2010) showing a mixed reproductive system of the fungus. The mixed reproductive system of pathogens such as R. solani AG-1 IA, that combines sexual and asexual reproduction, may facilitate the rapid emergence of host specialization, virulence and fungicide resistance (Milgroom, 1996;McDonald & Linde, 2002). The amount of inbreeding found within the sexual component of the reproduction of R. solani AG-1 IA most likely reflects the biology of the fungus in that particular environment, rather than a hint of homothallic or heterothallic mating systems (Anderson & Kohn, 1995;Kauserud & Schumacher, 2001).

Contemporary gene flow and historical migration
Evidence of isolation by distance was detected when pairwise genetic differentiation measuress RST/(1-RST) and geographical distances were compared (R 2 = 0.45; p ≤ 0.001, Figure1A), indicating open gene flow among regional populations within country and restricted gene flow across countries. Within China, significant differentiation was detected among the regional Northern -BJ, Southern HN-GX-JX and Western YU populations (RST values ranged from 0.28 to 0.36, p ≤ 0.00033; Table 2). In addition, the Southern population was differentiated from Central -WU and Western -YU populations (significant RST n IS F A I values ≈ 019) ( Table 2). No differentiation was detected among the other population pairs from China (Table  2; (Mantel, 1967) between pairs of 18 rice-infecting populations of Rhizoctonia solani AG-1 IA from China, Japan and the Philippines. B. Historical migration between pairs of populations of the pathogen under the best-fit migration model a,b . a For comparison of models of historical migration based on Bezier approximation scores to the marginal likelihood check Additional file 1. b Bayesian estimates of the migration parameters were calculated with MIGRATE-n v. 3.6.11, using a maximum likelihood test based on the Markov chain Monte Carlo method (Beerli & Felsenstein, 1999;Beerli, 2006;Beerli, 2009;Beerli & Palczewski, 2010). Colored horizontal arrows represent the mean Bayesian estimate of the number of migrants per generation (xNmDonor  Recipient) and the thin horizontal bars are the 95% credibility intervals (C.I.) for each estimate given its a posteriori distribution in parenthesis (Additional file 2). xNmDonor  Recipient is the number of immigrants per generation; where N is the population size, m is the fraction of the new immigrants of the population per generation, and x is an inheritance scalar and x= 2 for diploids. Theta (Θ) estimates and migration parameters estimates for all pairs of regional and country populations are available in Additional file 2. a Distances were computed as the sum of the squared size difference between two haplotypes determined based on nine polymorphic microsatellite loci (Slatkin, 1995); Shaded values indicate significance at p ≤ 0.00033, after Bonferroni correction (Bonferroni, 1935), based on 1,023 permutations.
No subdivision was detected among regional populations within Japan. All populations of R. solani AG-1 IA from Japan were differentiated from Southern China (HN, GX and JX) population (RST values ranging from 0.21 to 0.35, Table 2). With the exception of Eastern -MI and Western -CH and SH, most of the Japanese populations were differentiated from Western China -YU [RST values ranged from 0.13 to 0.24 (Table 2)]. On the other hand, no subdivision was detected among all the Japanese populations and Northern -BJ, Central -WU and Eastern -AN-ZH China.
The regional populations within the Philippines were not differentiated among themselves (Table 2). However, the four populations from the Phillipines were differentiated from most of the populations from China and Japan, with the exception of Southern Philipines and Western China populations, which were not significantly differentiated ( Table 2). The significant RST between pairs of populations from the Philippines and China varied from 0.12 to 0.64, while RST values between pairs of populations from the Philippines and Japan ranged from 0.31 to 0.61.
Because no subdivision was detected among the regional populations within Japan and within the Philippines, these populations were grouped into two country populations for further analyses of admixture and migration.
An assignment test implemented with Structure indicated that 17% of the multilocus genotypes ( Figure  2) sampled across the three country populations were admixed, i.e. immigrants to their given population. Western, Eastern and Central China were the populations with the higher admixture, ranging from 30 to 59% of the genotypes sampled in each given population, mostly associated with another Chinese population of the pathogen (Figure 2). This observation reinforces the levels of contemporary gene flow detected among the regional populations from China (Table 2). Japan was the population that had, in average, the highest number of admixed genotypes (N = 18.5), with the majority (N = 15.5) associated with a probable Central China origin. The admixture proportion in the Japanese population was 14%. No significant admixture was detected in the Philippines (Figure 2). Historical migration analyses supported unidirectional gene flow from Japan into the Northern, Central and Eastern China populations of the pathogen, as well as gene flow from Southern China into the population from the Philippines ( Figure 1B). On the other hand, historical migration among the Chinese populations of the pathogen indicated significant unidirectional migration to Central from both Eastern and Southern China populations ( Figure 1B). The historical exchange of migrants among all the other population pairs were negligible and not significantly different from zero.

B A
The null hypothesis of contemporary unrestricted gene flow over large spatial scales among the three country populations was rejected and isolation by distance was apparent. However, the unidirectional across-country historical migration detected suggested a man-mediated movement of inoculum. In Rhizoctonia diseases, gene flow can occur through the dispersal of fungal sclerotia. Sclerotia could be associated with dispersal over longer distances under special conditions, such as fields of paddy rice that are interconnected by an irrigation network, because the sclerotia can float in the irrigation water (Sumner, 1996). Finally, long distance dispersal from Japan to several Chinese regions and from Southern China to the Philippines could be due to the trading of contaminated rice seeds and planting materials over the past years (Sivalingam, Vishwakarma, & Singh, 2006;Singh et al., 2016). For instance, China has widely adopted the cultivation of the modern japonica rice cultivar Nongken 58 from Japan in the late 1950's, whose cultivation has spanned to 11 million ha from the 1960s to 1980s (Tang, Wei, & Javier, 2004;Tang, Ding, & Bonjean, 2010).
Another explanation for the long-distance dispersal of R. solani could have been mediated by migratory birds such as Chinese egret which breeds on mainland China rice fields and winters in Japan, the Philippines, and other Asian countries (Fasola, Galeotti, Dong, Dai, & Zhang, 2004;Wood et al., 2010).
The transport and long distance dispersal of inoculum facilitated by agricultural practices or migratory birds is an opportunity for pathogens to cross evolutionary boundaries, notwithstanding geographical or ecological barriers (Cunningham, Daszak, & Rodríguez, 2003). In our pathosystem, the transport of infested seeds is likely to have dispersed the pathogen across the three rice-cropping Asian countries.

Conclusion
Populations of R. solani AG-1 IA from rice are characterized by high levels of genetic diversity maintained by a mixed reproductive system and their potential for short and long distance gene flow mediated by infected seed dispersal through agricultural practices or migratory birds, possibly.