A platinum standard pan-genome resource that represents the population structure of Asian rice (2024)

Ethics statement

This work was approved by the University of Arizona (UA), the King Abdullah University of Science and Technology (KAUST), Huazhong Agricultural University (HZAU), the International Rice Research Institute (IRRI) and the International Center for Tropical Agriculture (CIAT). All methods used in this study were carried out following approved guidelines.

Population structure

We extracted 30 subsets of 100,000 randomly chosen SNPs out of the 3K-RG Core SNP set v0.4 (996,009 SNPs, available at https://snp-seek.irri.org/_download.zul). For each subset, we ran ADMIXTURE¹⁶ with the number of ancestral groups K ranging from 5 to 15. We then aligned the resulting Q matrices using CLUMPP software¹⁷. Since different runs at a given value of K often give rise to different refinements (splits) of the lower level grouping, we first clustered the runs for each K according to similarity of Q matrices using hierarchical clustering, thus obtaining several clusters of runs for each K. We discarded one-element clusters (outlier runs), and averaged the Q matrices within each remaining cluster. FigureS1 shows the admixture proportions taken from the averaged Q matrices of the final clusters for K = 5 to 15. The columns of these averaged Q matrices, representing admixture proportions for groups discovered in different runs, were then used to define the “K15” grouping. At K = 9, 12, and 13, the Q matrices converged to two different modes according to whether XI-1A or GJ-trop is split (these are labeled as K = 9.1, 12.1 and 13.1).

Group membership for each sample was defined by applying a threshold of 0.65 to admixture components. Samples with no admixture components exceeding 0.65 were classified as follows. If the sum of components for subpopulations within the major groups cA (circum-Aus), XI (Xian-indica), and GJ (Geng-japonica) was ≥0.65, the samples were classified as cA-adm (admixed within cA), XI-adm (admixed within XI) or GJ-adm (admixed within GJ), respectively, and the remaining samples were deemed ‘fully’ admixed. The newly defined groups were mostly align with the previous K = 9 grouping, or were refined and named accordingly (e.g. XI-1B1 and XI-1B2 are two new subgroups within XI-1B).

The phenogram shown in Fig.1 was constructed with DARwin v6 (http://darwin.cirad.fr/, unweighted Neighbor-joining) using the identity by state (IBS) distance matrix from Plink on the 4.8 M Filtered SNP set (available at https://snp-seek.irri.org/_download.zul). Colors were assigned to subpopulations based on K15 Admixture results. One entry, MH 63 (XI-adm) represents the admixed types among the XI group.

Sample selection, collection and nucleic acid preparation

To select accessions to represent the 12 subpopulations of Asian rice that lack high-quality reference genome assemblies, the following strategy was employed. The IBS distance matrix was used for a principal component analysis (PCA) analysis in R to generate 5 component axes. Then, for each of the 12 subpopulations, i.e. circum-Aus2 = cA2, circum-Basmati = cB, Geng-japonica (GJ) subtropical (GJ-subtrp), tropical1 (GJ-trop1) and tropical2 (GJ-trop2), and Xian-indica (XI) subpopulations XI-1B1, XI-1B2, XI-2A, XI-2B, XI-3A, XI-3B1, XI-3B2, the centroid of each group in the space spanned by first 5 principal components was determined from the eigenvectors, and the entry closest to the centroid for which seed was available was chosen as the representative for that subpopulation (Table1).

Single seed decent (SSD) seed from IR 64 and Azucena were obtained from the Rice Genetics and Genomics Laboratory, CIAT, in Cali, Colombia, and SSD seed from the remaining 10 accessions (Table1) were obtained from the International Rice Genebank, maintained by IRRI, Los Baños, Philippines. All seed were sown in potting soil and grown under standard greenhouse conditions at UA, Tucson, USA for 6 weeks at which point they were dark treated for 48-hours to reduce starch accumulation. Approximately 20–50 grams of young leaf tissue was then harvested from each accession and immediately flash frozen in liquid nitrogen before being stored at −80 °C prior to DNA extraction. High molecular weight genomic DNA was isolated using a modified CTAB procedure as previously described¹⁸. The quality of each extraction was checked by pulsed-field electrophoresis (CHEF) on 1% agarose gels for size and restriction enzyme digestibility, and quantified by Qubit fluorometry (Thermo Fisher Scientific, Waltham, MA).

Library construction and sequencing

Genomic DNA from all 12 accessions were sequenced using the PacBio single-molecule real-time (SMRT) platform, and the Illumina platform for genome size estimations and sequence polishing. High molecular weight (HMW) DNA from each accession was gently sheared into large fragments (i.e. 30 Kb–60 Kb) using 26-gauge needles and then end-repaired according to manufacturer’s instructions (Pacific Biosciences). Briefly, using a SMRTbell Express Template Prep Kit, blunt hairpins and sequencing adaptors were ligated to HMW DNA fragments, and DNA sequencing polymerases were bound to the SMRTbell templates. Size selection of large fragments (above 15 Kb) was performed using a BluePippin electrophoresis system (Sage Science). The libraries were quantified using a Qubit Fluorometer (Invitrogen, USA) and the insert mode size was determined using an Agilent fragment analyzer system with sizes ranging between 30 Kb–40 Kb. The libraries then were sequenced using SMRT Cell 1 M chemistry version 3.0 on a PacBio Sequel instrument. The number of long-reads generated per accession ranged from 2.01 million (LIMA::IRGC 81487-1) to 5.40 million (Azucena). The distribution of subreads is shown in Fig.S2 and the average lengths ranged from 10.58 Kb (Azucena) to 20.61 Kb (LIMA::IRGC 81487-1) (Table2). According to the estimated genome size of the IRGSP RefSeq, the average PacBio sequence coverage for each accession varied from 103x (LIMA::IRGC 81487-1) to 149x (IR 64) (Table2).

For Illumina short-read sequencing, HMW DNA from each accession was sheared to between 250–1000 bp, followed by library construction targeting 350 bp inserts following standard Illumina protocols (San Diego, CA, USA). Each library was 2 × 150 bp paired-end sequenced using an Illumina X-ten platform. Low-quality bases and paired reads with Illumina adaptor sequences were removed using Trimmomatic¹⁹. Quality control for each library data set was carried out with FastQC²⁰. Finally, between 36.52-Gb and 51.05-Gb of clean data for each accession was generated, and used for genome size estimation (TableS1) by Kmer analysis (Fig.S3) and the Genome Characteristics Estimation (GCE) program²¹.

Bionano optical genome maps

Bionano optical maps for each accession were generated as previously described²², except that ultra-HMW DNA isolation, from approximately 4 g of flash-frozen dark-treated (48 hour) leaf tissue per accession, was performed according to a modified version of the protocol described by Luo and Wing²³. Prior to labeling, agarose plugs were digested with agarase and the starch and debris removed by short rounds of centrifugation at 13,000 × g. DNA samples were further purified and concentrated by drop dialysis against TE Buffer. Data processing, optical map assembly, hybrid scaffold construction and visualization were performed using the Bionano Solve (Version 3.4) and Bionano Access (v12.5.0) software packages (https://bionanogenomics.com/).

De novo genome assembly

Genome assembly for each of the 12 genomes followed a five-step procedure as shown in (Fig.2):

Step 1: PacBio subreads were assembled de novo into contigs using three genome assembly programs: FALCON²⁴, MECAT2²⁵ and Canu1.5²⁶. The number of de novo assembled contigs obtained varied from 51 (e.g. NATEL BORO::IRGC 34749-1 and KETAN NANGKA::IRGC 19961-2) to 1,473 (CHAO MEO::IRGC 80273-1) for the 12 genomes (TableS2).

Step 2: Genome Puzzle Master (GPM) software²⁷ was used to merge the de novo assembled contigs from the three assemblers, using the high-quality O. sativa vg. indica cv. Minghui 63 reference genome sequence (MH63RS2)^13,14 as a guide. GPM is a semi-automated pipeline that is used to integrate logical relationship data (i.e. contigs from three assemblers for each accession) based on a reference guide. Contigs were merged in the ‘assemblyRun’ step, with default parameters (minOverlapSeqToSeq was set at 1 Kb and identitySeqToSeq was set at 99%). Redundant overlapping sequences were also removed for each assembled contig. In addition, we gave contiguous contigs a higher priority than ones with gaps to be retained in each assembly. After manual checking, editing, and redundancy removal, the number of contigs in each assembly ranged from 26 (NATEL BORO::IRGC 34749-1) to 588 (LIU XU::IRGC 109232-1) (TableS2).

Step 3: The sequence quality of each contig was then improved by “sequence polishing”: twice with PacBio long reads and once with Illumina short reads. Briefly, PacBio subreads were aligned to GPM edited contigs using the software blasr²⁸. All default parameters were used, except minimum align length, which was set to 500 bp. Secondly, the tool arrow as implemented in SMRTlink6.0 (Pacific Biosciences of California, Inc) was used for polishing the GPM edited contigs. The bwa-mem program²⁹ was then used for mapping short Illumina reads onto assembled contigs, and the tool pilon³⁰ was used for a final polishing step with default settings.

Step 4: The polished contigs for each accession were arranged into pseudomolecules using GPM, with MH63RS2^13,14 as the reference guide. The program blastn³¹ with a minimum alignment length of 1 Kb and an e-value < 1e⁻⁵ as the threshold was used to align the corrected contigs to the reference guide. In doing so, the corrected contigs were assigned to chromosomes, as well as ordered and orientated in the GPM assembly viewer function. The number of contigs after step 4 ranged from a minimum of 15 contigs (GOBOL SAIL (BALAM)::IRGC 26624-2) to a maximum of 104 contigs (IR 64) (Table3). The assembly size for the 12 accessions ranged from 376.86 Mb (CHAO MEO::IRGC 80273-1) to 393.74 Mb (KHAO YAI GUANG::IRGC 65972-1) (Table3) and the length of individual chromosome varied from 23.06 Mb (chromosome 9 of CHAO MEO::IRGC 80273-1) to 44.96 Mb (chromosome 1 of LIMA::IRGC 81487-1) (TableS4). The average N50 value was 23.10 Mb, with the highest and the lowest N50 values being 30.91 Mb (LIU XU::IRGC 109232-1) and 7.35 Mb (IR 64), respectively. The average number of gaps among the 12 new genome assemblies was 18, with 8 assemblies containing less than 10 gaps (Table3).

Step 5: To independently validate our assemblies, we generated and compared Bionano optical maps to each assembly. In total, 17 (Azucena) to 56 (LIU XU::IRGC 109232-1) Bionano optical contigs were constructed for all 12 rice accessions, which yielded contig N50 values of between 22.75 Mb (CHAO MEO::IRGC 80273-1) to 31.45 Mb (KHAO YAI GUANG::IRGC 65972-1) (TableS5). As shown in Figs.3 and S4, the chromosomes and/or chromosome arms of all 12 de novo assemblies were highly supported by these ultra-long optical maps. Although rare, a few discrepancies between the optical maps and genome assemblies could be found and are likely due to small errors and chimeras that were produced through both the optical map and sequence assembly pipelines¹⁵.

Following these five steps, we were able to produce 12 near-gap-free Oryza sativa platinum standard reference genome sequences (PSRefSeqs) that represent 12 of 15 subpopulations of cultivated Asian rice.

BUSCO evaluation

The Benchmarking Universal Single-Copy Orthologs (BUSCO3.0) software package³² was employed to evaluate the gene space completeness of the 12 genome assemblies. These genomes captured, on average, 97.9% of the BUSCO reference gene set, with a minimum of 95.7% (IR64) and a maximum of 98.6% (LARHA MUGAD::IRGC 52339-1 and KHAO YAI GUANG::IRGC 65972-1) (Table3).

Of note, when performing this analysis, we observed that on average 30 out of the 1,440 conserved BUSCO genes tested (https://www.orthodb.org/v9/index.html) were missing from each new assembly, 16 of which were not present in all 12, plus the IRGSP RefSeq-1.0, ZS 97, MH 63 and N 22 RefSeqs (Fig.S5). This result suggested that these 16 “conserved” genes may not exist in rice, or other cereal genomes, thereby artificially reducing the BUSCO gene space scores for our 12 assemblies. To test this hypothesis, we searched for all 16 genes missing in maize, which diverged from rice about 50 million years ago (MYA)^33,34,35. We found that 13 of the 16 genes in question could not be found in 3 high-quality recently published maize genome assemblies (Fig.S5) and therefore, concluded that 13 of the 16 “conserved” genes in the BUSCO database are not present in cereals, and should be excluded from our gene space analysis. Taking this into account, we recalculated the BUSCO gene space content for each of 12 assemblies and found that 10 of 12 assemblies captured more than 98% of the BUSCO gene set (Table3).

Transposable element (TE) prediction

To determine the pan-transposable element content of cultivated Asian rice, we analyzed the 12 new reference genomes, presented here, along with the MH 63, ZS 97, N 22 PacBio reference genomes. In addition, we also included a reanalysis of the IRGSP RefSeq-1.0, as it is conventionally considered the standard rice genome for which all comparisons are conducted.

A search for sequences similar to TEs was carried out using RepeatMasker³⁶, run under default parameters with the exception of the option: -no_is –nolow, and that an updated in-house version of the publicly available MSU_6.9.5 library³⁷, retrieved from https://github.com/oushujun/EDTA/blob/master/database/Rice_MSU7.fasta.std6.9.5.out, called “rice 7.0.0.liban” was used. The average TE content of this 16 genome data set was 47.66% with a minimum value of 46.07% in IRGSP RefSeq-1.0 and a maximum of 48.27% in KHAO YAI GUANG::IRGC 65972-1 (Table4). The major contribution to this fraction was composed of long terminal repeat retrotransposons (LTR-RTs, min: 23.55%, max: 27.27% and average: 25.96%) followed by DNA-TEs (min:14.87%, max, 16.18% and average: 15.26%). Long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs) were identified as on average 1.43% and 0.39% of the 16 genomes, respectively.

Structural Variants

Each genome assembly, as described above, was fragmented using the EMBOSS tool splitter³⁸ to create a 10x genome equivalent redundant set of 50 kb reads. These reads were then mapped onto every other genome assembly using the tool NGMLR³⁹. Finally, the software SVIM⁴⁰ was run under default parameters to parse the mapping output. Only insertions, deletions and tandem duplications up to a maximum length of 25 Kb were considered in this analysis.

The results of this analysis identified several thousand insertions and deletions whenever an assembly was compared to any other. Greater variability was found between varieties belonging to different major groups (e.g. Geng-japonica [GJ] vs. Xian-indica [XI] than occurred between those within these groups. The amount of genome sequences with structural variation between any two varieties ranged from 17.57 Mb to 41.54 Mb for those belonging to the indica (XI) varietal group (avg: 31.75 Mb) and from 18.55 Mb to 23.07 Mb (avg: 21.00 Mb) for those in the japonica (GJ) varietal group. When all 16 genomes are considered together, the range is between 17.57 Mb and 41.54 Mb, with an average value of 33.70 Mb (TableS6). The total unshared fraction collected out of all pairwise comparisons was composed for 89.89% by TE related sequences.