Scanning the Effects of Ethyl Methanesulfonate on the Whole Genome of Lotus japonicus Using Second-Generation Sequencing Analysis

Genetic structure can be altered by chemical mutagenesis, which is a common method applied in molecular biology and genetics. Second-generation sequencing provides a platform to reveal base alterations occurring in the whole genome due to mutagenesis. A model legume, Lotus japonicus ecotype Miyakojima, was chemically mutated with alkylating ethyl methanesulfonate (EMS) for the scanning of DNA lesions throughout the genome. Using second-generation sequencing, two individually mutated third-generation progeny (M3, named AM and AS) were sequenced and analyzed to identify single nucleotide polymorphisms and reveal the effects of EMS on nucleotide sequences in these mutant genomes. Single-nucleotide polymorphisms were found in every 208 kb (AS) and 202 kb (AM) with a bias mutation of G/C-to-A/T changes at low percentage. Most mutations were intergenic. The mutation spectrum of the genomes was comparable in their individual chromosomes; however, each mutated genome has unique alterations, which are useful to identify causal mutations for their phenotypic changes. The data obtained demonstrate that whole genomic sequencing is applicable as a high-throughput tool to investigate genomic changes due to mutagenesis. The identification of these single-point mutations will facilitate the identification of phenotypically causative mutations in EMS-mutated germplasm.

exposure in different organisms, including Arabidopsis thaliana (Greene et al. 2003;Till et al. 2011), Oryza sativa (Till et al. 2011), L. japonicus (Perry et al. 2009), Caenorhabditis elegans (Flibotte et al. 2010;Thompson et al. 2013), Solanum lycopersicum (Minoia et al. 2010), and Saccharomyces cerevisiae (Shiwa et al. 2012) at different rates. EMS also tends to produce random point mutations and induces a low level of chromosomal breaks and lethal effects (Greene et al. 2003). These effects provide a competent survival rate and allow subsequent analyses to be performed for both forward and reverse genetics.
The whole-genome sequence of Lotus japonicus ecotype Miyakojima (MG-20) is available, covering a total length of 315,073,275 bp (67% of the 472-Mb genome). A total of 91.3% of gene space is located in the genome sequence (Sato et al. 2008). The genome sequences of the chloroplast (150,519 bp) and mitochondrion (380,861 bp) also have been assembled, by Kato et al. (2000) and Kazakoff et al. (2012), respectively. Numerous bioinformatics materials on legume and nonlegume plants are also publicly available from various resources (Sato and Tabata 2005;Goodstein et al. 2012;Li et al. 2012). With the current high-throughput tools, genome sequencing can be performed at an affordable cost . Many programs also have been developed for analyzing sequencing data in silico, offering good strategies to employ our study.
In this study, we applied Illumina second-generation sequencing (2GS) to discover EMS effects and the mutation spectrum in the genome of model legume, L. japonicus ecotype MG-20. We randomly selected two mutagenized plant genotypes from an M3 population as our subjects that were deep-sequenced. Wild-type (WT) MG-20 also was resequenced and mapped to the reference genome from Sato et al. (2008) as a comparison to subtract natural variations and false positives. We aimed to scan the effects of EMS throughout the whole genome, regardless of the phenotypic characteristic, which resulted from the mutations in specific regions. We identified single-nucleotide polymorphisms (SNPs) and compared the base alterations that occurred between both the genomes. The data demonstrate how 2GS works as a high-throughput platform for rapidly identifying DNA changes in each EMS-induced genome. As an advantage over sequencing pooled mutants, scanning individual mutagenized genomes allows rapid analysis of the mutation spectrum and gives the actual picture of "corrupted" genetic structure. The output of this study also will provide information in the identification of genes mutated due to EMS mutagenesis and demonstrate the distribution of mutation is comparable in different germplasm of MG-20.

