Sugarcane (Saccharum spp.) has a complex genome because of its variable ploidy level, frequent aneuploidy and large genome size of approximately 10 Gb [1–5]. This crop is a member of the Poaceae family and the Andropogoneae tribe, which includes maize and sorghum [6, 7]. Modern sugarcane cultivars are the result of interspecific crosses between the domesticated species Saccharum officinarum L. (2n = 80) and the wild species S. spontaneum L. (2n = 40–120), followed by several backcrosses with S. officinarum [6, 8]. These cultivars have chromosome numbers ranging from 100 to 130, are vegetatively propagated, and result from the selection of populations derived from outcrossing heterozygous parents [1, 9]. Sugarcane has a very high photosynthetic efficiency and is a crop with major economic importance in many tropical and subtropical countries primarily because of its use in the production of sugar and bioethanol [10–12].

Polyploidy, an important driver of plant evolution in natural populations, has played a crucial role in the domestication of crops such as wheat, sugarcane, cotton and potato [13–16]. Sugarcane is predominantly an autopolyploid plant, and the understanding of its genome organization is limited [4, 7]. One possible way to increase knowledge of the genome organization of this species is by using genetic maps. High-resolution genetic linkage mapping may be used for quantitative trait loci (QTL) studies in mapping poulations and also such as a first step toward potential marker-assisted selection (MAS) in plants [17–22]. Several genetic linkage maps of sugarcane have been generated since a methodology based on single-dose markers (SDMs) was proposed by Wu et al. [23]. SDMs that segregate 1:1 and 3:1 in full-sib progenies (F₁ populations) [24] or 3:1 in populations created by selfing an individual are commonly used for constructing genetic maps in sugarcane [4, 25–35]. An integrated map of sugarcane with different types of molecular markers, such as microsatellites or single sequence repeats (SSRs) and target region amplification polymorphism (TRAP), extended the characterization of polymorphic variation throughout the entire genome [36–38]. However, in outcrossing heterozygous species such as sugarcane, for each segregating loci, different numbers of segregating alleles can exist, and a relative large number of markers is required to guarantee reasonable coverage of its genome [37, 39, 40].

Currently, high-throughput next-generation sequencing (NGS) technologies have provided new opportunities for discovering molecular markers, especially single nucleotide polymorphisms (SNPs), at the genome-wide level [41–43]. Some of these techniques, e.g., restriction-site associated DNA sequencing (RAD-seq) [44] and genotyping-by-sequencing (GBS) [42, 45, 46], employ a reduced genome representation that is achieved through restriction enzyme digestion, which could be particularly helpful for a complex genome such as that of a polyploid [47]. Moreover, these strategies can be used to species without reference sequence [48]. The GBS protocol has been widely used in a range of genetic studies in several species such as apple, barley, lettuce, switchgrass, maize, rice, wheat, and soybean [42, 45, 46, 49–55]. SNP datasets generated from GBS can be analyzed to detect associations between genotypes and phenotypes, perform diversity analyses, and construct genetic maps, among other applications [39, 56–58].

QTL mapping in sugarcane is a promising tool for characterizing the genetic architecture of several yield component traits of interest, such as sucrose yield, cane yield, stalk diameter, stalk height, stalk number, and stalk weight, as well as resistance to diseases, pests and abiotic stresses [10, 59–62]. Sugarcane is a semi-perennial crop with repeated measures data obtained for several harvests and locations, and QTL mapping studies are usually performed in two steps. First, adjusted phenotypic means are obtained; second, these means are searched for associations with molecular markers and/or along genetic maps [31, 32, 63]. Gazaffi et al. [62] proposed a method that considers an integrated genetic map in which QTL mapping is performed based on the advantages of the composite interval mapping (CIM) approach [64]. Briefly, a mapping model with three genetic effects is considered for genome scanning [62]. It is assumed that a QTL may also segregate in different patterns in progeny as a function of its genetic effects and of the linkage phase between markers and QTL alleles.

