Genome Sequencing-based Coverage Analyses Facilitate High-resolution Detection of Causal Deletions in Gamma-Irradiated Wheat Mutants

genomic gamma-irradiated mutants of Japanese elite cultivar “Kitahonami” and performed whole-genome resequencing of these mutant lines. The comparative analysis of depth-of-coverage between the wild-type and the mutants identied ~130 Mbp and 67.8 Mbp deletions in the mutants “30579” and “28511,” respectively. These deletions were tightly linked to the mutant phenotypes in the progeny populations generated by crosses between wild-type and the mutants. The simulation analyses showed that 2.5x depth-of-coverage was enough to detect large deletions in the gamma-irradiated hexaploid wheat mutants. These results indicate that the approach of short-read-based genome sequencings of the gamma-irradiated mutants is cost-effective in designing genetic markers for agriculturally benecial traits to breeding of wheat.

grain shapes and anthesis are generated by X-rays and gamma irradiation [1]. The mutant plants provide novel and agriculturally useful traits, contributing to wheat breeding as genetic resources. Generation frequency of mutant phenotypes and mutation types are dependent on irradiation dose [1,3]. Particularly, gamma irradiation often causes structural variations of translocations and inversions on chromosomes and large deletions with chromosomal breakages [1]. Based on the gamma irradiation feature of random occurrence of large deletions, radiation hybrid mapping is developed. Through the assessment of presence or absence of chromosome segments using molecular markers, radiation hybrid mapping allows building the high-resolution maps in wheat and its wild relatives [4,5]. In addition, causal genetic regions of the phenotypes of the gamma-irradiated wheat mutants are mapped using simple sequence repeat markers in segregating populations derived from crosses between wild-type and the mutants [6].
Development in genome sequencing technology allows to genome-wide genotyping and decoding genome sequences. In maize, the combined bulked segregant RNA-sequencing and exome sequencing approaches are applied to gamma-irradiated maize mutants, successfully identifying causative deletions of the mutant phenotypes [7]. On the other hand, since common wheat (Triticum aestivum L.) has a hexaploid species possessing the large and complicated genome with formula AABBDD, it was challenging to apply genome sequencing approach to gamma-irradiated wheat mutants. However, highquality genome sequence of the common wheat cultivar "Chinese Spring" (CS) has been released [8] and the wheat pangenome project provides ten reference-quality genome assemblies and ve scaffold assemblies of common wheat [9]. The presence of the reference genome and the reduction of cost of genome sequencing allow exploring genome-wide polymorphisms even in common wheat, unravelling evolutionary history among cultivated wheat and their wild relatives in high resolution [10]. Now, the application of genome sequencing to common wheat becomes more easily accessible. Genome sequencing of agronomically bene cial mutants generated by gamma irradiation could more e ciently identify causative genomic regions of the mutant phenotypes compared with the classical genotyping approaches.
We developed gamma-irradiated mutants of common wheat cultivar "Kitahonami." "Kitahonami" is a Japanese winter wheat elite cultivar bred from a cross between the cultivars "Kitami 72" and "Hokushin." It is adapted to the northern area of Japan and provides high-quality soft grains for the Japanese white salted noodle, udon [11]. If a mutant showing useful agronomic traits is found in the mutant population, it can be used for improving "Kitahonami" as well as other wheat cultivars. We have identi ed a grain hardness mutant and a pre-harvest sprouting (PHS)-tolerant mutant that were expected to be available for further improvement of wheat cultivars.
Grain hardness of wheat is one of the essential traits to decide grain our characteristics. A Hardness (Ha) locus on the short arm of chromosome 5D regulates the grain hardness of common wheat [12]. Allelic differences in two puroindoline genes, Pina-D1 and Pinb-D1, in Ha locus cause variations of grain hardness ranging from soft to hard in common wheat [13,14]. Reduction of gene expression of the Pina-D1 and Pinb-D1 genes by the RNAi-mediated gene silencing increases grain hardness in transgenic common wheat [15,16], implying that the existence of puroindoline genes is required for forming soft grain in common wheat.
The moist condition provided by precipitation causes PHS in common wheat, a phenomenon of immature seed germination in pre-harvest spikes. PHS gives devastating damages on yields and quality of wheat grains [17,18,19,20]. The beginning of the rainy season in Japan coincides with the wheat harvest time. PHS often occurs when the rainy season begins earlier than as usual. PHS resistance is one of the most essential wheat breeding goals in Japan [21,22]. PHS is a quantitative trait controlled by multiple genes and is essentially linked to seed dormancy. Several genes associated with seed dormancy and PHS including MOTHER OF FT AND TEL1 (MFT1), TaMKK3-A, Viviparous-B1 (Vp-B1), Qsd1, and Tamyb10, have been identi ed [23,24,25,26,27,28].
In this study, using the grain hardness mutant and the PHS-tolerant mutant, we developed a method to identify causal regions responsible for mutant phenotypes. At rst, whole-genome resequencing was applied to these mutants. To detect deleted genomic regions induced by gamma irradiation, we conducted a sliding window analysis of depth-of-coverage over chromosomes. Molecular markers to con rm the causal deletion were designed and applied to the F 2 populations generated by crosses between wild-type and the mutants. Finally, we evaluated how much depth-of-coverage and window size of sliding window analysis were optimal to detect large deletions in gamma-irradiated common wheat.
These analyses will provide the good approach to effective breeding methodology using mutants with useful agronomic trait.

