Whole genome de novo assemblies of three divergent strains of rice, Oryza sativa, document novel gene space of aus and indica

De novogenome assemblies and functional annotation

The three rice varieties were assembled using the ALLPATHS-LG whole genome assembler [47] using approximately 50× coverage of a 180 bp fragment library, approximately 30× coverage of a 2 kbp jumping library, and approximately 30× coverage of a 5 kbp jumping library (see Materials and methods). We selected this assembler based on its performance with these data compared with other assemblers and its high ranking in the Assemblathon I and II and GAGE evaluations [48]-[50]. The three assemblies were named Os-Nipponbare-Draft-CSHL-1.0, Os-IR64-Draft-CSHL-1.0, and Os-DJ123-Draft-CSHL-1.0, following nomenclature proposed by [51].

The assemblies were repeat-masked and annotated for protein-coding genes using the MAKER-P automated pipeline [52], combining both evidence-based and ab initio methods (Table S1 in Additional file 1). In addition to EST and full-length cDNA, we included as evidence the two published annotations of Nipponbare [51], and the published annotations of strains 93-11 and PA64s [28], thereby maximizing consistency and reducing bias of annotation across the three assemblies. Putative transposon-encoded genes were screened following analysis of InterPro domains (see Materials and methods), which flagged approximately 1% of initial gene calls in each of the three assemblies. Summary statistics for remaining genes are provided in Table 1 and in Table S2 in Additional file 1. Gene counts ranged from 37,758 (IR64) to 39,083 (Nipponbare), similar to the numbers reported by the Michigan State University (MSU) Rice Genome Annotation Project and Rice Annotation Project for the Os-Nipponbare-Reference-IRGSP-1.0 (39,102 and 35,681 respectively) [51]. Overall statistics for structural features, such as exons, introns, and coding regions were highly consistent between the three assemblies and with published annotations. For instance, average translated protein lengths compared across MSU, Rice Annotation Project, and the three de novo assemblies ranged from 280 to 288 amino acids (median values: 268 to 291 amino acids), suggesting that contiguity of the de novo assemblies did not limit ability to identify protein-coding genes. For each assembly, 61 to 62% of annotated loci possessed one or more InterPro domains and 77% showed homology to plant NCBI RefSeq genes.

Whole genome comparison to Nipponbare reference genomes

We evaluated the agreement between our de novo assemblies to the Nipponbare reference sequences using the GAGE assembly evaluation algorithm [50]. As expected, the de novo Nipponbare assembly very closely matches the reference Nipponbare sequence, with a 99.94% average identity and only 0.31% of the assembly not aligning to the reference (Tables 2, 3 and 4). Even at this very high agreement, there are several tens of thousands of small variations, and several hundred larger variations. These variations are a combination of true variations from our sample relative to the reference genome, of which we expect there to be few, and errors from ALLPATHS-LG when used with these libraries and coverage levels. Consequently, considering that the assembly has a 99.94% overall similarity, the upper-bound on the error rate of sequencing and assembling with ALLPATHS-LG is at most 0.06%.

The portions of the reference genome without any alignments from our Nipponbare assembly are scattered throughout the genome in 57,821 segments averaging 203 bp long. However, of this, only 301,525 bp are annotated to be within the coding sequence (CDS; 0.72% of the total CDS), and another 12,344 bp are annotated to be within non-coding exons. We further evaluated the unaligned regions by computing their read k-mer coverage from a sample of 400 million unassembled Nipponbare reads, and found the mean k-mer coverage of these regions exceeds 12,000×, while the mode k-mer coverage of the set is less than 100× (Figure S1 in Additional file 1). A full two-thirds (38,373/57,821) of these regions exceed 1,000× k-mer coverage, more than 10 times higher than unique segments of the genome. This implies the unassembled/unaligned regions are highly enriched for high copy repeats too complex to be assembled. In contrast, the genic regions are very well represented, suggesting it would be possible for a detailed analysis of the 'gene space' of the accessions from these assemblies.

Whole genome comparison to indicareference genomes

