{'fraction_non_alphanumeric': 0.03601266151211665, 'fraction_numerical': 0.021379274557625446, 'mean_word_length': 4.769046549917677, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An accurate and efficient method for large-scale SSR genotyping and applications

Abstract Accurate and efficient genotyping of simple sequence repeats (SSRs) constitutes the basis of SSRs as an effective genetic marker with various applications. However, the existing methods for SSR genotyping suffer from low sensitivity, low accuracy, low efficiency and high cost. In order to fully exploit the potential of SSRs as genetic marker, we developed a novel method for SSR genotyping, named as AmpSeq-SSR, which combines multiplexing polymerase chain reaction (PCR), targeted deep sequencing and comprehensive analysis. AmpSeq-SSR is able to genotype potentially more than a million SSRs at once using the current sequencing techniques. In the current study, we simultaneously genotyped 3105 SSRs in eight rice varieties, which were further validated experimentally. The results showed that the accuracies of AmpSeq-SSR were nearly 100 and 94% with a single base resolution for homozygous and heterozygous samples, respectively. To demonstrate the power of AmpSeq-SSR, we adopted it in two applications. The first was to construct discriminative fingerprints of the rice varieties using 3105 SSRs, which offer much greater discriminative power than the 48 SSRs commonly used for rice. The second was to map Xa21, a gene that confers persistent resistance to rice bacterial blight. We demonstrated that genome-scale fingerprints of an organism can be efficiently constructed and candidate genes, such as Xa21 in rice, can be accurately and efficiently mapped using an innovative strategy consisting of multiplexing PCR, targeted sequencing and computational analysis. While the work we present focused on rice, AmpSeq-SSR can be readily extended to animals and micro-organisms.


INTRODUCTION
Simple sequence repeats (SSRs), also known as short tandem repeats (STRs) or microsatellites, exist extensively in eukaryotic genomes. Most SSRs are non-coding and may affect gene expression, splicing, protein sequences and genome structures (1)(2)(3)(4)(5). The rate of SSR length mutations is estimated to be 10 −7 to 10 −3 per locus per generation in eukaryotes (6), which is much higher than the rate of ∼10 −9 for base mutations (7,8) and accounts for the high diversity of SSRs. Despite the high variability of SSR sequences, their flanking regions are often conserved within the same species, occasionally among closely related species (9,10) and even across species (11,12). SSRs have several advantages over other genetic variations, including co-dominance, high reproducibility and requiring a small amount of template DNA for experiment (13)(14)(15)(16). Importantly, the combination of the diversity of SSR sequences and the conservation of their flanking regions makes SSRs ideal genetic markers. Indeed, SSRs have been successfully adopted in various applications such as DNA fingerprinting, gene mapping, forensic analysis, marker assisted breeding and assessment of seed purity (16)(17)(18)(19)(20).
The most notable source of SSR diversity is the innate slippage of DNA polymerases during SSR replication (21-e88 Nucleic Acids Research, 2017, Vol. 45, No. 10 PAGE 2 OF 12 27). Such a slippage is also inherent in in vitro SSR polymerase chain reaction (PCR) amplification, which results in erroneous SSR alleles and makes accurate SSR genotyping difficult. Moreover, gel electrophoresis, the currently most popular approach to detecting SSR PCR products, is inaccurate due to its low resolution and is inefficient as it can only handle a small number of SSRs at a time (20). For example, only 28, 25 and 48 SSR loci are analyzed by gel electrophoresis and used to construct DNA fingerprints of jute (Corchorus spp.) (28), winter mushroom (Flammulina velutipes) (29) and pigeonpea (30), respectively. Such limited SSRs are not sufficient or robust to construct high-quality SSR fingerprints to discriminate close kinships. Whole genome resequencing can profile a large number of SSR loci at once (25,31,32). Nevertheless, since SSR sequences only constitute a small percentage of a whole genome, e.g. ∼3% of the human genome (33), whole genome resequencing dilutes the sequencing reads on SSRs with a huge number of other genomic reads, making it difficult to reach a required coverage of more than 10-to 100-folds (25) for accurate SSR genotyping with an acceptable cost. In addition, whole genome resequencing introduces additional problems to SSR genotyping, e.g. preferential amplification of specific SSR loci (34) and difficulties in data analysis with repeats in SSRs (35,36).
We developed a novel sequencing based SSR genotyping method by combining multiplex amplification of target SSRs, high-throughput sequencing of the amplicons (Ampli-Seq) and comprehensive computational and statistical analyses. For convenience, we call our new method AmpSeq-SSR. AmpSeq-SSR overcomes nearly all difficulties in the existing methods. In particular, AmpSeq-SSR is able to efficiently genotype a large number of SSR loci at once with ultra-deep coverages and has an accuracy close to 100% and a single-base resolution. We applied AmpSeq-SSR to genotyping a total of 3105 SSRs in eight rice varieties and validated all of them. We further demonstrated the utility and power of the new approach with two additional applications. The first is construction of rice fingerprints that contain 449 differential SSRs on average, which can accurately distinguish a variety under test. The second application is mapping of Xa21, a gene that confers a broad and persistent resistance against rice bacterial blight (BB), a devastating rice disease causing a substantial annual rice yield reduction worldwide. While the development of AmpSeq-SSR and our current work focus on rice, the new method is general and applicable to animal species and eukaryotic micro-organisms.

Rice plants used
The eight homozygous rice varieties (A-H in Supplementary Table S1) are commercially released varieties in China and are representatives of indica and japonica rice plants. The three pairs of nearly isogenic lines (NILs, I-N in Supplementary Table S1) share similar genetic backgrounds except the Xa21 gene and its linkage regions.

Target SSRs and design of multiplex primers
Forty-eight SSRs that are listed in the National Agricultural Standard of China (NY/T 1433-3014) and 3057 randomly selected SSRs from the Japanica reference genome (irgsp1.0) were chosen as target SSRs (Supplementary Table S2, the online file http://www.cse.wustl.edu/~zhang/ SSR ST2.pdf has the full list of primers). The service at https://ampliseq.com/ was used to design multiplex primers for 3105 target SSRs, which have amplicon lengths <250 bp on reference genome. The designed primers were synthesized by Thermo Company, USA.

Construction and high-throughput sequencing of Ampli-Seq libraries
We describe here the major steps for the construction and high-throughput sequencing of Ampli-Seq libraries. A stepby-step protocol is given in Supplementary Table S1 were also used to construct Illumina MiSeq Ampli-Seq libraries according to the user guide with modification. That is, additional PCR with 14 cycles were introduced for DNA amplification. The resulting libraries were sequenced on MiSeq system using the pairedend sequencing with a length of 2 × 300 bp. The library construction and sequencing for S5 and MiSeq systems were respectively performed in our lab and BestNovo Co., Ltd, China within a 30-day interval.

SSR genotyping by gel electrophoresis and Sanger sequencing
A 25 l PCR reaction system was used, which includes 0.05 M primer R and 0.05 M primer F, 10 ng template DNA and 12.5 l AmpliTaq Gold ® 360 Master Mix (4398876, Applied Biosystems™, USA). The PCR reaction procedure is 95 • C/5 min; 40 cycles of 95 • C/30 s and 60 • C/1 min; 72 • C/25 min. The amplicons were sent for Sanger sequencing (TsingKe Company, China) or were separated by 2% agarose gel electrophoresis, 6% polyacrylamide gel electrophoresis (PAGE) or capillary electrophoresis (CE). The CE was performed on ABI 3500 Genetic Analyzer (Applied Biosystems, USA) and the fragments were analyzed using GeneMapper softerware v4.1. The major and minor alleles of an SSR of CE were defined as the lengths of the largest and the second largest bands on electropherogram. The stutter ratio of an SSR was estimated as the ratio of the areas between the second allele and the major alleles of the SSR.
Processing of sequencing reads and SSR genotype calling. Here, we described the rationale and criteria for processing PAGE 3 OF 12 Nucleic Acids Research, 2017, Vol. 45, No. 10 e88 sequencing reads from homozygous varieties for SSR profiling. The related scripts were provided in Supplementary materials. The methods for heterozygous varieties were provided in the last section of 'Materials and Methods' section.
Reads were mapped to the reference genome by Bowtie2 (37) to determine which target SSRs they belonged to. However, many reads were left unmapped because of the variations in SSR lengths. To resolve this problem, sub-sequences of 40 bp were taken from the unmapped reads by a scheme of a sliding window of 40 bp. The sub-sequences were then mapped to the reference genome by BLASTN to determine potential SSRs they belonged to. To avoid the interference of SSR length variations on mapping, the SSR sequences on the reference genome were replaced with 'N'. After a read was located, the two boundaries of the SSR in the read were determined by aligning the read to the flanking sequences of the corresponding SSR on the reference genome (ref-SSR). Then, the tandem repeats of the ref-SSR motif between the two boundaries were determined as the SSR allele of the read.
The reads were tallied for each kind of allele at an SSR locus. The alleles with the largest and the second largest numbers of reads were designated as the major and minor alleles of the SSR locus, respectively. The major allele was also taken as the genotype because only a homozygous allele was expected for an SSR locus in the homozygous plants used in the current study. The ratio between the numbers of reads of the minor and major alleles was taken as the stutter ratio of the SSR locus.
To determine the amplicon length of a target SSR, the 20 bp upstream and downstream sequences of the amplicon on the reference genome were extracted and mapped to each sequencing read of the target SSR. The distance between the upstream and downstream sequences of each sequencing read was calculated and classified. The average distance in the class with the most reads plus 40 bp and the primer lengths was taken as the amplicon length of the target SSR, which corresponds to the amplicon length on electropherogram.

Estimation of AmpSeq-SSR accuracy and reproducibility
A valid SSR has a coverage of no less than a specified fold c and a stutter ratio no greater than a specified value s. For convenience, the specific genotyping condition of the valid SSRs is denoted as (c, s). A valid SSR is consistent if it has an identical genotype from the MiSeq and S5 systems, otherwise it is inconsistent. An inconsistent SSR is deemed to be incorrectly genotyped by the MiSeq or S5 system. The probability that a consistent SSR is the result of the same erroneous genotype from the two systems is negligible and therefore is ignored. Then, the genotyping errors of AmpSeq-SSR can be represented by half of the inconsistent SSRs. Therefore, the accuracy A(c,s) and the reproducibility R(c,s) of AmpSeq-SSR under a specific genotyping condition (c,s) are where M i and m i are the numbers of consistent and inconsistent SSRs in the ith variety, respectively. For a particular sequencing platform, the exact accuracy A(c,s) and reproducibility R(c,s) may slightly differ from the results from Equations (1) and (2) because different platforms have different error rates and sequence biases.

Improvement to the accuracy of AmpSeq-SSR by removing non-random error-prone SSRs
When an SSR is valid in two varieties by the two sequencing systems, we call the SSR comparable for the two varieties. If a comparable SSR is inconsistent for both of the two varieties, it is considered to be an error-prone SSR. Let n observed , n random and n non−random be the numbers of observed, random and non-random error-prone SSRs, respectively. Then, n observed = n random + n non−random .
We introduce a null hypothesis H 0 that n non−random = 0, i.e. all the observed error-prone SSRs occur randomly and the event that an SSR is inconsistent in one variety is independent of the event that the SSR is inconsistent in another variety. Let N ij be the number of comparable SSRs for the ith and jth varieties. Then, the total number of comparable in the current study). The probability of a random errorprone SSR can be estimated asq where p i and p j are the ratios between the inconsistent and comparable SSRs in the ith and jth varieties, respectively. The probability of having k random error-prone SSRs is binomially distributed, i.e. k ∼ Bi n(N,q). Then, the probability of having no less than n random error-prone SSRs in N comparable SSRs is Let 1 − α be the significance level for accepting H 0 (α = 1% in our study). When P(n = n observed ) > α, we accept H 0 . Otherwise, we reject it, meaning that non-random error-prone SSRs exist. When all non-random error-prone SSRs were identified and removed from the target SSRs, n observed = n random . Then, the inconsistent SSRs in the K varieties can be estimated as Kn random , where the value of n random is determined by the equations of P(n = n random ) > α and P(n = n random + 1) ≤ α. Following Equation (1) the accuracy of AmpSeq-SSR can be improved to: where M i is the number of consistent SSRs in the ith variety.
As an example, we used seven of the eight varieties to identify and remove the error-prone SSRs from the target SSRs, and calculated the improved accuracy of AmpSeq-SSR in where m ei is the number of the error-prone SSRs for the ith variety and one of the other seven varieties.

Estimation of the accuracy of Sanger sequencing
The MiSeq pair-end reads with identical forward and reverse sequences have negligible sequencing errors and have a quality comparable to Sanger sequencing. Therefore, the accuracy of SSR genotyping by Sanger sequencing can be estimated as the ratio between the reads for the correct (or major) allele of a valid SSR and all the reads of the SSRs.
where T i is the number of valid SSRs in the ith variety, c ik is the number of the reads which account for the major allele in the kth valid SSR of the ith variety and represents the times of correct SSR genotyping by Sanger sequencing, C ik is the total number of the reads of this SSR and represents all the times of SSR genotyping by Sanger sequencing.

Extension to heterozygous varieties and estimation of AmpSeq-SSR accuracy
To extend AmpSeq-SSR to heterozygous varieties and estimate its performance, we created a pseudo-heterozygous variety ij and generated a set of in silico sequencing data by randomly sampling and mixing 0.8 M reads from each of the ith and jth varieties. The genotypes from the ith and jth varieties were combined as the reference genotypes of the pseudo-heterozygous variety ij. To ensure the accuracy on the reference genotypes, we only considered those SSR loci with stutter ratios being lower than 0.5 and major alleles being covered by at least 50 reads in the ith and jth varieties. An error correction model was introduced to determine the homozygosis or heterozygosis of an SSR locus using the sequencing data from homozygous SSR loci following the principles outlined in (25). The model may include, for example, all the loci in rice varieties of the current study or the loci on the human X chromosome. In the current study, the sequencing data of variety G were used to build an error correction model to determine the heterozygosis of the SSR loci in pseudo-heterozygous varieties. The reads of all the homozygous SSR loci with the same genotypes were pooled to generate a set of in silico sequencing data of pseudo-SSR loci. For a pseudo-SSR locus of a specific genotype, the percentage of the reads for each observed allele among all the reads at this pseudo-SSR locus were used to estimate the probability that an allele would be detected at an SSR locus of this genotype. These probabilities were then used to calculate the probabilities of the six possible homozygous and heterozygous pseudo-genotypes from the three alleles with the most reads at an SSR locus of the pseudo-heterozygous variety ij. The homozygosis or heterozygosis of an SSR locus of the pseudo-heterozygous variety ij was determined by its pseudo-genotype with the highest probability.
The scripts in Supplementary materials were used to call for the alleles of an SSR locus in pseudo-heterozygous variety ij. The allele or the two alleles with the most reads were taken as the test genotypes for the homozygous or heterozygous SSR loci in pseudo-heterozygous variety ij, respectively. If a test genotype was the same as the reference genotype, the test genotype was taken as correct for the pseudoheterozygous variety. The accuracy of AmpSeq-SSR on the pseudo-heterozygous variety ij was then estimated as: where t i j and T i j were the numbers of the correct genotypes and reference genotypes in pseudo-heterozygous variety ij, respectively.

The AmpSeq-SSR method for SSR genotyping
We developed a novel method for large-scale SSR genotyping which is able to overcome the shortcomings of the existing methods ( Table 1). The new method, named as AmpSeq-SSR, has a higher resolution, can be more efficiently deployed, and is able to provide more accurate results that can also be compared from different experiments and/or methods. AmpSeq-SSR combines multiplexing PCR, targeted deep sequencing and comprehensive analysis ( Figure 1A). The current existing multiplex amplification techniques, such as that for whole exome sequencing, can amplify over 20 000 target loci in a single PCR reaction (38). Such a power of multiplex amplification was utilized in AmpSeq-SSR and 3105 SSRs were profiled in the current study. Since the genetic identities of SSR loci are mainly determined by their length variations, the SSR PCR products should be deeply sequenced on platforms (e.g. ABI S5 and Illumina MiSeq) that can produce long reads (e.g. 300 bp in the current study) to accommodate long SSR alleles. Sequencing reads were then mapped to target SSR loci to call SSR alleles. For homozygous samples, e.g. the rice varieties studied in the current study, only one allele was expected for a SSR locus and the allele with the most reads was taken as the SSR genotype (see 'Materials and Methods' section). For heterozygous samples, an additional step as described by (25) or (39) could be used to determine the SSR locus to be homologous or heterologous. Finally, the allele or the two alleles with the most reads were taken as the genotype of a homologous or heterologous SSR locus, respectively. Following the convention, an SSR genotype was represented as the amplicon length (see 'Materials and Methods' section), or the motif plus its repeat times, e.g. (AT)10.

Genotypes identified by AmpSeq-SSR are highly accurate
The SSR amplicons were designed to be <250 bp and could be sequenced through by MiSeq 2 × 300 bp paired-end sequencing mode from both forward and backward directions. On average, 1.38 million (M) (58.72%) of the 2.35 M paired-end reads for each of varieties A-H had identical forward and backward sequences and were expected to have negligible sequencing errors ( Table 2). As a result, 2427.75  Table 2). To assess the accuracy and reproducibility of AmpSeq-SSR, all of the 3105 target SSRs in varieties A-H were genotyped on the MiSeq and S5 sequencing systems by different labs (see 'Materials and Methods' section). An inconsistent SSR, which has distinct genotypes from the two sequencing systems (see 'Materials and Methods' section), was taken as a mistaken genotyping and used to calculate the accuracy of AmpSeq-SSR according to Equation (1). Even under a relaxed error-controlling criterion of the coverage no <10-folds and the stutter ratio no >0.5, AmpSeq-SSR had an accuracy of A (10, 0.5) = 99.73% when 10 584 of 10 641 SSRs had identical genotypes in the two sequencing systems (Table 3). With a lower stutter ratio of no >0.1 and a higher coverage of no <50-folds, AmpSeq-SSR had 100.00% accu-racy when all of the 4940 SSRs had identical genotypes in the two sequencing systems (Table 3).
A total of 21 pseudo-heterozygous varieties were generated (see 'Materials and Methods' section) from every pair of the 7 varieties A-H in Supplementary Table S1. By Equation (7), AmpSeq-SSR had an average accuracy of A ij = 94.47 ± 1.69% (Supplementary Table S3) on pseudoheterozygous varieties, which is comparable to the accuracy of SNP genotyping arrays on heterozygous varieties (40). The lower accuracy of AmpSeq-SSR on heterozygous varieties can be attributed to DNA polymerase slippage, especially when the real and slippage alleles had the same or similar lengths (25,39). It is noteworthy that AmpSeq-SSR should be more accurate on heterozygous varieties than other methods, such as whole genome sequencing (WSG) that is difficult to achieve a high sequencing coverage in order to statistically exclude errors from DNA polymerase slippage.  (25) and assumed that the average sequencing depth is 60-folds; 4 The greater number of the PCR cycles, the greater chance to incur DNA polymerase slippage; 5 The more PCR reaction times, the heavier workload and more time needed to perform PCR amplification; 6 Molecular ladder or reference sample; 7 For sequencing, the more highly the SSR loci are covered, the higher accuracy for SSR genotyping; 8 The more manual operation and the more detection times, the heavier workload and more time needed to detect and genotype SSR loci. We randomly selected three SSR loci in two varieties for validation of AmpSeq-SSR using Sanger sequencing, the gold standard for genotyping. All of the six genotypes of the chosen SSR loci were proven to be correct (see Figure 1B for an example). We also compared AmpSeq-SSR with CE, the currently most popular and accurate electrophoresis technique, on two randomly chosen SSRs ( Figure 1C and D). The lengths of the two SSRs detected by CE were respectively 5.53 and 4.47 bp shorter than those by AmpSeq-SSR ( Figure 1E), suggesting that the SSR genotypes from different batches of CE might be inconsistent and incomparable. The stutter ratios for CE were respectively 28.83-and 2.11-folds greater than that for AmpSeq-SSR ( Figure 1E), suggesting that the interference from CE on SSR genotyping was more serious than that from AmpSeq-SSR.