Plant materials
The Japanese elite cultivar "Kitahonami" is a soft winter wheat providing high-quality our for Japanese white salted noodles. To generate mutants, 250 Gy of gamma-rays were irradiated onto seeds of "Kitahonami." To screen grain hardness mutants, ~1,200 grains were randomly chosen from the M 2 mutants. The 26 M 2 grains showing amber-color grains were selected. Three individuals showing hard grain characteristics that were estimated by SKCS analysis were selected from the 26 individuals and were regarded as three lines. One line was selected from the three lines in M 4 . By repeating sel ng of the selected line, the M 7 mutant line, called "30579," was obtained. For the segregation analysis, F 2 progenies were generated from a cross between "Kitahonami" and "30579" (M 7 ).
To screen PHS-tolerant mutants, M 2 mutant seeds were randomly selected for M 3 mutants. In the eld of Kitami Agricultural Experiment Station, Kitami, Hokkaido, 5,000 M 4 mutants were grown. After owering, 1,500 spikes were harvested at the maturing stage, sprayed with water in the morning and evening, and incubated at 15 °C for seven days. Seven spikes with low germination and rooting were chosen as PHStolerant candidate lines. After M 5 generations, we performed PHS tolerance and seed dormancy tests on every generation, and selected candidate lines showing stronger PHS tolerance and seed dormancy than the wild-type. In M 7 , one candidate line showing stronger PHS tolerance and seed dormancy than the wild-type was selected and named "28511" mutant line. For the segregation analysis, F 2 progenies were generated from a cross between "Kitahonami" and "28511" (M 7 ).
Tests of grain hardness, PHS tolerance, seed dormancy, and ABA sensitivity Grain hardness was evaluated using SKCS 4100 (Perten, Stockholm, Sweden). The SKCS hardness index was obtained from crushing a sample of at least 50 kernels each from the wild-type and the grain hardness mutant "30579." To test PHS tolerance, spikes were detached during the maturity period after owering, incubated on trays at 15 ° C, sprayed with water in the morning and evening. After seven days, germination rate and rooting rate were calculated. For the seed dormancy test, 50 seeds 7 days after ripening were incubated in a petri dish at 10 °C or 15 °C. Germination rate was calculated based on the number of germinated seeds after three, ve, seven, or nine days.
To test ABA sensitivity, bioassay for ABA responsiveness were performed according to the protocol of Iehisa et al. [50]. Ten seeds of the wild-type or the PHS-tolerant mutant "28511" were treated with running water for six hours. The seeds were incubated on a wet lter paper in a petri dish at 4 °C under dark condition overnight and then incubated at 23 °C for 20 hours. Two petri dishes with lter paper were prepared. Six mL of 20 µM ABA was added to one of the petri dishes, and six mL of ultrapure water was added to the other one as a control. Five germinated seeds were put in each petri dish after measuring root length of the germinated seeds. Root length of each seed was measured 48 hours later. The elongated root length was de ned as the difference in root length after 0 hour and 48 hours. The suppression rate of root elongation was calculated by (elongated root length in the ultrapure water -(elongated root length in the ultrapure water -elongated root length in ABA))/elongated root length in the ultrapure water. Three independent bioassays were performed for the wild-type and the PHS-tolerant mutant "28511." Extraction of total DNA and genome sequencing Total DNA was extracted using the CTAB method from the leaves of wild-type and mutant plants of "Kitahonami." After RNase treatment, PCR-free DNA libraries for 150-bp paired-end DNA sequencing were constructed using TruSeq DNA library Preparation Kit (Illumina, San Diego, CA, USA). DNA sequencings of the wild-type was conducted on Ilumina HiSeq-X sequencer. DNA sequencing of the mutants was performed on NovaSeq sequencer. The sequencing reads were deposited in the database of the Komugi GSP (Genome Sequencing Program) (https://komugigsp.dna.affrc.go.jp/research/download/index.html).
Quality control, Alignments, and SNP calling Quality of the sequencing reads was checked using FASTQC ver. 0.11.7 [51]. Trimmomatic version 0.33 [52] was used to exclude short reads with an average minimum Phred quality score per four bp of less than 20 and a length of less than 100 bp. The ltered paired-end reads were aligned to the reference genome of CS version 1.0 [8] using BWA version 0.7.17 [53] with default options. Paired-end reads putatively generated by PCR duplications were removed using SAMtools version 1.9 [54]. The average number of aligned reads and genome coverage of short reads over the reference genome were calculated using BBmap version 37.77 [55]. SNP and indel calling was conducted using bcftools version 1.9 [56].
The command lines used for SNP and indel calling were described in Additional le 1: Table S3. Annotations of SNPs and indels were conducted using SnpEff version 4.3t [57].

Detection of large deletions induced by gamma irradiation
The depth-of-coverage for each genome position was estimated by SAMtools version 1.9. To detect deletions from the distribution of depth-of-coverage over chromosomes, moving average of depth-ofcoverage was calculated with a window size of three Mbp and a step size of one Mbp using Python version 3.7.1. The R package ggplot2 was used to visualize the distribution of the moving average of depth-of-coverage, SNP positions, and SNP density over chromosomes. Moving average of differences in depth-of-coverage (∆depth) between the wild-type and the mutants were also calculated. To estimate 95% con dence interval (q95) and 99% con dence interval (q99), the moving average of depth-ofcoverage was randomly extracted from each sample, and ∆depth between the samples was calculated. After this operation was repeated 4,000 times, the top 5% of ∆depth was designated as q95 and the top 1% of ∆depth as q99.
To evaluate how much depth-of-coverage is enough to detect large deletions induced by gamma irradiation, simulations was conducted by subsampling short reads using seqtk version 1. The command lines, packages, and scripts for the above analyses of deletions were described in Additional le 1: Table S3. The scripts of R and python are available from the Git-Hub repository (https://github.com/ShoyaKomura/WheatGammaIrradiationMutGenomeSeq).

Marker constructions for deletions and deleted genes
The borders of the large deletion in the PHS mutant of "Kitahonami" were estimated based on the distribution of the moving average of depth-of-coverage. The primer was constructed inside the deletion. Two primers were designed outside the deletion. The primers of Pina-D1 and Pinb-D1 described in Miki et al. [59] was used to con rm presence/absence of these genes in the wild-type, the grain hardness mutant, and their F 2 progenies. Sequence information of primers for markers and genes are listed in Additional le 1: Table S2. PCR condition of the markers and genes was as follows. After pre-denaturing of 94 °C for 2 minutes, 40 cycles of denaturing of 94 °C for 20 seconds, annealing of 60 °C for 30 seconds, and extension of 68 °C for 30 seconds were conducted, and then post-extension of 68 °C was performed for one minute. For PCR ampli cation of the co-dominant marker, the three primers were mixed before application of PCR.

Results
The hard grain mutant of Japanese elite wheat cultivar "Kitahonami" To con rm the capability of whole-genome resequencing approach to detect large deletions caused by gamma irradiation in the hexaploid wheat, we focused on grain hardness because grain hardness is regulated by Pina-D1 and Pinb-D1 genes, loss of which causes the grain hardness change from soft to hard [15,16]. Given that "Kitahonami" produces soft grains, the hard grain mutant was assumed to have mutations in the Pina-D1 and Pinb-D1 regions on the distal region of short arm of chromosome 5D. Based on the grain colour and the Single Kernel Characterization System (SKCS) value, we obtained the grain hardness mutant "30579" from the gamma-irradiated mutant population of "Kitahonami" (Fig. 1). Whereas the SKCS value of the wild-type is 24, the grain hardness of "30579" showed 102 of SKCS value (Table 1). Typically <40 and >70 of SKCS value is regarded as soft and hard grain, respectively. Therefore, the grains of mutant "30579" is hard. Grain protein content was also slightly elevated in the mutant (Table 1). Agronomic characteristics of the mutant "30579" were also compared with those of the wildtype. Grain weight (kg/10a) and 1000-grain weight of the mutant "30579" were low to those of the wildtype, indicating that the yield of the mutant was inferior to that of the wild-type. Grain protein content (%) 9.4 10.2 1 Maturity stage was examined in 2019-2020 season.
Comparative genomics of "Kitahonami" and its hard grain mutant "30579" revealed large deletions induced by gamma irradiation To identify the causal genome region of hard grain of the "Kitahonami" mutant "30579," genome resequencing of the wild-type and the mutant "30579" were performed. In total, 2.3 to 2.9 billion reads were obtained. After removing short reads with bad quality and PCR duplicates, we obtained 1.8 to 2.3 billion reads aligned to the reference genome of CS (Table 2). Average depth-of-coverage was 15.51 for the wild-type and 18.85 for the mutant. The aligned short reads covered over 94% of the reference genome sequences. The number of single-nucleotide polymorphisms (SNPs) and short indels between "Kitahonami" and CS were 30,597,871 and 1,328,442, while the number of SNPs and short indels between the mutant "30579" and CS were 32,254,934 and 1,713,770.  1 An average base quality score per 4bp >=20 The ltered reads rate = Total ltered reads / Total reads × 100 2 The aligned reads rate = Aligned reads / Total ltered reads × 100 3 The ltered reads rate = Aligned reads after removing PCR duplicates / Total ltered reads × 100 Page 9/32 4 The sets of read data were used for the simulations. The number in the paratheses indicates approximate depth-of-coverage.
SNP density distribution over chromosome for the mutant "30579" was almost identical to the wild-type (Fig. 2). The SNP density for the wild-type and the mutant "30579" was unevenly distributed over the chromosomes. In chromosomes 1A, 2A, 1B, 2B, 3B, 5B, and 6B, "Kitahonami" was genetically divergent from CS, while the chromosomes 3A, 4A, 6A, 4B, and 7B were genetically close in the two cultivars, particularly in their proximal region. D genome chromosomes showed less divergence than the other genomes between "Kitahonami" and CS. Less SNP density in the 420-450 Mbp position of chromosome 2A and the short arm of chromosome 5D was uniquely detected in the mutant "30579." Large deletions are often observed in gamma-irradiated mutants [29,30]. If a large deletion causes Pina-D1 and Pinb-D1 to be lost, a decrease of depth-of-coverage in the short arm of chromosome 5D should be uniquely observed in the mutant "30579." To identify genomic regions where depth-of-coverage uniquely decreased in the mutant, we analyzed moving average of depth-of-coverage per three Mbp over the chromosomes (Fig. 3). Depth-of-coverage in the mutant was close to zero in 420-450 Mbp position of chromosome 2A, around 90 Mbp position of chromosome 4B, and 0-130 Mbp position of chromosome 5D, indicating that these chromosomal regions were uniquely deleted from "30579" mutant. The genome regions of the mutant "30579" showing less SNP density than the wild-type (Fig. 2) were consistent with these deleted regions. The distribution of difference of depth-of-coverage (∆depth) per three Mbp between the wild-type and the mutant "30579" was also visualized over the chromosomes (Fig. 3). Peaks beyond the 99% con dence interval were observed at the uniquely deleted regions on chromosomes 2A, 4B, and 5B of the mutant "30579." Several sharp peaks of ∆depth over 99% con dence interval corresponded to regions with irregularly deep depth-of-coverage, where mapped reads were derived from repetitive sequences. Since the chromosome 5D region with the large deletion corresponded to the genome region containing Pina-D1 and Pinb-D1 [12], the hard grain mutant "30579" lost Pina-D1 and Pinb-D1.
Genotyping of the mapping population with indel markers con rmed causal regions of hard grains in mutant "30579" We developed a F 2 (n = 72) population of a cross between the wild-type and the mutant "30579" to evaluate segregation of grain hardness. The SKCS analysis was conducted to examine grain hardness of the seeds harvested from the F 2 population. Of the tested lines, 22, 42, and eight lines showed soft, intermediate, and hard grain phenotypes, respectively (Fig. 4a). To validate which deletions were linked to the hard grain phenotypes, indel markers were designed for each deletion (Fig. 4b). Genotyping of 22 soft grain lines and eight hard grain lines was performed using the indel markers. The indel marker of chromosome 5D was linked to the hard grain phenotypes, whereas the indel markers of chromosomes 2A and 4B were not (Fig. 4c). This result indicates that the large deletion on chromosome 5D caused "Kitahonami" to change from soft to hard.
Characteristics of the PHS-tolerant mutant of "Kitahonami" A PHS-tolerant mutant, called "28511," was selected from "Kitahonami" mutant population. The PHS tolerance test showed that the mutant "28511" was more tolerant to PHS than the wild-type (Fig. 5, Table  3). PHS tolerance of the mutant "28511" was comparable to that of "Kitakei 1831," used as a control variety for PHS tolerance. In the seed dormancy test, the mutant "28511" showed a lower germination rate than the wild-type under both 10 °C and 15 °C conditions in three seasons ( Table 3). The mutant "28511" had higher germination rate than "Kitakei 1831" in 2016-2017 season whereas had similar germination rate to "Kitakei 1831" in 2017-2018 season. Agronomic and quality characteristics of the mutant "28511" were also evaluated ( Table 4). Maturity stage, our yield, and our color of the mutant "28511" were almost the same as those of the wild-type. However, grain weight (kg/10a) and 1000-grain weight of the mutant "28511" were inferior to those of the wild-type. Moreover, the mutant "28511" showed higher our protein content than those of the wild-type. Whole genome resequencing of the PHS-tolerant mutant "28511" identi ed a large deletion on the long arm of chromosome 3B To identify the causal region of PHS tolerance of the mutant "28511," genome sequencings of the mutant were performed in the same method as the hard grain mutant "30579." Of the quali ed 2.1 billion reads, 99.9% were successfully aligned to the reference genome of CS (Table 2). Average depth-of-coverage for the mutant was 19.43 for the mutant. The number of SNPs and short indels between the mutant "28511" and CS were 35,368,664 and 1,865,274, respectively. Uneven distribution of SNP density over chromosome for the mutant "28511" was also observed as shown in that of the wild-type and the mutant "30579" (Additional le 2: Fig. S1).
To detect deletions in the mutant "28511," moving averages of depth-of-coverage per 3 Mbp and ∆depth per 3 Mbp were calculated over the chromosomes. At around 700 Mbp position of chromosome 3B, depth-of-coverage uniquely decreased in the mutant "28511" (Fig. 6). ∆depth also showed above 99% con dence interval. The reduced area extended to 67.8 Mbp, indicating that the mutant "28511" had a large deletion at the long arm of chromosome 3B. Such a remarkable reduction of depth-of-coverage was not observed in the other chromosomes.
Association between the large deletion at chromosome 3B and PHS tolerance To validate the large deletion at chromosome 3B, we constructed a co-dominant marker to detect the deletion (Fig. 7a, Additional le 1: Table S1). Primer pairs (pre-and post-deletion primers) were designed outside the deletion boundary. Another primer (in-deletion primer) was designed inside the deletion boundary. In wild-type, a 712 bp PCR fragment ampli ed by in-deletion primer and post-deletion primer was observed, whereas in the mutant, a 414 bp PCR fragment ampli ed by pre-deletion primer and postdeletion primer was observed (Fig. 7b). When the CS nulli-tetrasomic line of nulli-3B tetra-3D, the chromosome 3B of which is replaced with chromosome 3D, was used, no PCR fragment was detected, indicating that this marker was speci c to the deletion detected on chromosome 3B.
To con rm whether the large deletion was the causal genetic factor of PHS tolerance, we tested germination rate under the three conditions, three, seven, and nine days at 15 °C after seed sowing, for F 2 segregation population from a cross between the wild-type and the mutant "28511." Also, we examined the genotype for the F 2 population by using the co-dominant marker, whether the large deletion was associated with low germination rate. Under every condition, the wild-type genotype had signi cantly high germination rate, followed by the heterozygous genotype, and the mutant genotype had the lowest germination rate (Fig. 7c). This result indicates that a large deletion was associated with the low germination rate, which could facilitate PHS tolerance of the mutant "28511." Since ABA sensitivity rather than the ABA level in seeds regulates seed dormancy, wheat plants with low ABA sensitivity decrease seed dormancy, resulting in PHS [31]. To test the ABA sensitivity of the wild-type and the mutant "28511," suppression rate of germinated root elongation under exogenous ABA treatment was examined. No signi cant difference in the suppression between "Kitahonami" and the mutant "28511" was detected (Fig. 7d), implying that they have similar ABA sensitivity.
Vp-B1, associated with PHS tolerance [25,26], was located in the long arm of chromosome 3B where the large deletion was detected in the mutant "28511." A gene encoding GRAS family transcription factor was also found as a candidate gene of PHS tolerance in this region. Gibberellic acid (GA) promotes seed germination by inducing biosynthesis of α-amylase and protease in aleurone layers [32]. GRAS family transcription factor SCARECROW-LIKE 3 (SCL3) is a positive regulator of GA signaling in A. thaliana [33]. We designed two markers for each gene to test whether these genes were lacking in the mutant (Fig. 8a).
The markers for Vp-B1 (Vp-1B_2 and Vp-1B_3) showed no ampli cation in the mutant "28511" and the nulli-3B tetra-3D line. The markers for GRAS family transcription factor (GRAS-TF_2 and GRAS-TF_3) exhibited no ampli cation in the either the mutant "28511" or the nulli-3B tetra-3D line (Fig. 8b). These results con rmed that Vp-B1 and GRAS family transcription factor were deleted in the mutant "28511." Veri cation of detection ability in the gamma-irradiated wheat mutant with whole genome sequencing To reduce resequencing cost, it is essential to know how much depth-of-coverage is necessary to detect large deletions in a gamma-irradiated mutant genome. By subsampling short reads and adjusting average depth-of-coverage per genome, we evaluated the deletion detection power of the method based on depth-of-coverage, visualized over chromosomes. Four depth-of-coverage conditions, 2.5x, 5x, 10x, and 15x were tested using short reads from the PHS mutant "28511" and the wild-type (Fig. 9). The distributions of moving average of depth-of-coverage and ∆depth over the chromosomes were almost identical among all the conditions. The 67.8 Mbp deletion on chromosome 3B was signi cantly detected under all the depth-of-coverage conditions, implying that over 2.5x depth-of-coverage was enough to detect the large deletion in the gamma-irradiated wheat mutant.
In addition, to clarify how long deletions could be detected under the 5x depth-of-coverage conditions, simulation of deletion length was conducted by changing deletion length to 40 Mbp, 20 Mbp, 10 Mbp, 5 Mbp, 3 Mbp, and 1 Mbp (Fig. 10). When the deletion length was 1 Mbp or more, ∆depth was beyond 99% con dence interval. The length of the detectable deletion was found to depend on the window size of moving average. If window size was larger than the length of targeted deletion, the detectability of deletion decreased. For example, when the target deletion length was 3 Mbp, the estimated ∆depth for the 3 Mbp or 1 Mbp window size was over 99% con dence interval, but the estimated ∆depth for more than 3 Mbp window size was not. Another peak of ∆depth around 200 Mbp position on the chromosome.
This detected peak was caused by repeats derived from transposable elements. By con rming both distributions of depth-of-coverage and ∆depth, such irregular peaks can be distinguished from deletions.

Characterization of mutations detected in the gamma-irradiated wheat mutants
Nucleotide substitutions and small indels with less than 25 bp, which was the maximum length detected by indel calling of bcftools, were also detected in the gamma-irradiated wheat mutants (Table 5, Additional le 1: Table S2). Between the wild-type and the grain hardness mutant "30579," 2,412 SNPs and 329 small indels were detected. Between the wild-type and the PHS-tolerant mutant "28511," 2,715 SNPs and 266 small indels were detected. Gamma irradiation is assumed to cause these SNPs and indels. The ratio of transition and transversion (Ts/Tv ratio) between the wild-type and the mutants is lower than that between the wild-type and CS and is regarded as a natural variation. Over 98% of SNPs were detected in the intergenic regions. SNPs between the wild-type and the mutants were more randomly distributed over chromosomes compared with SNP density between wild-type "Kitahonami" and CS (Fig.   2, Fig. 11, Additional le 2: Fig. S1). The limited natural variations on chromosomes 3A, 4A, 5A, 6A, 4B, 7B, and all the D genome chromosomes were observed between these two cultivars, while the putative mutations induced by gamma irradiation covered these chromosomes.

Discussion
The usefulness of genome sequencing-based coverage analysis for deletions and duplications in the gamma-irradiated wheat mutants The present study showed that genome sequencing-based coverage analyses of the gamma-irradiated wheat mutants could e ciently detect large deletions. Combined with segregation analyses in progeny population, deletions responsible for mutant phenotypes were successfully identi ed. The hard grain characteristic needs the expression of both PINA and PINB. The con rmation of the deletion of Pina-D1 and Pinb-D1 at the Ha locus in the grain hardness mutant "30579" demonstrates the usefulness of the comparative analysis of depth-of-coverage between wild-type and the mutant. In addition to the detection of large deletions, the analysis of depth-of-coverage allows detecting large duplications in gammairradiated wheat mutants. A putative large duplication was detected on the short arm of chromosome 1B in "30579" (Fig. 3). Depth-of-coverage is used as a parameter for detecting copy number variations [34] and segmental duplications [35]. The approach based on depth-of-coverage ts the detection of deletions and duplications in the gamma-irradiated wheat mutants.
The simulation results showed that short sequencing read data corresponding to 5x depth-of-coverage over genome was enough to detect ve Mbp deletions visually. Furthermore, when ∆depth between wildtype and the mutant was used, even one Mbp deletion could be detected with over 95% con dential interval. The approach of ∆depth for the deletion detection is more powerful and provides evaluation criteria based on con dence interval. Detecting deletions in the gamma-irradiated mutants does not require high-coverage sequencing, reducing the cost of genome sequencing.

Characteristics of deletions and nucleotide mutations in the gamma-irradiated wheat mutants
The effectiveness of this approach could be unique to polyploid wheat with large genome. Double- The continuous existence of large deletions over generations could be associated with polyploidy. The hexaploid wheat has three homoeologous genes located on each of A, B, and D genomes. Even when one of the homoeologous genes are lost due to Mbp deletions, the other corresponding homoeologous genes can mitigate deleterious effects on viability. The aneuploid stocks are established in common wheat. The nullisomic lines, which lost one whole chromosome, are used for genetic research in wheat [38]. The deletion stocks of common wheat are also generated [39]. The large deletions in these stocks are stably inherited by the offspring. The large deletions above one Mbp generated by gamma irradiation are highly likely to be stably transmitted to the offspring due to compensation of the homoeologous chromosomes.
It is speculated that two double-strand breaks followed by removal of the large fragment and subsequent re-ligation at two separated ends of the same chromosome produce a large deletion [36,40]. Two double strand breaks could cause the large deletions inside the chromosomal ends in the tested gammairradiated mutants.
Genome sequencing of the gamma-irradiated wheat mutants also detected SNPs and small indels putatively induced by gamma irradiation. The random distribution of SNPs and small indels supports that the gamma irradiation induced these detected SNPs and indels (Fig. 11). In natural nucleotide variations, transition is more frequently observed than transversion due to differences in the ring structures of nucleobases and amino acid substitution frequency between transition and transversion, although transversion has twice as many nucleotide changes as transition. The Ts/Tv ratio of putative mutations introduced by gamma irradiation was lower than that of natural variations, although transition still occur more frequently than transversion (Table 5). Since mutations generated by gamma irradiation do not have evolutionary selective constraints on amino acid substitutions, the observed Ts/Tv ratio in the mutations generated by gamma irradiation could re ect only the differences in the ring structures of nucleobases. Given that the assembly size of the wheat reference genome is 14.5 Gbp [8] the number of SNP and indel per site are estimated to be ~2 x 10 -7 and ~2 x 10 -8 , respectively, although these values are potentially underestimated. Given that high impact SNPs and indels that potentially in uence phenotypes were rare in the "Kitahonami" mutants (Table 5), structural variations such as large deletions and duplications could be the main cause of phenotypic changes in the gamma-irradiated mutants of hexaploid wheat.
The causative gene of the PHS tolerance of the "Kitahonami" mutant "28511" Whole genome distribution of depth of coverage in the PHS tolerance mutant "28511" clari ed the 67.8 Mbp large deletion at chromosome 3B (Fig. 6). The co-dominant marker designed for detecting the 67.8 Mbp deletion revealed a statistically signi cant association between this deletion and PHS tolerance in the F 2 segregating population (Fig. 7c). Two PHS-related genes, GRAS family transcription factor (TraesCS3B01G441600) and Vp-1, are located in the deleted region. The deletion of these genes was con rmed by no PCR ampli cation of these genes from genomic DNA of the PHS-tolerant mutant.
GRAS family has the highly conserved C-terminal GRAS domain. A. thaliana has 33 genes encoding GRAS family transcription factor, which includes DELLA protein family involved in GA signaling in plants [33,41,42]. Rht-B1, which causes semi-dwar sm in wheat and led to the wheat green revolution, also belong to GRAS family [43]. TraesCS3B01G441600 is close to AtSCL1 and AtPAT1 and is separated from Rht-B1 and the gene encoding the positive regulator of GA signaling, AtSCL3 [33] (Additional le 2: Fig.  S2). TraesCS3B01G441600 can be classi ed into AtPAT1 subfamily [44]. AtPAT1 is involved in phytochrome light signaling rather than GA signaling. OsCIGR1 and OsCIGR2, rice GRAS family transcription factors belonging to AtPAT1 subfamily, are induced by exogenous GA in rice suspension culture but cannot be involved in mechanisms of α-amylase expression in aleurone layer [45]. Considering the function of genetically close homologs in A. thaliana and rice, TraesCS3B01G441600 is unlikely to be the causative gene of the PHS-tolerant mutant.
The Vp-1 positively contributes to seed maturation, dormancy, and ABA sensitivity in maize, A. thaliana and wheat [25,46,47,48], repressing germination. The loss of function of Vp-1 in maize and ABI3, the homolog of Vp-1, in A. thaliana arrests embryo maturation and decreases ABA sensitivity, resulting in precocious germination [46,47,48]. The "Kitahonami" mutant "28511" lost Vp-1 gene on chromosome 3B but increased PHS tolerance. In addition, the "Kitahonami" mutant "28511" did not reduce ABA sensitivity. These observations contradict those in diploid plants, maize, and A. thaliana. Given that hexaploid wheat has three homoeologous genes of Vp-A1, Vp-B1, and Vp-D1, encoding full-length protein [49], VP-A1 and VP-D1, could compensate for the VP-B1 function, preventing the reduction of ABA sensitivity and blocking precocious germination. To con rm this hypothesis, further experiments are needed.

Conclusion
In this study, we selected the grain hardness mutant "30579" and the PHS-tolerant mutant "28511" form gamma-irradiated mutants of Japanese elite cultivar "Kitahonami" and performed whole-genome resequencing of these mutant lines. The comparative analysis of depth-of-coverage between the wildtype and the mutants identi ed ~130 Mbp and 67.8 Mbp deletions in the mutants "30579" and "28511," respectively. These deletions were tightly linked to the mutant phenotypes in the progeny populations generated by crosses between wild-type and the mutants. The simulation analyses showed that 2.5x depth-of-coverage was enough to detect large deletions in the gamma-irradiated hexaploid wheat mutants. These results indicate that the approach of short-read-based genome sequencings of the gamma-irradiated mutants is cost-effective in designing genetic markers for agriculturally bene cial traits to breeding of wheat. Morphology of the grain hardness mutant "30579" Grains of the mutant "30579" (left) and wild-type "Kitahonami" (right).

Figure 2
Page 24/32 Distribution of SNPs between wheat cultivars "Kitahonami" and "Chinese Spring" or between the grain hardness mutant "30579" and "Chinese Spring" over the chromosomes. SNP density was estimated by counting number of SNPs per 10 Mbp. SNP density is displayed as the gradation color from white to black. The higher the density, the blacker it is. In the regions with an asterisk, the mutant shows less SNP density compared with the wild-type.

Figure 3
Depth of coverage of wild-type "Kitahonami" and the grain hardness mutant "30579" over chromosomes Distribution of moving average of depth-of-coverage and difference of depth-of-coverage (∆depth) between wild-type "Kitahonami" and the mutant "30579" was visualized over chromosomes. q95 and q99 indicate 95% and 99% con dence intervals, respectively. Window size is 3 Mbp, and step size is 1 Mbp.
An asterisk indicates a region with a large deletion. # indicates a large duplication.  Pre-harvest sprouting tolerance of "Kitahonami" mutant "28511" Eight spikes of wild-type "Kitahonami" (A) and the mutant "28511" (B) are shown.

Figure 6
Depth of coverage of wild-type "Kitahonami" and the pre-harvest sprouting mutant "28511" over chromosomes Distribution of moving average of depth-of-coverage and difference of depth-of-coverage (∆depth) between wild-type "Kitahonami" and the mutant "28511" was visualized over chromosomes.
q95 and q99 indicate 95% and 99% con dence intervals, respectively. Window size is 3 Mbp, and step size is 1 Mbp. An asterisk indicates a region with a large deletion.

Figure 7
Associations between PHS tolerance and the deletion on chromosome 5D in the F2 population (a) A diagram of co-dominant markers to detect the deletion on chromosome 5D. A 712 bp PCR fragment is ampli ed by "in deletion" primer and "post deletion" primer is expected in wild-type Kitahonami (WT). A 414 bp PCR fragment is ampli ed by "pre deletion" primer and "post deletion" primer in the PHS tolerance  Relationship between PHS tolerance and presence/absence of the candidate genes on chromosome 3D.
(a) Primer positions of Viviparous-B1 and the GRAS transcription factor gene to con rm presence/absence of these genes in "Chinese Spring" (CS) deletion lines at chromosome 3B (Endo and Gill, 1996). (b) The two primer sets were designed to amplify two fragments for each gene. PCR ampli cations of Viviparous-1B and the GRAS transcription factor gene in wild-type "Kitahonami" (WT), the PHS-tolerant mutant "28511", CS, the CS 3BL ditelosomic line (CS ditelo 3BL), and the CS nullitetrasomic line of nulli-3B tetra-3D (CS N3BT3D). Three independent tests of seed dormancy were conducted.

Figure 9
The tests of deletion detection power under different depth of coverage conditions Distribution of moving average of depth-of-coverage and difference of depth-of-coverage (∆depth) between wild-type "Kitahonami" (WT) and the PHS-tolerant mutant "28511" on chromosome 3B under four depth coverage conditions: 15x, 10x, 5x, and 2.5x. are simulated. q95 and q99 indicate 95% and 99% con dence intervals, respectively. The number of reads of 15x, 10x, 5x, and 2.5x was described in Table 2.

Figure 10
The