Plant materials
A total of 4920 seeds of MG-20 were treated with 0.5% (v/v) EMS and grown as described in Biswas et al. (2009). One hundred MG-20 seeds were soaked in sterile water as a germination control. After 3 wk, the phenotypes of the plants that survived (67.3%) were observed to examine physical effects on the plant growth. Physiological tests also were performed as reported by Biswas et al. (2009) on abscisic acid (ABA) insensitivity and by Chan et al. (2013) on ethylene insensitivity. In this study, two mutated germplasms (M3) were isolated from the ABA assay (called AM and AS) for sequencing to screen EMS effects on their genome sequences and compare them with the resequenced WT MG-20 genome. In short, AM is a homozygous dominant, ABAinsensitive mutant (confirmed by stability in segregating families), whereas AS is the WT ABA phenotype segregant that was generated from a self-regeneration of a heterozygous ABA-insensitive mutant. Thus, AM and AS should carry the same spectrum of SNPs due to EMS mutagenesis.

DNA extraction
Genomic DNA was extracted from each individual plant after 1 mo of growth. The cetyl trimethylammonium bromide (CTAB) extraction method was adapted from Stewart and Via (1993) with minor modification. Plant tissues (about 1 g) were ground to powder in liquid nitrogen before adding 1 mL of CTAB extraction buffer. The mixture was incubated at 65°for 30260 min. Five-hundred microliters of the mixture was transferred into a new 2-mL tube and 500 mL of CTAB buffer was added to each tube. Both tubes were incubated as before. Chloroform purification was performed followed by isopropanol precipitation and washing with 70% ethanol. The nucleic acid containing pellet was air dried and dissolved in 100 mL of sterile water and subjected to RNase treatment (Sambrook and Russel 2001). Extracted genomic DNA was assessed using a spectrophotometer and agarose gel electrophoresis.

Sequencing and SNP identification
Whole-genome paired-end, 100-bp, short-sequence reads (.10· coverage) for AM, AS, and WT MG-20 were generated using the Illumina Genome Analyzer IIx according to the manufacturer's instructions. These three datasets were then mapped to the MG-20 reference genome (www.kazusa.or.jp/lotus/) using SOAP2 v2.21 with the option -r0 to retain only uniquely mapping read pairs (Li et al. 2009). SNPs were called using SGSautoSNP 2.001 (Lorenc et al. 2012) wherein AM, AS, and WT were referred as different cultivars. To avoid false-positive output, only homozygous SNPs were selected for further analysis. The three datasets of paired reads have been deposited into the National Center for Biotechnology Information Short Read Archive database under accession: SRX719550 (AM), SRX729747 (AS), and SRX131060 (WT). All custom scripts are available on request by E-mail at dave. edwards@uq.edu.au.

SNP analysis
The resequenced WT genome was used for comparison to identify variants or base changes in the AS and AM genome sequences. The frequency of transitions and transversions generated in both mutated germplasms was also calculated. SNPs were categorized using SnpEff 3.0j (Cingolani et al. 2012) according to their effect on L. japonicus MG-202annotated genes (Sato et al. 2008).
Selection of genes involved in ABA perception and signaling pathways A number of ABA candidate genes also were selected to test whether a SNP is located in their sequences. As a preliminary study, a total of 32 genes reported to be involved in ABA signaling were selected as candidate genes. These genes commonly are reported in the ABA-gene interaction in ABA perception and signaling pathways (Wasilewska et al. 2008;Cutler et al. 2010;Kim et al. 2010).