Improvement to the accuracy of AmpSeq-SSR
The accuracy and reproducibility of AmpSeq-SSR can be improved on SSRs with adequate sequencing coverages and lower stutter ratios (see 'Materials and Methods' section, Table 3). The stutter ratios of SSRs in varieties B-H were significantly decreased by 50.88% on average after removing the target SSRs with stutter ratios >0.2 in variety A (Sup- plementary Table S4). The SSRs with moderate GC contents of 40-70% had relatively uniform coverages of 0.2-to 5.0-folds ( Figure 1F). Lower coverages on regions of low GC content have also been observed in WSG of sugar beet (41) and transcriptome sequencing (42). Therefore, the phenomenon of low coverages on low GC content regions may be due to PCR during library construction and sequencing since it is not limited to special amplicon types, sequencing platforms or specific species. We called an SSR error-prone when it was inconsistent in two varieties. An error-prone SSR can occur randomly or non-randomly. Random error-prone SSRs are evidently unavoidable. A non-random error-prone SSR in one variety tends to be error-prone in another variety. Therefore, nonrandom error-prone SSRs, if exist, can be identified from the genotyping data of the tested varieties using a statistical analysis (see 'Materials and Methods' section). Following Equation (3), the probability of no non-random errorprone SSRs for the eight rice varieties in our current study is P (n = n observed ) = 1.29 × 10 −14 ≤ α = 1% (Supplementary Table S5). Therefore, non-random error-prone SSRs exist and the accuracy of AmpSeq-SSR can be improved by removing them from the target SSRs (see 'Materials and Methods' section). Following Equation (4), the accuracy of AmpSeq-SSR can be improved to A(10, 0.5) no−random = 100.00% when all non-random error-prone SSRs are removed, which can be realized when a sufficient number of varieties were analyzed. As an example, we used seven of the eight varieties to identify as many as 70.18% error-prone SSRs (Supplementary Table S6). Following Equation (5), the accuracy of AmpSeq-SSR in the eighth variety was actually improved to A(10, 0.5) improved = 99.92%. (For the calculation parameters, referred to Supplementary Table S6).
Most of the target SSRs were randomly selected from the reference genome (see 'Materials and Methods' section) and their distributions did not skew toward short SSRs (Figure 1G) or long motifs ( Figure 1H), which can be more accurately genotyped as suggested in (25) and our analysis (data not shown). Therefore, the high accuracy of AmpSeq-SSR did not come from the over-representativeness of short SSRs or long motifs.

