Draft Genomes of Two Artocarpus Plants, Jackfruit (A. heterophyllus) and Breadfruit (A. altilis)

Two of the most economically important plants in the Artocarpus genus are jackfruit (A. heterophyllus Lam.) and breadfruit (A. altilis (Parkinson) Fosberg). Both species are long-lived trees that have been cultivated for thousands of years in their native regions. Today they are grown throughout tropical to subtropical areas as an important source of starch and other valuable nutrients. There are hundreds of breadfruit varieties that are native to Oceania, of which the most commonly distributed types are seedless triploids. Jackfruit is likely native to the Western Ghats of India and produces one of the largest tree-borne fruit structures (reaching up to 45 kg). To-date, there is limited genomic information for these two economically important species. Here, we generated 273 Gb and 227 Gb of raw data from jackfruit and breadfruit, respectively. The high-quality reads from jackfruit were assembled into 162,440 scaffolds totaling 982 Mb with 35,858 genes. Similarly, the breadfruit reads were assembled into 180,971 scaffolds totaling 833 Mb with 34,010 genes. A total of 2822 and 2034 expanded gene families were found in jackfruit and breadfruit, respectively, enriched in pathways including starch and sucrose metabolism, photosynthesis, and others. The copy number of several starch synthesis-related genes were found to be increased in jackfruit and breadfruit compared to closely-related species, and the tissue-specific expression might imply their sugar-rich and starch-rich characteristics. Overall, the publication of high-quality genomes for jackfruit and breadfruit provides information about their specific composition and the underlying genes involved in sugar and starch metabolism.