This study is the first to report on the development and application of GBS for mapping studies in sugarcane. Our objectives were to (i) establish a pipeline for identifying SNPs and insertion-deletion (indels) from GBS data in a sugarcane F₁ population, (ii) construct the first GBS-based integrated genetic map with additional SSR and TRAP markers in this bi-parental mapping population, (iii) identify QTLs related to yield component traits based on the integrated genetic map, and (iv) perform annotation of the sequences that originated the associated markers with mapped QTLs to search putative candidate genes that may be involved in yield traits in sugarcane. We discuss these results in the context of where GBS is likely to be most useful in sugarcane crop development.

The simultaneous identification and genotyping of SNPs and indels is possible because of important recent advances in sequencing [41–50]. GBS is the preferred high-throughput genotyping method for plants with some level of genetic complexity; this method involves complexity reduction and multiplex sequencing to produce high-quality polymorphism data at a relatively low cost per sample [101]. Using four pseudo-references to discover GBS-based markers, we obtained more markers suitable for linkage analysis (Table 3) than any other previously published study on sugarcane mapping. The strategies adopted for the discovery of GBS-based markers allowed us to relate the sugarcane markers to sorghum chromosomes [73] and to potential genetic regions sampled from the methyl-filtered sugarcane genome [72], RNA-seq sugarcane transcriptome [66] and SUCEST project sequences [65].

The highest alignment (87.94%) of the 3.1 million resulting tags against the methyl-filtered sugarcane genome also contains most of the high-quality SDMs (Table 2). This great alignment rate may be related to the greater amount of scaffolds for this reference compared with the other three references and to the fact that the PstI enzyme used for library formation is sensitive to DNA methylation [102]; thus, more polymorphic sites are expected from the methyl-filtered genome. Additionally, GBS-based markers had more markers mapped to parent SP80-3280 than to parent RB835486 (Table 4). This result can be explained by some factors: a) the possible presence of more PstI enzyme restriction sites in cultivar SP80-3280 than in RB835486, leading to more polymorphic sites in the first cultivar, as shown for barley cultivars by Liu et al. [54]; b) also could be due to more similarity between the genome of SP80-3280 with the reference; or c) due to different methylation effects between the two cultivars. Furthermore, inconsistencies in the number of sites sequenced per sample [102] and in the number of reads per site [103, 104], in addition to the filtering steps applied to the GBS libraries to obtain the markers, can influence the observed result. Other factors that can influence these results are the quality and quantity of the biological replicates used for GBS-based marker calling. Because the sequencing of the samples can present failures that will be included in the downstream process, better dosage and ploidy level estimates for each marker in the SuperMASSA software can be hampered.

The analysis of the loci with high coverage that were filtered for missing data after analysis using SuperMASSA software showed that for the ploidy levels under consideration, the number of loci varied within each ploidy class (Fig. 1), suggesting that the number of chromosomes within the HGs is not constant in sugarcane, as reported previously [39, 105]. As stated by Garcia et al. [39], technique artifacts yielding either strong bias or too much noise should explain marker misclassification, i.e., loci not included in the 6-14 expected ploidy range. In fact, the graphical Bayesian model used in the analyses benefits smaller ploidies due to parsimony when the skewed clusters are confounded or, conversely, favors a higher number of clusters by attempting to explain a diffuse scatterplot [106]. In addition, Garcia et al. [39] hypothesized that poor-quality data can also be generated by biological events such as copy number variations or paralogous regions. We used the same ad hoc criteria as Garcia et al. [39] to classify each locus into one quality category based on the a posteriori probabilities for each ploidy category A as represented by 60.4% of loci of all ploidies on average (Fig. 1), which is smaller than the 77.6% of Sequenom-based data that were previously studied [39]. The GBS read count data worked slightly poorer because of their broad genome coverage and eventual technique artifacts. Despite this reduction, a large part of the loci could be further exploited for mapping purposes.