Highly discriminative fingerprints derived via AmpSeq-SSR
The discriminative power of fingerprints on rice varieties is proportional to the number of SSRs included in the finger-prints. In this study, we employed as many as 3105 SSRs for fingerprinting rice varieties A-H, which are much more than the 48 SSRs widely used for rice identification currently adopted in the National Agricultural Standard of China (Standard No. NY/T 1433-3014) (Figure 2A). The power of the fingerprints is also positively correlated with the discernable differences between the SSR genotypes of two rice varieties. The resolution of SSR genotypes derived from AmpSeq-SSR is single base and therefore any subtle difference between SSR genotypes can be clearly observed. As a result, 449.71 SSRs on average were identified by AmpSeq-SSR to possess differential genotypes between two varieties ( Figure 2B and Supplementary Table S7). To our knowledge, 449.71 differential SSRs between two fingerprints are the largest ever reported, suggesting the great power of AmpSeq-SSR to unambiguously distinguish any rice variety under test.
Although electrophoresis can discern the differences of SSR amplicon lengths, it cannot distinguish base changes or base differences. Among the 449.71 differential SSRs, 33.68 (7.07%) had different bases but the same amplicon lengths between two varieties so that they were deemed to be mistaken as identical on electropherograms. For example, amplicon AMPL1141969 has the same amplicon length but distinct SSR genotypes in varieties F and G, which could be clearly identified by AmpSeq-SSR ( Figure 2C) but not by electrophoresis ( Figure 2D). Furthermore, when the differences of amplicon lengths were below the resolution, they might also be indiscernible on electropherograms (e.g. Figure 2E). Among the 449.71 differential SSRs, 221.32 (51.72%) had amplicon length differences no >5 bp (Figure 2B and Supplementary Table S7), which was the resolution of CE ( Figure 1E), the most accurate electrophoresis technique.