Phenotypic effects
All control seeds successfully germinated and were well-developed. EMS-treated seeds had a germination rate of 67.3% (3313/4920). Among well-developed EMS-treated plants, 76 M1 plants showed abnormal phenotypes after 3 wk of growth. Almost 1% of these plants (30/3313) showed an albino phenotype, in which yellow sectors were observed. Some plants also had pale green patches (0.6%, 23/3313), early branching (0.2%, 7/3313) or a looped base (0.39%, 13/3313). Two other plants showed either unusual leaf shape or early flowering. Vivipary also was detected on a pod among this population. In this study, the impaired phenotypes of L. japonicus indicated the successful treatment of mutagenesis using 0.5% EMS, and therefore, the population could be utilized for subsequent analyses. We chose mutants at the third generation (M3) as subjects to ensure the mutation is stable and fixed (Serrat et al. 2014).
Sequencing and read mapping output A total of 32, 965,291, 34,020,296, and 25,737,274 paired reads were generated from the WT, AS, and AM genomes, respectively (Table 1). At least 23% of reads from the sequenced genomes mapped to the reference. The percentage of mapped read pairs was 28.17%, 30.28%, and 23.59% from WT, AS, and AM, respectively, which resulted in more than 19· genome coverage. As a result, the genome coverage was sufficient to be applied for SNP calling between mutants and WT reads.

Frequency of mutation
After read mapping, homozygous SNPs were predicted, to identify SNPs that are unique in the mutant genome as opposed to WT. As a result, the frequency of mutation could be observed in assembled sequences of all chromosomes (Table 2). Chromosomes 1 and 6 have the longest and shortest length of assembled sequences, respectively. Meanwhile, the unmapped contigs cover approximately 32.9 Mb of total assembled length. Our SNP calling showed that Chromosome 1 has the greatest number of SNPs with one mutation per 170 kb and 165 kb in the AS and AM genomes, respectively. Mutation frequency was the lowest in Chromosome 6, which had less than 6% of total SNP in both genomes. The change rates were one mutation per 222 kb for AS and 217 kb for AM. In total, mutation frequency in AM and AS was nearly identical 1490 SNPs and 1447 respectively. This frequency is reflected in the change rate of mutation in the AM genome (one homozygous mutation in every 202 kb) and AS genome (one homozygous mutation in every 208 kb).
The frequency of transitions and transversions generated was analyzed (Table 3). In our mutants, EMS has generated 64.7% (AS) and 62.3% (AM) transitions as a percent of total mutation. Both mutants showed a bias to G/C-to-A/T transitions, which were the most frequent mutations that occurred (45.0% in AS and 34.9% in AM). Transversion mutation was also detected, but at lower percentages as listed in Table 3. The lowest percentage of mutation was C/G-to-G/C changes in AS and A/T-to-T/A changes in AM.

Distribution of transition and transversion mutations
Based on the mutation frequency, how each mutation type was located in individual chromosomes can be observed by calculating the percentage of mutation type in relation to the SNP total in each chromosome ( Figure 1). Both genomes comprised of C/T or G/A transitions as the most frequent mutation type in their chromosomes or unmapped regions. In the AS genome, G/A transitions were the highest in Chromosomes 1, 2, 5, and 6. Chromosomes 3 and 4 and unmapped regions had C/T transitions as the highest percentage of mutation type. The lowest percentage of mutation type was C/G transversions which were present the least in all individual chromosomes and unmapped regions of the AS mutant. The distribution of transitions and transversions in the AM genome was relatively similar to the AS genome. The greatest percentage of mutation type in each chromosome of the AM genome was C/T transitions except Chromosome 3. Unmapped regions and Chromosome 3 of AM have G/A transitions as the greatest percentage of mutation type. Meanwhile, the lowest percentage of mutation type was detected to be either A/T or C/G transversions in the AM genome.

Distribution of SNPs across mutagenized genomes
We determined how many SNPs occurred every 1000 kb to show the distribution of mutations across our mutagenized genomes (Figure 2). Regardless of mutation types, the distribution of SNPs is unique between AS and AM when comparing the same chromosome. SNPs were randomly located along the genomes with no specific chromosome position being particularly abundant or lacking in SNPs for both genomes. However, Chromosomes 1 and 2 were highly "corrupted" in their arms. Meanwhile, some chromosomes (Chromosomes 3, 4, and 5 for AS; Chromosome 3 for AM) have a high peak of SNPs toward their center. The density of mutations in every 1000 kb was relatively apparent in Chromosomes 1 and 2 of AS and AM. Meanwhile, Chromosome 6 had less dense mutations for both genomes.

