FASTX-Toolkit (version 0.0.14)² was used to pre-process the raw reads. First, base pairs with low-quality scores (<20) were trimmed, and the resultant reads <50 bp were removed. Then, the processed reads were mapped to a synthetic reference (see Design of MNP marker and multiplex PCR primer subsection) using Bowtie2 (version 2.1.0).

For each read that covered the entire genomic region of an MNP locus, the sequence located within the locus was used as its allele. To reduce errors caused by incorrect read alignment, any variation within consecutive mismatches with reference (mismatch interval ≤2 bp), repeat sequences (SSR), indels, and their flanking 2 bp of an allele were ignored. In addition, loci with at least one typed read were considered detected loci, and those with >100 typed reads were deemed valid MNP loci.

Valid MNP loci were genotyped (Supplementary Figure S2). Briefly, all alleles identified in a valid MNP locus were referred to as candidate alleles of the locus in a microorganism strain. The candidate allele with the most supported reads was the major allele of the MNP locus in the strain if it was supported by at least 50% of the mapped reads on the locus; otherwise, this locus was excluded from the following analysis.

The remaining candidate alleles were a mixture of genuine minor alleles at lower frequencies and artifacts caused by PCR or sequencing errors. To distinguish genuine alleles from false ones, candidate alleles were first subjected to strand bias filters: (1) if the strand bias was greater than 10-fold or the differences between the strand bias of the candidate allele or that of the major allele were more than five-fold, the candidate allele was removed; (2) Fisher’s exact tests were applied to determine the strand bias, and candidate alleles with multiple tests corrected p-values <0.05 were abandoned.

Subsequently, a statistical model was applied. A candidate allele with c supporting reads was primarily assumed to be the product of the major allele solely due to amplicon sequencing errors. Under this hypothesis, the number of supporting reads of a candidate allele follows a binomial distribution with parameters d (marker sequencing depth) and e (amplicon sequencing error rate). If the candidate allele had an excessively larger number of supporting reads than expected, namely, if P(X≥c) was significantly low, the null hypothesis was overridden. When multiple candidate alleles were tested, the p-values were corrected for multiple testing, and candidate alleles with an FDR <0.5% were deemed true minor alleles.

Amplicon sequencing errors supposedly occur independently across reads; thus, the probability of multiple errors occurring is lower within a read sequence than that of a single error. Therefore, the parameter e used in the statistical model depends on the number of SNPs (n) between the candidate allele and major allele. The minor alleles detected at the homozygous MNP loci can be regarded as products of technical error. In this study, the parameter e for each n, denoted as e (n), had the highest ratio of reads that an erroneous allele bearing n SNPs could take up at a homogeneous MNP locus. Microorganisms typically live in populations, and monoclonal organisms are often composed of individuals of multiple generations. Therefore, the homozygous MNP loci in microorganisms could not be validated. Note that rice is a diploid plant, and F₁ hybrids are produced by hybridization. The MNP loci identified as homozygous in both parents of hybrids are expected to be homozygous in the hybrid. Any minor alleles identified at these MNP loci can be attributed to amplicon sequencing errors. Therefore, we utilized the minor alleles detected at the expected homozygous MNP loci in the individual leaves of rice plants by amplicon sequencing to estimate the parameters [e (n1)] and MNPs [e (n1)]. To minimize sampling errors, only loci with depth >1,000 were considered. Based on the 930 MNP markers of the 16 individual rice plants, e (n = 1) and e (n > 1) were estimated to be 1.03 and 0.0994%, respectively.

We designed the MGV-Seq method to profile MGVs. The MGV-Seq procedure began with the selection of MNP markers (Fang et al., 2021) using the (re)sequencing data of representative strains from a species (see Materials and methods section). Briefly, the genome database of target species was constructed by collecting genome assemblies of conspecific strains from a public database, such as National Center for Biotechnology Information. Then, the genome assembly of a representative strain was used as the reference genome to call SNPs in the other genome assemblies, and a sliding window was used to scan those genomes assemblies to screen windows with at least two SNPs. The length of the sliding window and step between windows is determined based on the length of the amplicon that can be processed by the sequencer used. The windows with high discriminative power within species and uniform distribution in the reference genome were selected as candidates of MNP markers (see Materials and methods section). Primers for the multiplex amplification of the candidate MNP markers were designed using the available tools and synthesized after adding adaptors adapted to the sequencer used. Using the synthesized primers and a suitable kit for amplicon sequencing, a small number of representative strains was tested to screen MNP loci that could be detected in all strains to form the final MNP marker set and primer pools of MGV-Seq.