Mapping of Xa21 gene by AmpSeq-SSR
Rice BB is a devastating disease that causes a significant annual rice yield reduction worldwide and Xa21 is one of the most effective genes to control BB (43,44). Through back cross of over six generations, we introduced Xa21 from the donor variety IRBB21 into three receptor parental varieties, IRBB24, 9311 and D62B, and developed their respective NILs, R24, R11 and R62 (Supplementary Table S1). Each parental variety and its NIL share similar genomic background except for the target gene Xa21 and its linked region (hereinafter referred to as 'target region'). Used as genetic markers, the SSRs applicable to Xa21 mapping should be in the target regions of the NILs, have distinct genotypes between a receptor parent and its NIL but the identical genotype among the three NILs, R24, R11 and R62.
The 3105 target SSRs that were used for fingerprinting were adopted to map Xa21. Among them, 9, 14 and 55 SSRs have different genotypes between IRBB24 and R24, 9311 and R11, and D62B and R62, respectively ( Figure 3A and Supplementary Table S8). One of these SSRs, in amplification AMPL1562757 (hereinafter referred to as 'SSR21'), has the same genotype of (CGC)7 in the three NILs (Figure 3A and B). Therefore, using AmpSeq-SSR Xa21 can be easily mapped to the region near SSR21, which is about 0.6 Mbp from the actual locus of Xa21 (45). This gene mapping process was done within 24 h ( Figure 1A), which is significantly faster than any conventional method. However, PAGE electrophoresis failed to detect the key differ-ences of SSR21 genotypes between the receptor parents and their NILs ( Figure 3C). A close inspection of the Ampli-Seq reads revealed that the length differences of the SSR21 genotypes were compensated by non-SSR sequences in amplicons ( Figure 3D).