Finally, the comparison between the IR64 and DJ123 assemblies shows that they differ from each other nearly as much as either one differs from the Nipponbare reference sequence. These results suggest that the aus genome harbors a greater amount of novel variation than previously recognized. It also highlights the value of taking an unbiased, de novo assembly approach when evaluating genomic variation among varieties and subpopulations to capture genome-specific variations.

Pan-genome analysis

We next evaluated the 'pan-genome' of the three de novo assemblies to identify sequences that were conserved across the genomes as well as sequences specific to just one genome (see Materials and methods). Using the whole genome alignment information, we classified each base of each genome as being specific to that genome (unaligned to either other genome), or shared by one or both genomes. The majority of the assembled sequences (approximately 302 Mbp per genome) and exonic sequences (approximately 55.5 Mbp per genome), were shared among the three genomes, although 4.8 Mbp to 8.2 Mbp (423 kbp to 930 kbp exonic) were found to be genome-specific (Figure 2A). Since a gene sequence may be partially shared or partially genome-specific, we assigned each gene to the sector on the Venn diagram for which the majority of the exonic bases were assigned over all transcripts associated with each gene. For example, if 90% of a gene is shared among all three genomes, but 10% is genome-specific, we would assign it to the center (fully shared) sector under the majority rule. This will not necessarily characterize changes in gene function if critical protein domains are shared/unshared, but highlights the major trends between the lineages and discovers 297 to 786 genome-specific loci.

Using the same k-mer analysis techniques we applied for the reference analysis, we further classified the genome-specific bases as being unique or repetitive, using a threshold of 100× average k-mer coverage to classify unique sequences. From this, we identified only 1.2 Mbp to 1.5 Mbp of non-repetitive sequence specific to each genome, meaning that most of the genome-specific bases were actually repetitive (Table 8). Since repetitive sequences are also the most likely to be unassembled, as observed in our comparison to the reference genomes, we further examined the genome-specific exonic bases and refined our initial estimates to 555 kbp to 760 kbp of non-repetitive, genome-specific sequences intersecting annotated genes by at least 100 bp (Table 9). Note these segments may include flanking promoter and other regulatory regions in addition to the exons themselves. From this catalog, we selected 10 of the largest regions in each of the genomes for PCR validation, and were able to confirm the computational analysis with 100% success (Figure 3; Table S4a-c in Additional file 1).

Genome	Genome-specific	Regions	Mean ± standard deviation	Maximum size
Nipponbare	1,574,801 bp	2,250	699 ± 904	8,463
IR64	1,336,650 bp	1,702	785 ± 1066	10,400
DJ123	1,263,681 bp	1,569	805 ± 1,058	8,960

Identified sequences must be at least 100 bp long with no alignments to the other genomes, average between 10× and 100× k-mer coverage in that genome, and average below 10× k-mer coverage using the reads from the other two genomes.

Genome	Genome-specific	Regions	Mean ± standard deviation	Maximum size
Nipponbare	760,064 bp	779	975 ± 1,197	8,463
IR64	637,470 bp	583	1,093 ± 1,430	10,400
DJ123	555,507 bp	492	1,129 ± 1,409	8,960

Regions identified to be specific to a given accession using the criterion used in Table 8 and that intersect an annotated gene region by at least 100 bp.

For Nipponbare and IR64, we determined the positions of the non-repetitive segments along the different reference chromosomes, and found the segments were broadly distributed. For Nipponbare, we could localize 2,208 of the genome-specific regions, and found that one region occurred, on average, every 162 ± 362 kbp, following an approximately exponential distribution (data not shown). For IR64, we could localize 1,074 of the genome-specific regions, and found one region occurred, on average, every 338 ± 752 kbp, also from an approximately exponential distribution. The distributions suggest that the genome-specific bases are not highly localized, as an exponential distribution in spacing can occur if there is a uniform probability distribution of a site occurring at any position at random.