The presence of repeat elements, paralogs, and incomplete or inaccurate reference genome sequences can create ambiguities in GBS-based marker calling [107]. After we selected the loci that were classified into category A and ploidies ranging from 6 to 14 (Table 2), we continued on to the final steps of quality control and redundancy analyses that showed a low redundancy considering simultaneously all four references. Aitken et al. [35] presented the first sugarcane genetic map with DArT markers and did not remove any redundant markers. Sugarcane has a large and complex genome, and a low level of redundancy is important for showing the true coverage of the genome. Heslot et al. [52] showed that DArT markers were significantly more redundant than GBS markers, and they suggested that GBS markers were significantly more evenly distributed across the wheat genome. These authors also concluded that GBS is the marker platform of choice for further diversity analyses and genomic selection.

The integrated genetic map of sugarcane obtained in this paper presents improvements in comparison with previous works. Here, the number of markers used for linkage analysis is more than twice the number of markers used for the development of the largest map [35]. Furthermore, this genetic map of sugarcane is the first to use a high-throughput approach for genotyping. Co-dominant biallelic markers can segregate in a 1:2:1 fashion (‘B3’ cross type), which is even more informative for map integration purposes. The previously published genetic maps had molecular markers treated as dominant even when are potentially co-dominant, that is, all clusters that have at least one copy of the allele will collapse into a single cluster [37]. The mapping population was formed by a cross between polyploid heterozygous parents, and for each segregating loci, there could be different numbers of segregating alleles and different dosages that are potentially expressed. Thus, accessing the dosage information of the SDMs with a segregation pattern of 1:2:1 was important; the double SDMs with a segregation pattern of 3:1 (‘C8’ type) were also important. This information was used to construct an integrated genetic map for sugarcane that increased the genome coverage.

The 223 LGs obtained here had a cumulative map length of 3,682.04 cM and an average marker density of 3.70. The number of LGs exceeds the number of chromosomes of modern sugarcane cultivars, which can range from 100 to 130 [1, 9], and 56 LGs showed lengths shorter than 2 cM. This result indicates the presence of gaps and that the map is not yet well saturated. In 2007, Oliveira and collaborators [38] claimed that because there is a constraint to discarding markers in multiples doses, i.e., duplexes of monoparental origin, triplex or higher multiplex markers, gaps are evidently expected; the same discussion should be applied to this study. The number of unlinked markers (87.07%) is higher than that obtained in other sugarcane maps [4, 31, 35, 37, 38, 81, 108–110] and reflects the highly stringent criteria used to construct an integrated genetic map of sugarcane that is reliable for performing QTL mapping analysis.

To increase the understanding of the genetic architecture of sugarcane, a necessary requirement is the availability of good genetics maps, i.e., maps with a high density of markers and with high coverage of the genome [111–113]. The complexity of the sugarcane genome, the cost of generating a large number of markers, and the absence of a statistical genetic model that could consider other segregation ratios beyond 1:1, 1:2:1 and 3:1 have limited the development of high-density genetic maps. These limitations have delayed practical applications of genomic tools in sugarcane, in contrast to other crops that have already advanced to MAS and genomic selection. Sugarcane still does not have its genome completely sequenced, and sorghum genome is widely recognized as a reference genome for comparative analysis with sugarcane [67]. The origin of modern sugarcane cultivars raises issues that are related to not only the extent and nature of the divergence of the sugarcane and sorghum genomes but also the relations (in terms of meiosis and dosage) among homo(eo)logous loci [68]. Differences in chromosome structures between the ancestor species and pairing behavior in modern cultivars suggest that the hybrid monoploid number is likely to be greater than 10 in sugarcane hybrids [38, 114]. Probably because of aneuploidy, an unequal number of chromosomes in each HG is likely to occur; this inequality was reflected in the 18 HGs with differences in genome coverage. Moreover, translocation events may have occurred in sugarcane between regions equivalent to sorghum, as discussed previously [27, 28, 115–118]. Although these comparative studies proposed hypotheses about the evolutionary aspects of sugarcane and sorghum, the results showed variations that were primarily derived from the low resolution of the genetic maps used and to the coverage of the sugarcane genome. In addition, it is important to highlight that advances in the assembly of polyploid genomes will enable the use of the full sugarcane genome as a reference in the future [119].