DISCUSSION
As a classic molecular marker, SSRs have been employed in broad applications, such as fingerprinting, gene mapping, forensic analysis and variety identification (16)(17)(18)(19)(20). Compared to the existing methods, AmpSeq-SSR is high applicable to those applications thanks to its high accuracy, high efficiency and low cost.

AmpSeq-SSR is more accurate than the existing methods for SSR genotyping
Accuracy is always the primary consideration for any genotyping method. The current study demonstrated that AmpSeq-SSR had a much higher accuracy than the existing electrophoresis-based methods for SSR genotyping ( Figures 1E, 2B-E, 3C and D; Supplementary Table S7). Except gel electrophoresis, Sanger sequencing is occasionally used for SSR genotyping. By double sequencing of the same sequences, AmpSeq-SSR could detect and avoid most of the sequencing errors to a level comparable with Sanger sequencing, i.e. ∼10 −6 per base. As a by-product, the accuracy of SSR genotyping by Sanger sequencing could be estimated with AmpSeq-SSR reads by Equation (6) as A(10, 0.5) Sanger = 83.89%, which is much lower than the AmpSeq-SSR accuracy of A(10, 0.5) = 99.73%. The whole genome resequencing can also be used for SSR genotyping (25,31,32). However, for a high accuracy that AmpSeq-SSR has, the target SSRs need to be covered for ∼500-folds (Table 2), which requires coverage of the whole genome for thousands of folds, resulting in an unacceptable cost.
By overcoming most of the shortcomings in the existing methods for SSR genotyping, the accuracy of AmpSeq-SSR reached A(10, 0.5) = 99.73% (Table 3), which can be further improved to A(10, 0.5) improved = 99.92% (for the parameters used, refer to Supplementary Table S6) or even potentially A(10, 0.5) no−random = 100.00%. The high accuracy makes AmpSeq-SSR capable to accurately identify a few or even one distinct gene from the plant varieties developed by transgene, backcrossing and mutation. However, a few distinct genes can be elusive within a background of hundreds of uncertain SSRs on electrophoresis profiles (Figure 2B and Supplementary Table S7). Moreover, AmpSeq-SSR has potential advantages for legal applications such as forensic analysis or granting rights under plant variety protection (PVP), for which any mistake may have serious consequences.