Genome-specific loci, as well as those shared between two genomes but not the third, exhibited shorter CDSs and greater novelty compared with genes shared among all three genomes (Figure 2B). For example, loci common to all genomes had a median coding length of 888 bp compared with median values ranging from 483 to 654 bp for the genome-specific gene sets. Likewise, the core, fully shared set of genes averaged 4.9 exons per transcript compared with a range of 2.9 to 3.0 genes per transcript amongst genome-specific genes. A smaller fraction of genome-specific loci contained InterPro domains compared with the core set (40% versus 63%), and fewer showed homology to plant RefSeq proteins (57% versus 79%). However, artifacts of inaccurate annotation may contribute to this trend [53], so we investigated if these differences were negatively influenced by assembly quality, especially if genome-specific genes tended to terminate in scaffold gaps more frequently than core genes. We observe a modest effect, and genes shared by all three strains have a median distance of 12 to 14 kbp (5′ and 3′ flanking distances), whereas genome-specific genes have a median distance of 8 to 10 kbp (Table S7 in Additional file 1). Only 7 to 12% of the genes have a flanking distance of less than 1 kbp, suggesting this is not a major factor in our analysis. Indeed, our results are similar to studies in yeast, Drosophila, and vertebrates that have found that novel and recently evolved genes tend to encode smaller proteins than conserved or ancient genes [53]-[55].

To characterize potential function of genome-specific genes we further examined genes with annotated InterPro domains. Notably, genes with domains related to disease resistance were the most prevalent type among genome-specific genes. For example, 12% of genes specific to IR64 possessed the NB-ARC motif (IPR002182), the central nucleotide-binding domain of plant R-genes. This domain, and others associated with R-genes, also prevailed among the DJ123-specific and Nipponbare-specific gene sets, accounting for 9% and 5% of genes, respectively. In contrast, only 0.35% of genes shared among all three genomes encode the NB-ARC domain. Genes shared between just two genomes showed intermediate frequencies of disease resistance genes (1.5 to 2.5%). Similar distributions were seen for genes classified with the Gene Ontology (GO) term 'defense response' (GO:0006952). These results are consistent with Ding et al. [14], who showed high levels of 'genome asymmetry' among R genes when comparing the Nipponbare and 93-11 reference assemblies. A large diversity of other protein domain classes, such as those associated with receptor and non-receptor protein kinases, transcription factors, metabolic enzymes, proteases, and transporters, were also found in the genome-specific gene sets. A complete listing of putative strain-specific genes, their InterPro domains, GO terms, and summary of homology search results are provided in Additional file 2. We anticipate these findings will greatly enhance the ongoing 3,000 rice genomes project [36] and other resequencing projects that had previously focused on single nucleotide variations relative to the Nipponbare reference.

Detailed regions

We chose four agronomically relevant regions of the rice genome that were previously reported to harbor differences among the three varieties or subpopulations to illustrate the utility of these high quality whole genome assemblies for understanding the variation in genome structure underlying salient phenotypic variants.

Sub1 locus

The Submergence 1 (Sub1) locus on chromosome 9 is a major QTL determining submergence tolerance in rice [33]. The Sub1 locus is a cluster of three genes encoding putative ethylene response factors. Sub1B and Sub1C are present in all rice accessions tested to date, while Sub1A may be present or absent. Originally identified in the aus accession FR13A, Sub1A appears to be found only within the Indica varietal group [33]. Sub1A has two alleles: Sub1A-1 is found in submergence-tolerant varieties, while Sub1A-2 is found in intolerant varieties. A haplotype survey in O. sativa varieties also identified nine Sub1B and seven Sub1C alleles [33].

