Identification of Genome-Wide Variations among Three Elite Restorer Lines for Hybrid-Rice

Rice restorer lines play an important role in three-line hybrid rice production. Previous research based on molecular tagging has suggested that the restorer lines used widely today have narrow genetic backgrounds. However, patterns of genetic variation at a genome-wide scale in these restorer lines remain largely unknown. The present study performed re-sequencing and genome-wide variation analysis of three important representative restorer lines, namely, IR24, MH63, and SH527, using the Solexa sequencing technology. With the genomic sequence of the Indica cultivar 9311 as the reference, the following genetic features were identified: 267,383 single-nucleotide polymorphisms (SNPs), 52,847 insertion/deletion polymorphisms (InDels), and 3,286 structural variations (SVs) in the genome of IR24; 288,764 SNPs, 59,658 InDels, and 3,226 SVs in MH63; and 259,862 SNPs, 55,500 InDels, and 3,127 SVs in SH527. Variations between samples were also determined by comparative analysis of authentic collections of SNPs, InDels, and SVs, and were functionally annotated. Furthermore, variations in several important genes were also surveyed by alignment analysis in these lines. Our results suggest that genetic variations among these lines, although far lower than those reported in the landrace population, are greater than expected, indicating a complicated genetic basis for the phenotypic diversity of the restorer lines. Identification of genome-wide variation and pattern analysis among the restorer lines will facilitate future genetic studies and the molecular improvement of hybrid rice.


Introduction
As the main staple food for more than half of the world's population, rice (Oryza sativa L.) is one of the most important food crops. In 1973, the field production of Indica hybrid rice succeeded when Chinese rice breeders completed the three-line breeding system [1]. A land area of approximately 130,000 hm 2 was soon developed for hybrid rice cultivation, greatly increasing rice yield in China. In the three-line breeding system, the cytoplasmic male sterility (CMS) line is crossed with the restorer line to produce the F 1 hybrid rice, and with the maintainer line for self-reproduction. The restorer line is widely considered as being key to further improve the resistance, yield, quality, and heterosis of hybrid rice [1,2]. IR24, an elite rice variety introduced in China by the International Rice Research Institute, was the most common restorer line used during the 1970s until the early 1980s. MH63, which was developed from a cross between IR30 and Gui630 [1], is thus far the most widely used restorer line in China. Its popularity can be attributed to the fact of being a coparent of ShanYou63, the largest hybrid rice acreage that has created substantial economic and social benefits. SH527 is a heavy-panicle restorer line bred in the 1990s [3]. More than 40 new elite hybrid rice varieties have been bred using SH527 as the male parent, among which 5 were chosen for super hybrid rice development. At present, many hybrid rice varieties generated from SH527 are widely grown in China. IR24, MH63, and SH527 thus represent the first-, second-, and third-generation restorer lines, respectively, of the three-line breeding system. Although they are all significant backbone parents at different stages of hybrid rice development, their field performances and combining abilities differ considerably. Further research on the genetic diversity of these lines, which might be related to their varying performances, can improve our understanding of restorer lines and promote improved restorer line selection and super hybrid rice breeding.
The genomic sequences of the Japonica cultivar Nipponbare [4] and the Indica cultivar 9311 [5] were recently released. The availability of high-throughput sequencing technology not only increases sequencing throughput but also allows for simultaneous sequencing of a large number of samples [6,7] in addition to decreasing time and cost. These merits open the door to highthroughput re-sequencing and genotyping of various rice strains. A genetic map with a resolution of recombination breakpoints within an average of 40 kb were previously constructed for ,150 rice recombinant inbred lines by utilizing whole-genome re-sequencing data generated using the Illumina Genome Analyzer [8]. Six elite maize inbred lines, including the parents of the most productive commercial hybrid in China, were recently re-sequenced and more than 1,000,000 SNPs, 30,000 indel polymorphisms and 101 low-sequence-diversity chromosomal intervals were uncovered in the maize genome [9]. Huang et al. [10] identified approximately 3.6 million single-nucleotide polymorphisms (SNPs) by sequencing 517 rice landraces and constructed a high-density haplotype map of the rice genome. Moreover, they pioneered genome-wide association studies for 14 agronomic traits of the O. sativa indica subspecies. Molecular marker screening has suggested narrow genetic backgrounds for rice restorer lines [11,12], which play a vital role in hybrid rice production. However, the current lack of information on genetic variation over the entire genome has limited further research into this topic. In the present study, we conducted re-sequencing and genome-wide variation analysis of IR24, MH63, and SH527 using the Solexa sequencing technology. Identification of genome-wide single-nucleotide polymorphisms (SNPs), insertion/deletion polymorphisms (InDels), and structural variations (SVs) as well as pattern analysis among these lines has the potential to provide valuable resources for future genetic studies and the molecular improvement of hybrid rice.

