Maize (Zea mays ssp. mays L., NCBI:txid381124) is one of the most important crops grown worldwide for food, forage, and biofuel, with an annual production of more than 1 billion tons [1]. Owing to the rapid human population growth and economic demand, maize has been predicted to account for 45% of the total cereal demand by 2050 [2]. In addition, it is an important model organism for fundamental research in genetics and genomics [3].

Because of its importance in crop science, genetics and genomics, several reference genomes of common maize inbred lines used in breeding have been released since 2009 [4–8]. However, comparative genomic analyses have found that maize genomes exhibit high levels of genetic diversity among different inbred lines [1, 7, 9]. Meanwhile, accumulating studies have suggested that one or a few reference genomes cannot fully represent the genetic diversity of a species [7, 10, 11].

The maize cultivar Dan340 is an excellent backbone inbred line of the LvDa Red Cob Group that has several desirable characteristics, such as disease resistance, lodging resistance, high combining ability, and wide adaptability. More than 50 maize hybrid breeds have been derived from Dan340 since 2000, and their planting area has reached 19 million ha. It is considered that Dan340 originated from a landrace in China and exhibits significant genetic differences from other maize germplasms that represent the most important core maize germplasms in China [12]. Therefore, Dan340 could serve as a model inbred line for the genetic dissection of desirable agronomic traits, combining ability, heterosis, and breeding history.

In the present study, we constructed a high-quality chromosome-level reference genome for Dan340 by combining PacBio long HiFi sequencing reads, Illumina short reads, and chromosomal conformational capture (Hi-C) sequencing reads. The completeness and continuity of the resulting genome are comparable with those of other important maize inbred lines: B73 [4], Mo17 [7], SK [13], PH207 [5], and HZS [8]. Furthermore, comparative genomic analyses were performed between Dan340 and other maize lines. Genes and gene families specific to Dan340 were identified. In addition, large numbers of structural variations between Dan340 and other maize inbred lines were detected. The assembly and annotation of this genome will increase our understanding of the intraspecific genomic diversity in maize and provide a novel resource for maize breeding improvements.

The inbred line Dan340 (Figure 1) was selected for genome sequencing and assembly because it is an elite maize cultivar that plays an important role in maize breeding and genetic research. The plants were grown at 25 °C in a greenhouse of the Beijing Academy of Agriculture and Forestry Sciences, Beijing, China. Fresh and tender leaves were harvested from the best-growing individual, immediately frozen in liquid nitrogen, and then preserved at −80 °C in the laboratory prior to DNA extraction. Genomic DNA was extracted from the leaf tissue of a single plant using the DNAsecure Plant Kit (Tiangen Biotech Co., Ltd., Beijing, China). To ensure that the DNA extracts were useable for all types of genomic libraries, their quality and quantity were evaluated using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and electrophoresis on a 0.8% agarose gel, respectively.

In recent years, third-generation DNA sequencing technologies have undergone rapid technological innovation and are now widely used in genome assembly. In this study, PacBio circular consensus sequencing (CCS) libraries were prepared using the SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, CA, USA; Ref. No. 101-685-400), following the manufacturer’s protocols, and they were subsequently sequenced on the PacBio sequel II platform (Pacific Biosciences, RRID:SCR_017990). As a result, 63.53 Gb (approximately 27× coverage) of HiFi reads was generated and used for the genome assembly.

In addition, one Illumina paired-end sequencing library, with an insert size of 350 bp, was generated using the NEB Next Ultra DNA Library Prep Kit (NEB, Ipswich, MA, USA) following the manufacturer’s protocol and then sequenced using an Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA, RRID:SCR_016385) at the Novogene Bioinformatics Institute, Beijing, China. Approximately 80.66 Gb (∼34×) of Illumina sequencing data were obtained.