Coding and noncoding effects
Here, we used SnpEff to predict the effect of mutations on coding regions of annotated genes of L. japonicus for both mutants, as shown in Figure 3. Our results showed that the greatest number of SNPs (34%) was predicted to be located in intergenic regions, followed by downstream and upstream regions of predicted genes (27% each). Only 7% and 3% of SNPs were located in exon and intron sequences, respectively. Intragenic regions also were predicted to have low EMS effect. The mutation percentage was very low at splice site donor and acceptors. In our mutagenized genomes, EMS effects did not contribute to a high fraction of nonsynonymous and synonymous changes (5% and 2%, respectively).

Candidate genes
As preliminary research, we selected 32 notorious ABA genes (Table  4) to determine whether any SNPs were located in their sequences. Orthologous genes of L. japonicus were identified based on A. thaliana and Glycine max sequences. Six were found to have orthologous sequences, which were positioned in unmapped regions of the MG-20 genome. The rest of the orthologs were found in the assembled chromosomes. SNPs were identified in three candidate loci, which were ABI2, ABI3, and ABI4. ABI2 and ABI4 were mutated in both mutants, whereas SNP was only found in ABI3 of the AS mutant.

DISCUSSION
The effect of EMS mutagenesis initially can be observed through phenotypic changes of mutagenized plants, and numerous reports n have been published on the effects. The abnormal phenotypic variance due to EMS exposure also has been reported previously in L. japonicus (Perry et al. 2003) and other plants such as Glycine max (soybean; Carroll et al. 1986 in L. japonicus mutagenized by EMS. Response to high temperature also was affected in EMS-induced Arabidopsis mutants impaired in ABA and salicylic acid syntheses (Zhu et al. 2012). Sequencing data often are used to mine SNPs or polymorphisms among different cultivars for identification of traits and allelic variations in crops Edwards et al. 2013;Zander et al. 2013). The high-throughput technology also is applied in discovering causative mutations by pooling backcrossed segregant populations to increase SNP frequency and reveal mutated regions (Ashelford et al. 2011;Mokry et al. 2011;Hartwig et al. 2012). Here, we individually sequenced two selected mutant genomes and used 2GS data to compare mutation distribution between these mutants.
To our knowledge, this is the first report on the application of SGSautoSNP for discovery of SNPs induced by EMS mutagenesis. This tool uses assembled reads to identify homozygous SNPs without using the reference genome, which is only used for mapping reads. SGSAutoSNP was used to identify only homozygous SNPs because the identification of heterozygote SNPs leads to a high number of false positives, which would lead to an overrepresentation of the mutation frequency in the genome (Lorenc et al. 2012). We also managed to rapidly identify the mutation spectrum occurring in the mutant genomes using the output data.
Our read mapping has been successfully performed before SNP identification in each genome. However, a low percent of reads was mapped to our reference. A high percentage of paired reads was reported to map to the reference in other species such as Caenorhabditis elegans (average 92%; Zuryn et al. 2010) and A. thaliana (average 73%; Austin et al. 2011). In our case, only an average of 27% paired reads mapped. Similar results were observed in other species using the same SNP prediction software [for example in canola (Dalton-Morgan et al. 2014) and wheat (Lai et al. 2015)]. The number of reads mapped relies on the parameters set up during the mapping procedure. We only selected paired reads that mapped and ignored reads that aligned as single reads and only used reads that mapped to a single unique location to generate more accurate read mapping (Li et al. 2009;Lorenc et al. 2012). The assembled length of pseudomolecules in our MG-20 reference covers 67% of estimated genome size (Sato et al. 2008), which also influences read mapping output. In addition, reads that aligned to multiple positions also were removed to increase SNP calling accuracy and avoid false positives (Lorenc et al. 2012;Shiwa et al. 2012). These factors have reduced the number of reads mapped in this study.
Chromosome 3 is the largest among the six chromosomes of the MG-20 genome. However, the assembled length of Chromosome 1 is the longest, and presumably more complete, based on the genome assembly (Sato et al. 2008), allowing for better read mapping and a greater rate of SNP identification. This reflects the greatest number of SNPs predicted in Chromosome 1. Chromosome 6 has the lowest number of identified SNPs, because it is the shortest chromosome in assembled length (Sato et al. 2008). In addition, Chromosome 6 might contain more repetitive sequences and therefore less unique reads mapped, leading to a lower rate of SNPs being identified. Furthermore, the different rates of base changes between chromosomes may demonstrate the different capacity levels of each chromosome in tolerating mutation impact. There is insufficient evidence to reach clear conclusions; however, the presence of selection pressure, such as clustering of housekeeping genes on chromosomes, could contribute to inability to tolerate mutations. Additionally, not all SNPs will be identified across the genome due to repetitive regions and the stringency of read mapping, reflecting SNP/mutation density across chromosomes.
A different number of SNPs between AS and AM was expected, because it is impossible to obtain an exact total of SNPs in different mutated genomes. The mutation rate of our mutants was 1/208 kb and 1/202 kb, which were greater than previously reported in L. japonicus (1/502 kb; Perry et al. 2009). Different rates of mutation also were reported in various plants mutagenized by EMS, as summarized by Sikora et al. (2011). Although only homozygous SNPs n n A/C-to-C/A 12.7 15.3 G/T-to-T/G 12.6 13.3 A/T-to-T/A 6.9 4.5 C/G-to-G/C 3.1 4.6 A high frequency of transition mutation was observed as expected were taken into account in this study, the frequency of mutation was quite high. If heterozygous mutations were called, the frequency will be greater. Here, the frequency of mutation demonstrates the effectiveness of EMS to produce high mutagenesis as reported by many researchers (Dube et al. 2011;Shiwa et al. 2012;Serrat et al. 2014). The mutation spectrums in our germplasms were not consistent with a previous report on TILLING work of mutagenized L. japonicus Gifu (Perry et al. 2009) wherein 97.6% of base changes were G/C-to-A/T mutation upon EMS mutagenesis. Similar to TILLed Arabidopsis, a high frequency of G/C-to-A/T transition (99%) also was identified (Greene et al. 2003). A lower rate of G/C-to-A/T base changes was observed in other plants, including barley (70%; Caldwell et al. 2004), rice (70%; Till et al. 2007), and tomato (60%; Minoia et al. 2010), showing a possibility of the presence of various base changes upon EMS exposure. However, transition mutations could not be denied as the most common base mutation in the EMS sphere as reported previously (Lawley and Martin  1975; Sikora et al. 2011;Shiwa et al. 2012). In the TLLING analysis, a number of interesting genes with known sequences were selected for the identification of base changes in pooled mutated genomes (McCallum et al. 2000). In this study, we adopted a high-throughput technology (2GS) to scan homozygous SNPs that are present throughout the individual mutated genomes compared with our resequenced WT MG-20. This provided a wider range of EMS effects in a genome. Our approach also neglected natural variation between the reference and the resequenced MG-20 genomes. Our reference was used merely for read mapping and has not been counted during SNP calling (Lorenc et al. 2012).
Although the percentage of transversion mutations was relatively low in both genomes, they should not be omitted in a mutation study because they potentially are causal mutations for phenotypes of interest. Taking into account only homozygous mutations, our data also show that the presence of specific base changes is comparatively distributed among individual chromosomes and not abundantly located in a specific chromosome for both alleles. Comparing these two genomes, we found that transition and transversion mutations were present nearly consistently in the same chromosomes. These results reveal that the frequency of EMS-induced transition and transversion mutations is comparable between different individual mutants that were derived from the same mutagenized population.
As mentioned previously, the total assembled length affects the SNP total identified in each chromosome and consequently, their distributions. The longer the assembled length, the denser the mutation distribution was found to be along the chromosome. The mutation distribution also was affected by the occurrence of assembled sequences at specific regions (www.kazusa.or.jp/lotus/). Chromosome 6 is the best example, in which low SNP peaks were detected and the mutations distributed broadly. Read mapping is very difficult in centromeric regions, and therefore the density in these regions has not been investigated in this study. Distribution of mutations was scattered due to the effect of EMS as a mutagen that causes a random mutation (Lawley and Martin 1975;Greene et al. 2003;Sikora et al. 2011;Tagu et al. 2014). A random distribution also was reported by Thompson et al. (2013) in the genome of single C. elegans from different strains, which were mutagenized by EMS. They propagated mutagenized worms through 10 generations to obtain stable mutants. However, this is time-consuming for legumes, which have a longer generation period. Nevertheless, our 2GS data of third-generation mutants also could provide the literal effect of EMS mutagenesis throughout the individual genomes. EMS induces base changes or nucleotide substitution, which consequently alter codon sequences, leading to either nonsynonymous or synonymous effects. In genetics studies, nonsynonymous change is a favorable mutation effect because it gives a clue on which gene may be associated with a specific phenotype (Ng and Henikoff 2006). Our results showed that EMS has arbitrary and broad effects of mutation on codon sequences. Mutated coding regions could represent nonsynonymous changes, which are useful to identify gene of interest for desirable phenotypes.
To extend the implication of our data, preliminary work has been employed to specifically detect the presence of mutations in our genes of interest. A total of 32 genes reported to be involved in ABA signaling were selected as candidate genes. These genes are commonly reported in the ABA2gene interaction in ABA perception and signaling pathways (Wasilewska et al. 2008;Cutler et al. 2010;Kim et al. 2010). We chose to use ABA candidate genes because our mutants Figure 3 Mutation effects on codon sequences by type and region in our mutant genomes. Single-nucleotide polymorphisms were observed highly located in intergenic regions, upstream, and downstream parts of the annotated genes. Only a small percentage of nonsynonymous changes was predicted.
were isolated from a mutagenized population effected in response to ABA. We did identify base changes in several candidate genes, indicating a low rate of EMS effects on these sequences. In addition, identified SNPs were not only in the ABA-insensitive AM mutant, indicating they were background mutations. On the other hand, comparing mutations between AM and AS mutants would provide clues on the causal mutation of mutant phenotypes. However, other candidate genes also can be considered because ABA is involved in a wide range of plant growth and responses to environmental stresses. Its effects are varied in regulating the events of plant physiology, which require gene interaction and/or cross talks with other hormones (as reviewed in Kermode 2005;Fujita et al. 2006;Rock et al. 2010). Because this research worked on an individual genome of mutants, another strategy has been developed to sequence pooled DNA from different mutants to remove the effect of background mutation and increase the likelihood of identifying the causative gene.
n In conclusion, this scanning revealed a detailed effect of EMS mutation on the whole genome of an individual mutant. Our results presented an overview of point mutations that occurred in the genome of mutants, which were usually pooled to identify SNPs. Here, EMS has produced a number of abnormal plants in our mutant population of L. japonicus. Our 2GS data also revealed how EMS efficiently mutates genomic sequences in an individual mutagenized plant. As expected, EMS created a random spectrum of mutation across the whole genome of our mutants and biased to G/C-to-A/T changes. However, other transition and transversion mutations also were identified with quite apparent fractions. Calling only homozygous SNPs has put a high confidence in identified base changes occurred. Mutation distribution apparently is distinct between mutated genomes, which derived from the same mutagenised population. Effects of SNPs on coding and noncoding regions could be manipulated to identify a causal mutation of a phenotype of interest. Our next question will be which gene is responsible for our abnormal phenotype of interest. This 2GS data will be further analyzed to discover the causative mutated gene.

ACKNOWLEDGMENTS
We thank the Centre for Integrative Legume Research, the University of Queensland for giving the opportunity to run this work. Special thanks to Dr Stephen Kazakoff, Queensland Centre for Medical Genomics, the University of Queensland for the technical support. A scholarship from the Ministry of Higher Education, Malaysia, in financially supporting our student also is acknowledged.