A mosaic monoploid reference sequence for the highly complex genome of sugarcane

Sugarcane (Saccharum spp.) is a major crop for sugar and bioenergy production. Its highly polyploid, aneuploid, heterozygous, and interspecific genome poses major challenges for producing a reference sequence. We exploited colinearity with sorghum to produce a BAC-based monoploid genome sequence of sugarcane. A minimum tiling path of 4660 sugarcane BAC that best covers the gene-rich part of the sorghum genome was selected based on whole-genome profiling, sequenced, and assembled in a 382-Mb single tiling path of a high-quality sequence. A total of 25,316 protein-coding gene models are predicted, 17% of which display no colinearity with their sorghum orthologs. We show that the two species, S. officinarum and S. spontaneum, involved in modern cultivars differ by their transposable elements and by a few large chromosomal rearrangements, explaining their distinct genome size and distinct basic chromosome numbers while also suggesting that polyploidization arose in both lineages after their divergence.

S ugarcane-a member of the Poaceae family-produces 80% of the world's sugar and has recently become a primary crop for biofuel production 1 . Cultivars are vegetatively propagated, while sugarcane breeding is still essentially focused on conventional methods since sugarcane genetics knowledge has lagged behind that of other major crops 2 . This is due to the interspecific, polyploid, and aneuploid nature of modern sugarcane cultivar genomes, a complexity that exceeds that of most if not all other crops. Modern sugarcane cultivars are derived from a few interspecific hybridizations performed a century ago between S. officinarum and S. spontaneum, two highly polyploid species. S. officinarum (2n = 8× = 80, x = 10) 3,4 has a high sugar content and is believed to have been domesticated around 8000 years ago from the wild species S. robustum (mainly 2n = 60, 80, and up to 200) 5,6 . S. spontaneum is a wild species with various cytotypes and many aneuploid forms (2n = 5× = 40 to 16× = 128; x = 8) 3,7 and was used to incorporate disease resistance, vigor, and adaptability 4 . Recovery of high sugar production was achieved by backcrossing the first hybrid to S. officinarum and then enhanced via 2n chromosome transmission by the female S. officinarum 8,9 . The resulting cultivars are highly heterozygous, aneuploid, and have 100-130 chromosomes, most of which are derived from S. officinarum, 10-20% from S. spontaneum, and 10% from interspecific recombinants 10,11 (Fig. 1).
In the last 20 years, molecular geneticists have achieved significant advances in developing molecular resources and enhancing the overall understanding of the sugarcane genome. Several geneticmaps have been produced (reviewed in ref. 12 ) mainly based on single-dose markers that are the most informative markers in this high polyploidy context 13 . None of these maps were saturated, but they revealed that chromosome assortment in modern cultivars results from general polysomy with some cases of preferential pairing [14][15][16][17][18][19][20][21] , while the meiosis of modern sugarcane cultivars mainly involves bivalent pairing 8,22,23 . Meanwhile, comparative mapping with other Poaceae species revealed extensive genome-wide colinearity with sorghum, which thus became a model for sugarcane 15,[24][25][26][27] . QTL studies have been performed to identify genomic regions involved in important agronomic traits, generally revealing minor effects [28][29][30][31] , except for the presence of a few major resistance genes 20,[32][33][34] . New approaches are now being developed, such as genome-wide association (GWA) studies [35][36][37] and genomic selection 38 , both of which exploit the high level of linkage disequilibrium observed within sugarcane cultivars 39,40 . Genomic resources such as large expressed sequence tag (EST) libraries 41 and BAC libraries 42,43 have been developed. A sugarcane reference genome would greatly help to achieve further progress on molecular genetics and functional genomics geared toward assisting in the development of improved cultivars. However, sequencing such a highly complex genome poses challenges that have not been addressed in any prior sequencing project due to the high polyploidy, aneuploidy, and interspecific structure of this genome, with a complete set of hom(oe)ologs predicted to range from 10 to 12 copies 44,45 . This polyploidy results in a total genome size of about 10 Gb for sugarcane cultivars, while the monoploid genome size is about 800-900 Mb, close to that of sorghum (750 Mb) 46 .
In this context, an attractive initial strategy would be to assemble the gene-rich, recombinationally active part (i.e., euchromatin) of a monoploid (i.e., basic) sugarcane genome. Sorghum and rice genome sequencing showed that euchromatin is gene-rich, repeat-poor, and accounts for most of recombination, while being of similar physical size in sorghum and rice (250-300 Mb) despite their overall genome-size difference (450 Mb for rice and 750 Mb for sorghum) 47 . Comparative analysis among sugarcane hom(oe)ologous chromosome segments (BAC) showed very high gene colinearity and sequence conservation 48,49 . One haplotype sequence would thus provide a good reference for the other haplotypes (=hom(oe)ologous chromosomes). In addition, these studies revealed a high level of microcolinearity between sorghum and sugarcane [47][48][49][50][51] . Sorghum could be used as a guide to identify a minimum tiling path (MTP) of BACs covering the euchromatin of a monoploid genome of sugarcane.
While new sequencing and assembling methods tailored for complex genomes continue to be developed, we deployed this original strategy to produce a first sugarcane reference genome sequence assembly.
This international initiative focused on the R570 cultivar, which has the best-characterized sugarcane genome to date 10,17,18,21,24,44,45 . This cultivar was obtained by CERF (now eRcane) in Réunion (http://www.ercane.re/). It has broad adaptability and is used as a parent stock in many breeding stations worldwide, while also being commercially very successful in Réunion, Mauritius, and Guadeloupe. The genome of this cultivar has about 115 chromosomes, including 10% of whole chromosomes derived from S. spontaneum and 10% from S. officinarum/ S. spontaneum recombinant chromosomes, with the remaining being whole chromosomes from S. officinarum 10 . Genetic maps of this cultivar were constructed 17,18,21 (http://tropgenedb.cirad. fr/tropgene/JSP/index.jsp) and aligned with sorghum 24 . A BAC library of 103,296 clones 43 representing 14× monoploid genome coverage and 1.3× total genome coverage of this cultivar 44 is available and has been used by the research community for comparative structural analyses 26,[47][48][49]52,53 .
Here, we implemented this strategy using WGP TM technology 54 to align sugarcane BACs from the R570 cultivar with the sorghum genome and identify a MTP of 4660 sugarcane BACs originating from distinct hom(oe)ologous chromosomes that we sequenced, assembled, and annotated. We developed a SNPbased genetic map for the R570 cultivar and performed global recombinants Fig. 1 Schematic representation of the genome of a typical modern sugarcane cultivar. Each bar represents a chromosome, in orange or red when originating from S. officinarum or S. spontaneum, respectively. Chromosomes aligned on the same row are hom(oe)ologous chromosomes (HG). The key characteristics of this genome are the high polyploidy, aneuploidy, bispecific origin of the chromosomes, the existence of structural differences between chromosomes of the two origins, and the presence of interspecific chromosome recombinants genome-wise comparisons between S. officinarum, S. spontaneum, and sorghum. This first sugarcane sequence represents the gene-rich part of the monoploid sugarcane genome.

Results
A MTP for the gene-rich part of the genome. The overall strategy that we used to produce a sugarcane monoploid reference sequence assembly is represented in Fig. 2 Table 2). A sorghum orthologous gene could be found for 20,809 (82%) of these sugarcane-predicted genes among the 27,532 sorghum protein-coding gene models predicted on the sorghum genome. Many of the remaining 4507 gene models were short monoexonic predictions with no hits in the Uniprot/Trembl sequence database (http://www.uniprot.org) and no Interpro domains and likely corresponded to overpredictions. Among the 1267 gene models that had at least one Interpro domain, the most frequently found domains corresponded to proteins involved in signaling or defense responses such as kinases, receptor-like kinases, and NB-ARC proteins. Proteins involved in protein degradation by the proteasome were also found, specifically F-box proteins and proteins with an SKP1/BTB/POZ domain (Supplementary Table 3). Many of these proteins are encoded by genes that are known to evolve rapidly 55,56 .
Among a set of 9871 conserved grass genes, 79% (7791) were captured by the STP.
Gene colinearity conservation at the BAC level (i.e.,~100 kb) was analyzed by comparing BACs of the sugarcane STP to the corresponding orthologous regions in sorghum. Among the 20,809 sugarcane-predicted gene models for which an ortholog could be found in the sorghum sequence, 83% were located within the corresponding orthologous sorghum segments, thus leaving 17% of sugarcane-predicted genes in the non-collinear position in sorghum. Since we selected the BACs to be sequenced based on their global alignments with sorghum, this level of noncolinearity could be higher in more rearranged regions between sugarcane and sorghum, or among sugarcane hom(oe)ologs.
Finally, based on a sequence comparison of 13,069 pairs of orthologous genes, we evaluated the sugarcane/sorghum divergence at 8.5 My (median ks =~0.111), which is in the upper range of previous estimations [48][49][50][51] .
Differential transposable element (TE) amplification among Saccharum species. TEs were annotated in the STP and their proportions estimated in WGS data from the R570 cultivar, one S. spontaneum accession and one S. officinarum accession ( Table 2). They represented 43% of the STP as compared to at least 52% for the whole R570 cultivar genome in agreement with the STP focusing on the gene-rich euchromatin fraction of the genome as opposed to the TE-rich heterochromatin. The most represented TE families in both S. officinarum and S. spontaneum are the LTR retrotransposon families Gypsy Chromovirus and Copia Maximus/SIRE. The insertion dates of the complete LTR retrotransposons displayed an L shape, revealing continuous transposition activity over the last few million years    Fig. 1). TE represented a larger part of the S. officinarum genome compared to the S. spontaneum genome (~54 vs.~47%), in agreement with their distinct basic genome sizes of~950 vs.~800 Mbp 46,57 . The larger basic genome of S. officinarum mainly resulted from a more intense amplification of the Gypsy Chromovirus and Copia Maximus/SIRE families in this species' lineage after its divergence from S. spontaneum around 1.5-3.5 Mya 48,53 . The majority of TE families are less abundant in S. spontaneum compared to S. officinarum, with the exception of two Gypsy families, Ogre_TAT and Athila that have been amplified more intensively in S. spontaneum after the divergence of these two species.
Chromosome structural variations among Saccharum species. A SNP-based genetic map was built using 186 individuals from a self-progeny of cultivar R570 genotyped by DArTseq. A total of 12,468 single-dose SNPs were assembled into 132 cosegregation groups (CGs) that included 5-624 DArTseq markers, with an average of 94 markers.
The origin of CGs was then tentatively assigned to ancestral species (S. officinarum vs. S. spontaneum) based on 5377 SNPs differentiating S. officinarum/S. robustum vs. S. spontaneum accessions within a panel of 34 accessions genotyped with DArTseq. Most CGs (74%) were assigned to S. officinarum and 26% were assigned to S. spontaneum or to recombinant chromosomes between these two species (Fig. 5). CGs assigned to S. officinarum were generally much smaller (i.e., contained fewer markers) than CGs assigned to S. spontaneum, as also observed in previous sugarcane-mapping studies 18,58 . This is explained by the fact that the S. spontaneum genome is less redundant than the S. officinarum genome in sugarcane cultivars and thus bears more single-dose markers amenable to genetic mapping.
Synteny conservation between sugarcane CGs and sorghum chromosomes was analyzed using 5406 DArTseq markers from the R570 genetic map that could be aligned with the sorghum sequence. Most sugarcane CGs (125, 95%) had a large majority of markers aligned with one single sorghum chromosome (Fig. 5a, b, Supplementary Fig.2), while for seven CGs a large majority of markers aligned with two distinct sorghum chromosomes (Fig. 5c). Sugarcane CGs were assembled in hom(oe)ology groups (HGs, one HG = one basic chromosome set) based on these syntenic relationships ( Fig. 5 and Supplementary Figs. 2, 3).
Interestingly, all seven CGs that aligned with two sorghum chromosomes were assigned to S. spontaneum chromosomes or to recombinant chromosomes between S. spontaneum and S. officinarum (Fig. 5). These results suggested that a few large chromosome rearrangements differentiate S. spontaneum on one side and S. officinarum and sorghum on the other. In S. spontaneum, one HG was orthologous to sorghum chromosome 2 and one arm of sorghum chromosome 8, whereas another HG was orthologous to sorghum chromosome 9 and the other arm of sorghum chromosome 8. A third S. spontaneum HG was orthologous to sorghum chromosome 6 and one arm of sorghum chromosome 5. The R570 SNP map did not feature any S. spontaneum CG aligned with the second arm of sorghum chromosome 5. However, preliminary SNP mapping data on another sugarcane clone (MQ 76-53, Supplementary Fig. 3) revealed a few CGs orthologous to sorghum chromosome 7 and part of sorghum chromosome 5 (Fig. 5). In addition, Aitken et al. 15 also observed a few CGs in cultivar Q165 orthologous to chromosomes 5 and 7 of sorghum. This suggests that one S. spontaneum basic chromosome could be orthologous to sorghum chromosome 7 and one arm of sorghum chromosome 5. This configuration may have been lacking in the R570 genetic map, probably because our map is far from saturation, and CGs only partially cover the corresponding chromosomes. In addition, one Table 2 Transposable element (TE) contents in the single tiling path (STP) and in whole-genome sequencing data from the R570 cultivar, one S. officinarum (S. OFF) and one S. spontaneum (S. SPONT) R570  large CG assigned to HG08 in cultivar R570 apparently resulted from an interspecific recombination event between S. spontaneum and S. officinarum chromosomes: this recombination could have removed the portion orthologous to sorghum chromosome 7 from the original S. spontaneum chromosome (see the arrow in Fig. 5b). A similar situation was observed for one interspecific recombinant CG assigned to HG07 (see arrow in Fig. 5b). Previous mapping studies have suggested a few large chromosome structural variations between sugarcane and sorghum with some sugarcane CGs having segments orthologous to two sorghum chromosomes 15,24,59 . While our genetic map remains unsaturated, our relatively high number of markers with assignments to S. officinarum vs. S. spontaneum enabled us for the first time to describe chromosome rearrangements involving two sets of three ancestral chromosomes that were each rearranged in two chromosomes (Fig. 5). These rearrangements corroborated and explained a variation in basic chromosome numbers in S. officinarum vs. S. spontaneum, i.e., x = 10 and x = 8, respectively 3,60 . The coexistence of distinct chromosome organizations in the genome of modern cultivars most likely has an impact on their meiosis and may be a source of the aneuploidy observed in modern cultivars. A basic chromosome number of x = 10 is found in both sorghum and S. officinarum and is recognized as being ancestral in Saccharinae-Sorghinae 50 . These rearrangements must have arisen in the S. spontaneum lineage after its divergence from the lineage of S. officinarum and its recognized wild ancestor S. robustum. The mechanisms involved will require further investigation but could imply two successive reciprocal translocations, as reported in Brassicaceae and Poaceae 61,62 . This scenario also suggests that polyploidization occurred independently in S. officinarum/S. robustum and S. spontaneum lineages after their divergence around 1.5 to 3.5 Mya 48,53 in discordance with previous hypotheses 50 .

