High-quality genome assembly
We assembled a highly continuous and complete genome sequence of Heinz 1706 using an integrated genome sequencing approach that combined 131.78 Gb (168.52×) of SMRT data, 226.97 Gb (290.24×) of BioNano data, 140.52 Gb (179.70×) of Hi-C data, and 50.93 Gb (61.53×) of Illumina short-read data (Supplementary Table S1). The PacBio long reads with an N50 read length of 32.82 kb were assembled with CANU software[20], generating a 875.21-Mb genome with a contig N50 of 17.83 Mb (Table 1). To reduce fragmentation and fill in gaps, BioNano data and Hi-C data were used to assist with scaffold construction using Aigner and Assembler[21], HERA[22], and Juicer[23] software. A Hi-C-based physical heatmap comprising 12 groups was generated (Supplementary Fig. S1) and used to create 12 pseudo-chromosomes that anchor ~790.59 Mb of the genome and harbor 97.61% (33,562) of the predicted protein-coding genes. The genome assembly was polished with Illumina short reads for error homozygous SNPs or indels using Pilon software[24]. As a result, we generated a 799.09-Mb genome assembly, SLT1.0 (Fig. 1, Table 1).
Table 1
Genome assembly and annotation of SLT1.0
|
|
SLT1.0
|
SL4.0
|
SL3.0
|
Genome assembly (Mb)
|
799.09
|
782.52
|
828.08
|
Non-N bases
|
797,955,212
|
782,475,302
|
746,357,470
|
Number of gaps
|
210
|
286
|
22,700
|
Number of total contigs
|
1,615
|
504
|
-
|
Longest contig length (Mb)
|
47.16
|
26.29
|
-
|
N50 of contigs (Mb)
|
17.83
|
6.01
|
-
|
Number of unplaced contigs
|
112
|
176
|
4,374
|
Unplaced contigs sequence length (Mb)
|
8.50
|
9.64
|
20.85
|
Number of genes
|
34,384
|
34,075
|
35,768
|
Percentage of gene length in genome (%)
|
16.21
|
15.56
|
17.33
|
Mean gene length (bp)
|
3,766.53
|
3,572.44
|
4,011.09
|
Gene density (per Mb)
|
43.03
|
43.55
|
43.19
|
Mean coding sequence length (bp)
|
223.02
|
228.01
|
219.97
|
Mean exon length (bp)
|
310.11
|
275.03
|
308.36
|
Mean intron length (bp)
|
270.41
|
606.69
|
632.38
|
Masked repeat sequence length (Mb)
|
558.49
|
546.95
|
507.14
|
Repeats percentage of genome size (%)
|
69.89
|
69.90
|
61.24
|
The conserved genes from the Benchmarking Universal Single-Copy Orthologs (BUSCO) gene set[25] were used to gauge the accuracy and completeness of the SLT1.0 assembly. The results showed that the SLT1.0 assembly contained 97.70% complete genes and 0.30% fragmented genes. The value of the LTR Assembly Index (LAI) was 12.41, which was consistent with that of the previously released SL4.0 tomato reference genome (LAI 12.54). More than 99.88% of the genome assembly had greater than one-fold coverage with Illumina short reads. All these evidences demonstrated the high continuity and completeness of the SLT1.0 genome assembly.
High-quality genome annotation
Except for ab initio prediction and protein-homology-based prediction, we also used transcriptome data, including the bulked RNA-seq data with a mapping rate of 99.73%, and previously-released RNA-seq data from various tissues[3] with a mapping rate of 97.97%, to facilitate gene annotation of the assembled genome. In total, we predicted 34,384 protein-coding genes with an average length of 3,766.53 bp and 6.55 exons per gene in the SLT1.0 genome (Table 1 and Supplementary Table S2). Gene completeness was estimated to be 98.20% based on the BUSCO gene set[25], and the protein-coding genes were unevenly distributed along the chromosomes (Fig. 1). Comparative analysis showed that 234 genes in the SLT1.0 genome corresponded to 488 genes in the SL4.0 genome (Supplementary Table S3). Gene collinearity analysis identified 33 collinear gene blocks between the SLT1.0 and SL4.0 genomes, harboring 28,892 (84.03%) and 28,389 (83.30%) homologous genes, respectively (Fig. 2A). Some unplaced contigs in the SL4.0 genome were successfully assigned to chromosomes in the SLT1.0 genome. These results highlight the high accuracy and completeness of the SLT1.0 genome assembly and gene models.
A comprehensive analysis of the genome sequences identified 965 collinear chromosomal blocks between the SLT1.0 and SL4.0 genomes. These blocks contained 32,922 and 32,554 genes, accounting for 95.75% and 95.54% of the SLT1.0 and SL4.0 genomes, respectively. However, we detected a 2.76 Mb inversion from 39.17 to 41.93 Mb on chromosome 2 of the SLT1.0 genome (Fig. 2B). The continuous interaction signals on the Hi-C heatmap, as well as PCR and Sanger sequencing, showed that this region was not misassembled (Fig. 2B-C, Supplementary Table S4). This result indicated that heterozygous variation may exist in the previously reported Heinz 1706 accession.
Transposable element analysis
A total of 524.84 Mb of repetitive sequences were identified, accounting for 65.66% of the SLT1.0 genome assembly, which was similar to that reported in the SL4.0 genome (508.89 Mb, 65.03%) (Supplementary Table S5). Among these repetitive sequences, long terminal repeats (LTRs) were the predominant TE family, covering 50.25% (401.60 Mb) of the genome. Gypsy-type LTRs (344.52 Mb) were the most common subfamily and six times more abundant than Copia-type LTRs (50.09 Mb). We used a combination of methods, including LTR-FINDER[26], LTR-Harvest[27], and LTR-Retriever[28], to identify intact LTRs. A total of 3,220 LTRs were detected in the SLT1.0 genome assembly, including 1,553 Gypsy-type LTRs and 1,346 Copia-type LTRs. The estimated insertion time of the LTR retrotransposons showed that Gypsy and Copia-type LTRs had a recent and similar burst 0.60-1.00 million years ago (Mya) (Fig. 3A), and were enriched far from coding genes (Fig. 3B). These results indicated that the burst of Gypsy-type LTRs may be the major driving force for the expansion of the tomato genome.
To identify the centromere regions, we detected the top 12 TE subfamilies, including 11 Gypsy and one unknown-type subfamilies, which together comprised over 15.47% of the genome (Fig. 3C). The density of these TE subfamilies along all the chromosomes showed that only the Unknown-type rnd-1_family-4 subfamily (1.65% of the genome) was enriched near centromeres but absent from the rest of the genome (Fig. 3D, Supplementary Fig. S2). In addition, we found that 65.21% of the unanchored Contig/Scaffold sequence length comprised highly repetitive regions. Overall, we predicted 12 potential centromeric regions ranging from 1.90 to 6.90 Mb on the 12 chromosomes.
Comparison of the SLT1.0 and S. pimpinellifolium LA2093 genomes
Structural variations (SVs) between wild and cultivated species can cause many phenotypic differences in domestication traits such as fruit weight and quality[29]. Based on protein homologies between the SLT1.0 and LA2093 genomes, we found that 23,544 genes (68.47%) in the SLT1.0 genome had one-to-one collinear relationships with 23,474 genes (65.64%) in the LA2093 genome (Fig. 4A). In addition, genome collinearity analysis showed that syntenic genomic blocks occupied 95.63% of the SLT1.0 genome and 96.67% of the LA2093 genome, respectively. We also identified 6,647 SVs (more than 1 kb in length) between the SLT1.0 and LA2093 genomes, including 3,054 (45.95%) SVs in 2,862 genes (Fig. 4B). GO analysis showed that these genes were significantly enriched in the function of oxidation-reduction process, photosynthetic electron transport chain and proton-transporting ATP synthase complex (Supplementary Fig. S3). We also identified 4,493,889 SNPs and 2,459,597 indels between the two genomes (Fig. 4B), including 418,844 SNPs and 245,310 indels located in 29,862 genes. We noted that 45,229 nonsynonymous SNPs resided in 18,178 genes and 9,148 frameshift indels in 1,559 genes, including 7,788 located in domestication regions[30]. They were significantly enriched in macromolecular complex, pigment metabolic process, nutrient reservoir activity, and intracellular organelle parts (Fig. 4C), suggesting these genes may have contributed to disease resistance and fruit traits during tomato domestication.
High-density genetic map construction
A high-quality tomato reference genome can provide new insights into the genetic basis of important agronomic traits. Here, we constructed a RIL7 population of 247 progenies derived from a cross between the cultivated line OH88119 (64.27×) and the wild line PI128216 (42.35×). Resequencing of these progenies generated 1.7 Tb of data with an average depth of 6.91× and 97.98% coverage of the SLT1.0 genome. After aligning the reads to the genome, we identified 4,739,716 SNPs between the parental lines, 2,818,901 (59.35%) of which were genotyped in the RIL7 population. To construct a high-density tomato genetic map, we employed a Hidden Markov Model (HMM) to infer all recombination events in the RIL population. A total of 17,726 recombination events were detected across the whole genome, with one recombination event at an average interval of 3.78 kb (Fig. 5A). We found that 1,477 crossovers per chromosome ranged from 1.01 kb to 7.92 Mb in size and that the recombination rate varied across the genome with an average of 11.22 cM/Mb. Furthermore, recombination rate increased with distance from the centromeric regions on chromosomes 5 and 9, but this trend was not evident on the other chromosomes (Fig. 1).
We defined 35,836 raw bins (an average physical length of 22.06 kb) across the entire RIL population and constructed a high-quality set of 17,741 bin markers with a physical length of greater than 1 kb and a minor allele frequency of greater than 0.05. Based on recombination rate, we created a genetic linkage map of 13,973 unique bin markers anchoring twelve linkage groups that corresponded to the twelve chromosomes (Fig. 5B). This genetic map spanned 786.81 Mb, representing more than 98.46% of the assembled SLT1.0 reference genome. Taken together, these analyses produced a high-quality genetic map for genetic research on important agronomic traits in tomato.
QTL mapping of important agronomic traits in tomato
The RIL population exhibited diversity for several important agronomic traits. The cultivated line OH88119 has a larger fruit and seed size, both of which are key traits associated with human selection in tomato. To gauge the accuracy of the genetic bin map, we performed an association study on three agronomic traits (seed size, leaf architecture, and trichomes) and identified five significantly associated QTLs in the SLT1.0 genome (Fig. 5C, Supplementary Fig. S4). Plant glandular trichomes can produce secondary metabolites that defend against herbivores and pathogens[31, 32]. Previous studies have shown that a 618-bp H gene encoding a C2H2 zinc finger protein regulates the formation of multicellular trichomes in tomato leaves[33]. Here, the leaf trichome locus was localized to a smaller interval of 2.82 kb that included two bins on chromosome 10 (Supplementary Fig. S4), suggesting the high quality of the genetic bin map.
The seed is an important plant reproductive organ, but larger seeds affect fruit taste, reducing the economic value of the berry crop. We used the genetic bin map to identify two seed size QTLs above the LOD threshold associated with seed length (Fig. 5C). These major-effect QTLs were resolved into a 1.32-Mb region with 40 bins ranging from 5.75 to 7.07 Mb on chromosome 4 (sl4.1) and a 0.74-Mb region with 36 bins ranging from 4.08 to 4.82 on chromosome 11 (sl11.1). Intriguingly, we found that sl11.1 on chromosome 11 was located in domestication sweep DS169[30], indicating that it may have played an essential role during human selection. We identified 11,709 parental SNPs in this region, including 303 nonsynonymous variations in 87 genes. Among these SNPs, 132 nonsynonymous variations in 65 genes were also present in the RIL population. After exploring their functions, we identified six candidate genes for seed length, including a TIFY domain protein, a lipid droplet-associated protein (LDAP)-interacting protein (LDIP), and three small auxin-up RNA (SAUR) proteins (Fig. 5D). In particular, one gene encoding an LDIP, SlT11G005670, appeared to be a strong candidate for seed length. Based on phylogenetic analysis of homologous genes, we found that the LDIP gene influences lipid droplet size and neutral lipid homeostasis in Arabidopsis seeds[34] (Fig. 5E).