An improved haplotype resolved genome reveals more rice genes

Nipponbare is a japonica rice cultivar that has been widely used as the standard reference genotype for rice^[1]. The rice (Nipponbare) genome was one of the first crop genomes to be sequenced more than 20 years ago^[2]. The 1^st sequence of the rice genome was completed in 2002 and was a major milestone in the field of plant genomics by the International Rice Genome Sequencing Project, 2005^[3]. These international collaborative efforts provided the first genome of a crop plant. The Nipponbare genome assembly contained gaps, primarily due to repetitive DNA sequences. In 2005, these gaps were estimated to be approximately 18.1 Mb in total, with the majority originating from centromeres and telomere regions. Sequencing technological advancements and ongoing research efforts, have improved the rice genome sequence over time^[4,5]. Thorough efforts were made to improve the quality of the Nipponbare reference genome assembly in 2013, resulting in greatly enhanced accuracy of cDNA sequences and gene annotation, while it remained incomplete^[5]. In the human genome, recent significant strides have been made in assembling and characterization the previously unexplored 8% of the human genome, especially including telomere sequences^[6].

Reference genome assemblies often contain gaps, especially regions with repetitive sequences, termed the 'dark side' of the genome^[7,8]. New sequence technology allows improved assembly quality, with less gaps, leading to a more complete and accurate representation of the genome. The achievement of a higher quality and more complete reference genome will provide new insights into genomics and breeding, supporting pan-genome studies and genome wide association studies^[9]. Recently many other Oryza genomes have been sequenced and assembled, including indica and wild rice species^[9−11] Most recently the Nipponbare genome sequence gaps and telomere sequence were addressed^[12,13]. Despite these advancements, a fully haplotype resolved assembly has not been reported. In this study, we have used PacBio HiFi reads to produce a more accurate genome sequence assembly. The novel genome assembly is not only almost 11.3 Mb longer than the IRGSP-1.0 reference but also exhibits improvements in all chromosomes (Fig. 1), including the addition of telomeric regions in all chromosomes (T2T), with the addition of fully resolved haplotypes (haplotype 1 and haplotype 2, telomere-to-telomere-T2T) (Tables 1, 2 & Supplemental Tables S1−S10).

Figure 1. 

Comparison of the UQ Nipponbare genome with previously published reference genome assembly IRGSP.1.0 Nipponbare. (a), (b) For UQ Nipponbare all the chromosome sizes are larger and most include telomeres as compared to IRGSP.1.0 Nipponbare. (c) Whole genome dot plot of UQ Nipponbare genome vs IRGSP.1.0 Nipponbare.

Comparative analysis of annotations of the new genome (UQ Nipponbare) and the IRGSP-1.0 reference revealed the presence of 3,050 additional genes, for which more than 95% had supporting transcript evidence (Supplemental Fig. S1).

These findings underscore the potential of new sequencing technologies to significantly augment reference genomes, potentially leading to more comprehensive genetic information. These results also suggest that applying advanced sequencing technologies to other established genomes may yield similar benefits, potentially enhancing our knowledge of these species. This study highlights the continuous evolution of genomics and underscore the importance of staying at the forefront of sequencing technologies for the accurate representation of complex genomes.

PacBio HiFi reads and Hi-C reads were used to generate a contig assembly with Hifiasm^[14] producing a haplotype phased assembly. The contig level assembly produced single contigs for nine chromosomes, while the remaining three chromosomes were each covered by two contigs each. Hi-C data were employed to hierarchically cluster the assembled contigs into 12 pseudo-chromosomes, by using the YaHS scaffolding tool^[15]. The T2T assembly had a single scaffold for each of the 12 pseudo-chromosomes and was larger in size than the corresponding IRGSP.10 genome (Fig. 1). The results of BUSCO analysis showed that the collapsed assembly covered 99.3% of the universal single copy genes with an N50 of 30.7 Mb (Supplemental Tables S1−S3). The UQ Nipponbare collapsed genome assembly is larger in size compared to the IRGSP.1.0 Nipponbare reference genome assembly. In Fig. 2, additional non-aligning regions of each chromosome in the UQ Nipponbare collapsed assembly are highlighted, along with the structural variants in the comparison of the previously published IRGSP1.0 Nipponbare reference genome assembly.

Figure 2. 

Sequence collinearity and structural variants, including inversions, translocations, duplications, and non-aligning regions, were analysed between the UQ Nipponbare genome assembly and the IRGSP-1.0-Nipponbare genome assembly. The two assemblies were aligned using minimap2, and the resulting BAM file was indexed with samtools. Detection of structural variations between these two genomes was performed using SyRI^[26−28]. The non-aligning regions of chromosomes 1 and 11 are highlighted in the bottom section of the figure.