Factors that made AmpSeq-SSR accurate
The slippage of DNA polymerases during PCR amplification is the major source of error for SSR genotyping (21)(22)(23)(24)(25)(26)(27). The more PCR cycles, the higher chance for slippage errors to accumulate. Because visual inspection requires sufficiently bright bands, gel electrophoresis for SSR genotyping needs more PCR cycles than AmpSeq-SSR, i.e. 40 versus 16 PCR cycles in this study, resulting in more serious slippage (i.e. higher stutter ratios) for gel electrophoresis ( Figure 1E). Note that a sequencing library can be constructed using a PCR-free protocol (25), so that polymerase slippage during library construction can be avoided for AmpSeq-SSR.
Duplex sequencing can improve the accuracy of SSR genotyping by reducing sequencing errors. In addition to inherent and unavoidable sequencing error, the PCR cycle number (n) and sequencing coverage (N) are primary factors contributing to the genotyping error probability (P). An SSR locus cannot be genotyped correctly when slippage reads take up more than 50% of the total coverage of the locus. Therefore, the error probability can be estimated as: where r is the slippage probability of a read during a single PCR cycle, and R is the percentage of slippage reads in the final PCR product. The PCR cycle number (n) smaller than the electrophoresis based methods and the sequencing coverage (N) greater than WSG lead to a potentially higher accuracy of AmpSeq-SSR for SSR genotyping. The SSR often slips for a motif (2-4 bases) (25), making it difficult to distinguish the allele from slippage on electropherogram. Without slippage, accurately comparing SSR genotypes on electropherogram might still be nontrivial. For example, the SSRs in amplicon AMPL1141969 of variety F had negligible stutter ratio of no >0.01 (Supplementary Table S2) but migrated differently between two neighboring wells of the same gel ( Figure  2D).

