Human-Activity-Mediated Gene Flows Are Involved in the Establishment of Magnaporthe oryzae Population Structures: A Population Genomic Analysis

1. Introduction

Rice blast is a devastating disease that occurs in rice-growing areas worldwide, it is caused by the filamentous ascomycete fungus Magnaporthe oryzae and poses a serious threat to global food security [1]. M. oryzae can rapidly overcome the resistance genes in rice and can coexist with resistant varieties within a few years of their initial deployment in rice agrosystems [2,3]. The pathogen infects rice throughout the growth period and also has a wide host range that includes more than 50 cultivated and wild monocot plants, such as rice (Oryza sativa), barley (Hordeum vulgare), wheat (Triticum aestivum), finger millet (Eleusine coracana), goosegrass (Eleusine indica), perennial ryegrass (Lolium perenne) and more [4]. The asexual reproduction is prevalent in most rice fields, resulting in local populations of M. oryzae often exhibiting only one mating type[5,6,7,8,9].

It is necessary to better understand the features of the genetic structure of M. oryzae in order to provide the genetic basis for developing sustainable and effective prevention and control strategies for rice blast disease [10]. Genetic diversity can reflect the survival ability and adaptive potential of natural populations in the face of rapidly changing biotic and abiotic backgrounds [11]. For the plant pathogen M. oryzae, genetic diversity is usually estimated using molecular markers, including microsatellite markers, or DNA fingerprinting via amplified fragment length polymorphism. These methods have been applied in several countries and continents, including China [12,13,14], India [15,16], Thailand [17], the Philippines [18], Africa [19], Europe [20] and America [21], and have revealed regional diversity variations underling local adaptations [22]. The establishment of genetic structures can be driven by geographic isolation and ecological factors. In a previous population genomic analysis, three main genetic clades were identified in M. oryzae, which were associated with the distribution of mating types [23,24]. In a worldwide population structure analysis of M. oryzae, multiple endemic and pandemic lineages were identified, which were distributed in specific rice-growing areas [25]. Based on the amplified fragment length polymorphism, very few genetic differences have been found between the geographically distant populations of M. oryzae in Iran and Uruguay, although evidence of gene flow has been observed [26]. As a broad host pathogen, host specialization has also been reported to be a non-negligible factor in the genetic differentiation of M. oryzae [27,28]. However, the limited genetic markers have hindered population genetic analysis from accurately identifying the comprehensive ensemble population structures of rice-infecting lineages due to insufficient information regarding the incomplete divergent populations. In particular, there remains a scarcity of studies addressing the spatiotemporal dynamics of genetic diversity and population structures of M. oryzae through large-scale genomic analysis.

Gene flow, or migration, implies the movement of genetic material among spatially or temporally separated populations and acts as a major driving force for organisms to establish population structures. It can accelerate novel variations through gene recombination, migration to found new populations and takeovers of other local populations. Island [29] and stepping-stone models [30] have indicated that the migration of one or more subpopulations can decrease the genetic correlation with geographical distance among sexually reproducing species. For plant pathogens, spatial movement occurs in many forms, such as short-distance transfer through rainwater and long-distance global transfer, which counteracts the disadvantage of immobility for large-scale epidemics. Population structures are generally established by genetic variations and natural factors, as well as anthropogenic changes, especially among crop pathogens that are closely linked to human survival. Therefore, human movement significantly influences the migration of plant pathogens [31,32].

In agricultural ecosystems, human movement can aid the extensive spread of plant pathogens by breaking inherent mobility limitations or geographic barriers, such as rivers and mountains, that can restrict population expansion, thereby contributing to the complexity of population structures. Therefore, the globalization of agricultural products, largely driven by human activities, has led plant pathogens to evolve toward metapopulation formation [31]. For example, the transportation of infected seeds is one of that ways that M. oryzae has migrated, demystifying the close genetic relationships that have arisen between geographically separate populations [33]. Through the genome sequence analysis of global M. oryzae populations over different time periods, a sexual recombination signature was detected in the Southeast Asian endemic lineage, suggesting the occurrence of gene flows among geographic population distributions [25]. Based on microsatellite markers, M. oryzae has been found to have weak geographic structures in three islands groups with limited natural migration, which was induced by the transportation of infected seeds around the Philippines [18]. Based on amplified fragment length polymorphism analysis, frequent gene flows have been discovered in East African populations [34], as well as between the different provinces in Korea [35]. This has also occurred with other plant viruses; for example, turnip mosaic potyvirus spread from west to east regions in Eurasia due to historical trade arteries, such as the Silk Road [36].