For the two phased haplotypes (T2T) of UQ Nipponbare; haplotype 1 covered 98.9% of the single copy orthologs with an N50 of 30.6 Mb, whilst haplotype 2 covered 94.9% single copy orthologs with an N50 of 29.3 Mb (Supplemental Tables S4−S9 & Supplemental Figs S2 & S3). The haplotype 1 chromosomes were larger than the haplotype 2 chromosomes (Fig. 3). This first haplotype resolved Nipponbare genome incorporated 3,050 new genes compared to IRGSP-1.0, and is expected to be a valuable and significant resource for rice researchers, and these additional genes had a wide range of functions (Supplemental Table S11). Of these additional genes, 58 genes fell in new regions that were missing in the IRGSP genome, but most genes were in regions that were not new due to the improved accuracy of sequencing.

Figure 3. 

Sequence collinearity and structural variants, including inversions, translocations, duplications, and non-aligning regions, were examined between the UQ Nipponbare haplotype 1 genome assembly and the UQ Nipponbare haplotype 2 genome assembly. The same method used in Fig. 2 was used for this analysis.

The CTAB method^[19] was used to extract DNA from young leaves of a rice (Oryza sativa cv Nipponbare) plant grown in a glasshouse at the University of Queensland (Australia). The high quality DNA extracted was sequenced using a PacBio (Pacific Biosciences) Sequel II platform to produce HiFi sequences.

Genome assembly

Approximately 54.9 Gb of HiFi reads were obtained. HiC reads (59.6 Gb) were downloaded from the NCBI Sequence Read Archive database (SRR6470741). De novo haplotype-resolved assembly of these reads was performed using hifiasm with parameters '--write-ec --write-paf -l0'^[14]. The contig assemblies were scaffolded using the YaHS tool^[15]. QUAST and BUSCO were used to evaluate the quality and completeness of a genome assembly^[16,17]. A telomere identification toolkit (tidk) was used to search for tandem repeats of the telomeric sequence 'TTAGGG' and 'TAAACCC' and the exact location in the Nipponbare collapse, haplotype-1 and haplotype-2 assembly (https://github.com/tolkit/telomeric-identifier).

Genome annotation

Repetitive DNA sequences were obtained from a Oryza repeat database^[18] and used to mask the genome with the Repeatmasker soft masking option^[19]. Protein sequences of Viridiplantae from OrthoDB v.11^[20] and RNA-sequencing (RNA-seq) reads from the NCBI Sequence Read Archive database (SRR23560402, SRR23560417, SRR23560416, SRR23560419, SRR23560418, SRR23560409, SRR23107175, SRR23107177, SRR23107178, SRR8051554, SRR7974062, SRR8051550) were obtained. Quality and adapter trimmed RNA-seq reads were aligned to the masked genomes using HISAT2^[21]. Annotation of protein-coding genes in Nipponbare was conducted using a combination of homology-based prediction, de novo prediction, and transcriptome-based prediction methods using Braker^[22]. BUSCO was used to assess the genome annotation completeness. The Large Gap Mapping tool (length fraction; 0.9 similarity fraction; 0.9) of CLC was used to identify the new genes with the comparison of IRGSP-1.0 Nipponbare genes (CLC Genomics Workbench 23.0.05, QIAGEN, USA, www.clcbio.com) and further transcript evidence for these new genes was estimated (Supplemental Table S12). The functional annotation of the identified additional genes was performed using OmicsBox 2.2.4^[23]. CDS sequences were subjected to a BLASTX analysis with a specific e-value of 1.0E-10 against the non-redundant protein sequences database, utilizing Viridiplantae taxonomy. Subsequently, the CDS sequences were processed through InterProScan, and GO terms were extracted for all matches acquired via the BLAST search, employing Gene Ontology mapping with the Blast2GO annotation tool. The annotations generated from InterProScan and Blast2GO were then harmonized by merging the respective GO terms (Supplemental Table S11).

The authors confirm contribution to the paper as follows: experiment interpretation: Abdullah M; project supervision: Furtado A, Masouleh AK, Henry RJ; data analysis: Abdullah M, Furtado A, Masouleh AK; draft manuscript preparation and data interpretation: Okemo P; study conception: Henry RJ. All authors approved the final version of the manuscript.

The whole genome sequence data reported in this paper have been deposited in the Genome Warehouse in National Genomics Data Center^[24,25], Beijing Institute of Genomics, Chinese Academy of Sciences / China National Center for Bioinformation, under accession number GWHEQBP00000000 that is publicly accessible at https://ngdc.cncb.ac.cn/gwh.