Results
Field performances of the restorer lines and their hybrid descendants IR24, MH63, and SH527 (Fig. 1A) are considered hybrid rice core restorer lines because of the large number of elite commercial hybrid rice cultivars and useful restorer lines bred and generated from them. Based on their cross genealogies, MH63 and SH527 were both indirectly generated from IR24 (Fig. 1B), indicating that these three lines originate from the same restoring genes. We examined the field performances of these lines by selfing ( Table 1). Performances of the hybrid rice made by crossing these three lines with six other widely used CMS lines, namely, G3A, Zhongjiu A, II-32A, G46A, 92A, and Chuangu A, were also examined ( Table 1). No obvious differences were found in the yield components of MH63 and IR24 except for plant height, while the hybrid rice of MH63 was significantly different from that of IR24 in growth period, plant height, panicles per plant and seed setting rate. Between SH527 and IR24, significant differences were detected in plant height, panicles per plant and 1000-grains weight. Significant differences between their hybrid rice were also detected in growth period, plant height, seed setting rate and 1000-grains weight. In general, from the breeding stage of IR24, MH63 to SH527, combinations of these changes lead to an apparent yield increase for hybrid rice, although no obvious yield differences were found in the restorer lines themselves. Since the yield increase was evaluated on the average performance of hybrid rice generated from these three restorer lines with several common CMS lines, the yield increase of hybrid rice reflect an obvious genetic improvement of the restorer lines, possibly by improving the combining ability of the restorer lines.

