Current status and impending progress for cassava structural genomics.

KEY MESSAGE: We demystify recent advances in genome assemblies for the heterozygous staple crop cassava (Manihot esculenta), and highlight key cassava genomic resources. Cassava, Manihot esculenta Crantz, is a crop of societal and agricultural importance in tropical regions around the world. Genomics provides a platform for accelerated improvement of cassava's nutritional and agronomic traits, as well as for illuminating aspects of cassava's history including its path towards domestication. The highly heterozygous nature of the cassava genome is widely recognized. However, the full extent and context of this heterozygosity has been difficult to reveal because of technological limitations within genome sequencing. Only recently, with several new long-read sequencing technologies coming online, has the genomics community been able to tackle some similarly difficult genomes. In light of these recent advances, we provide this review to document the current status of the cassava genome and genomic resources and provide a perspective on what to look forward to in the coming years.

Introduction

Manihot esculenta Crantz, referred to as cassava, yuca, or manioc, among other common names, is an important staple food source for over a billion people in tropical regions around the globe (Lebot, 2019). This crop develops large storage roots that are harvested for starch. Cassava leaves are also eaten but are considered a relatively minor source of calories. Cultivated cassava (M. esculenta ssp. esculenta) was domesticated in the Amazon, likely from M. esculenta ssp. flabellifolia (Pohl) (Olsen 2004; Léotard et al. 2009; Olsen and Schaal 1999; Olsen and Schaal 2001; Hillocks et al. 2002). While M. glaziovii Muell. Arg., the Ceará or India rubber tree, is sometimes erroneously referred to as “wild cassava,” genome sequencing has confirmed that it is a distinct species that diverged from cassava approximately 3 million years ago (Bredeson et al. 2016). However, some cassava accessions contain genomic regions derived from M. glaziovii, the result of an interbreeding program in the 1930s intended to confer disease resistance (Jennings 1976; Bredeson et al. 2016).

In contrast to the majority of staple food sources that are planted each year from botanical seed, such as cereals and legumes, cassava is propagated clonally through stem cuttings. A relevant consequence of this agricultural practice is that it circumvents sexual reproduction. Among sexual organisms, recombination that occurs during meiosis is a major driver of diversification that fuels evolution. Thus, in clonally propagated plants like cassava, we might expect limited genetic diversity and a slower rate of evolution that ultimately yields a weak crop. However, there are thousands of distinct varieties of cassava and this crop is prized for agriculturally important phenotypes including tolerance to abiotic and biotic stresses and an ability to produce high yields without external inputs such as irrigation and chemical fertilizer (Food and Organization 2013). Cassava breeding is difficult, in part because of high levels of heterozygosity, asynchronous flowering, and low seed set (Halsey et al. 2008). However, cassava is able to outcross and so through determined breeding efforts, significant increases in cassava yields have been observed over the mid and late 1900s (Ceballos et al. 2020); Fig. 1).

Fig. 1

Global cassava yields show potential for improvement. a Cassava global median yields trends (1961–2017). Natural cubic spline smoothed trendline (blue) and standard errors (shaded ribbon). Cassava yields rose significantly since the 1980s, possibly due to improvements in germplasm and breeding. More recently however, there has been a plateau in yields. b Large disparity in global yields around cassava producing regions also suggests there is still potential for large scale gains. Yields in each cassava producing country plotted relative to the maximum produced in 2017 (32 Tons/ha, Lao People’s Democratic Republic) (Source: FAOSTAT, December 2019)

For an in depth review of cassava agricultural attributes, we point readers to existing, recent reviews (Food and Organization 2013; Fanou et al. 2018; Bart and Taylor 2017; Lin et al. 2019). In this review, we focus on recent advances within the cassava genomics community, highlighting publicly available resources and specific aspects of cassava that make genomics research on this crop difficult but fascinating. We compare and contrast different sequencing and genomics technologies and refer to examples from other organisms as a guide for future efforts on cassava. This review is targeted to a broad audience interested in generating and using genomics data for cassava research. Further, we hope this review will be useful and interesting to researchers outside the cassava community who are interested in heterozygous genomes.

Cassava Reference Genome

The year 2019 marked the ten-year anniversary of the first public release of a cassava reference genome assembly, making this a fitting time to review the history of its development as a public resource. Cassava was one of the first “orphan” crop genomes to be sequenced, and the quality of the cassava genome sequence has steadily improved over the past decade, taking advantage of ongoing developments in genome sequencing and bioinformatics (Table 1). While plant chromosomes are tens to hundreds of megabases (Mb) in length, genome sequencers are only able to “read” shorter sequences (hundreds to thousands of base pairs, depending on technology). These reads are computationally combined, or “assembled,” in order to reconstruct chromosomal sequences. Here we briefly summarize the development of the cassava genome.

Table 1

Cassava AM560-2 reference genome assemblies

	v4.1	v5.1	v6.1	v7.1
Release	2009	2014	2016	2019
Primary sequence technology	454	454	Illumina	PacBio
Primary scaffolding data	454 mate pair	Composite genetic map	Illumina mate pair, fosmid, genetic maps	Illumina mate pair, fosmid, genetic maps
Total scaffold length	533 Mb	534 Mb	582 Mb	669 Mb
Total contig length	419 Mb	419 Mb	496 Mb	667 Mb
Number of chromosomal scaffolds	0	18	18	18
Bases in chromosomes	0	305 Mb	444 Mb	606 Mb
Contig N50 length (contiguity)	11 kb	11 kb	27 kb	693 kb
Annotated genes	30,666	30,666	33,033	33,849