In the IR64 and DJ123 de novo assemblies reported here the Sub1 locus lies within a single scaffold and haplotypes can be easily identified (Table S5 in Additional file 1). In the IR64 assembly the Sub1A gene is present as the Sub1A-2 allele, previously identified in submergence-intolerant accessions including IR64 [33]. For the Sub1B and C genes, IR64 carries the alleles Sub1B-1 and Sub1C-3, as reported [33]. Sub1A is absent from the DJ123 assembly, suggesting that this aus variety is not submergence tolerant. DJ123 carries a novel Sub1B allele (Sub1B-10), and the previously identified Sub1C-6 allele. In the Nipponbare assembly, Sub1B and Sub1C lie within a single scaffold, and the alleles identified are in agreement with published results [33]. Nipponbare is not submergence tolerant and the Sub1A gene is absent in Nipponbare according to previous reports. Our de novo assembly is unresolved in the region that corresponds to the Sub1A gene, but a k-mer analysis using the methods and data applied above clearly shows a lack of coverage in the DJ123 and Nipponbare sequencing reads across the locus except for high copy repeats dispersed in the sequence (Figure 4, top and bottom). Conversely, the k-mer coverage of the IR64 assembly is uniformly at the single-copy coverage level (approximately 100×), except for a small number of localized gaps in coverage, corresponding to SNPs in IR64 relative to the reference Sub1A sequence, and the high frequency repeats (Figure 4, middle). In contrast, the k-mer coverage across Sub1B and Sub1C is consistently approximately 100×, except for isolated sharp gaps corresponding to variations relative to the reference sequences (Figures S5 and S6 in Additional file 1).

Pup1 region

Phosphorus uptake1 (Pup1) is a major rice QTL associated with tolerance to phosphorus deficiency in soils [60],[61]. The Pup1 locus is a large, 90 kb region originally identified in Kasalath, an aus variety that is tolerant to phosphorus deficiency, but is absent in phosphorus starvation-intolerant varieties, including Nipponbare [62]. A gene encoding a protein kinase, Pstol1, located within the 90 kb indel, is responsible for the P-uptake efficiency phenotype [41].

Of the three de novo assemblies reported here, the 90 kb indel is absent from both Nipponbare and IR64, but a large portion of it, including the Pstol1 gene, is present in the aus variety DJ123 (Figure 5; Table S5 in Additional file 1). Although it is at least partially present, the region of the 90 kb indel described in Kasalath could not be fully resolved in our DJ123 assembly. This suggests that the 90 kb indel may be truncated and/or rearranged in some aus varieties. Interestingly, as shown in Figure 5, the Kasalath reference sequence contains unresolved gaps flanking regions of very high k-mer coverage; therefore, longer reads may be necessary to assemble this region with confidence. The Pstol1 gene sequence is complete in DJ123, and shows six SNPs relative to the Kasalath sequence (also apparent as abrupt drops in coverage in the k-mer coverage plot; Figures S9 and S10 in Additional file 1). These SNPs were confirmed via Sanger sequencing on genomic DNA (Figure S11 in Additional file 1). One of these SNPs introduces a premature stop codon, resulting in a protein that is only 136 amino acids long (the intact PSTOL1 protein is 324 amino acids) and, therefore, presumably non-functional.

Plant material

DNA sequencing

The DNA sequencing was performed in the Cold Spring Harbor Laboratory Genome Center using Illumina HiSeq 2000 instruments. For each of the three varieties, three libraries were sequenced following the requirements and recommendations of the ALLPATHS-LG whole genome assembler: (1) a 180 bp fragment library sequenced as 2 × 100 bp reads; (2) an approximately 2 kbp jumping library sequenced as 2 × 50 bp reads; and (3) an approximately 5 kbp jumping library sequenced as 2 × 50 bp reads.

For the 180-bp overlap library the sample was mechanically fragmented by using the Covaris S2 System and then prepared based on the New England Biolabs NEBNext Illumina library protocol and ligated to standard Illumina paired-end adapters. To maximize sample throughput the samples were size-selected in 50-bp windows between 290 and 310 bp using the Caliper XT instrument. Each library was PCR enriched for 12 cycles and quantified using the Bioanalyzer.

For the jumping libraries, the Illumina mate-pair library protocol was used. The DNA was fragmented into 2 kb and 5 kb segments. We again used the Covaris S2 System using programs that we developed in the lab. The fragmented DNA was then end-repaired with biotin-labeled dNTPs. The labeled fragments are circularized and fragmented again into 400 bp pieces. Fragments with the biotin labels are enriched, end-repaired, and ligated with adapters used for downstream processes. Each library was PCR enriched for 18 cycles and size-selected for 350 to 650 bp fragments. The final library consists of fragments made up of two DNA segments that were originally separated by approximately 2 kbp or approximately 5 kb. Each of the libraries was sequenced to 30× to 80× sequence coverage, as recommended by the assembler.