Discussion
The sequence reported here is the first assembly of the sugarcane genome with 382 Mb of sequence in 3965 high-quality contigs. This was achieved by isolating BAC clones corresponding to a monoploid genome of sugarcane to overcome assembly difficulties linked to the high ploidy, aneuploidy, and heterozygosity of sugarcane cultivars. This sequence assembly is proposed as a reference for the euchromatin gene-rich part of the genome, known to account for most of the recombination, and thus representing the most useful part for breeding. We confirmed the overall gene synteny conservation and colinearity between sugarcane and sorghum, while identifying 17% of genes predicted to be non-collinear in the surveyed regions. We showed that the two species involved in modern cultivars differ for LTR retrotransposon amplification after their divergence, and that the larger basic genome size of S. officinarum compared to S. spontaneum is mainly due to the amplification of two families, i.e., the Gypsy RLG_Chromovirus and Copia RLC_Maximus/SIRE. In addition, we showed that these two species differed by two rearrangements, each involving three sets of ancestral chromosomes, which explains their distinct basic chromosome numbers. Two chromosome organizations thus coexist in the genome of modern cultivars. We propose a chronologic sequence for these rearrangements and polyploidization in the Saccharum genus.
This monoploid reference sequence provides an essential genome template for aligning sequencing data, such as genotyping by sequencing, WGS, and RNA-Seq data to explore hom(oe) ologous allelic variation and perform genetic (e.g., QTL and GWAS) and genomic studies in cultivars and sugarcane germplasm. These studies, as before, will require methods (e.g., SNP calling) adapted to the high polyploid context of sugarcane. This reference will also facilitate the identification of candidate genes/ loci in regions found to be associated with targeted phenotypic traits by permitting the location of these regions within the annotated genome. Finally, it will provide an essential framework to help current whole-genome sequencing initiatives meet the challenge of assembling such a complex genome.