Next, MGV-Seq can be used to detect MNP markers in samples from pure culture to uncultured diseased host samples (Figures 1A,B). The PCR amplification products of each sample were ligated with a unique DNA barcode. Finally, the barcoded PCR products were mixed to form a library for high-throughput multiplex sequencing (Figure 1C). A customized computational tool was developed for MNP genotyping (Figure 1C). Briefly, we first assigned sequencing reads to samples based on DNA barcodes and mapped the reads to the reference genome. Next, all reads of an MNP locus were tallied to call for candidate alleles. The candidate allele supported by >50% of the total reads on a locus was defined as the major allele of the locus. All candidate minor alleles for an MNP locus were assumed to originate from the major allele due to amplicon sequencing errors. The true minor alleles of the locus were determined by strand bias filtering and statistical modeling (Figure 1C, see Materials and methods section). Finally, the MNP genotypes of each sample were constructed by integrating the alleles of all detected MNP loci.

A procedure for applying MGV-Seq to homogeneous and heterogeneous strain determination, strain identification, strain authentication, and monitoring strain microevolution was proposed (Figure 1D). The heterogeneity of the tested samples was determined based on the number of minor alleles detected in the MNP genotypes. By comparing the MNP genotypes of the tested sample with those of the conspecific, parental, or standard strains, MGV-Seq can be further used for identifying the dominant strain in the tested sample, monitoring the microevolution of daughter strain lines, and authenticating the strain to be used (see Materials and methods section).

This study used the phytopathogen Xanthomonas oryzae (X. oryzae) as a model microorganism. A set of 213 MNP markers was developed (Supplementary Table S1; see Materials and methods section) based on 458 genome assemblies of X. oryzae strains downloaded from the National Center for Biotechnology Information (NCBI), including 427 Xoo and 19 Xoc genomes (Supplementary Table S2). Widely used P1–P10 Xoo strains from IRRI (Quibod et al., 2016) and daughter lines from three different laboratories in Beijing, Wuhan, and Hainan, China, were collected (Supplementary Figure S1). The P1–P10 strains were short-handed for PXO61, PXO86, PXO79, PXO71, PXO112, PXO99, PXO145, PXO280, PXO339, and PXO341 in this study. Biological samples with the same name but from different laboratories are referred to as homonyms in this study, while DNA samples sequenced at the same run are referred to as a sequencing batch. Repeats of the same DNA tested in the same and different sequencing batches are called technical and batch repeats, respectively.

The dominant strain within each sample was identified by comparing the major alleles on the detected MNP loci with the reference genotypes of the conspecific strains. Therefore, the reproducibility of MGV-Seq for strain identification can be estimated using the reproducibility of MGV-Seq for detecting major alleles. A total of 4,971 MNP pairs among the 36 technical repeats and 4,997 MNP pairs between the 36 batch repeats were generated. The major alleles in all MNP pairs were 100% reproducible (Figure 2A and Supplementary Tables S3, S4). Reproducible major alleles were regarded as correct because there is only a tiny chance that major alleles of one MNP locus, detected by two repeats, were consistent but incorrect. Therefore, MGV-Seq demonstrated 100% reproducibility and accuracy for the technical and batch repeats for major allele detection. These 12 strains included three homonyms of the four strains. Pairwise comparison of homonymous strains was performed to examine the reproducibility of MGV-Seq for detecting differences across strain homonyms. A total 4 × (32) = 12 pairs of homonymous strains were formed and each pair performed 4×4 = 16 repeats from the four replicates of each strain across the two batches. The reproducibility of MGV-Seq for detecting differences across strain homonyms was 100% (Supplementary Table S5).