High efficiency of AmpSeq-SSR extends its applications
We simultaneously amplified 3105 SSRs in eight rice varieties in a single PCR reaction and sequenced them in a single sequencing run. At present, such a combination of multiplex amplification and high-throughput sequencing can examine even more SSRs. Consider the MiSeq sequencing platform, for example, which can produce M = 25 million reads in a single run. When covered by m = 10 reads, an SSR can be genotyped with almost 100% accuracy (Table  3). Therefore, AmpSeq-SSR using the MiSeq platform can accurately genotype M/m = 2.5 million SSRs at once. In contrast, gel electrophoresis based SSR genotyping individually amplifies and separately examines each SSR, and thus is labor intensive and inefficient. Furthermore, thousands of SSR PCR reactions require a significant amount of template DNAs, e.g. an excessive ∼31 g for 3105 SSRs.
The ability to genotype a large number of SSRs is the key to the success of AmpSeq-SSR in many applications. The 3105 target SSRs used by AmpSeq-SSR for fingerprinting are critical for the great discriminative power of the resulting fingerprints. Another good example is the mapping of Xa21 gene. Because of the successive backcrossing and selection, the target region of Xa21 gene became rather small, i.e. ∼2 Mbp, as estimated in our previous study (46). This small region reduces the chance to have usable SSRs for gene mapping. On average, 449.71 (29.61%) SSRs of the 1518.61 comparable SSRs on the ∼400 Mbp rice genome had differential genotypes between two fingerprints (Supplementary Table S8). Therefore, 7.59 (2 × 1518.61/400) SSRs are expected to exist in the target region. The probability of having at least one differential SSR genotype between the three pairs of NILs within the target region is [1 − (1 − 29.61%) 7.59 ] 3 = 80.54%, suggesting a great chance for AmpSeq-SSR to successfully map a target gene.

