Gene prediction was conducted through a combination of homology-based prediction, ab initio prediction and transcriptome-based prediction methods. Protein repertoires of plants including A. comosus10, Amborella trichopoda45, Arabidopsis thaliana (phytozomev10), Brachypodium distachyon (phytozomev10), D. officinale9, O. sativa (phytozomev10), P. equestris7, Vitus vinifera (phytozomev10), Sorghum bicolor (phytozomev10) and Zea mays (phytozomev10) were downloaded and mapped to the G. elata genome using TBLASTN (E-value ≤ 1e−5). The BLAST hits were conjoined by Solar software46. GeneWise (version 2.4.1)47 was used to predict the exact gene structure of the corresponding genomic region on each BLAST hit. Homology predictions were denoted as “Homology-set”. RNA-seq data derived from protocorm, juvenile tuber, immature tuber, mature tuber, and scape (Fig. 1a) were assembled by Trinity (version 2.0) 41. The Trinity assembly included 183,515 contigs with an average length of 592 bp. These assembled sequences were aligned against the G. elata genome by PASA (Program to Assemble Spliced Alignment)48. Valid transcript alignments were clustered based on genome mapping location and assembled into gene structures. Gene models created by PASA were denoted as PASA-T-set (PASA Trinity set). Besides, RNA-seq reads were directly mapped to the genome using Tophat (version 2.0.8)49 to identify putative exon regions and splice junctions; Cufflinks (version 2.1.1) was then used to assemble the mapped reads into gene models (Cufflinks-set). Augustus (version 2.5.5)50, GeneID (version)51, GeneScan (version 1.0)52, GlimmerHMM (version 3.0.1)53, and SNAP (version)54 were also used to predict coding regions in the repeat-masked genome. Of these, Augustus, SNAP and GlimmerHMM were trained by PASA-H-set gene models. Gene models generated from all the methods were integrated by EvidenceModeler (EVM)48. Weights for each type of evidence were set as follows: PASA-T-set > Homology-set > Cufflinks-set > Augustus > GeneID = SNAP = GlimmerHMM = GeneScan. The gene models were further updated by PASA2 to generate UTRs, alternative splicing variation information, which generated 26,872 gene models. Gene models only supported by ab initio evidence were filtered out. To reduce the possibility of missing and poorly annotated genes, we invested additional effort in annotating some gene families that could be missed by automated genome annotation, such as NBS-encoding genes. In total, 1943 protein sequences containing an NB-ARC domain were searched against the G. elata genome using TBLASN with a threshold of 1e−5. All BLAST hits in the genome, together with 5000 bp flanking regions on both sides, were annotated by the GeneWise program. The resulting predictions were surveyed to verify whether they encoded NBS or LRR motifs using Pfam. We also focused on other genes, such as those related to photosynthesis, and transporter, and these were manually annotated through a combination of BLAST search and motif verification. Ultimately, a comprehensive non-redundant reference gene set was produced that contained 18,969 protein-coding gene models. Functional annotation of the protein-coding genes was carried out using BLASTP (E-value cut-off 1e−05) against two integrated protein sequencing databases, SwissProt and TrEMBL55. Protein domains were annotated by searching against InterPro (Version 5.16)56 and Pfam (Version 3.0) database57, using InterProScan (version 4.8) and HMMER (version 3.1b1) (http://hmmer.janelia.org), respectively. The GO terms for genes were obtained from the corresponding InterPro or Pfam entry. The pathways in which the genes might be involved were assigned by BLAST against the KEGG databases (release 20150831)58 with the E-value cut-off of 1e−05.