The genetic maps and field data obtained through designed experiments are required for QTL mapping studies. For sugarcane, multiple harvest-location trials may be used to infer the genetic architecture of quantitative traits. However, this inference makes the data analysis more complex and challenging because of the interactions that it generates, e.g., genotype by environment interactions. To solve this problem, a mixed model approach has been used to obtain highly accurate genetic estimates [10, 31, 120], and for segregating populations, these results will be the input for QTL mapping.

In this study, QTL mapping was performed by applying the statistical model proposed by Gazaffi et al. [62], which extends the CIM [64] for a full-sib progeny. The primary advantage of CIM is that it is more precise and effective at mapping QTLs in comparison with single-marker analysis and IM, especially when QTLs are present outside the mapping window [64]. The results obtained from the CIM method are usually comparable to those obtained from multi-QTL analysis if a high-density genetic map is employed to better represent the number of loci underlying the quantitative traits [79]. Several QTL mapping studies in sugarcane have been published [31, 34, 61, 81, 121–136]. The comparison between the results may be biased by several issues, such as the different rates of polymorphisms in parents, the number of progeny, the evaluation methodologies for phenotypic traits, the methodologies used for QTL detection, genetic map contruction, and experimental design, among others. For example, Pastina et al. [31] worked with a population of 100 individuals from a cross between cultivars SP80-180 and SP80-4966 to construct an integrated genetic map that was 2,468.14 cM in length. These researchers used IM to test presence of putative QTLs and a multi-QTL model to declare QTLs and found 46 QTLs. There were 13 mapped QTLs for the tonnes of cane per hectare (TCH), 14 for sugar content in tonnes of sucrose per hectare (TSH), 11 for FIB and eight for POL%C. Singh et al. [34] studied the progeny of 207 individuals derived from a cross between cultivars Co86011 and CoH70, and they constructed two separate genetic maps, one for each parent. Through the CIM model, these researchers found 31 QTLs, with seven for BRIX and four for stalk number (SN). Thus, specific objectives must be taken into consideration for QTL mapping in sugarcane.

For QTL mapping in this study, we performed a permutation test to obtain the threshold for declaring significant QTLs [97]. The CIM model was a useful tool once it was able to identify regions with next QTL considering Araras and Ipaussu over the harvests (three years of evaluations). Seven QTLs were identified, being that a region located in LG4 at 43.32 cM showed QTLs for BRIX (B1-B3) and POL%C (P1-P2). The marker associated with the QTLs was identical for both traits, and the region that gave rise to this marker could be evaluated for future applications by sugarcane breeding programs. POL%C and BRIX are correlated traits [10], and although the commercial cultivars used as parents of the mapping population presented a small contrast in terms of phenotypic averages, especially for sucrose content, these results show that a combination of different alleles in each parent segregates and contributes to the observed variation in progeny. Furthermore, this result was expected because the parents of the mapping population are cultivars that were improved primarily to increase the sucrose content.

The percentage of phenotypic variation explained by each QTL ranged from 2.71% (FIB1) to 9.19% (B1) for Araras and from 5.38% (SD1) to 8.09% (P1) for Ipaussu. The sampling of the genome in single doses requires that QTL also segregate as single doses. Furthermore, the use of improved parents that have close phenotypic averages and that have some level of fixed alleles for traits of interest could decrease the chances of detecting QTLs with high rates of explained phenotypic variation. Nevertheless, the QTLs described here can be considered reliable because they have all taken into account the phenotypic average of three harvests and because it was identified QTLs in same position into LG4 for both locations.