High reproducibility of AmpSeq-SSR extends its applications
Reproducibility is important for a new technique. Based on Equation (2), AmpSeq-SSR had a reproducibility close to 100% (Supplementary Table S9), ensuring it to produce highly consistent and comparable results under various conditions. Therefore, AmpSeq-SSR can be used to collaboratively construct fingerprint database without sharing the original biological resources (e.g. specific rice varieties), which are often taken as trade secrets or national strategic resources. More importantly, a fingerprint from AmpSeq-SSR can be freely and accurately compared with all records in a database to determine its distinctness from all the existing varieties, which is the legal precondition for a PVP grant. However, the fingerprints from electrophoresis can be accurately compared only when they are constructed in parallel, making it nearly infeasible to compare one variety against thousands of existing varieties.

AmpSeq-SSR is economical
The cost of AmpSeq-SSR is low, i.e. ∼$0.015 per SSR in this study, which is more affordable than other techniques (47). The high consensus of SSR genotypes between MiSeq and S5 (Table 3) indicated that the cost could be further decreased by the more economical single-end sequencing strategy.

The advantages of AmpSeq-SSR over the existing methods
AmpSeq-SSR has comparative advantages over the existing methods for large-scale SSR genotyping, including accuracy, resolution, throughput, efficacy and comparability of results from different experiments and/or methods (Table 1). Besides, comparing with WGS, AmpSeq-SSR is able to avoid the interference of homologous sequences on target SSRs by designing multiplex primers that do not anneal to homologous SSRs, which would be difficult if not impossible to realize for WGS. To overcome this shortcoming, WGS reads have to be discarded when the regions flanking SSRs are not perfectly matched to the reference, e.g. that appeared in (25), resulting in a genotyping failure for SSRs with variations in their flanking regions. In order to correct genotyping errors, a WGS-based strategy needs to resort to an ultra-high-throughput sequencing to deeply sequence target SSRs (Table 3). On the other hand, ultra-high-throughput sequencing typically produces short reads, producing a skewed distribution of reads toward short SSRs. AmpSeq-SSR has evident advantages and shows no obvious shortcomings over the other methods ( Table 1), suggesting that AmpSeq-SSR is the most preferable method for large-scale SSR genotyping.

ACCESSION NUMBERS
All sequencing data of the eight rice varieties produced and analyzed in the current study have been deposited into NCBI Sequence Read Archive (SRA) under accession number of SRP089769.

AVAILABILITY
The package of scripts that implement the AmpSeq-SSR method is freely available at a public repository https://github.com/SystemsBiologyOfJianghanUniversity/ AmpSeq-SSR.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.