Cassava genomics has emphasized the “whole genome shotgun” approach, in which a large number of random, overlapping genome fragments (“reads”) are sequenced. Since many reads are redundantly derived from each genomic locus, such reads can be identified bioinformatically and combined to reconstruct the underlying genome sequence. For outbred species, however, a single individual has two haplotypes, inherited from its maternal and paternal parents (Fig. 2b). In M. esculenta, haplotypes typically differ by ~0.7% at the single nucleotide level, which may include short insertions and deletions. More dramatic differences (e.g., large insertions and structural rearrangements) may also be present but have not yet been comprehensively investigated. Thus, the task of assembly is greatly simplified if a homozygous genotype can be used as the reference, so that the two haplotypes can be collapsed into a single haploid representation of the genome sequence (Fig. 2a). Accordingly, inbred lines were used for the first plant genome sequences, Arabidopsis and rice (Arabidopsis Genome Initiative 2000; International Rice Genome Sequencing Project 2005).

Fig. 2

Reference genome assembly strategies. a Generating a haploid representation (reference assembly) of a diploid inbred genome is relatively straightforward. Due to homozygosity, sequence reads from the two haplotypes assemble together. b The heterozygosity present in a diploid outbred genome means that sequences from maternal and paternal haplotypes (blue and gold) will tend to assemble separately. In this case, to generate a haploid reference assembly, researchers can either combine maternal and paternal contigs into a haploid representation for each chromosome (haplotype-mosaic reference assembly), or they can try to fully assemble the maternal and paternal chromosomes, choosing one or the other to represent each chromosome in the reference (haplotype-phased reference assembly). Gray, assembly gaps

To minimize the bioinformatic challenges of heterozygosity in cassava, the reference genotype was chosen to be the inbred accession AM560-2. This accession was generated by Hernan Ceballos at the International Center for Tropical Agriculture (CIAT) by repeatedly selfing MCol1505, a Colombian cassava with good agronomic performance and excellent cooking quality (H. Ceballos, personal communication). As the product of three generations of selfing (S3), an estimated 94% of the AM560-2 genome is homozygous (Bredeson et al. 2016). In the reference genome sequence, the remaining heterozygous 6% of the genome (distributed across 10 long blocks larger than 1 Mb) is represented as an alternating mosaic of the two haplotypes.

In typical usage, a reference genome serves as a unifying platform for studying the gene content, gene expression, and genetic variation of a species, and provides common reference coordinates for labeling genes and genetic variants. Since the chromosomes of an outbred, sexual species are continually recombining with one another, there is no “special” genotype, and any haploid reference is generally as good as any other for the purpose of developing genetic markers, aligning genetic data from other individuals, and studying gene expression through methods such as RNA-seq. Thus the AM560-2 reference has general utility as an organizing platform, and has even been used to analyze sequence from other species like M. glaziovii (Bredeson et al. 2016). The primary drawback to having a single reference arises from insertion-deletion-duplication variation at the scale of hundreds to thousands of bases. If the reference genotype has a deletion relative to other members of the species, or if other members have an insertion, the corresponding sequence data from those individuals may not align to the reference, since it does not include the relevant locus. This drawback can be addressed by expanding the concept of a reference genome sequence to a “pan-genome” that seeks to represent all the common structural variants in a species (Gao et al. 2019; Gordon et al. 2017; Li et al. 2014; Zhao et al. 2018; Hirsch et al. 2014).

The AM560-2 cassava reference genome has seen four major releases, evolving along with developments in genome sequencing technologies and algorithms. This evolution is summarized in Table 1, which shows quantum improvements based on several metrics for genome assembly quality. The ideal assembly would capture the entirety of every chromosome, from telomere to telomere, without any gaps. In practice, transposable elements and other repetitive sequences scattered throughout the genome are difficult to resolve. Assembling repetitive sequence is especially problematic in the extended pericentromeric regions of cassava (Bredeson et al. 2016). Long and accurate reads are useful for overcoming these challenges, since small differences between repeated units, or variation in flanking sequence, can be used to bioinformatically assign each distinct repetitive copy to a unique position in the genome sequence. However, in the absence of a read that passes through a repetitive or otherwise complex region, current practice is to link the sequence on either side of the repetitive region, placing a string of “N”s in between to signify the unknown sequence. This is known as a gap in the assembly. Depending on the linking or “scaffolding” technology used to span gaps, it is often possible to estimate gap sizes. For a detailed review of contemporary approaches to plant genome assembly, see (Michael and VanBuren 2020).

Contiguous sequences free of gaps are called “contigs,” and sequences containing gaps are referred to as “scaffolds” or “super-contigs.” A measure of the typical length of a sequence is the “N50 length,” which is the weighted median contig length of assembled sequence, that is, the median number of bases between gaps in the assembly (Lander et al. 2001). Thus, it is a measure of assembled sequence contiguity. Table 1 shows that the typical length of contiguous sequence between gaps has grown by a factor of 60 between the 2009 and 2019 releases and is now nearly 700 kilobases (kb). Similarly, the total contig length has grown by 50%, reflecting the greater coverage of repetitive regions within the latest reference genome. The number of predicted genes has largely remained stable at ~ 30–33k, although this metric obscures the considerable improvements in the quality of individual gene predictions, which are now bolstered by extensive transcriptomic data developed over the past decade, as well as interspecific comparisons that take advantage of sequence similarity between encoded proteins across plant species. The stability of the gene number, despite dramatic improvements in contiguity and total sequence, implies that the major improvements have come in resolving repetitive sequences. This interpretation is supported by the anticorrelation between genes and repeats, with genes clustered near the distal ends of chromosomes and repetitive elements clustered in extended pericentromeric regions (Fig. 3).