One Hi-C library was constructed using young leaves following previously published procedures [14], with slight modifications outlined in our published protocol [15] (Figure 2). In brief, approximately 5-g leaf samples from seedlings were cut into minute pieces and cross-linked using a 4% formaldehyde solution at room temperature in a vacuum for 30 min. Each sample was mixed with excess 2.5 M glycine for 5 min to quench the cross-linking reaction and then placed on ice for 15 min. The cross-linked DNA was extracted and then digested for 12 h with 20 units of DpnII restriction enzyme (NEB, Ipswich, MA, USA, Catalog #R0543S) at 37 °C. Next, the resuspended mixture was incubated at 62 °C for 20 min to inactivate the restriction enzyme. The sticky ends of the digested fragments were biotinylated and proximity ligated to form enriched ligation junctions and then ultrasonically sheared to a size of 200–600 bp. The biotin-labelled DNA fragments were pulled down and ligated with Illumina paired-end adapters, and then amplified by PCR to produce the Hi-C sequencing library. The library was sequenced using an Illumina HiSeq X Ten platform with 2 × 150 bp paired-end reads. After removing low-quality sequences and trimming adapter sequences, 304.37 Gb (approximately 130×) of clean data were generated and used for the genome assembly.

To obtain a high-quality genome assembly of Dan340, we employed both PacBio HiFi reads and Illumina short reads, with scaffolding informed by high-throughput Hi-C. The assembly was performed in a stepwise fashion. First, a de novo assembly of the long CCS reads generated from PacBio single-molecule real-time (SMRT) sequencing was performed using Hifiasm [16] (RRID:SCR_021069). A total of two SMRT cells produced 4,073,418 subreads, with an average length of 15,598 bp and a read N50 of 15,715 bp. The generation of HiFi reads and adapter trimming was performed using PacBio SMRTLink (Version 8.0) [17] with default parameters, followed by the deduplication of reads using pbmarkdup (Version 0.2.0) [18], as recommended by PacBio. Next, HiFi reads were aligned to each other and assembled into genomic contigs using Hifiasm [16] with default parameters. Next, the primary contigs (p-contigs) were polished using Quiver [19] by aligning the SMRT reads. Then, Pilon [20] (RRID:SCR_014731) was used to perform the second round of error correction using the short paired-end reads generated by the Illumina Hiseq platforms. Subsequently, the Purge Haplotigs pipeline [21] was used to remove redundant sequences formed due to heterozygosity. The draft genome assembly was 2348.68 Mb; it reached a high level of continuity and a contig N50 length of 45.11 Mb.

To reduce Hi-C reads having a bias due to experimental artefacts, we removed the following read types using HiCUP [22] (RRID:SCR_005569): (a) reads with ≥10% unidentified nucleotides (N); (b) reads with >10 nt aligned to the adapter, allowing ≤10% mismatches; and (c) reads with >50% bases having a Phred quality <5. Next, the filtered Hi-C reads were aligned against the contig assemblies using BWA (Version 0.7.8, RRID:SCR_010910) [23]. Reads were excluded from subsequent analyses if they did not align within 500 bp of a restriction site or did not uniquely map. Also, the number of Hi-C read pairs linking each scaffold pair was tabulated. ALLHiC (Version 0.8.12) [24] was used in simple diploid mode to scaffold the genome and optimize the ordering and orientation of each clustered group, producing a chromosome-level assembly. The Juicebox Assembly Tools (Version 1.9.8, RRID:SCR_021172) [25] were used to visualize and manually correct the large-scale inversions and translocations to obtain the final pseudo-chromosomes (Figure 3). Finally, 2315 scaffolds (representing 91.30% of the total length) were anchored to 10 chromosomes (Figure 4).

The final assembly of the Dan340 genome was 2348.72 Mb, including 2738 contigs and 2315 scaffolds, with N50 of 41.49 Mb and 215.35 Mb, respectively (Table 1).

We assessed the quality of the assembly using several independent methods. First, the short reads obtained from the Illumina sequencing data were aligned to the final assembly using BWA [26]. Our results showed that the percent of reads mapped to the reference genome was 97.48%. Second, a total of 248 conservative genes existing in six eukaryotic model organisms were selected to form the core gene library for the Core Eukaryotic Genes Mapping Approach (CEGMA) [27] (RRID:SCR_015055) evaluation. To evaluate its integrity, our assembled Dan340 genome was aligned to this core gene library using TBLASTN (RRID:SCR_011822) [28], GeneWise (Version 2.2.0, RRID:SCR_015054) [29], and the GeneID tools (Version 1.4 RRID:SCR_021639) [30]. Our results showed that 238 complete (95.97%) and 243 partial (97.98%) genes were detected in our assembly. Third, the completeness was assessed using the benchmarking universal single-copy orthologs (BUSCO) [31] (RRID:SCR_015008). The final assembly was tested against BUSCO (v.3) with embryophyta_odb10 database [32], which includes 1614 conserved core genes. Our results showed that 98.08% (1583), 1.11% (18), and 0.81% (13) of the plant single-copy orthologs were present in the assembled Dan340 genome as complete, fragmented, or missing genes, respectively. Fourth, the long-terminal repeat (LTR) Assembly Index (LAI) metric was used to evaluate the assembly continuity in Dan340 and three other maize genomes (B73, Mo17, and SK; Figure 5). The intact LTR retrotransposons were identified in the four genomes using LTRharvest (Version 1.6.1, RRID:SCR_018970) [33], LTR_Finder (Version 1.07, RRID:SCR_015247) [34], and LTR_retriever (Version 2.9.0, RRID:SCR_017623) [35]. The LAI pipeline was executed using the following parameter settings: -t 20 -intact genome.fasta.pass.list -all genome.ltr.fasta.out. Our Dan340 genome had an LAI score of 25.13, which was relatively high among the four maize genomes compared in this study. B73, Mo17, and SK scored 24.94, 24.45, and 27.12, respectively (Figure 5 and Table 2). A higher LAI score indicates a more complete genome assembly because more intact LTR retrotransposons are identified, as was the case of our Dan340 genome. Furthermore, whole-genome sequence alignments of Dan340 to the genomes of the other three maize inbred lines demonstrated that our assembly has highly collinear relationships with other published maize genomes (Figure 6). Taken together, our assessment results suggest that the Dan340 genome assembly is of high quality.

Repeat sequences of the Dan340 genome were annotated using both ab initio and homolog-based search methods. For the ab initio prediction, RepeatModeler (Version 1.0.8, RRID:SCR_015027) [36], RepeatScout (Version 1.0.5, RRID:SCR_014653) [37], and LTR_Finder [34] were used to discover transposable elements (TEs) and to build a TEs library. An integrated TEs library and a known repeat library (Repbase Version 15.02, homolog-based, RRID:SCR_021169) were subjected to RepeatMasker (Version 3.3.0 RRID:SCR_012954) [38] to predict the TEs. For the homolog-based predictions, RepeatProteinMask was performed to detect the TEs in our genome by comparing it against a TE protein database. Tandem repeats were ascertained in the genome using Tandem Repeats Finder (Version 4.07b, RRID:SCR_022193) [39]. As a result, 1723.99 Mb of repeat sequences were identified, accounting for 73.40% of the genome size. Among these repeat sequences, 1555.57 Mb were predicted to be long-terminal repeat (LTR) retrotransposons, and 44.53 Mb were predicted to be DNA transposons, accounting for 66.23% and 1.60% of the genome, respectively. Furthermore, among the LTR retrotransposons, the Gypsy and Copia superfamilies comprised 23.81% and 12.75% of the genome, respectively. Thus, retrotransposons accounted for a large proportion of the Dan340 genome, which was consistent with the genomic characteristics of other maize inbred lines (Table 2).

All repetitive regions except the tandem repeats were soft-masked for protein-coding gene annotations. Five ab initio gene prediction programs, Augustus (Version 3.0.2, RRID:SCR_008417) [40–42], GENSCAN (Version 1.0, RRID:SCR_013362) [43], GeneID [30], GlimmerHMM (Version 3.0.2, RRID:SCR_002654) [44], and SNAP (Version 2013-02-16, RRID:SCR_007936) [45], were used to predict genes. In addition, the protein sequences of five homologous species (Sorghum bicolor, Setaria italica, Hordeum vulgare, Triticum aestivum, and Oryza sativa) were downloaded from Ensembl and NCBI. Homologous sequences were aligned against the genome using TBLASTN (E-value 1 × 10−5). GeneWise [29] was employed to predict gene models based on the sequence alignment results.

For the RNA-seq predictions, fresh samples of six tissues (stem, endosperm, embryo, bract, silk, and ear tip) were collected. The total RNA was extracted from each sample using an RNAprep Pure Plant Kit (Tiangen Biotech Co., Ltd., Beijing, China). The isolated, purified RNA, having fragment lengths of approximately 300 bp, was the template for constructing a cDNA library. The NEBNext Ultra RNA Library Prep Kit from Illumina (New England Biolabs, Ipswich, MA, USA) was used to construct the cDNA library following the manufacturer’s instructions. The sequencing was performed on an Illumina HiSeq X Ten platform, and 150-bp paired-end reads were generated. Raw reads were trimmed by removing the adapter sequences, reads with more than 5% of unknown base calls (N), and low-quality bases (base quality less than 5). Clean paired-end reads were aligned to the genome using TopHat (Version 2.0.13, RRID:SCR_013035) [46] to identify exon regions and splice positions. The alignment results were then used as input for Cufflinks (Version 2.1.1, RRID:SCR_014597) [47] to assemble the transcripts into the gene models. In addition, RNA-seq data were assembled using Trinity (Version 2.1.1, RRID:SCR_013048) [48], creating several pseudo-ESTs (short for expressed sequence tags). These pseudo-ESTs were also mapped to the assembled genome using BLAT [49] (RRID:SCR_011919), and gene models were predicted using PASA [50] (RRID:SCR_014656). A weighted and non-redundant gene set was generated using EVidenceModeler (EVM, Version 1.1.1, RRID:SCR_014659) [51], which merged all the gene models predicted by the above three approaches. Finally, PASA was used to adjust the gene models generated by EVM. As a result, 39,733 protein-coding genes were annotated in our final set. To better understand gene functions, we used all our 39,733 protein-coding genes as queries against public protein databases, including NCBI non-redundant protein sequences, Swiss-Prot, Protein family, Kyoto Encyclopedia of Genes and Genomes (KEGG), InterPro, and Gene Ontology (GO). In total, 39,646 genes (99.8%) were annotated using these databases, and 24,402 (61.41%) were supported by RNA-seq data. Furthermore, the number of genes, the gene length distribution, and the exon length distribution were all comparable to those of other maize inbred lines and common crop species (Table 3).

Transfer RNA (tRNA) genes were predicted using tRNAscan-SE (Version 1.4, RRID:SCR_010835) [52] with the default parameters. Ribosomal RNAs (rRNAs) were annotated based on their homology levels with the rRNAs of several species of higher plants using BLASTN with an E-value of 1 × 10−5. The microRNA (miRNA) and small nuclear RNA (snRNA) fragments were identified by searching the Rfam database (Version 11.0, RRID:SCR_007891) using Infernal (Version 1.1, RRID:SCR_011809) [53, 54]. Finally, 4547 miRNAs, 5963 tRNAs, 63,564 rRNAs, and 1422 snRNAs were identified, with average lengths of 126.79, 75.25, 309.47, and 132.10 bp, respectively (Table 4).

We applied the OrthoMCL pipeline [55] to identify orthologous gene families among the four maize inbred lines, including Dan340, B73, Mo17, and SK. The longest protein from each gene was selected, and the proteins with a length of less than 30 amino acids were removed. Subsequently, pairwise sequence similarities between all input protein sequences were calculated using BLASTP [56] (RRID:SCR_001010) with an E value cut-off of 1 × 10−5. Markov clustering (MCL) of the resulting similarity matrix was used to define the ortholog cluster structure of the proteins, using an inflation value (-I) of 1.5 (default setting of OrthoMCL).

Next, comparative analyses were performed among Dan340, B73, Mo17, and SK (Figure 7A). The genes from the Dan340 genome and those from B73, Mo17 and SK were clustered into 27,654 gene families. Of these, 15,690 families were shared among the four maize inbred lines, representing a core set of genes across these maize genomes. We found 1806 genes from 359 gene families that were specific to Dan340, of which many had functional GO annotations related to “protein phosphorylation”, “single-organism catabolic process”, and “pheromone binding” (Figure 7B). Using the KEGG functional enrichment, the most enriched pathways of the Dan340-specific genes were “antifolate resistance”, “epithelial cell signaling in Helicobacter pylori infection”, and “pentose and glucuronate interconversions” (Figure 7C).

In addition, OrthoMCL was used to identify the core and dispensable gene sets based on gene families. The gene families that were shared among the four inbred lines were defined as core gene families. Furthermore, gene families shared among three inbred lines, between two inbred lines, and those only present in one inbred line (private gene families) are also displayed in Figure 7D.

To investigate the genetic and structural variations between Dan340 and other maize inbred lines, we first aligned the other three genomes to the Dan340 reference genome based on MUMmer (V 4.0.0 beta2) [57] (RRID:SCR_018171) with parameter “—mum -g 1000 -c 90 -l 40”. Then the alignment files were filtered to generate 1-to-1 mapping by delta-filter with parameters:“-m -i 90 -l 100”. The genomes of B73 and Mo17 were downloaded from MaizeGDB [58], and the genome of SK was obtained from the National Genomics Data Center [59]. Next, the output of Nucmer was analyzed using SyRI [60] with default parameters to identify variation. On the basis of the above pipeline, we obtained structural variation sets and generated into the vcf file. We also used PBSV (Version 2.2.2) [61] to investigate the genetic and structural variations (details and outputs available via the GigaDB entry [62]).

The high-quality Dan340 reference genome allowed us to identify large SVs in different maize inbred lines. By aligning the genome of B73 to the Dan340 genome, we identified 36,363 structural variations (longer than 500 bp) between the two representative maize genomes, including 15,923 insertions, 16,173 deletions, 141 inversions, and 4,126 duplications (Table 5). Furthermore, the structural variations presented in Mo17 and SK were also detected in this study (Tables 6 and 7). The dataset generated by PBSV is available in GigaDB [62]. These datasets provide abundant variation resources for future molecular improvements and breeding in maize.