The common QTLs between locations may be regarded as potential regions to search genes that are involved in controlling quantitative traits. In addition, a strong marker-QTL association detected in full-sib progenies also could be an impact on crop improvement via clonal propagation because probability of crossover between the marker and the QTL is low [34]. Therefore, the similarity analysis and annotation of sequences that originated the markers with mapped QTLs are important for identified putative candidate genes in sugarcane, although the QTLs regions are relatively large and an uncertain number of genes may be involved with the evaluated traits. Sugarcane has high genetic complexity and its genome still does not was completely sequenced [67], whereas inferences about putative candidate genes could be contribute for new insights and open new fronts of research to mining and validation of genes of interest.

For the BRIX trait, we can highlight the similarity of the marker SCSFAM1074E10_287, located in QTL B2, with extended synaptotagmin-1-like, which is a member of a membrane-trafficking protein family that is characterized by an N-terminal transmembrane region, a linker of variable size, and two C-terminal C2 domains in tandem [137]. C2 domains, identified as a conserved sequence motif in protein kinase C [138], are autonomously folded protein modules that generally act as calcium (Ca²⁺) and phospholipid-binding domains and that were shown to represent autonomously folded Ca²⁺-binding domains in synaptotagmins [139]. In addition, Ca²⁺ acts as a second messenger in the signal transduction pathways of hormones and environmental stimuli (touch, wind, chilling, light, and elicitors) [140], and several proteins that are involved in photosynthesis depend on Ca²⁺ [141]. In sugarcane, Papini-Terzi et al. [142] identified differentially expressed genes in genotypes contrasting for sucrose content and showed that among these genes, nine were associated with calcium signaling and a calcium-dependent protein kinase (ScCDPK-27). Further researches are needed to expand the inference found in this study for pathways and regulatory networks of sugarcane in order to relate with sucrose biosynthesis and, consequentely, with BRIX trait.

For the FIB trait, a QTL (FIB1) showed different similarities for each of the two adjacent markers and we can highlight the similarity of the marker mf125302_409 with zinc finger protein CONSTANS-LIKE 15. The CONSTANS (CO) protein is a zinc finger transcription factor that contain two conserved domains (i.e., a B-box zinc finger domain and a CCT [CO, CO-like, TOC1] domain), located in the region near the amino- and carboxy-terminus, respectively [143, 144]. The CO proteins play a central role in the photoperiod pathway of Arabidopsis by mediating the circadian clock and floral integrators via positive regulation of FLOWERING LOCUS T expression [145–147]. In cotton (Gossypium spp.), which is the most important natural source of fiber for the textile industry, the CO5 protein (CONSTANS-LIKE 5) was up-regulated for cell wall modification and developing fibers in MD52ne, a near-isogenic line. Furthermore, using an F₂ population derived from a cross between MD52ne and MD90ne, stable QTLs for bundle fiber strength and fiber length were found, and CO5 was present on the QTL region for fiber length [148]. Although there is an indicative of relation between the candidate gene and FIB trait, advanced studies should be conducted to validate and prove the real function and effects into phenotypic expression in sugarcane.

Our understanding of the genetic architecture of triats of interest in sugarcane is increasing with the development of new analytical methods. The estimation of ploidy and allelic dosage through markers generated by GBS as well as the inclusion of these markers in an integrated genetic map of sugarcane were first observed in this study, and these markers showed great potential for QTL mapping. The CIM approach that provided additive and dominance effects and estimated the segregation patterns for all mapped QTLs was efficient for detecting possible stable QTLs among the evaluated locations. The verification of possible candidate genes for mapped QTLs as a preliminary analysis showed importance for new insights into the comprehensive relations between phenotypes and genotypes. It is still necessary to develop statistical approaches to enable the inclusion of markers at multiple doses to enhance the coverage by linking the SDMs that are pulverized by the genome. Moreover, QTL mapping with markers in multiple doses must be considered a major step in the understanding of regions that control quantitative traits in polyploid organisms and perhaps permit the verification of the allelic expression of phenotypic traits in the future.