The sensitivity and specificity of MGV-Seq were estimated based on its performance in detecting mixed samples. A mixed sample was characterized by the presence of minor alleles belonging to MNP genotypes of the different strains. MGV-Seq was tested on 14 single clones of WH-P8 designated as CP8-1–CP8-14, eight single clones of WH-P6, designated as CP6-1–CP6-8, and a series of mixtures by spiking the DNA of one WH-P8 clone into one WH-P6 clone at increasing reference ratios of 0.1, 0.3, 0.5, 0.7, 1, 3, 5, and 7% w/w, at the first sequencing batch (Supplementary Table S3). The genotyping results of the single clones and mixtures were analyzed and compared using the MGV-Seq-adapted computational pipeline developed by us and LoFreq, a previously reported ultra-sensitive SNP-based variation detecting tool (Wilm et al., 2012). Minor alleles detected in the mixture belonging to either WH-P6 minor alleles or WH-P8 major alleles were called true positive minor alleles (TP), while others were called false positive minor alleles (FP). The sensitivity and specificity of MGV-Seq were referred to as the lowest ratio of spiking DNA in the detected mixture and the percentage of TP to the sum of TP and FP, respectively. MGV-Seq reported 20–44.4% TP for 0.1–7% reference ratios of mixtures and achieved a perfect (100%) specificity (TP/(TP+FP)) at all mixture ratios (Figures 2B,C). In comparison, the percentages of TP detected by LoFreq for the 0.1–7% reference ratios ranged from 12–51.7% (Figure 2B). The results showed that the sensitivity of MGV-Seq and LoFreq in the detection of mixed samples was as low as 0.1%; however, LoFreq had lower specificity for six of the eight reference ratios tested, for example, 92.3% specificity at 0.5% mixture ratio (Figure 2C).

For technical comparison, whole genomes of the 10 Xoo strains from IRRI were sequenced (Supplementary Table S6). A total of 1,400 MNP loci detected by MGV-Seq in the 10 Xoo strains (IRRI-P1-IRRI-P10 at the second sequencing batch) were covered by the WGS data, and the major alleles were 100% verified by the WGS data (Supplementary Table S4), suggesting the genotyping reliability of MGV-Seq.

Genetic homogeneity or heterogeneity of the strain was determined based on the minor alleles detected at the MNP loci of each strain. We applied MGV-Seq to detect 40 laboratory-adapted Xoo strains corresponding to 10 distinct parental strains each with four homonyms at the second sequencing batch. The four homonymous strains were P1–P10 from the IRRI and daughter lines preserved in the Beijing, Wuhan, and Hainan laboratories, designated as IRRI-P1–IRRI-P10, BJ-P1–BJ-P10, WH-P1–WH-P10, and HN-P1–HN-P10, respectively (Supplementary Figure S1). On average, 1,676,528 reads mapped to 203 MNP loci per strain, with an average coverage of 8,245 folds per marker (Supplementary Table S3). No minor alleles were detected at the MNP loci in four of the 40 Xoo strains, which were considered homogeneous strains. MGV-Seq identified 254 minor alleles in the remaining 36 Xoo strains (Figure 3A), and the number of minor alleles in each strain ranged from 1–37 (Figure 3B), indicating that 90% of the tested laboratory-adapted strains were heterogeneous. Moreover, 31.9% (81 of 254) of minor alleles had a frequency <0.9% (Figure 3D; Supplementary Table S7), which was considered to be the lower limit of the SNP genotyping error rates (Bayer et al., 2017, Fang et al., 2017, Li, 2016, Sasaki et al., 2018). Therefore, accurately capturing low-frequency minor alleles using SNP genotyping-based methods is difficult.

We then identified these 40 strains by conducting a genetic cluster analysis based on the MNP genotypes of the 40 strains, the previously published genomes of P1–P10, named Pub-P1–Pub-P10 (Quibod et al., 2016, Salzberg et al., 2008, Song et al., 2023). and the genome of BLS254, representative strain of Xoc (Feng et al., 2015, Hajri et al., 2012, Lee et al., 2005, Quibod et al., 2016, Salzberg et al., 2008). As illustrated by the cluster dendrogram, 39 of the 40 strains were Xoo pathovars, but one (HN-P5) was distant from all analyzed Xoo strains but closer with Xoc strains, which was identified as Xoc (Figure 4A). Alarmingly, BJ-P2 and HN-P2 were not clustered with their parental strains but were closest to the P6 strains, WH-P3 and HN-P3 were closest to the P9 strains, WH-P4 and WH-P10 were closest to some P2 strains, and HN-P5 were not clustered with any Xoo strains, respectively, indicating that there was a deviation between the identified and claimed identities of these strains. Pairwise comparison of P1–P10 homonymous strains from the four laboratories was performed to examine MGV across strain homonyms. Of the total 10 × (42) = 60 pairs, 17 (28.3%) had distinct major alleles with ratios of 14.0–76.2% (Figure 4B; Supplementary Table S8), suggesting that MGV was prevalent among these strain homonyms. An ultra-high ratio of major allele differences existed between HN-P5 and other P5 strains, confirming the MGV-Seq genotyping result that HN-P5 belonged to Xoc rather than Xoo.

