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

As the human population grows from 7.8 billion to 10 billion over the next 30 years, breeders must do everything possible to create crops that are highly productive and nutritious, while simultaneously having less of an environmental footprint. Rice will play a critical role in meeting this demand and thus, knowledge of the full repertoire of genetic diversity that exists in germplasm banks across the globe is required. To meet this demand, we describe the generation, validation and preliminary analyses of transposable element and long-range structural variation content of 12 near-gap-free reference genome sequences (RefSeqs) from representatives of 12 of 15 subpopulations of cultivated Asian rice. When combined with 4 existing RefSeqs, that represent the 3 remaining rice subpopulations and the largest admixed population, this collection of 16 Platinum Standard RefSeqs (PSRefSeq) can be used as a template to map resequencing data to detect virtually all standing natural variation that exists in the pan-genome of cultivated Asian rice.


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 (Table 1).
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 (Table 1) 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 18 . 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). www.nature.com/scientificdata www.nature.com/scientificdata/ 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 www.nature.com/scientificdata www.nature.com/scientificdata/ 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) ( Table 2). 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) ( Table 2).
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 19 . Quality control for each library data set was carried out with FastQC 20 . Finally, between 36.52-Gb and 51.05-Gb of clean data for each accession was generated, and used for genome size estimation (Table S1) by Kmer analysis (Fig. S3) and the Genome Characteristics Estimation (GCE) program 21 .
Bionano optical genome maps. Bionano optical maps for each accession were generated as previously described 22 , 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 23 . 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 24 , MECAT2 25 and Canu1.5 26 . 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 (Table S2). www.nature.com/scientificdata www.nature.com/scientificdata/ Step 2: Genome Puzzle Master (GPM) software 27 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) (Table S2).
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 28 . 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 29 was then used for mapping short Illumina reads onto assembled contigs, and the tool pilon 30 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 31 with a minimum alignment length of 1 Kb and an e-value < 1e −5 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) ( Table 3). 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) ( Table 3) Table S4). 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 (Table 3).
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) ( Table S5). 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 www.nature.com/scientificdata www.nature.com/scientificdata/ 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 15 .
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.
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 (Table 3).

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 36 , 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 37 , 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 ( Table 4). 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.  www.nature.com/scientificdata www.nature.com/scientificdata/ Structural Variants. Each genome assembly, as described above, was fragmented using the EMBOSS tool splitter 38 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 39 . Finally, the software SVIM 40 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  (Table S6). The total unshared fraction collected out of all pairwise comparisons was composed for 89.89% by TE related sequences.

Technical Validation
Assembly accuracy. Bionano optical maps were generated and used to validate all 12 genome assemblies.

code availability
The population re-analysis of 3K-RG dataset and 12 genome assemblies were obtained using several publicly available software packages. To allow researchers to precisely repeat any steps, the settings and the parameters used are provided below: Population structure: ADMIXTURE was run with default options. The R scripts for further population structure analysis, including setting up CLUMPP files, can be found in Github repository https://github.com/dchebotarov/Q-aggr.