Libraries were sequenced on one or more lanes of an Illumina HiSeq 2000 using paired-end 50- or 100-bp runs. Image processing and base calling were performed as the runs progressed with Illumina’s Real Time Analysis (RTA) software. The binary base call files were streamed to a shared Linux server for further processing. The Illumina Casava pipeline (v1.8) was used to process the binary files to fastq files containing the base-called reads and per base quality scores. Only reads passing the standard Illumina quality filter were included in the output files.

Genome assembly

The ALLPATHS-LG version R41348 assembly algorithm was used for the assemblies. It consists of five major phases: (1) pre-assembly error correction, (2) merging of the overlapping fragment reads into extended reads, (3) constructing the unipath graph from the k-mers present in the reads, (4) scaffolding the unipaths with the jumping libraries, and (5) gap closing. To complete the five phases, the algorithm requires an overlapping pair fragment library and at least one jumping library, although the authors recommend at least two jumping libraries of approximately 2 kbp and approximately 5 kbp or larger. We assembled each of the genomes using approximately 50× coverage of the fragment library and approximately 30× coverage of each of the two jumping libraries using the recommended parameters, except we lowered the MIN_CONTIG size to 300 bp from the default 1,000 bp. This parameter controls the minimum contig size to be used for scaffolding, and our previous testing determined this change leads to (modestly) improved contig and scaffold statistics.

We also evaluated using SOAPdenovo2 [66] and SGA [67] for the assemblies (Table S3 in Additional file 1), using the same fragment, 2 kbp, and 5 kbp libraries but both assemblers had substantially worse contiguity statistics under a variety of parameter settings. For SOAPdenovo2, we corrected the reads using the Quake error correction algorithm [68], and then ran seven assemblies with the de Bruijn graph k-mer size set to k = 31 through k = 45 (odd values only, as required). In every attempt the scaffold N50 size was below 10 kbp compared with >200 kbp for our best ALLPATHS-LG assembly. For SGA, we evaluated four assemblies with the string graph minimum overlap length of k = 71 through k = 77 (odd values only, as required), but the scaffold N50 size was below 15 kbp in every attempt. We hypothesize that ALLPATHS-LG achieved superior results because the algorithm automatically measures many of the properties of the sequencing data, and could therefore self-adjust the various cutoffs used by the algorithm for error correction, contigging, and scaffolding.

Applying nomenclature proposed by [51], we have named these assemblies to convey accession, quality, origin, and iteration as follows: Os-Nipponbare-Draft-CSHL-1.0, Os-IR64-Draft-CSHL-1.0, Os-DJ123-Draft-CSHL-1.0.

Genome annotation

Repeat elements were masked using RepeatMasker [69] with a rice repeat library available from the Arizona Genome Institute. Protein-coding genes were annotated using MAKER-P version 2.30, installed on the Texas Advanced Computer Center Lonestar cluster and provisioned through an iPlant Collaborative allocation [52],[70]-[72]. Sequence evidence used as input for MAKER-P included Oryza expressed sequences (EST, cDNA, and mRNA) downloaded from the National Center for Biotechnology Information (NCBI), and annotated coding and protein sequences available for Nipponbare (IRGSP1.0 and MSU release 7) [51], 93-11 [28], and PA64s (Table S1 in Additional file 1). Ab initio gene predictions made using FGENESH [73] were incorporated exogenously into the MAKER-P pipeline using the pred_gff parameter. The SNAP [74] ab initio predictor was run within MAKER-P using the O.sativa.hmm parameter provided with SNAP. To annotate protein domain structure and assign GO terms we used InterProScan 5 software [75], available within the iPlant Discovery Environment [76]. Among resulting InterPro domains we curated 21 as being associated with transposon-encoded genes and screened out MAKER-P annotations with these domains (IPR000477, IPR001207, IPR001584, IPR002559, IPR004242, IPR004252, IPR004264, IPR004330, IPR004332, IPR005063, IPR005162, IPR006912, IPR007321, IPR013103, IPR013242, IPR014736, IPR015401, IPR018289, IPR026103, IPR026960, IPR027806). To identify homologies we conducted BLASTP alignment to the plants subsection of NCBI RefSeq (release 63), using an e-value threshold of 1e-10.