The significant degree of genetic discrepancy in the reference strain across daughter lines raised concerns about the origin of the genetic differences. Thus, we analyzed the distribution of distinct major alleles among the daughter lines and different strains. Different Xoo strains diverged for much longer than daughter lines and thus should have accumulated more mutations. Distinct major alleles among the daughter lines (14.0–76.2%, Supplementary Table S8) were greater than those among different strains (0–23.5%, Supplementary Table S9), and the distribution of distinct major alleles in the daughter lines was discontinuous (Figure 4C). We further analyzed the distribution of minor alleles among the 40 strains. Eleven (4.3%) of the 254 minor alleles appeared only once and differed from the major alleles at any locus, suggesting that spontaneous mutations may have caused them. The remaining 243 (95.7%) minor alleles were the same as the major alleles of the other 1–39 Xoo strains tested together (Figure 3C; Supplementary Table S7). These results strongly suggested that human error, such as sample mislabeling or contamination, rather than spontaneous mutations, might be the primary factor underlying the genetic discrepancy in the examined Xoo daughter lines captured by MGV-Seq. Genotypes of a randomly selected MNP locus in the strain homonyms were successfully validated using Sanger sequencing (Figure 4D).

A large number of MGV of Xoo homonyms suggested that these strains might have distinct pathogenicity, one of the most important phenotypic features of Xoo. We investigated the pathogenicity of P2, P3, P4, P5, and P10 homonyms from Wuhan, Beijing, and Hainan on the susceptible rice variety IR24 and a near-isogenic line IRBB2 carrying the Xoo resistance gene Xa2 (Figure 5; Supplementary Table S10). Most of the 10 pairs of strain homonyms with distinct major alleles induced significantly different disease reactions in IR24 and IRBB2 rice varieties, showing the effect of MGV on pathogenicity reproducibility. One (WH-P3 vs. HN-P3) out of the five pairs of Xoo homonyms without distinct major alleles had a difference in lesion lengths (LL) of leaves, which may be due to the effect of minor alleles on pathogenicity reproducibility.

We used MGV-Seq to detect X. oryzae strains in rice leaves—with similar leaf blight symptoms—collected from the same paddy field in 2019 and 2020. MGV-Seq reported 185 and 144 valid MNP loci in the two rice leaf samples, respectively. And 7 of the 185 valid MNP loci and 47 of the 144 valid MNP loci had minor alleles (Figure 6A). The number of minor alleles of each MNP loci was one, with frequencies ranging from 0.5 to 19.8% (Figures 6B,C; Supplementary Table S11), suggesting the presence of low-frequency genetic variation in these samples. We identified the strains in the two rice leaf samples by comparing the major alleles of MNP loci detected with the 458 genomes in the reference database of X. oryzae. To illustrate the genetic distance between the rice sample genotypes and the reference strains, we plotted a clustering tree (Figure 6B). For the clarity and readability of the clustering tree, the 10 strains, with genomes closest to that of P1-P10, in the X. oryzae reference database based on the calculated average nucleotide identity, and the X-representative strains were selected as the reference strains for plotting the clustering tree. All the results presented in the heatmap (Supplementary Figure S3) and the clustering tree showed that the dominant strain in the 2019 sample was closest to the PXO211 strain of Xoo, whereas that in the 2020 sample was closest to the L8 strain of Xoc (Figure 6B; Supplementary Table S12).

The object of strain authentication is to determine whether the homonymous strains used at different periods (from the same laboratory or among different laboratories) are genetically consistent with their standard strain. The feasibility of MGV-Seq-based strain authentication was validated via the authentication of P6 and P8 homonymous lines. We took IRRI-P6 and IRRI-P8 as the standard strains, respectively. The homonymous lines of P6 included BJ-P6, HN-P6, WH-P6 and eight WH-P6 clones (CP6-1 – CP6-8) and that of P8 included BJ-P8, HN-P8, WH-P8 and 14 WH-P8 clones (CP8-1 – CP8-14). A pairwise comparison of the MNP major alleles of WH-P6, WH-P8 and their single-colony strains with IRRI-P6 and IRRI-P8, respectively, revealed that all the daughter lines showed 100% major allele similarity (MAS) with their standard strains IRRI-P6 and IRRI-P8, demonstrating that they were authenticated properly (Supplementary Table S13). Indeed, pathogenicity identification showed no significant differences among the three homogeneous WH-P6 clones (Figure 5; Supplementary Table S10).

Overall, these results demonstrated that MGV-Seq can identify microorganisms at the strain level. Additionally, it can accurately detect low-abundance variations in cultured strains and uncultured host-contaminated samples (directly). Thus, MGV-Seq was readily available for daughter-line microevolution monitoring and strain authentication.