Fig. 3

Repeats, genes, and recombination frequency in the AM560-2 v7 cassava genome. Repeat density (light blue lines), gene count (blue lines), and recombination rate (gold lines) are plotted. Genic regions are anticorrelated with repetitive regions (Y-axis). Regions with low recombination frequency tend to co-occur with areas of high repeat density, thus, these hard-to-assemble regions also tend not to benefit from scaffolding information provided by a genetic map. Repeat density is measured as the fraction of bases that are annotated as repetitive in 1 Mb sliding windows sampled every 100 kb along the AM560-2 v7 chromosomes. The gene count was also taken with 1 Mb sliding windows every 100 kb. Recombination rate is measured as the number of recombinations per 1 Mb sliding window (100 kb step) using the first derivative of a natural cubic spline-smoothed fit line to the ICGMC 2014 framework map anchored to the v7 genome sequence. The marker positions of the framework map are plotted with vertical black ticks below the X-axis

A major development in the evolution of the cassava reference genome was the construction of a dense genetic map, which allowed assembled sequences to be organized into chromosome-scale sequences (these are sometimes referred to as “pseudomolecules”) (ICGMC 2014). Smaller scaffolds that are not situated into chromosome-scale scaffolds are often called “unplaced.”

The first cassava reference genome, released in 2009 as version 4.1 (v4.1), was constructed using the Roche 454 sequence technology, in a collaboration led by Steve Rounsley, Dan Rokhsar, Chinnappa Kodira, and Tim Harkins, with support from the Bill and Melinda Gates Foundation, the U.S. Department of Energy Joint Genome Institute, and 454 Genomics. Although fragmented, it assembled the vast majority of gene space (Table 1) (Prochnik et al. 2012) and provided a first look at the repeat composition of the cassava genome. Version 5 (v5) improved upon the v4 assembly by elevating it to chromosome scale, ordering and orienting 57% of the v4 assembled sequences into n = 18 chromosomes by taking advantage of a dense composite genetic map ((ICGMC 2014; De Carvalho and Guerra 2002). With the advent of new sequencing technologies, the version 6 (v6) assembly of AM560-2 started afresh from deep (120×) Illumina short-read sequence assembled de novo into contigs. These contigs were ordered and oriented into scaffolds with Illumina mate pairs and fosmid ends, and then to chromosome scale using genetic maps (Bredeson et al. 2016). V6 represented an advance over v5 by capturing 18% more total contig sequence, increasing overall sequence contiguity, and incorporating 45% more sequence into chromosomal scaffolds. It has served as the cassava community’s primary reference platform for genomic analyses since its release (Bredeson et al. 2016; Kuon et al. 2019; Wolfe et al. 2019; Kayondo et al. 2018; Nzuki et al. 2017; Ramu et al. 2017; Andrade et al. 2019).

The version 7 (v7) reference assembly was released on Phytozome in 2019. This latest assembly incorporated PacBio continuous long read (CLR) sequences, most longer than 10 kb. Long reads more effectively traverse genomic repeats, allowing them to be captured in assembled contigs more completely. In contrast, the assembly of short reads, e.g. Illumina, often leads to the fragmentation or collapse of repetitive sequences, as their sequence similarity with many other loci in the genome renders them refractory to assembly. Raw CLR data suffer from an order-of-magnitude higher error rate than Illumina, but modern computational methods have largely overcome this issue (Chin et al. 2013). The Canu assembler (Koren et al. 2017) was used to construct de novo contigs from error-corrected CLR data 4 kb or longer, constituting 34× depth coverage, which were scaffolded using an approach similar to that of v6. The use of long reads facilitated an increase of over 100 Mb of included sequence, and over 25× increase in contiguity. We note that while the current 669 Mbp reference assembly is somewhat shorter than the 750 Mbp haploid genome size estimated for other genotypes by flow cytometry (Kuon et al. 2019), it is close to the size estimated directly from shotgun sequencing. Any missing sequences are likely to be confined to highly repetitive regions, since genic regions are essentially completely assembled within v7. The v7 assembly is available at Phytozome and Cassavabase, as discussed below.

Version 8 (v8), currently being assembled at UC Berkeley, will combine High-throughput Chromatin Conformation Capture (Hi-C) data with the contigs and long-range linking information used to assemble v7. Hi-C captures DNA-DNA contact pairs between folded DNA sequences within the cell nucleus, and two closely linked loci along the same linear chromosome are more likely to be in contact than two distant loci (Lieberman-Aiden et al. 2009). Thus, Hi-C data can be used to place DNA sequences relative to one another during genome assembly (Burton et al. 2013; Putnam et al. 2016; Dudchenko et al. 2017). For cassava, Hi-C is useful for incorporating additional contigs/scaffolds into chromosomes, and refining contig/scaffold orders and orientations, especially in areas of low recombination rate and for small scaffolds, where the genetic map does not have much resolution (Fig. 3).

Researchers who have investigations underway using a given version of a genome assembly may feel frustrated to hear that a new one is being released. Must all the analyses be redone using the new assembly? For many purposes (e.g., QTL mapping, or RNAseq analysis) the cassava v6 reference is satisfactory. There are, however, genomic regions that are fragmented or even absent in the earlier assemblies, so it is in general worth examining loci of interest in the newer assembly. It’s important to note that a given gene may be annotated or named differently in different genome assembly versions due to improvements in the completeness and correctness of the underlying assembly, the complement of expression- and homology-based evidence supporting it, and the computational algorithms used to evaluate and integrate that evidence for modeling that gene.

Additional cassava genome assemblies

While the AM560-2 reference genome is a starting point for cassava genomics, increasingly the sequencing of additional cassava diversity is beginning to provide a broader perspective on cassava genetic variation. Wang and colleagues reported the assembly of two cassava varieties, W14 and KU50 (Wang et al. 2014). The former was reported to be M. esculenta ssp. flabellifolia (Pohl) Cif and is referred to as a “wild cassava” or “wild ancestor”. Subsequent analysis using the sequence reads deposited by Wang et al., however, showed that the plant sequenced as W14 does not belong to M. esculenta, but is more closely related to M. glaziovii (Bredeson et al. 2016). Including this accession in a phylogenetic analysis derived from a broader selection of Manihot could help determine whether it is an M. glaziovii or a member of another Manihot species. KU50 is an improved cassava variety popular in Southeast Asia (Ceballos et al. 2020). The W14 and KU50 assemblies were constructed from Illumina whole genome shotgun short-insert and mate-pair libraries, 454 sequencing of BAC clones, and long-range linkages from the AM560-2 reference genome. The W14 and KU50 assemblies are important contributions to our understanding of Manihot genomic diversity but are less complete than the subsequently released AM560-2 reference v6 and v7. These later sequences are therefore recommended as the current best references for mapping short read datasets for genotyping purposes.

The advent of long-read sequencing technologies makes it possible to tackle outbred genomes and, in principle, to separately assemble each haplotype with the aim of recovering all 2n = 36 chromosomes of a diploid cassava. Important first steps in this direction were taken by Kuon and collaborators (Kuon et al. 2019). These authors sequenced the African varieties TME3 and 60444, selected for their resistance and susceptibility, respectively, to cassava mosaic disease (CMD). The genomes were assembled with the Canu v. 1.4 assembler, IrisView, and HiRise software using a combination of PacBio CLR reads (70× coverage), Bionano optical mapping, and Hi-C technologies. The resulting chromosomal scaffolds were deposited in NCBI (Table 2).

Table 2

Accessing key cassava genomic resources

Resource	Portal	url
AM560-2 v6.1	Phytozome	https://phytozome-next.jgi.doe.gov/info/Mesculenta_v6_1
AM560-2 v7.1	Phytozome	https://phytozome-next.jgi.doe.gov/info/Mesculenta_v7_1
W14 assembly (unresolved Manihot)	NCBI	https://www.ncbi.nlm.nih.gov/assembly/GCA_000737105.1
KU50 assembly	NCBI	https://www.ncbi.nlm.nih.gov/assembly/GCA_000737115.1
TME3 assembly	NCBI	https://www.ncbi.nlm.nih.gov/assembly/GCA_003957995.1
60444 assembly	NCBI	https://www.ncbi.nlm.nih.gov/assembly/GCA_003957885.1
HapMapI SNV calls (download)	Phytozome	https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Mesculenta, under Mesculenta → v6.1 → diversity
Genomic and phenotypic data and breeding metadata	Cassavabase	https://cassavabase.org/
ICGMC genetic map viewer	Cassavabase	https://cassavabase.org/cview/map.pl?map_id=2
HapMapII sequence and SNP data	Cassavabase	ftp://ftp.cassavabase.org/HapMapII/
HapMapII SNPs (browser)	Cassavabase	https://cassavabase.org/jbrowse_cassavabase/?data=data%2Fjson%2Fhapmap_variants
Bart Lab Cassava Atlas	none	http://shiny.danforthcenter.org/cassava_atlas/

A major challenge for “diploid” genome assembly of heterozygous genomes is separating haplotypes whose sequences can be very similar. The ultimate goal of “diploid” genome assembly is the generation of complete phased sequences for both haplotypes, i.e., all 2N chromosomes of a diploid. For TME3, Kuon et al. make considerable progress towards this goal, and report the assembly of a “primary” haplotype totaling 732 Mbp of contig sequence, close to the genome size estimated by flow cytometry. They also assemble over 200 Mb of alternative haplotype sequence, providing direct access to allelic variation in this important cultivar, with similar results obtained for 60444. To attempt to link these contigs into phased chromosomal haplotypes, Kuon et al. generated high coverage optical maps and HiC chromatin data for both accessions, a first for cassava. Kuon et al. describe two assemblies derived from these data. The first, a “primary” assembly deposited in NCBI (Table 2), combined contigs with optical map and chromatin data to obtain a composite assembly that mixed both primary contigs and alternate haplotype contigs. This assembly contains 300 Mb of “gap” sequence represented by strings of “N”s, where the length of the gap is estimated using the optical map. Comparison of the deposited assemblies from Kuon et al. with the AM560-2 reference shows that some gaps appear to arise from the artifactual interaction of optical map and HiC scaffolding (Bredeson, unpublished). For example, haplotypes may be juxtaposed, with gaps from the bionano assembly of each haplotype retained. So, a sequence that looks like A1-B1-C1-D1 and A2-B2-C2-D2 in the two haplotypes may be represented as A1-NN-C1-D2-NN-B2 in the primary assembly, creating both artifactual gaps as well a local mis-orderings. Thus, in many cases the “gaps” in this sequence are not truly missing but may be misplaced nearby.

To remove these artifacts, Kuon et al. performed additional processing to remove residual duplicate sequences and gaps using tools such as Purge Haplotigs (Roach et al. 2018). The resulting final assemblies represent mosaics of the two haplotypes of each accession, rather than complete haplotype-resolved diploid assemblies, but still capture the genetic structure of these two important cultivars and appear to capture additional repetitive sequence not found in the AM560-2 reference. We note, however, that the final pseudomolecules described by Kuon et al. total 796 Mb for TME3 and 854 Mb for 60444, are longer than the estimated 750 Mb from flow cytometry (Kuon et al. 2019), presumably due to undetected residual duplication. These assemblies supported the analysis of gene families that had differentially expanded and diverged among TME3, 60444, and AM560, including genes involved in antiviral and other disease response processes. Developing methods for accurately and efficiently performing diploid assembly is an area of ongoing and rapid development in bioinformatics, and we expect that future “diploid” assemblies of cassava and other outbred crops will become commonplace in the next few years. These assemblies will capture both haplotypes at high contiguity and with phase accuracy (i.e., even distant segments of a given maternal or paternal chromosome will be on the same reported pseudomolecule sequence). These will be particularly useful for detailed examination of specific loci (e.g., disease resistance, quality traits) and the genetic variation at such loci across germplasm, to guide both breeding and engineering of improved cassava, as well as the exploration of larger structural variation across cassava and related species.

As more and more genomic variation is discovered in cassava, it will be a challenge for users to take advantage of these resources. For example, a diploid genome sequence is not ideal for mapping resequencing data from other accessions, since reads from the target accession will not have a unique locus to map to but can align more-or-less equally well to either haplotype of the diploid genome sequence. We therefore expect that a high-quality haploid reference will continue to be important as an organizing principle for this flood of genome data. Even in genomic regions with structural variants (e.g., insertions, duplications, inversions) not found in the reference, it will be important to locate these on the reference genome and with regard to markers on the genetic map, providing a common coordinate system. Following the lead of the human genome we anticipate the development of collections of alternate haplotypes (as done for example for the highly polymorphic human HLA region) and ultimately a graphical “pan-genome” representation of the collection of diploid genomes that can be directly queried by users. Such tools are still in development but given the genomic sophistication of the cassava research and breeding community we anticipate their rapid adoption.

Resequencing and Cassava Genetic Diversity

In addition to de novo assembled genomes, many cassava varieties have been ‘re-sequenced’ using short read sequencing technologies. Here we use re-sequencing to mean that sequence was generated with a uniformly-random sampling approach from across the genome of a cultivar of interest (whole genome shotgun [WGS]), but not assembled de novo. Such shotgun sequences can be aligned to the reference to learn about the genome of the re-sequenced variety. Comparisons between re-sequenced individuals can characterize family-level relatedness and larger-scale population structure. In the HapMapI study, Bredeson et al. performed whole genome resequencing of 53 M. esculenta, mostly cultivated cassavas, and five other Manihot (Bredeson et al. 2016). It can be challenging to track down reliable information for a given cassava variety. The HapMapI team compiled provenance/pedigree, descriptive, and/or phenotype information, as available, for the accessions they re-sequenced or genotyped. This information is included in Supplemental File 1, along with identifiers for obtaining the raw sequence from the NCBI Sequence Read Archive (SRA). The HapMapII project incorporated the WGS data from HapMapI and added more diversity, bringing the total to 241 total accessions, comprising 203 cultivated cassavas plus hybrids and wild relatives (Ramu et al. 2017). The raw sequence data generated for HapMapII are available for download on Cassavabase, along with the single nucleotide polymorphism (SNP) data for all 241 accessions (See Table 2).

High throughput genotyping technologies

In the cassava genomics era, there has been a shift from using individual or sparse markers, such as simple sequence repeats (SSRs), to generating and using sets of markers derived from SNPs, which are denser in the genome. SNP-based approaches have included EST-derived marker genotyping arrays, such as KBiosciences KASPar arrays and a 1536 SNP Illumina GoldenGate array (Ferguson et al. 2012a, b; Rabbi et al. 2012; Ferguson et al. 2019), and reduced representation approaches, such as genotyping-by-sequencing (GBS) (Elshire et al. 2011; Hamblin and Rabbi 2014). GBS has been applied broadly in the cassava genomics community (Rabbi et al. 2014a, b; ICGMC 2014; Bredeson et al. 2016; Wolfe et al. 2016a, b; Rabbi et al. 2017; Kayondo et al. 2018; Esuma et al. 2016). Today, researchers can perform GBS in-house or use a fee-for-service GBS facility; DArTseq is a similar technology offered as fee-for-service (Diversity Arrays Technology Pty, Canberra, Australia). Both GBS and DArTseq rely on the same principle of using restriction enzyme digestion paired with Illumina sequencing to sample a reduced―yet broadly distributed and repeatable―fraction of the genome. SNPs are identified, and genotypes called, from the resulting sequence reads. DArTseq service includes this analysis, returning reports with genotype calls for each individual across a set of loci. Several computational approaches for analyzing GBS data have been used in cassava, e.g. TASSEL/TASSEL-GBS (Bradbury et al. 2007; Rabbi et al. 2014a, b; Glaubitz et al. 2014; Ramu et al. 2017), and the pipeline implemented by the International Cassava Genetic Map Consortium (ICGMC 2014).

The HapMapI project used GBS to genotype 268 African cassavas, and the SRA identifiers are listed in the paper’s Supplementary File 1 (Bredeson et al. 2016). These data were demultiplexed (sorted by sample) and the GBS-facilitating sequences removed, so that the deposited sequence reads contain only the cassava sequence.

Genetic maps

Cassava researchers have been making genetic maps for over 20 years (Fregene et al. 1997). Since 2012, a series of papers have demonstrated the utility of SNP markers for genetic mapping in cassava, especially when combined with phenotype data to map important traits (Soto et al. 2015; ICGMC 2014; Rabbi et al. 2014a, b; Rabbi et al. 2012). The International Cassava Genetic Map Consortium used GBS data from 10 mapping populations to construct a dense composite reference map, comprising 22,403 genetic markers organized into the expected 18 linkage groups (ICGMC 2014). This composite map has ongoing utility, including for improving or evaluating new genome assemblies (ICGMC 2014; Kuon et al. 2019; Bredeson et al. 2016). From version 5 onward, the chromosomes of the AM560-2 reference assemblies have been numbered according to the composite map, e.g. chromosome 7 corresponds to linkage group VII.

Omics datasets and data access

Many transcriptomics datasets have been published for cassava. Collectively, these datasets describe gene expression patterns across abiotic and biotic stress (Rauwane et al. 2018; Muñoz-Bodnar et al. 2014; Amuge et al. 2017; Zhang et al. 2015; Ding et al. 2016; Vanderschuren et al. 2014; Cohn et al. 2014, 2016) and across different cassava tissue types (Wilson et al. 2017). These transcriptomics datasets have also fueled a significant number of manuscripts characterizing diverse classes of proteins (Hu et al. 2016a, b; Liao et al. 2016; Wei et al. 2016; Ou et al. 2018; Li et al. 2017; Fan et al. 2016; Hu et al. 2015). As of November 17th, 2019, there were 463 Manihot esculenta transcriptomic samples available through the SRA run selector tool; a subset of these are organized into 23 BioProjects. A high-quality reference genome is also valuable for proteomics research. For example, Vanderschuren et. al. analyzed proteomic changes that occur during post physiological deterioration (PPD) (Vanderschuren et al. 2014), Wang et al. reported proteomics analysis of cassava root formation (Wang et al. 2016) and multiple groups have characterized the epigenetic profile of cassava by mapping bisulfite sequencing (BS-seq) data onto the cassava reference genomes (Wang et al. 2015; Zou et al. 2017). Application of these diverse omics datasets is described below in the “Applications of genomics” section.

Making genomics data accessible and easily visualized to researchers is challenging. Genomics data are famously large and this presents challenges for storage and processing. In addition, most platforms, once constructed, are rigid and can only be used for the intended analysis types. Currently, the cassava community has several points of entry for accessing genomics data. First, the National Center for Biotechnology Information (NCBI) hosts genomics data from all organisms and provides a suite of tools for query and analysis (described within the handbook here: (The NCBI Handbook 2013)). Phytozome is a plant specific public repository for comparative genomics funded through the Department of Energy (https://phytozome-next.jgi.doe.gov/) and has previously been reviewed (Goodstein et al. 2012). M. esculenta references v6.1 and v7.1 and their annotations are available in Phytozome as browsers and for download. The v6.1 browser also contains tracks for single nucleotide variants (SNVs) and indels identified in the 60 Manihot analyzed in the HapMapI project (Bredeson et al. 2016), and methylation and RNAseq data from accession TME7 (Wang et al. 2015). See Table 2 for how to access these data in Phytozome.

In addition to the traditional NCBI and Phytozome portals, new platforms dedicated to the cassava scientific community have been developed. Cassavabase is a genomic-based breeding database. As such, it contains comprehensive management and analysis tools for genomic, genotypic, phenotypic (cassava trait ontology) data and breeding metadata. The system implements a “digital ecosystem”, which allows all breeding processes to be tracked in the database and all breeding data to be collected digitally. Cassavabase manages accessions, pedigrees and seed lots, calculation of field layout randomizations, collection of phenotypic data using tablets and other methods such as drones, collection of crossing data, as well as DNA samples for genotyping. Breeders use it to design their field layouts, barcode the field using built-in barcode tools, collect phenotypic and genotypic data from the plants, and plan and track their crossing experiments. Cassavabase can be used to implement genomic selection (GS) via the SolGS tool (Tecle et al. 2014): phenotypic and genotypic data can be correlated and phenotypic predictions generated on lines that have been genotyped, but not phenotyped. By leveraging previous phenotyping rather than conducting it each time, these genomic predictions accelerate the breeding cycle by a significant factor, and thus increase genetic gain per unit time. Cassavabase currently contains over 30 years’ worth of experimental data, across 250 locations in a dozen countries, including 258 phenotypic traits. Researchers use Cassavabase to explore or analyze various genomic resources by using tools such as BLAST or by identifying new trait-associated genome regions using tools such as GWAS. Cassavabase contains user friendly applications and visualization tools for genomic resources, such as a pedigree viewer, the ICGMC genetic map viewer, and a browser showing the SNPs identified by the HapMapII project (see Table 2). High (GBS) and low density (ie: Kompetitive allele specific PCR, KASP) genotyping assay data are available for more than 30,000 accessions. Users who wish to submit newly re-sequenced resources information may do so using the “Sequencing Status” feature on each germplasm page. All aforementioned functionalities are publicly available and can be accessed by creating a free user account.

For genomics analyses and visualization other than those described above, most scientists rely on lab specific tools and while these methods are reported during publication, this system does not facilitate subsequent queries of the data by future researchers, especially those without a strong computational skillset. The makers of Rstudio (Copyright 2017 Rstudio Inc.) developed a way to turn a now commonly used, open-source data processing language, R (R Core Team 2014), into interactive web applications. They call this development R Shiny (Chang et al. 2017)). This development makes it easy for intermediate and expert R users to create web-ready applications for data queries and visualization. Shiny, just like all R packages, is a large collection of tools and functions that can be pieced together to create a fully functioning website. Wilson et al. used this platform to create a web-based interface, Bart Lab Cassava Atlas, for exploring tissue specific gene expression patterns for cassava (Wilson et al. 2017). Recently, they added data from a previous paper (Cohn et al. 2014), to this tool and in the future, it may be possible to further expand this platform with the additional datasets described above.

Applications of genomics

When a reference genome is published, the authors may draw attention to specific aspects of the genome or use the newly available resource for a specific purpose. In the case of the cassava reference genome, the first scientific manuscript describing the genome assembly accompanied v4.1 (Prochnik et al. 2012). This manuscript communicated the details of this resource to the growing cassava community and was quickly incorporated into additional research efforts. Many papers have used the reference genome to investigate a specific class of proteins. For example, Lozano et al. (2015) used cassava gene annotations from v4.1 (Prochnik et al. 2012) to catalog the NBS-LRR disease resistance genes, and a subsequent chromosome-scale assembly, v5.0, (ICGMC 2014) to show that this class of genes are arranged in clusters in the genome. Similarly, WRKY transcription factors, potassium ion (KUP) transporters (Ou et al. 2018), calcium-dependent protein kinases (Hu et al. 2016a, b), bZIP and ERF transcription factors (Hu et al. 2016a, b; Fan et al. 2016), to name a few, have been the foci of recent papers. This resource has also been used for gene discovery. For example, Cohn and Bart et al. used the reference genome to identify cassava genes upregulated during attack by a bacterial pathogen (Cohn et al. 2014). In another transcriptomics focused study, Amuge et al. investigated genes that were differentially expressed between susceptible and tolerant cassava varieties during attack by Ugandan Cassava Brown Streak Virus (Amuge et al. 2017). The M. esculenta genome is paleotetraploid, i.e., is the result of an ancient genome duplication (Bredeson et al. 2016). Thus, in conducting these types of analyses, researchers should be aware that the cassava genome may contain two copies (paralogs) of a gene of interest. For example, paralogs CYP79D1 and CYP79D2 both code for the enzyme CYP79D, a member of the cyanogenic glucoside biosynthetic pathway (Andersen et al. 2000).

A major goal of the HapMapI project was to characterize global cassava diversity. The authors leveraged their reference genome assembly v6.1 as the platform for cataloging SNVs in 60 resequenced Manihot accessions, including diverse cassava varieties, putative progenitor species, wild relatives; and 268 African cassavas (Bredeson et al. 2016). Major findings included introgressed M. glaziovii segments in many cassava accessions, particularly in improved varieties, and widespread relatedness among cultivated cassavas. In evaluating relatedness among the cassavas via identity by descent, they showed that some cassavas considered to be different accessions are clones of one another, notably TME3, TME14, and TME7; and TME204 and TME419. Here, this means there is no more parsimonious way to explain their high degree of relatedness (sharing > 98% of variable sites relative to the reference) than that they are related via clonal propagation. Measured differences must be due to genotyping error or somatic mutations acquired over time since cloning. Phenotypic differences between clonally related accessions, then, may be due to somatic mutation or epigenetic changes.

The HapMapII project used reference assembly v6.1 with the WGS data from 241 accessions to reveal that cultivated cassava varieties harbor a substantial burden of deleterious mutations (genetic load) maintained in the heterozygous state (Ramu et al. 2017). Subsequent QTL and GWAS studies showed that loci associated with resistance to cassava brown streak disease (CBSD) overlap with regions of M. glaziovii introgression (Nzuki et al. 2017; Kayondo et al. 2018). Rabbi et al. showed that a region associated with both dry matter content and carotenoid content coincides with an introgressed region on chromosome 1, which also has elevated linkage disequilibrium (LD), indicating reduced recombination frequency (Rabbi et al. 2017). Leveraging the WGS data from HapMapI/II, and GBS data from several cassava breeding programs, Wolfe and colleagues identified introgressed M. glaziovii segments in over 2700 cassavas, with American and African populations showing on average 2% to 4.2% of M. glaziovii introgression, genome wide. A reduced level of recombination was observed in the introgressed regions on chromosomes 1 and 4. When heterozygous, the introgressed segments are associated with preferred traits such as CBSD resistance and higher dry matter content, but over progressive rounds of genomic selection homozygotes are selected against. Carrying these haplotypes in the heterozygous state, thus also plays a role in the maintenance of genetic load (Wolfe et al. 2019).

Genomics tools have also helped narrow the genomic location of important loci, such as those that confer resistance to Cassava Mosaic Disease, and the major loci for root carotenoid and dry matter content (Wolfe, Rabbi, et al. 2016a, b; Rabbi et al. a, b, 2017). In addition to gene discovery, high quality genomics resources are invaluable for new genome editing technologies. The design of CRISPR-Cas9 guide sequences, for example, relies on a good quality genome assembly and accurate annotation of gene models. Where available, accession-specific re-sequencing data may be used to ensure that the target sequence matches the guide and to predict potential off-targets. In a proof of principle study for genome editing in cassava, Odipio and colleagues (Odipio et al. 2017) disabled Phytoene desaturase, resulting in the expected albino phenotype (Odipio et al. 2017). This technology is also being deployed to create agriculturally important phenotypes. For example, Gomez and Lin et al. showed that disabling nCBP-1 and nCBP-2 conferred resistance to CBSD (Gomez et al. 2019). A team from ETH Zurich generated waxy (amylose-free) cassava by mutating the GBSS gene, and reduced amylose cassava by editing PTST1 (Bull et al. 2018). In the first report of a knock-in at an endogenous cassava locus, Hummel and colleagues used genome editing to generate glyphosate tolerant plants by swapping the promoter and changing two amino acids of the EPSPS gene (Hummel et al. 2018).

Notably, the rate of published research articles that use existing cassava -omics datasets is increasing overtime, demonstrating the need for such resources and reflecting both an increase in the number of researchers interested in cassava and an increased comfort with using these resources among the cassava community. Cassava has become a model for how to bring an orphan crop into the genomics era. We look forward to further democratization of both the generation and the usage of cassava genomic resources.

Future perspectives

The term “finished” is hardly ever applied to genome assemblies, and rightly so. With the exception of the simplest viral and bacterial genomes, even the highest quality genome sequences are incomplete and a constant focus for further improvement. With this in mind, the cassava community eagerly awaits new genome assemblies including the release of the AM560-2 v8 assembly. This assembly should be the most complete and accurate version available, particularly in centromeric and distal subtelomeric regions.

In addition to increased contiguity of genomic sequence, we can look forward to fully resolved haplotypes for heterozygous genomes such as cassava. Advancements in sequencing technologies have reduced the costs and improved the capacity for sequencing and assembling genomes. As more non-model plants enter the genomic era, it is becoming evident that the utility of haploid genome assemblies is limited in some cases. This is especially problematic in plants with high heterozygosity and strong inbreeding depression, where the generation of homozygous lines for sequencing may be impossible and with limited relevance to the natural state of these genomes. These haploid assembly strategies collapse both haplotypes into a single sequence representing a pseudo-haploid reference (Fig. 2b). During this process of collapsing the haplotypes, assembly errors may be introduced, including the generation of variants which may not be present in either haplotype (Church et al. 2015).

In an effort to overcome these challenges, recent advancements in genome-phasing algorithms have been developed. The trio binning (Koren et al. 2018) approach is the current gold standard for genome phasing. While effective and accurate, this approach relies on the ability to sequence a trio, including two parents and their progeny. This is prohibitive in some cases, such as sequencing farmer preferred cassava lines, where the parents are often unknown. Furthermore, due to the need to sequence long reads from the sample of interest, as well as short reads from the parental lines, this method is more resource input intensive. Using the information from the sequencing of the parents, this approach bins maternal and paternal sequence reads of the progeny prior to their assembly. Thus, each haplotype is then assembled independently and accurately.

In lieu of parental material, other approaches have been developed for phased assembly. The FALCON/FALCON-Unzip and Hifiasm tools attempt to separate phases during the assembly process by identifying and extracting regions containing polymorphisms and structural variants (Chin et al. 2016; Cheng et al. 2020). The output of the FALCON/FALCON-Unzip algorithms is an assembled primary contig alongside “haplotigs” representing a corresponding heterozygous assembly for a specific region. These algorithms, however, have difficulty in resolving regions of high haplotypic divergence (> 4%) (e.g. Padgitt-Cobb et al. 2019) and thus assemble both haplotypes in the primary assembly, similarly to that described in (Kuon et al. 2019). Other orthogonal tools have been developed to facilitate the haplotig purging process (Roach et al. 2018; Guan et al. 2020) and can help resolve the haplotypes. However, this de novo approach does not fully resolve phase switching errors, in which sequences from different haplotypes are assembled adjacent to one another. By utilizing information from long range genomic sequencing (Hi-C), a new algorithm (FALCON-Phase, (Kronenberg et al. 2018)) attempts to resolve these errors and can be applied iteratively on both contigs and scaffolds, thus allowing resolution of chromosome level phases. The end result is two assembly files which represent an approximate depiction of the two haplotypes. This phasing approach has been successfully demonstrated for the human genome, as well as several heterozygous animal species and an Arabidopsis F1 hybrid (Kronenberg et al. 2018). It thus might also be applicable for genome assembly of plant species with high heterozygosity, such as cassava.

Structural variations are likely important to cassava genomic studies and breeding but are not easily detected using short read resequencing. With the availability of long read technologies such as PacBio and BioNano, it may now be possible to survey cassava structural variations by sequencing and assembling many accessions, as is being done in other crops such as tomato (Alonge et al. 2020). Lack of concordance between physical positions of two haplotypes is caused by large structural variations and this can complicate downstream analyses. For example, genes might not directly align between the haplotypes. This can also cause difficulty when comparing mapping results from GBS for example, as SNPs will not land in exact physical positions on both phases. New browsers need to be designed to allow users to explore this type of structural variation (Garrison et al. 2018). The ability to resolve this information is crucial for accurate understanding of genetics in these species and further identification of important genes and variants which may have been masked in a collapsed haploid assembly. Cassava’s high genetic load and clonal propagation have maintained the highly heterozygous nature of its genome for generations. While the technologies and approaches described in this review have provided valuable new insights into cassava structural genomics, in the near term, researchers may be best served by combining these genomics approaches with more standard techniques such as sequence capture and sequencing of bacterial artificial chromosome (BAC) clones to fully resolve particularly important genomic locations. In conclusion, by employing old and new genome phasing strategies we are poised to make valuable gains in understanding of cassava evolution, domestication and breeding.