Genome sequencing and variation identification
The genotypes of IR24, MH63, and SH527 were determined with approximately 10-fold coverage by genome sequencing using the Solexa sequencing technology. According to the protocol, three DNA libraries were constructed and 12.48 G bases were generated (raw sequence data obtained have been deposited in the NCBI Short Read Archive with accession number SRP006823). The alignment of reads was used to build consensus genome sequences for each rice accession. Furthermore, approximately 10.78 G high-quality raw databases were aligned with the reference sequence of cultivar 9311 using SOAPaligner [13] (http://soap.genomics.org.cn/). In total, an effective depth of 306 coverage was achieved, with an average of 106 for each restorer line ( Table 2). The resulting consensus sequence of each rice accession covered approximately 84.8% of the reference genome (84%-85.99%), indicating a close relationship between the samples and cultivar 9311.
SNPs, InDels, and SVs were then examined with SOAPsnp11 and SOAPsv using a conservative quality filter pipeline [14], yielding 267,383 SNPs from the genome of IR24, 288,764 SNPs from that of MH63, and 259,862 SNPs from that of SH527 (Table 3, http://rice.sicau.edu.cn/re-sequencing/variation/9311. rar). These outcomes resulted in a non-redundant collection of 568,787 SNPs after excluding the shared SNPs of each sample by synteny analysis ( Fig. 2A-2C). In total, 100,095 InDels ranging from 1 to 5 bp in length and 5,561 SVs across the whole genome were identified. Because of inherent relationship between the samples, the overall genome diversity among these re-sequenced elite restorer lines was much lower than that reported for a more diverse population [10], which is also in accordance with the close relationship among the three lines revealed by genealogy analysis. A phylogenetic tree [15] was constructed using several authentic collections of SNPs. An extremely closed genetic relationship was observed between sequencing samples, and a relatively distant relationship was observed between samples and the reference (Fig. 2D), which is consistent with a previously reported result of low genome diversity among rice restorer lines [11,12]. The frequencies of SNPs, InDels, and SVs for each sample were plotted at a 100 kb sliding window with a step size of 50 kb along each chromosome. SNP/InDel/SV frequency was defined as the corresponding number of SNPs/InDels/SVs divided by the number of nucleotides within the 100 kb interval, excluding the uncovered nucleotides. Each sample was compared with the corresponding intervals to identify regions that showed nonrandom variation frequencies. In total, 227/936 SNP high/low regions, 298/889 InDel high/low regions, and 188/1899 SV high/low regions were identified between IR24 and MH63; 339/ 914 SNP high/low regions, 440/1,030 InDel high/low regions, and 267/2,052 SV high/low regions were identified between IR24 and SH527; and 507/825 SNP high/low regions, 523/1,266 InDel high/low regions, and 235/2,684 SV high/low regions were identified between MH63 and SH527. Out of these, 135/450 SNP high/low regions, 229/297 InDel high/low regions, and 87/1,058 SV high/low regions were found to be identical among the three restorer lines (Figs. 3 and 4).

Variations between samples
As differences between the samples (i.e., not between the samples and the reference) may reflect the genetic improvement of the recent restorer lines (such as SH527 and MH63) from older lines (such as IR24), an analysis of the variations and their distributions among the samples was performed.  Table 4. Furthermore, only 10 different SNPs and 12 different InDels (allelic pleomorphic loci with different nucleotides in each line) were identified by the variation consensus comparative analysis of the three sequenced lines, although large numbers of shared SNPs and InDels were found ( Table 5).
The SNPs in coding regions were analyzed to gain further insights into the potential functional effects of the detected SNPs (Table 6). Between IR24 and MH63, 13,160 shared SNPs, of which 2,290 were synonymous coding sequences (Syn CDS) and 2,902 were non-synonymous coding sequences (Non-syn CDS), and 291 different SNPs, of which 54 were Syn CDS and 99 were Non-syn CDS, were found. Between IR24 and SH527, 14,473 shared SNPs (2,522 Syn CDS and 3,366 Non-syn CDS) and 594 different SNPs (94 Syn CDS and 138 Non-syn CDS) were found in  (Table 7). Different CDS-located InDels were not detected. Three hundred thirty-one large-effect SNPs that were expected to affect the integrity of encoded proteins were also identified. These included changes introduced by premature termination codons (premature termination; 238 SNPs), elimination of translation initiation sites (ATG change; 11 SNPs), and replacement of nonsense with sense codons (stop change; 82 SNPs). Of these large-effect SNPs, only 10 SNPs (2 ATG changes, 5 premature terminations, and 3 stop changes) were observed from the different SNPs; the rest were from the shared SNPs (Table 8).
GO and PFAM analyses were further carried out for the shared and different SNPs (InDels) in genes between samples to explore gene functions. In both the shared and different SNPs (InDels), the top GOs were protein kinase activity, nucleic acid binding, protein binding, DNA binding, and catalytic activity (Fig. 5 and 6). Genes coding for leucine-rich repeats and NB-ARC domains were found to have a significantly higher ratio of nonsynonymous-tosynonymous SNPs than average. As these domains are common in proteins that mediate disease resistance in plants, our finding is consistent with these proteins being particularly diverse due to pathogen pressure.

Variation analysis on important rice genes
Several important rice genes related to yield, quality, resistance, and development processes were subjected to molecular cloning and functional analysis. Natural variations among the genes, which might explain the phenotypic differences of the sequenced sample, were then evaluated. A large number of SNPs (Table 9) were detected both in the DNA sequence and in the coding regions of genes related to disease/insect resistance, such as Pib [16], Xa1 [17], Pi9 [18], Xa21 [19], Xa26 [20] and Bph14 [21]. Although found to have many SNPs, genes related to rice developmental processes, yield, and quality, such as ALK [22], qSW5 [23], GS3 [24], Gn1a [25], HTD2 [26], GW2 [27] and EUI1 [28], had rare or no variations in the coding regions, which might explain the functional conservation. In addition, only a few InDels (or none in some cases) were found in the coding regions ( Table 10), suggesting that SNPs, not InDels, effectively contribute to functional variation of the genes. When compared to the 9311 sequence, a number of SNPs were found both in the DNA sequence (,60) and in the coding regions (,40) of Rf1a [29], a possible allelic gene for Rf4 [30], which is the major restoring gene of the WA-CMS line, while the sequence difference in this gene between the sequencing samples was limited. These variations may account for the differences between the sequenced samples (restorer lines) and the reference cultivar 9311(non-restorer lines) in terms of their restoring ability.

Discussion
In the present study, we conducted re-sequencing and genomewide variation analysis of three famous representative restorer lines, namely IR24, MH63, and SH527, with the aim of uncover genetic variation at a genome-wide scale by using the Solexa sequencing technology. Identification of genome-wide SNPs, InDels, and SVs, as well as pattern analysis of restorer lines can provide valuable resources for future genetic studies and the molecular improvement of hybrid rice.
We firstly used the 9311 [5] and Nipponbare [4] sequence as the reference genome, respectively. The genome size of 9311 is 374,545,499, of which the effective size is 359,401,158 (excluding the N bases in the reference). On the other hand, the genome size of Nipponbare is 382,150,945, of which the effective size is 372,089,805. When the Nipponbare genome was used as the reference, the number of SNPs detected was noticeably higher (data not shown). However, quality of original sequence data such as mapped bases, sequencing depth, and coverage decreased, rendering the SNP data less reliable. Given that genetic variations between restorer lines, not the japonica and indica rice varieties, underlie the mechanism of their phenotypic differences, the 9311 genome sequence was then used as the only reference for detecting SNPs, InDels, and SVs, and for assembling the consensus sequence to exclude the large amount of background variations   that account for differences between the japonica and indica rice varieties.
Interestingly, approximately 76,000, 71,000, and 76,000 heterozygous SNPs in IR24, MH63, and SH527, respectively, were identified throughout the whole rice genome, leading to an estimated heterozygosity rate of approximately 1.98-2.0610 24 , which is lower than that for other species, such as pandas [31] and humans [32]. The heterozygosity rate showed, to some extent, an un-purified genetic background of the sequenced rice varieties and indicated that the rice restorer lines still have high genetic variability, supporting the sporadic phenotypic variability of individuals observed within a rice line, even it is strictly selfpollinated. Thus we may speculate that, besides spontaneous mutations, genomic heterozygosity might also play a role in phenotypic variations. These results might also suggest that selfpollinated plants have the potential to maintain a relatively high heterozygosity rate. More plant lines should be studied to confirm this idea.
Here we report variations over the whole genome among elite rice restorer lines. Our results indicate that genetic variations among these lines, although far lower than those reported for a more diverse landrace population [10], are greater than expected, indicating a complicated genetic basis for the phenotypic diversity of the restorer lines. Although several candidate genes have been proposed to account for the varying performances of rice lines and selected for functional analysis, further analysis of more restorer lines is necessary to better understand the mechanism by which restorer lines are improved by breeding. Furthermore, several follow-up steps can be taken to pinpoint candidate genes that may contribute to phenotypic diversity in rice cultivars. This study therefore lays the groundwork for long-term efforts to uncover genes and alleles important for cultivar improvement in rice restorer lines.

Sampling
Seedlings of IR24, MH63, and SH527 and six other widely used CMS lines, namely, G3A, Zhongjiu A, II-32A, G46A, 92A, and  Chuangu A, were planted in the experimental field of the Rice Research Institute, Sichuan Agricultural University, Wenjiang. When they reached the flowering stage, these three restorer lines were crossed with the six CMS lines to obtain the F 1 hybrid rice. The three elite restorer lines, together with the F 1 hybrid rice, were then planted in the following year for phenotypic evaluation and field test. All the restorer lines and the F 1 hybrid rice were planted across 20 lines, with three replicates totaling 12 plants in each line. Eight middle plants of the 10 middle lines were surveyed, and data were recorded for statistical analysis. To compare the field performances of these elite restorer lines, we used a One-way ANOVA and LSD's test of DPS Software (http://www.chinadps. net/index.htm). To compare the contribution of restorer lines to their hybrids' field performances, we used a Two-way ANOVA and LSD's test of DPS Software (http://www.chinadps.net/index. htm) [33].

DNA isolation and genome sequencing
Total genomic DNA was extracted from the leaf tissues of one individual for each line using a DNeasy Plant Mini Kit (Qiagen). The DNA of each line was then randomly fragmented. After electrophoresis, DNA fragments of the desired length were gelpurified. Adapter ligation and DNA cluster preparation were performed and subjected to Solexa sequencing.

Read mapping
The raw pair-end (PE) sequencing reads were aligned to the 9311 reference genome sequence using SOAPaligner [13] under the following conditions: if an original read cannot be aligned to the reference sequence, the first nucleotide from the 59 end and two nucleotides from the 39 end will be deleted and then realigned to the reference. If the alignment still cannot be achieved, two more nucleotides from the 39 end will be deleted. The procedure was repeated until the alignment was available or the read was less than 27 bp long. Average sequencing depth and coverage were calculated using the alignment results.

Assembly of consensus sequences and SNP/InDel detection
Based on the alignment results, and taking into consideration the analysis of data characters, sequencing quality, and other factors influencing the experiments, a Bayesian model was applied to calculate the probability of genotypes with the actual data. The genotype with the highest probability was selected as the genotype of the sequencing individual at a specific locus, and a quality value was designated accordingly to reveal the accuracy of the genotype. Polymorphic loci against the reference sequence were selected from the consensus sequence and then filtered under certain requirements (e.g., the quality value must be greater than 20 and the result must be supported by at least two reads) using SOAPsnp [14]. Mapped reads that satisfied the PE requirements and   contained alignment gaps at one end were also used to detect the short InDels. The maximum gap length allowed in the alignments was 5 bp. Gaps that were supported by at least three gapped PE reads were extracted in InDel calling.

SV detection
According to the principle of PE sequencing, under normal situations, one read of PE should be aligned to the forward sequence and another should be aligned to the reverse. The distance between the two aligned positions at the reference should be in accordance with the insert size. Thus, the alignment of the two paired reads to the genome is regarded to be of normal direction and appropriate span. If the direction or span of the alignments of the two paired reads is different from that expected, then the region might have SVs. Abnormal PE alignments observed in our analysis were further analyzed by clustering and compared with previously defined SVs. In this manner, the SVs were detected using SOAPsv [14], with support from at least three abnormal PE reads. Currently, the types of SVs that can be detected include deletion, replication, reversion, and transposition, among others.

SNP annotation
The localization of SNPs in coding regions, noncoding regions, start codons, stop codons, and splice sites were based on the annotation of gene models provided by the Rice Genome Sequencing Project of 9311 [34]. The characterization of synonymous or non-synonymous status of SNPs within the CDS was conducted using Genewise version 30 [35]. The GO/PFAM annotation data were further used to functionally annotate each gene [36].

Variation frequency distribution
The frequencies of SNPs, InDels, and SVs for each sample were plotted over a 100 kb sliding window with a step size of 50 kb along each chromosome to explore the genomic distribution of DNA polymorphism in these lines [37]. The scanned regions were defined as high-or low-variation frequency regions if variation rates were higher than 4 fold or lower than 1/20th of the average rate over the whole genome (ARG), respectively. The deviation ratio (DR) of samples in a given window was first calculated as the sum of the ratio of each sample that deviated from the average rate, then the ARG was defined as the arithmetic average of all the windows across chromosomes.

Variations between samples
The SNPs/InDels/SVs detected for each individual line were further compared between samples to identify the shared and unique SNP/InDel loci. Only those loci for which at least one effective sequence read was mapped for every individual were selected for comparison. A phylogenetic tree was constructed using the MEGA4 software [15] based on these data on SNPs.