Methods
Sugarcane BAC selection through whole-genome profiling. WGP of 20,736 BACs from the R570 sugarcane BAC library 43 was performed as described by Van Oeveren et al 54 Table 1). Then 22,299 WGP tags that aligned at multiple positions on sorghum but in the same 300-kb window as the anchored BACs were used as additional information to refine BAC positions. Another 577 already-sequenced BAC clones were anchored onto the sorghum genome using BLAST 63 . Finally, among the 12,309 anchored BACs, an MTP of 4660 BAC clones was selected based on visual inspection of their position on the sorghum genome and of BAC clone overlaps.
BAC sequencing and assembly. An initial set of 577 BAC clones were sequenced individually, including 245 BAC clones reported by De Setta et al. 52 , Garsmeur et al. 49 , and Jannoo et al. 48 , and 332 BAC clones sequenced by DOE-JGI using 2 × 250-bp paired-end single-indexed Illumina and assembled with Phrap (Green, http://www.phrap.org/phredphrap/phrap.html). For 4083 BAC clones, DNAs were extracted individually and then pools of BAC DNA were obtained. PACBIO RS II was used for sequencing with a targeted 100× depth per BAC. Half of the BACs were sequenced by KeyGene in pools of 24 clones and the other half by DOE-JGI in pools of 96 clones. Assembly was performed using HGAP3 (version 2.3.0) followed by consensus-sequence calling with Quiver (version 2.1). For each sequenced pool, WGP tags corresponding to the pooled BACs were aligned with the contigs obtained. BLAST alignments of WGP tags with identity and coverage equal to or greater than 99% were used to assign contigs to the corresponding BAC clones, resulting in the assembly of 3958 BAC sequences. Whole BAC sequence contigs were aligned to the sorghum sequence using the BLAST algorithm with the best hits being retained to produce Fig. 4.
To construct STP, contigs (larger than or equal to 20 kb) from the 4535 sequenced BAC clones were aligned with the sorghum-annotated nucleotide gene set (version 3.1) using BLAT (https://github.com/dib-lab/ged-docs/wiki/ BLAT, parameters -noHead -extendThroughN -q=rna -t=dna). Sorghum gene alignments were used to inform the process of reducing the initial MTP in an STP representing a single copy of sugarcane gene space. BAC segments sharing at least two sorghum genes were considered as redundant overlapping regions. These BAC clones were aligned with one another using Gepard (http://cube.univie.ac.at/ gepard). Overlapping BAC regions sharing identical sorghum genes were then collapsed into a single path, thus maximizing the number of retained genes and minimizing sequence loss. Sorghum gene positions were finally used to order, orient, and join the BAC contigs using 10,000 Ns to form 10 super-scaffolds.
Ancestral origin of sugarcane BAC clones. To identify the origin, S. officinarum vs. S. spontaneum, of sugarcane BAC sequences, Illumina WGS data from one S. officinarum accession (LA Purple) and one S. spontaneum accession (SES234B) were aligned separately to BAC sequences using BWA-MEM. The cumulative read coverage of each BAC clone was calculated and compared using SAMtools (min depth = 5 reads, min quality = 20) and bedtools genomecov (-bga). BAC contigs were then classified if there was at least 10% difference in the total coverage of the BAC clone by S. officinarum or S. spontaneum reads. Among the 4535 BAC contigs sequenced, 86% could be assigned with 78% to S. officinarum and 22% to S. spontaneum. The remaining BAC contigs (>50 kb) were analyzed using speciesspecific 24-bp kmers (24 mers) extracted from (i) repetitive regions (read depth > 2 × SD of mean), and (ii) well-covered regions (read depth within 1 SD of the mean) of the previous set of S. officinarum and S. spontaneum BAC contigs assigned to S. officinarum and S. spontaneum. These kmers were used to mask BAC contig sequences, and the ancestral origin was determined based on the majority of sequences masked by the species-specific kmers. Finally, a total of 4281 BAC contigs were classified, representing 360 Mbp (94.4%) out of the 382 Mbp of the STP, with 272 Mb (71%) assigned to S. officinarum and 88 Mb (23%) to S. spontaneum.
R570 RNA-Seq data. RNA was extracted from R570 cultivar leaves, roots, and stems using a method adapted from Bugos et al 64 . Three cDNA libraries (one per tissue) were built using the TruSeqStranded mRNA sample preparation kit (Illumina) and subsequently paired-end sequenced (2 × 125 bp) on an Illumina HiSeq 2500 system. A total of 354 million reads were produced and 63, 54, and 57 million paired-read clusters were built for leaves, roots, and stems, respectively, and assembled using Trinity in 340,409 transcripts larger than 300 bp (https://github. com/trinityrnaseq/trinityrnaseq/wiki).
Sugarcane genes annotation. Structural annotation was performed using EuGene 4.2a (http://eugene.toulouse.inra.fr), which combines evidence from alternative splice site detection, sequence model training, genome masking, transcript mapping, and protein homology to define confident gene models. The R570 cultivar RNA transcripts, together with transcripts from six sugarcane genotypes produced by Cardoso-Silva et al. 65 , were aligned on the sugarcane STP using GMAP (https:// github.com/juliangehring/GMAP-GSNAP/) and best-scoring hits were kept. The resulting alignments were integrated as sugarcane transcript evidence in the prediction process. Protein sequences from the reference sorghum proteome (Phytozome, release 3.1.1) were aligned with the sugarcane STP using BLASTX to obtain protein homology evidence.
Gene function annotation was assigned based on sequence and domain conservation. Protein sequences were aligned with SwissProt, TrEMBL, and NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-05051-5 ARTICLE NATURE COMMUNICATIONS | (2018) 9:2638 | DOI: 10.1038/s41467-018-05051-5 | www.nature.com/naturecommunications sorghum datasets by BLASTP using an E-value threshold of 1e-10. Best-hit BLAST results were then used to define the gene functions. Moreover, InterProScan = (https://www.ebi.ac.uk/interpro/search/sequence-search) was performed to annotate protein domains, extending the annotation to Gene Ontology terms associated with these protein domains. The predicted gene models were finally screened to remove models that could correspond to TEs using CENSOR (https:// www.girinst.org/censor/) 66 resulting in a total of 25,316 sugarcane genes annotated on the STP.
TEs annotation. Sets of 2 million Illumina reads (150 bp) were retrieved from WGS data for the R570 cultivar, one S. officinarum accession (LA Purple) and one S. spontaneum accession (SES 234). A similar dataset was created from the sugarcane STP through generation of in silico Illumina reads with WGSIM (https:// github.com/lh3/wgsim). These four datasets were used to perform a graph-based clustering analysis of repetitive sequences using the RepeatExplorer pipeline 67 . A total of 179 clusters were obtained, manually inspected using the .htlm summary output file generated by the pipeline, and clusters corresponding to chloroplast or containing less than 3000 Illumina reads were not considered for further analyses. The 135 remaining clusters were manually annotated using a library of 58 previously annotated sugarcane elements 52 and using similarity searches against Repbase (https://www.girinst.org/repbase/), a protein domain database of plant mobile elements (included in RepeatExplorer). For each annotated cluster, the WGS read numbers and percentage were calculated in datasets corresponding to STP, R570, S. officinarum, and S. spontaneum (Table 2).
TE annotation of the sugarcane STP was performed with the TEannot pipeline from the REPET package v2.5 (https://urgi.versailles.inra.fr/Tools/REPET), which requires a reference library of TEs. This latter reference library was built through clustering of the following TE sequences: -The 58 previously annotated sugarcane TEs.
The percentage of TE in the STP was calculated based on the sequence proportion covered by each TE type ( Table 2).
The timing of insertion of full-length LTR retrotransposons was determined by estimating the amount of divergence between LTR sequences. The two LTRs of each retrotransposon were aligned using MAFFT (https://mafft.cbrc.jp/alignment/ software/), and the genetic divergence was estimated using the Kimura 2 parameter method. The insertion dates were estimated using a substitution rate of 1.3e-08 per site per year 68 .
DArTseq genotyping. A mapping population encompassing 186 individuals from a self-progeny of cultivar R570 and a panel of 34 accessions belonging to the Saccharum genus were genotyped through DArTseq at Diversity Arrays Technology Pty Ltd., Australia. Two adapters corresponding to two different restriction enzyme overhangs were used; a PstI-compatible adapter, including a barcode region for multiplexing and a reverse NspI-compatible adapter. Only fragments that contained both adapters (PstI-NspI) were used to build sequencing libraries (Illumina, USA). The resulting libraries (two copies per sample) were sequenced in batches of 94 per lane on Illumina Hiseq2500.
Construction of the R570 cultivar SNP-based genetic map. Identification of high-confidence single-dose SNP markers: DArTseq data from the 186 individuals from the self-progeny of R570 were analyzed to identify high-confidence singledose SNP markers useful for genetic mapping. The analysis was based on a de novo approach involving the assembly of a pseudomolecule that was further used as a template for read mapping. This pseudomolecule was assembled by clustering all DArTseq reads using CD-Hit-EST 69 with the following parameters: -c 0.9, -al 0.9, -p 1, and -M 0. Clusters were then randomly assembled to form the pseudomolecule, with each cluster being separated by a stretch of 1000 "N". A pipeline (available at https://github.com/SouthGreenPlatform/VcfHunter/) was developed to perform the following complete process: (1) DArTseq reads were aligned on the pseudomolecule using BWA-MEM (http://bio-bwa.sourceforge.net); (2) a local realignment of reads around indels was performed with the GATK IndelRealigner (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/ org_broadinstitute_gatk_tools_walkers_indels_IndelRealigner.php); (3) metrics were generated at each single-nucleotide position within realignments (count number of A, C, G, and T at each position) with Bam-readcount (https://github. com/genome/bam-readcount); (4) SNP positions on the pseudomolecule and allele count were converted in a standard variant-calling format file (VCF); and (5) the VCF was then screened to determine the genotype (homozygous vs. heterozygous) of each individual and at each position. Several criteria were taken into account to improve the genotype determination accuracy. First, a 30× minimum and 1000× maximum read depth was required and only diallelic polymorphic markers were considered. For heterozygous genotypes, the minor allele then had to be represented by at least two reads with a frequency of at least 4%. Ambiguous genotypes with a minor allele frequency ranging from 1 to 4% were converted as missing data (uncalled genotype). Homozygous genotypes were determined when no variant was observed or with a frequency of <1% (considered as sequencing error). (6) The segregation ratio of single-dose polymorphic markers within the population was then evaluated using a χ 2 test (p = 0.05). SNP markers displaying an expected 3:1 ratio and with less than 20% missing data were considered as high-confidence single-dose SNPs. A total of 13,062 single-dose SNPs were identified.
Construction of the R570 cultivar SNP-based genetic map and alignment with the sorghum genome: The genetic map was built with the 13,062 single-dose SNP markers using JoinMap 4.1 70 with a LOD score threshold of 10 and a maximum recombinant fraction of 0.35. CGs containing at least five markers were kept. Markers were ordered with MSTmap 71 with standard parameters.
DArTseq markers from sugarcane CGs were aligned with the sorghum genome using BLASTN and 5406 DArTSeq markers with hit identity and coverage equal to or greater than 95% were kept. Orthologous relationships between sugarcane CGs and sorghum chromosomes were then established based on the proportion of DArTseq markers from one sugarcane CG aligned with each of the sorghum chromosomes. A minimum of five markers, representing at least 10% of aligned markers for a CG was used. For 125 CGs, a large majority of markers (mean = 88%) aligned with a single sorghum chromosome. For seven CGs, a large majority of markers (mean = 81%) aligned with two sorghum chromosomes. The resulting genetic map with markers aligned with the sorghum chromosomes was constructed using Circos 72  , were aligned on the same pseudomolecule developed for the analysis of DArTseq genetic mapping data. A raw VCF file comprising the SNP positions and allele counts was created and filtered, as described above for identifying high-confidence SNP markers (pipeline, steps 1-5). A SNP was considered as originating from S. officinarum if present in at least two S. officinarum or S. robustum accessions and absent from all S. spontaneum accessions. Reciprocally, a SNP that was present in at least two S. spontaneum accessions and absent from all officinarum and robustum accessions was considered as originating from S. spontaneum. This enabled the determination of the ancestral origin of 5377 SNP markers from the R570 genetic map.
Sugarcane-sorghum comparative analysis. Identification of orthologous gene pairs: Sorghum gene models from the v2.1 DOE-JGI annotation set (http:// phytozome.jgi.doe.gov/) were filtered to remove those corresponding to TE-like genes using CENSOR, as done previously for sugarcane genes. For genes with several predicted alternative transcripts, we only kept the longest one. The resulting 27,532 sorghum genes were aligned with the 25,316 annotated sugarcane genes using BLASTP, and orthologous gene pairs with identity equal to or greater than 70% and coverage equal to or greater than 40% were retained. Since both genomes were annotated with different pipelines, some genes that had been predicted in one genome could have been missed in the other genome. We checked if genes for which no orthologs were found using BLASTP could be retrieved in the genomic sequences. These genes were thus aligned with genomes using TBLASTN and the matching regions were realigned using EXONERATE 73 (model protein2genome: bestfit' and cutoff identity equal to or greater than 70%). A second complementary approach was performed using BLAT.
Microcolinearity of sugarcane and sorghum genes: BAC sequences from the sugarcane STP were aligned with BLASTN on the sorghum genome to identify orthologous genomic regions between the two genomes. Then the gene content was compared for each orthologous region. Sugarcane/sorghum orthologous gene pairs located in the expected orthologous regions were considered as colinear genes, while orthologous gene pairs located in regions that differed from those expected were considered as non-collinear genes.
Divergence timing: Protein sequences from sugarcane gene models were aligned with the sorghum proteome using BLASTP and best bidirectional mutual hits were retained. The 17,704 resulting orthologous gene pairs were realigned using CLUSTALW (https://www.ebi.ac.uk/Tools/msa/clustalw2/). These alignments were used to guide nucleic coding sequence alignments with PAL2NAL (http://www. bork.embl.de/pal2nal/). Synonymous substitution rates (Ks) were calculated using the Nei-Gojobori method implemented in the PAML package. All of these steps were automated with an available Python script (synonymous_calc.py, https:// github.com/tanghaibao/bio-pipeline/tree/master/synonymous_calculation). Divergence time (T) was then estimated using T = MKs/2k, where MKs represents the median Ks (0.1268) calculated from 17,704 orthologous gene pairs and k is the substitution rate of 6.5 × 10 −9 per site and per year 75 .
Code availability. The pipeline developed to perform the complete process of SNP mapping, calling, and filtering is available at https://github.com/ SouthGreenPlatform/VcfHunter/.