Introduction
The family Moraceae contains at least 39 genera and approximately 1100 species [1][2][3]. Species diversity of the family is primarily centered in the tropics with variation in inflorescence structures, pollination forms, breeding systems, and growth forms [2]. Within the Moraceae family, the genus Artocarpus is comprised of approximately 70 species [2,4]. The most recent evidence indicates that Borneo was the center of diversification of the Artocarpus genus and that species diversified throughout South and Southeast Asia [2]. All members of the genus have unisexual flowers and produce exudate from laticifers. Inflorescences consist of up to thousands of tiny flowers, tightly packed and condensed on a receptacle [2]. In most species, the perianths of adjacent female flowers are partially to completely fused together and develop into a highly specialized multiple fruit called a syncarp, which is formed by the enlargement of the entire female head. Syncarps of different species range in size from a few centimeters in diameter to over half a meter long in the case of jackfruit [2,5]. Many Artocarpus species are important food sources for forest fauna, and about a dozen species are important crops in the regions where they are from [2,6].
Jackfruit (A. heterophyllus Lam.) is thought to have originated in in the Western Ghats of India and is cultivated as an important food source across the tropics. It is monecious, and thought to be pollinated by gall midges [7]. In some areas it is propagated mainly by seeds [8], however, clonal propagation via grafting is increasing in areas where it is grown for commercial use [9]. On average it contains more than 100 seeds per fruit with viability of less than a month [10,11]. The male flowers are tiny and clustered on an oblong receptacle, typically 5-10 cm in length. Limited studies exist on the range of cultivated varieties of jackfruit, but they are often grouped into two main types, varieties with edible fleshy perianth tissue (often referred to as "flakes") that are either (a) small, fibrous, soft, and spongy, or (b) larger, less sweet crisp fruit [10,12]. The latter type is often more commercially important.
Breadfruit (A. altilis [Parkinson] Fosberg) is most likely derived from the progenitor species A. camansi Blanco, which is native to New Guinea [10,13]. As humans migrated and colonized the islands of Remote Oceania, indigenous people selected and cultivated varieties from the wild ancestor over thousands of years [13], giving rise to hundreds of cultivated varieties [10,[13][14][15]. Cultivated varieties were traditionally propagated clonally by root cuttings but can now be commercially propagated by tissue culture [16,17]. Among the hundreds of varieties, some are diploid (2n = 2x =~56) and may produce seeds, while other varieties are seedless triploids (3n = 2x =~84), and still others are of hybrid origin with another species, A. mariannensis Trécul [13,[18][19][20]. A small subset of the triploid diversity is what has been introduced outside of Oceania [19,21].
To diversify the global food supply, enhance agricultural productivity, and eradicate malnutrition, it is necessary to focus on crop improvement of plants that are utilized in rural societies as a local source of nutrition and sustenance. This study is part of the African Orphan Crops Consortium (AOCC), an international public-private partnership. A goal of this global initiative is to sequence, assemble, and annotate the genomes of 101 traditional African food crops [22,23]. Both breadfruit and jackfruit are nutritious [24][25][26][27] and have the potential to increase food security, especially in tropical areas. Until now limited genomic information has been available for the Artocarpus genus as a whole. Microsatellite markers have been used to characterize cultivars and wild relatives of breadfruit [8,19,21,28], jackfruit [29], and other Artocarpus crop species [6,30,31]. Additionally, an assembled and annotated reference transcriptome of A. altilis has been generated [20]. Twenty-four transcriptomes of breadfruit and its wild relatives revealed signals of positive selection that may have resulted from local adaptation or natural selection [20]. Finally, a low coverage whole genome sequence has been published for A. camansi [32], but full genome sequences for jackfruit and breadfruit are still not available. Here, we report high-quality annotated draft genome sequences for both jackfruit and breadfruit. The results help explain their energy-dense fruit composition and the underlying genes involved in sugar and starch metabolism.

Sample Collection, NGS Library Construction, and Sequencing
Genomic DNA was extracted from fresh leaves of A. heterophyllus (ICRAFF 11314) and A. altilis (ICRAFF 11315), grown at the World AgroForestry (ICRAF) campus in Kenya, using a modified CTAB method [33]. Extracted DNA was used to construct four paired-end libraries (170, 350, 500, and 800 bp) and four mate-pair libraries (2, 6, 10, and 20 Kb) following the standard protocols provided by Illumina (San Diego, CA, USA). Subsequently, the sequencing was performed on a HiSeq 2000 platform (Illumina, San Diego, CA, USA) using a whole genome shotgun sequencing strategy. To improve the data quality, the poor quality reads were filtered using SOAPfilter (v2.2) [34] with the following parameters: (1) low-quality bases with Q = below 7 and 15 for A. altilis and A. heterophyllus, respectively were trimmed from forward and reverse reads; (2) reads with ≥30% low quality bases (quality score ≤ 15); (3) reads with ≥10% uncalled ("N") bases; (4) reads with adapter contamination or PCR duplicates; (5) reads with undersized insert sizes were discarded. Finally, more than 100 high-quality reads were obtained for each species (see Supplementary Materials: Table S1).
For transcriptome sequencing, the RNA was extracted from different tissues of A. altilis (various leaf stages, leaf buds, and roots) and A. heterophyllus (various leaf stages, leaf bud, stem, bark, roots, germinated seed, and seedling). The RNA was extracted using the PureLink RNA Mini Kit (Thermo Fisher Scientific, Carlsbad, CA, USA) according to the manufacturer's instructions. For each sample, RNA libraries were constructed by following the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) manual, and were then sequenced on the Illumina HiSeq 2500 platform (paired-end, 100-bp reads), generating more than 47 Gb of sequence data for each species. Data were then filtered using a similar criterion as used to filter DNA NGS data, with a slight modification: (1) reads with ≥10% low-quality bases (quality score ≤ 15) were removed; and (2) reads with ≥5% uncalled ("N") bases were removed (see Supplementary Materials: Table S2). All of the transcriptome data were compiled, and the combined version was used to check the completeness of the whole genome sequence assembly.

Evaluation of Genome Size
Clean reads of the paired-end libraries (170, 250, and 500 bp) were used to estimate the genome size by k-mer frequency distribution and heterozygosity analysis. The genome size was estimated based on the following formula: where N represents the number of used reads, L represents the read length, K represents the k-mer value in the analysis, and K_depth refers to the location of the main peak in the distribution curve [35]. The heterozygosity was evaluated by the GCE software [36].

De Novo Genome Assembly
The de novo genome assembly tool, Platanus (Platanus, RRID: SCR_015531) [37], was used to construct the contigs and scaffolds in three steps: contig assembling, scaffolding and gap closing. In contig assembling, paired-end libraries ranging from 170 to 800 bp were used with the parameters "-d 0.5 -K 39 -u 0.1 -m 300". In the scaffolding step, paired-end and mate-pair information were used to with parameters "-u 0.1". Lastly, the paired-end reads were used for gap closing using GapCloser version 1.12 (GapCloser, RRID: SCR_015026) [34], with the parameters "-l 150 -t 32 -p 31".

Genome Assembly Evaluation
The genome assembly completeness was assessed using BUSCO (Benchmarking Universal Single-Copy Orthologues), version 3.0.1 (BUSCO, RRID: SCR_015008) [38]. Next, the unigenes generated by Bridger software [39] from the transcriptome data of each species were aligned to the assembled genomes using BLAT (BLAT, RRID: SCR_011919) [40] with default parameters. In order to confirm the accuracy of the assembly, some of the paired-end libraries (170, 250, and 350 bp) were aligned to the assembled genomes, and the sequencing coverage was calculated using SOAPaligner, version 2.21 (SOAPaligner/soap2, RRID: SCR_005503) [41].
We calculated the GC content and average depth with 10 kb non-overlapping windows. The distribution of GC content indicated a relative pure single genome without contamination or GC bias (Supplementary Materials: Figure S3). Moreover, the GC content of A. altilis and A. heterophyllus genomes were also compared with three rosids species (Fragaria vesca, Malus domestica, and Morus notabilis).

Repeat Annotation
Repetitive sequences were identified by using RepeatMasker (version 4-0-5) [42], with a combined library consisting of the Repbase library and a custom library obtained through careful self-training. The custom library was composed of three parts: the MITE, LTR, and an extensive library, which were constructed as described below. First of all, the library of miniature inverted-repeat transposable elements (MITEs) was created by annotation using MITE-hunter [43] with default parameters. Secondly, the library of long terminal repeats (LTR) was constructed using LTRharvest [44] integrated in Genometools (version 1.5.8) [45] with parameters "-minlenltr 100 -maxlenltr 6000 -mindistltr 1500 -maxdistltr 25,000 -mintsd 5 -maxtsd 5 -similar 90 -vic 10" to detect LTR candidates with lengths of 1.5 kb to 25 kb, including two terminal repeats ranging from 100 bp to 6000 bp with ≥ 85% similarity. In order to improve the quality of the LTR library, we used several strategies to filter the candidates. As intact PPT (poly purine tract) or PBS (primer binding site) was necessary to define LTR, we subsequently used LTRdigest [46] with a eukaryotic tRNA library [47] to identify these features, and then removed the elements without appropriate PPT or PBS locations. Subsequently, to remove contamination, like local gene clusters and tandem local repeats, 50 bp of flanking sequences on both sides of each LTR candidate were aligned using MUSCLE (MUSCLE, RRID:SCR_011812) [48] with default parameters: if the identity ≥ 60%, the candidate was taken as a false positive and removed. LTR candidates which nested with other types of elements were also removed. Exemplars for the LTR library were extracted from the filtered candidates using a cutoff of 80% identity in 90% of sequence length. Furthermore, the regions annotated as LTRs and MITEs in the genome were masked, and then put into RepeatModeler version 1-0-8 RepeatModeler, RRID: SCR_015027) to predict other repetitive sequences for the extensive library.
Finally, the MITE, LTR and extensive libraries were integrated into the custom library, which was combined with the Repbase library and then taken as the input for RepeatMasker to identify and classify repetitive elements genome wide.

Gene Prediction
Repetitive regions of the genome were masked before gene prediction. Based on the RNA, homologous and de novo prediction evidence, the protein-coding genes were identified using the MAKER-P pipeline (version 2.31) [49]. For RNA evidence, the clean transcriptome reads were assembled into inchworms using Trinity version 2.0.6 [50], and then fed to MAKER-P as EST evidence. For homologous evidence, the protein sequences from four relative species in the rosids (F. vesca, M. domestica, M. notabilis, Prunus persica, Ziziphus jujuba) were downloaded and provided as protein evidence.
For de novo prediction evidence, a series of training attempts were made to optimize different ab initio gene predictors. At first, a set of transcripts were generated by a genome-guided approach using Trinity with parameters "-full_cleanup -jaccard_clip -genome_guided_max_intron 10,000 -min_contig_length 200". The transcripts were then mapped back to the genome using PASA (version 2.0.2) [51] and a set of gene models with real gene characteristics (e.g., size and number of exons/introns per gene, features of spicing sites) was generated. The complete gene models were picked for training Augustus [52]. Genemark-ES (version 4.21) [53] was self-trained with default parameters. The first round of MAKER-P was run based on the evidence above with default parameters except "est2genome" and "protein2genome" set to "1", yielding only RNA-and protein-supported gene models. SNAP [54] was then trained with these gene models. Default parameters were used to run the second and final round of MAKER-P, producing final gene models.

Ks-Distribution Analysis
The coding sequences and annotations for Morus notabilis were downloaded from the NCBI, reference RefSeq assembly accession GCF_000414095.1 [66]. The coding sequences and annotations for Ziziphus jujube [67] were downloaded from the Plaza4 database [68]. The headers of the fasta files, as well as the ninth columns of the gff3 files were edited to make the datasets compatible with the software packages used for downstream analysis.
Ks-distribution analyses were performed, using the wgd-package [69]. For each species, the paranome was obtained by performing an all-against-all BlastP [70], with MCL clustering [71]. Codon multiple sequence alignment was done using MUSCLE [48]. Ks-distributions were constructed using codeml from the PAML4 package [72] and Fast-Tree [73] for inferring phylogenetic trees used in the node weighting procedure, other software used by the wgd. Thereafter, i-ADHoRe [74] was used to get anchor-point distributions and produce dot-plots. Lastly, Gausian mixture modes were fitted using 1-5 components.

One vs. One Synteny
One-vs.-one synteny analysis was performed for pairs of the above-mentioned species, using the "work-flow 2" script that is part of the wgd-package [69]. In order to compare A. altilis, A. heterophyllus, we use MCScanX to detect the synteny and collinearity, and MCscan (Python version) for the visualization.

Phylogenetic Analysis and Divergence Time Estimation
We identified 486 single-copy genes in the nine species, and subsequently used them to build the phylogenetic tree. Coding DNA sequence (CDS) alignments of each single-copy family were constructed following the protein sequence alignment with MUSCLE (MUSCLE, RRID: SCR_011812) [48]. The aligned CDS sequences of each species were then concatenated to a super gene sequence. The phylogenetic tree was constructed with PhyML-3.0 (PhyML, RRID: SCR_014629) [78] with the HKY85+ Gamma substitution model on extracted four-fold degenerate sites. Divergence time was calculated using the Bayesian relaxed molecular clock approach using MCMCTREE in PAML (PAML, RRID: SCR_014932) [72], based on the published calibration times (divergence between Arabidopsis thaliana and Rosales was 108-109 Mya, divergence between P. mume and P. persica was 24-72 Mya) [66]. The divergence time between M. notabilis and Artocarpus was predicted to be 61.8 (54.1-76.0) Mya ( Figure 1A). Subsequently, to study gene gain and loss, CAFE (CAFE, RRID:SCR_005983) [79] was employed to estimate the universal gene birth and death rate λ (lambda) under a random birth and death model with the maximum likelihood method. The results for each branch of the phylogenetic tree were estimated ( Figure 1A). Enrichment analysis on GO and the pathway of genes in expanded families in the Artocarpus lineage were also calculated.

Genome Sequencing and Assembly
A total of eight libraries were constructed including four short-insert libraries (170 bp, 350 bp, 500 bp, and 800 bp) and four mate-pair libraries (2 kb, 5 kb, 10 kb, and 20 kb) for Illumina Hiseq2000 sequencing. In total, 273 Gb and 227 Gb of raw data was generated from A. heterophyllus and A. altilis, respectively (Supplementary Materials: Table S1). We used the GCE software to evaluate the heterozygosity, and the results showed that the heterozygous ratio is 1.13% and 0.911% for A. altilis and A. heterophyllus, respectively. The K-mer distributions of A. altilis and A. heterophyllus showed two distinct peaks (Supplementary Materials: Figures S1 and S2), the first peak was the heterozygous peak, the second peak was the homozygous peak, where the second peak was confirmed as the main one for each of the species. Based on K-mer frequency methods [36], the A. heterophyllus and A. altilis genomes were estimated to be 1005 MB and 812 Mb, respectively (Supplementary Materials: Figure S1, Table S3). The genome sizes of A. altilis and A. heterophyllus were relatively close to the genome size of species in the genus Artocarpus based on existing data in the 1C-values database of 1.2 pg.
Using the SOAPdenovo2 program [41], all of the A. heterophyllus high-quality reads were assembled into 108,267 scaffolds, totaling 982 Mb ( Table 1). The N50s of contigs and scaffolds were 27 kb and 548 kb, with the longest being 255 kb and 3.1 Mb, respectively (Table 1, Supplementary Materials: Figure  S3). Similarly, for A. altilis, the N50s of contigs and scaffolds were 17 kb and 1.5 Mb with the longest being 174 kb and 7.4 Mb respectively (Table 1, Supplementary Materials: Figure S3). These results indicate the high quality of the assemblies for both species. The GC content of the A. heterophyllus and A. altilis genomes were 32.9% and 32.3%, respectively. The GC depth graphs and distributions indicated there was no contamination in the genome assemblies (Supplementary Materials: Figure S4). Evaluation of the quality and completeness of the draft genome assembly was done with the Benchmarking Universal Single-Copy Orthologs (BUSCO) datasets [38]. Of the total of 1440 BUSCO ortholog groups searched in the A. heterophyllus assembly, 932 (64.7%) BUSCO genes were "complete single-copy", 437 (30.3%) were "complete duplicated", 15 (1%) were "fragmented", and 56 (4%) were "missing" (Table 2). Similarly, in A. altilis, 988 (68.6%) BUSCO genes were "complete single-copy", 383 (26.6%) were "complete duplicated", 14 (~1%) were "fragmented", and 55 (3.8%) were "missing" ( Table 2), suggesting that the quality of the genome assembly is high. From the 1440 core Embryophyta genes, 1371 (95.2%) and 1369 (95.1%) were identified in the A. altilis and A. heterophyllus assemblies, respectively ( Table 2). We observed a significant difference in the number of duplicated core genes in A. altilis and A. heterophyllus (Table 2), which might be ascribed to the genome duplication in these species. The results also indicated that the assembly covered more than 90% of the expressed unigenes, (Table 3). As expected, after the comparative GC content analysis the close peak positions showed A. altilis, A. heterophyllus, and M. notabilis are closer than other species in GC content (Supplementary Materials: Figure S5).

Gene Annotation
A combination of de novo and homology-based methods (using transcript data as evidence) were used to identify repeat sequences. We found that up to 51.01% of the A. heterophyllus and 52.04% of the A. altilis assembled sequences were repeat sequences, comprised mostly of transposable elements and tandem repeats. Interestingly, the amounts of these elements were higher than what is observed in orange (20%, 367 Mb) [80], peach (29.6%, 265 Mb) [81], pineapple (38.3%, 526 Mb) [82], and others (Table 4). This is consistent with the finding that larger fruit tree genomes often retained higher percentages of repetitive elements compared to the smaller fruit tree genomes [83]. Among the repetitive sequences, 37.0% and 46.0% were of the long terminal repeat (LTR) type, respectively (Table 4), indicating LTRs are the most abundant transposable elements in A. heterophyllus and A. altilis genomes. Using a comprehensive annotation strategy, we annotated a total of 35,858 A. heterophyllus genes and 34,010 A. altilis genes (Table 5). This was close to the number of genes (39,282) predicted in Dimocarpus longan, an exotic round to oval Asian fruit [83]. The average A. heterophyllus gene length was 3472 bp, the average length of the coding sequence (CDS) was 1241 bp, and the average number of exons per gene was 5.5 (Table 5, Supplementary Materials: Figure S6). We predicted a total of 466 rRNA, 159 miRNA, 1554 snRNA genes and 713 tRNA in A. altilis; and a total of 2706 rRNA, 168 miRNA, 1005 snRNA genes and 689 tRNA in A. heterophyllus (Table 6) (Table 7). BUSCO evaluation showed that more than 89% of 1440 core genes were complete, suggesting an acceptable gene annotation for A. altilis and A. heterophyllus genomes (Supplementary Materials: Table S4)

Gene Family Evolution and Comparison
Orthologous clustering analysis was conducted with the A. altilis and A. heterophyllus genomes following comparison with seven other plant genomes: A. thaliana, F. vesca, M. domestica, M. notabilis, P. mume, P. persica, and Z. jujuba. A Venn diagram shows that A. altilis, A. heterophyllus, A. thaliana, M. notabilis, and Z. jujuba contain a core set of 9462 gene families in common; there were 1028 orthologous families shared by three Moraceae species; while 329 gene families containing 515 genes were specific to A. altilis; and 420 gene families containing 907 genes were specific to A. heterophyllus ( Figure 1C).
Phylogenetic analysis showed that A. heterophyllus and A. altilis were more closely related to mulberry than to Jujube ( Figure 1A), further supporting a previous phylogeny of Artocarpus [2]. CAFE [79] was used to identify gene families that had potentially undergone expansion or contraction. We found a total of 2822 expanded gene families and 1497 contracted families in A. heterophyllus, as well as 2034 expanded and 1800 contracted families in A. altilis ( Figure 1A). The genes in the expanded and contracted families were assigned to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [84]. The A. heterophyllus-expanded gene families were remarkably enriched in metabolism related pathways/functions, including starch and sucrose metabolism (ko00500, p = 0.003), glycan degradation (ko00511, p = 0.007), glycolysis/gluconeogenesis (ko00010, p = 0.016), and others (Supplementary Materials: Table S7). KEGG enrichment analysis of A. altilis revealed that pathways associated with photosynthesis, such as carbon fixation in photosynthetic organisms (ko00710, p = 0.017), other types of O-glycan biosynthesis (ko00514, p = 0.018), and photosynthesis (ko00195, p = 0.006) were particularly enriched (Supplementary Materials: Table S7).
Collinearity between the largest orthologous scaffolds were determined ( Figure S7B) and the result indicates conserved shared synteny. In order to determine whether there is any evidence for whole genome duplications in A. heterophyllus and A. altilis, the distance-transversion rates at four-fold degenerate sites (4DTv) was calculated ( Figure 1D, Supplementary Materials: Figure S7A). Two 4DTv values that peaked at 0.07 and 0.08 for orthologs between A. heterophyllus, and between A. altilis, respectively, which highlighted the recent whole-genome duplication of these two species. The results of the Ks distributions mostly corroborate the findings of the 4DTv analysis. The results suggest that the whole genome duplication event was shared by A. altilis and A. heterophylus. Their divergence is recent, as suggested by the overlap of their WGD peaks (Figure 2), meaning that they have equal substitution, duplication, and loss rates. Thus, for further analysis (one-vs.-one synteny with the close relatives M. notabilis and Z. jujube), only A. altilis was used. These results suggest that the Artocarpus genome duplication event occurred after divergence from the common ancestor they share with M. notabilis (Figure 3), thus, between 62 and 10 MYA. A total of 33,614 gene pairs were identified in 2694 syntenic blocks.

Gene Family Expansion and Tissue Specific Expression of Starch Synthesis-Related Genes
The copy number of starch synthesis related genes were compared between A. heterophyllus, A. altilis, closely related species, as well as some other starch-rich plant species (Figure 4). We observed a remarkable copy number expansion of the UGD1 gene in A. heterophyllus compared with the other species. The enzyme encoded by UGD1, catalyzes the conversion of Glucose-1-P into UDP-GlcA, thereby stalling the starch synthesis process [85] (Figure 4). Interestingly, the tissue-specific expression pattern of UGD1 contrasts with other starch synthesis genes in A. heterophyllus ( Figure 5A). For instance, in A. heterophyllus there is a suppression of UDPG transcription in the stem, while the other starch biosynthesis genes are activated. However, differential expression of UGD1 was not shown in A. altilis ( Figure 5B). This unusual expression pattern of UGD1, as well as the gene copy number expansion, might lead to the failure of starch accumulation in A. heterophyllus rather than A. altilis. However, this needs to be further validated by real-time qPCR for confirmation of the tissue-specific expression. For the GO enrichment, expansion of gene families were related to small molecule binding or single organism signaling (Supplementary Materials: Table S6) in A. altilis. Moreover, there were some expansion of gene families related to molecule binding, reproductive process, and cellular response to stimulus in A. heterophyllus. Gene families belonging to expanded pathways in A. altilis were mainly related to plant-pathogen interaction, lysine biosynthesis or photosynthesis. In contrast, the gene families that were expanded in A. heterophyllus belonged to pathways involving secondary metabolite biosynthesis, phenylpropanoid biosynthesis, and fatty acid metabolism. In contrast, the biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and fatty acid metabolism were enriched in the expanded gene families in A. heterophyllus (Supplementary Materials: Table S7).

Conclusions
Here, we report the characterization of the genomes of jackfruit (A. heterophyllus) and breadfruit (A. altilis). The publication of these high-quality draft genomes and annotations may provide plant breeders and other researchers with useful information regarding trait biology and their subsequent improvement. In particular, we highlight genes unique to A. heterophyllus and A. altilis due to their high sugar and starch content, respectively, which are desirable characteristics in these edible plants. The information provided in the draft genome annotations can be used to accelerate genetic improvement of these crops. The availability of these genomes on the AOCC ORCAE platform (https://bioinformatics.psb.ugent.be/orcae/aocc) will enable various stakeholders to access and improve the annotations of these genomes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/1/27/s1, Figure S1. K-mer (K = 17) analysis of the two genomes, Figure S2. Distribution of sequencing depth of the assembly data, Figure S3. Distribution of the length and number of the scaffold in two species, Figure S4. The distribution of GC content, Figure S5. Comparison of GC content across closely related species, Figure S6. Statistics of gene models in A. altilis, A. heterophyllus, F. vesca, M. domestica, M. notabilis, Prunus persica, and Ziziphus jujube, Figure S7. The collinearity between two species, Table S1. Statistics of the raw and clean data of DNA sequencing, Table S2. Summary statistics of the transcriptome data, Table S3. Estimation of the genome size based on K-mer statistics, Table S4. BUSCO evaluation of the annotated protein-coding genes in A. altilis and A. heterophyllus, Table  S5. Analysis of gene families of different species, Table S6. Enriched GO terms (level 3) of genes in families with expansion, Table S7. Enriched pathways of genes in families with expansion.  We also thank Arthur Zwanepoel for his insights and technical assistance. This work is part of 10KP project.

Conflicts of Interest:
The authors declare that they have no competing interests. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.