In this study, we first collected the published whole-genome sequences of 189 rice-infecting M. oryzae isolates from five continents and 22 re-sequenced genomes from China. Then, we analyzed the genetic diversity and population structure characteristics to evaluate the population divergence driven by spatiotemporal changes. We found that a M. oryzae population from Asia exhibited the highest genetic diversity and could be divided into three divergent clades. The diffusion route of M. oryzae followed human activity, suggesting that the Asian population served as the genetic pool for rice-growing regions worldwide, with China being the migration origin within Asia. Our study provides a detailed understanding of the spatiotemporal dynamics of the genetic diversity and population structures, which could be useful for developing cultivars with different resistant genes and could contribute new insights into disease management.

2. Materials and Methods

2.1. Isolate Sampling, DNA Preparation, Genome Sequencing and Published Genome Collection

We collected 22 rice-infecting M. oryzae isolates from rice fields in China. The sampling locations are described in Table S1. Genomic DNA was extracted from fresh mycelium cultures after single-spore isolation, as previous described [37]. An Illumina paired-end DNA sample prep kit was employed to construct the Illumina fragment libraries and an Illumina HiSeq2500 instrument was used for the sequencing.

A total of 189 M. oryzae genomes have been sequenced from rice, which we have downloaded from the NCBI database using the ncbi-genome-download script (https://github.com/kblin/ncbi-genome-download). The genome information is also listed in Table S1.

2.2. Identification of Single-Nucleotide Polymorphism

To obtain high-quality single-nucleotide polymorphism sites from the sequencing reads, the raw data were first trimmed using Trimmomatic software [38]. The eighth version genome of 70-15 was taken as the reference genome and the paired-end reads were aligned against it using the BWA v0.7.17 as the default parameters [39,40]. Then, single-nucleotide polymorphism (SNP) mining was implemented using the Genome Analysis Toolkit (GATK v4.1.4.1), which compared the reference genome to the standard settings [41]. After that, VCFtools was used to filter the biallelic site with minor allele count smaller than 3 and missing data lower than 20% [42].

For the genomes downloaded from the NCBI database, genome-to-genome SNP identification was performed using the -maxmatch and -c 100 options in MUMMer [43] against the reference genome.

2.3. Characterization of Genetic Diversity and Population Genetic Structures

To estimate the genetic diversity of M. oryzae populations that were divided by spatiotemporal boundaries, the nucleotide diversity of the whole-genome SNPs was calculated using the “PopGenome” v2.1.6 package with a 20 kb sliding window and 2 kb step size [44]. The statistical significance between each pair of geographical populations was estimated using Student’s t-test (two-tailed). A SNP-based neighbor-joining phylogenetic tree was reconstructed according to the bitwise distance with 1000 bootstrap replicates using the “poppr” packages in R [45]. A principal component analysis (PCA) of the genetic structures of global populations from Asia, Europe, North America and South America was implemented on the matrix of binary allele sites using the “adegenet” package [46].

To establish the population structure compositions, the SNP dataset was pruned based on linkage disequilibrium (LD) with a value of r2  = 0.5 using Plink v.1.9 (https://www.cog-genomics.org/plink/). To evaluate the individual ancestry components and admixture proportions, a population admixture analysis was performed using the ADMIXTURE v1.3.0 program [47]. In ADMIXTURE, 5-fold and 10-fold cross-validation procedures with K values between 2 and 10 were implemented to identify the optimal number of population clusters (K) with the lowest cross-validation error values. The admixture proportions of each sample from the geographic populations were visualized using PopHelper [48].

2.4. Phylogeographic Analysis of M. oryzae

To detect the reticulation events of the populations, such as hybridization, horizontal gene transfer and recombination, the phylogenetic networks were reconstructed using the SplitsTree algorithm [49]. In order to assess the robustness of the tree topology, 1000 bootstrap replicates were performed using the “Neighbor-Net” and “Uncorrected p-distance” parameters.

To ascertain the phylogeographical dispersal routes (including migration routes and directionality) of M. oryzae across the world and in Asia (four regions), we employed the asymmetric discrete trait substitution model with the Bayesian stochastic search variable selection approach [50] using BEAST 1.10.4 [51]. Based on the Bayesian information criterion, the best nucleotide substitution model was estimated using jModeltest [52] after converting the SNP dataset into the Phylip format using vcf2phylip v2.0 [53]. The posterior distributions of the parameters were estimated using two independent Markov chain Monte Carlo (MCMC) methods for 100 million generations, with sampling every 20,000 generations. Tracer v1.7 was used to check model convergence (effective sample size >200) after excluding 10% of the initial samples as burn-in [54]. Using the SpreaD3 v0.9.6 software, a Bayes factor analysis (BF >3) was performed to identify diffusion routes that were strongly supported by the data [55].

3. Results

3.1. Spatiotemporal Dynamics of Genetic Diversity in M. oryzae Populations

After the SNP mining and filtering, a final SNP dataset of 605,369 whole-genome SNPs was obtained from the samples, which was applied for the evaluation of the genetic diversity of each geographical population of M. oryzae. A nucleotide diversity statistical analysis revealed the variations among the populations in Asia, Europe, Africa, North America and South America, indicating an average nucleotide diversity index (Pi) of 553.57–774.30 for each analysis window (Figure 1). Additionally, according to the t-test statistic between each of the pairwise populations, the M. oryzae samples from Asia had the highest nucleotide diversity (774.30). By contrast, the populations in Europe had the lowest nucleotide diversity (553.57), while those in Africa, South America and North America had moderate diversity.

3.2. Genetic Structures of Worldwide Geographic Populations

A neighbor-joining (NJ) tree was constructed for the global populations of M. oryzae based on the whole-genome SNPs, which showed that the worldwide samples formed three distinct sub-clades (Figure 2A). Out of these clades, the samples from Asia were distributed into three sub-clades. However, for the populations from Europe, North America, South America and Africa, although there were a few isolates that were localized in other clades, the majority of the individual isolates were assigned to a single clade. Therefore, there was no obvious segregation of individuals based on geographical location. The results of the PCA corresponded to the phylogenetic tree topology, indicating that individuals were grouped to three clusters with incomplete separation among the geographic populations (Figure 2B).

Furthermore, population structure composition analysis using ADMIXTURE revealed the individual ancestry and admixture proportions within each population. The optimum cluster number was K = 3 (Table S2), as it had the lowest cross-validation error and was thus selected as the number of potential ancestors for the worldwide M. oryzae populations (Figure 2C). Among them, the genetic compositions of the populations from Asia, Africa, Europe and North America contained three ancestral lineages, although the individuals were dominated by a single ancestral lineage. However, some samples contained two ancestral lineages. Notably, the European populations consisted of three independent ancestral lineages, without any mixed multilineage isolates. In addition, the South American populations were descended from two independent lineages.

3.3. Migration Patterns of M. oryzae

The genetic differentiation between populations (FST) ranged from 0.088 to 0.237 (Figure 3A, Table S3) with the 20 Kb window length and 2 Kb steps, indicating moderate genetic differentiation. This result implies the likelihood of M. oryzae migration occurring among the worldwide geographic populations. At the same time, the results from the phylogenetic network reconstruction exhibit the reticulate relationships within the three clades, suggesting the occurrence of spatial migration (Figure 3B). In addition, the migration pathways of this pathogen that were inferred from the phylogeographic analysis revealed that migration events could be detected from the spatial diffusion of M. oryzae across the world and Asia. The results also show that five major migration events occurred across five continents. Of these, Asia was the major source gene pool, which then flowed to other continents (Figure 4A). Interestingly, Europe served as a pivotal transfer station, facilitating gene flow to North America. In Asia, five major migration events were detected among China (including Taiwan and the mainland), Japan and South Korea (Figure 4B). Mainland China acted as the primary source of gene migration, initially spreading to Japan and South Korea and subsequently extending to Taiwan via Japan. Additonally, bidirectional gene flows were observed between South Korea and Japan.

3.4. Demographic History of M. oryzae Populations

To investigate the demographic history of M. oryzae, effective population sizes were calculated using a pairwise sequentially Markovian coalescent (PSMC) model. As previously reported, rice-infecting isolates of M. oryzae evolved from grass-infecting isolates around 1000 years ago [23]. Therefore, the effective population sizes were predicted for the last 800 years (Figure 5). Our finding indicated that M. oryzae underwent a striking population expansion across five continents approximately 100 years ago.

4. Discussion

For fungal pathogens of plants in agroecosystems, genetic diversity can reflect their capacity to adapt and survive under biotic and abiotic stresses [56,57]. This diversity is closely linked to their population origin and evolution, influencing their formation and maintenance [58,59]. In this study, Asian populations of M. oryzae exhibited the highest nucleotide diversity based on whole-genome polymorphism, in contrast to populations from Africa, Europe, North America and South America. This suggests that Asian populations could serve as potential reservoirs of M. oryzae variants containing more novel genotypes. Furthermore, among the three Asian countries studied, including Taiwan island and mainland of China, Japan and South Korea, the highest nucleotide diversity was observed in Chinese populations. This suggests that China serves as a center of genetic diversity of M. oryzae in Asia, as indicated by the microsatellite markers [24]. Based on its history of rice domestication, China is also deemed to be the center of origin for rice [60,61], meaning that M. oryzae would have accumulated more abundant genetic resources for adaptation and expansion across its host range.

For M. oryzae populations, the temporal dynamics of genetic diversity are often closely related to rice cultivation. According to the US Department of Agriculture (https://www.ers.usda.gov/data-products/rice-yearbook/), worldwide rice cultivation areas have expanded rapidly, increasing by 11% from 1960 to 2019. The main rice-planting acreage and production regions are concentrated in Asia, where M. oryzae populations have maintained a high level of genetic diversity. Factors contributing to this diversity include genomic diversification, genetic instability and high variation rates, such as the presence of abundant variable number tandem repeats (VNTRs) and transposable elements (TEs) [62,63,64]. The rapidly decreasing genetic diversity observed in the 2020s could potentially be attributed to insufficient sampling from this period. In addition, the frequent host-shifting and host range expansion of this pathogen in crops and weeds have historically contributed to its sustained high genetic diversity [37,65,66], along with the gene flow events between cereal- and grass-specific lineages [28]. However, the continually growing genetic diversity of M. oryzae poses a significant threat to rice production, as these pathogens with diversified genetic resources could more rapidly overcome resistant genes in rice. While this genetic diversity could provide important new insights into the development and deployment of resistant rice varieties in practice. Under these scenarios, our findings underscore the need for greater attention to the effective monitoring and control of rice blast disease.

However, past hybridization events have imprinted on individual genome samples. In this study, the phylogenetic reticulate network and population differentiation indices indicated that incompatible events occurred within worldwide geographic populations, reflecting their harbored admixture status. In addition, three main clades were identified in worldwide populations from the genetic structure analysis, with each clade containing individuals from different continents, suggesting that geographical location could not sufficiently account for the population genetic differentiation of M. oryzae. Although geographical isolation plays an important role in shaping genetic divergence and leading to the geographic differentiation of many plant pathogen populations [67,68], the genetic differentiation of M. oryzae appears to be more vulnerable to its mating types [23]. While both MAT1-1 and MAT1-2 mating types have been identified in China [69], the individuals from China are distributed across three divergent clades. The population differentiation indices, as represented by FST values, demonstrate that the African populations are genetically distant from other populations, aligning with previous findings that African populations maintain relatively distant genetic relationships with Asian populations and are predominantly distributed in MAT1-2 [19,70].

In terms of the population structures of M. oryzae, populations with only one mating type occupied the principal position, indicating that clonal lineages dominated by asexual reproduction were prevalent in natural populations. The species can be geographically separated into a series of genetically less closely related subpopulations, according to population communication decreases [71] or gradual subpopulation increases following extinction and recolonization events. However, migration mediated by human activity and the transportation of infected seed [33] has made plant pathogen populations deviate from their evolutionary routes through recombination. In this study, an incomplete separation signature was also detected in the SplitsTree analysis for the reticular structures in the clades, which contained samples from Asia, Africa, Europe and North America. Additionally, five migration events among five continents were identified, consistent with findings reported by Tharreau et al. [33], indicating a global distribution of virulent genotypes. In this study, Asia was considered as a gene pool, which agreed with the result that Southeast Asia was recognized as a center of origin and diversity based on microsatellite markers [72]. This result is further corroborated by previous genetic diversity studies, which consistently revealed a higher genetic diversity within Asian populations [22,72]. Additionally, the Asian populations exhibited the highest nucleotide diversity in this study, underscoring Asia’s role as a center of genetic diversity and a source of migration. This finding aligns with the putative origins of rice in southern China and northeast India [60,61,73,74]. As a seed-borne pathogen, M. oryzae can infect two major subspecies: Oryzae sativa subsp. Japonica and Indica. The Japonica subspecies originated from rice domestication and then diverged into temperate and tropical Japonica subspecies, which progressively led to the formation of the Indica subspecies. This implies the possibility of pathogen migration following the domestication and introduction of rice crops [73,75,76]. Within Asia, migration events were observed in South Korea, Japan and China. Given that China is considered as the global origin of M. oryzae, it likely served as a foundational source for migration within Asia. Our results suggest that migration pathways originated from mainland China, initially extending to Japan and subsequently to Taiwan. These migration pathways have mainly been attributed to rice cultivars, as Japonica rice is preferred in both Japan and Taiwan. Furthermore, the limited trade and agricultural exchange between Taiwan and mainland China in recent decades could also have contributed to the limited gene flow. Similarly, a reciprocal gene flow was observed between populations in Japan and South Korea, potentially influenced by frequent interstate trade between these two countries. A similar pattern of results has been observed in Fusarium head blight pathogens, where agricultural practices and human migration have been critical driving forces in shaping the population genetic structure of Fusarium asiaticum [77].

In summary, gene flows at different geographical scales have exerted a significant influence on the establishment of the genetic structure of M. oryzae. Understanding the spatiotemporal dynamics of population compositions following the success of introduced fungal pathogens is essential for combatting plant disease pandemics [78,79].