Whole genome comparisons

We used the MUMmer [77] whole genome alignment package and the GAGE assembly comparison scripts to compare the de novo assemblies to the reference Nipponbare and Indica genomes. Briefly, we aligned the assemblies to the genomes using nucmer using sensitive alignment settings (-c 65 -l 30 -banded -D 5). For base level accuracy evaluations, we used the GAGE assembly comparison script, which further refines the alignments by computing the best set of one-to-one alignments between the two genomes using the dynamic programming algorithm delta-filter. This algorithm weighs the length of the alignments and their percentage identity to select one-to-one non-redundant alignments. This effectively discards spurious repetitive alignments from consideration, allowing us to focus on the meaningful differences between the genomes. Finally, the evaluation algorithm uses dnadiff to scan the remaining, non-repetitive alignments to summarize the agreement between the sequences, including characterizing the nature of any non-aligning bases as substitutions, small indels, or other larger structural variations. To characterize the unaligned regions of the reference genome, we converted the whole genome alignments into BED format. For this we did not exclude repetitive alignments, so that we could focus on novel sequence instead of copy number differences. We used BEDTools [78] to intersect the unaligned segments with the reference annotation, and summarized the size distributions of the unaligned segments using AMOS [79].

K-mer analysis

To evaluate the repeat composition, we selected a random sample of 400 million unassembled reads from each of the three genomes and used Jellyfish [80] to count the number of occurrences of all length 21 k-mers in each read set. Length 21 was selected to be sufficiently long so that the expected number of occurrences of a random k-mer was below 1, but short enough to be robust to sequencing errors. The modes of the 3 k-mer frequency distributions, excluding erroneous k-mers that occurred less than 10 times, were 60× (Nipponbare), 64× (DJ123), and 73× (IR64) drawn from an approximately negative binomial distribution (Figure S1 in Additional file 1). These values correspond to the average k-mer coverage for single copy, non-repetitive regions of the genome. See Kelly et al. [68] for a discussion of k-mer frequencies. We then used the AMOS program kmer-cov-plot [79] to report the kmer coverage along the two reference genomes using the three databases of read k-mer frequencies. Unlike read alignments, which may be sensitive to repeats and variations, evaluating k-mer coverage is very robust to determine repetitive content [81],[82]. Single nucleotide variants are also readily apparent in these plots as abrupt gaps in coverage kilobase pairs long, while indels will be present as longer gaps in coverage [83].

Pan-genome analysis

The pan-genome analysis followed the reference-based analysis above, using nucmer to align the genomes to each other, BEDTools to find the genome-specific and shared regions of the genomes, and the jellyfish/AMOS k-mer analysis as described above to classify unique and repetitive sequences. We also used BEDTools to intersect the genome-specific/shared regions against their respective annotations to determine how the exonic bases were shared across the genomes. We summarized the genome-specific/shared exonic bases into gene counts by counting the total number of shared or specific exonic bases across all possible transcripts for a gene, and assigned the gene to the sector of the Venn diagram with the most bases associated with it. For the purposes of the Venn diagram (Figure 2A), wherever possible, the Nipponbare base or gene counts were used, followed by the values from IR64, and then followed by the DJ123 specific values, although the values were all largely consistent.

PCR and sequencing validation of specific regions

The same algorithms and parameters as the pan-genome analysis were also used to characterize the specific regions identified in the paper. PCR and/or sequencing validation were performed on genomic DNA extracted from tissue collected from independently grown plants obtained from the same seed source used for Illumina sequencing. Genomic DNA was extracted from young leaf tissue using the Qiagen Plant DNeasy Mini kit. Primers used for validation of 10 of the longest genome-specific sequences from each rice line, and of the S5 and Pup1 loci, are listed in Table S4a-c in Additional file 1. Sanger sequencing was performed at the Biotechnology Resource Center at Cornell University.

Data access