Genetic complexity of miscanthus cell wall composition and biomass quality for biofuels

Miscanthus sinensis is a high yielding perennial grass species with great potential as a bioenergy feedstock. One of the challenges that currently impedes commercial cellulosic biofuel production is the technical difficulty to efficiently convert lignocellulosic biomass into biofuel. The development of feedstocks with better biomass quality will improve conversion efficiency and the sustainability of the value-chain. Progress in the genetic improvement of biomass quality may be substantially expedited by the development of genetic markers associated to quality traits, which can be used in a marker-assisted selection program. To this end, a mapping population was developed by crossing two parents of contrasting cell wall composition. The performance of 182 F1 offspring individuals along with the parents was evaluated in a field trial with a randomized block design with three replicates. Plants were phenotyped for cell wall composition and conversion efficiency characters in the second and third growth season after establishment. A new SNP-based genetic map for M. sinensis was built using a genotyping-by-sequencing (GBS) approach, which resulted in 464 short-sequence uniparental markers that formed 16 linkage groups in the male map and 17 linkage groups in the female map. A total of 86 QTLs for a variety of biomass quality characteristics were identified, 20 of which were detected in both growth seasons. Twenty QTLs were directly associated to different conversion efficiency characters. Marker sequences were aligned to the sorghum reference genome to facilitate cross-species comparisons. Analyses revealed that for some traits previously identified QTLs in sorghum occurred in homologous regions on the same chromosome. In this work we report for the first time the genetic mapping of cell wall composition and bioconversion traits in the bioenergy crop miscanthus. These results are a first step towards the development of marker-assisted selection programs in miscanthus to improve biomass quality and facilitate its use as feedstock for biofuel production.


Background
Miscanthus is a perennial C4 grass capable of producing high biomass yields in temperate climates [1]. It is a crop characterized by high resource-use efficiency owing to its early spring emergence and long vegetative phase, as well as its rhizomatous growing habit, which allows the recycling of nutrients between growing seasons [2][3][4].
These characteristics make miscanthus an interesting lignocellulose feedstock for the production of cellulosic biofuels [5].
So far, M. × giganteus is the only species of the genus Miscanthus that is commercially exploited for biomass production [6,7]. M. × giganteus (2n = 3x = 57) is derived from a natural cross between the diploid M. sinensis (2n = 2 × = 38) and the Japanese allotetraploid species M. ogiformis (2n = 4 × = 76), which is often erroneously referred to as tetraploid M. sacchariflorus [8,9]. Its success is mainly due to its high productivity. In a quantitative review of biomass yields of M. × giganteus across 100 diverse field trial locations, the average dry matter yield was 22 t ha −1 yr −1 [10]. However, the genetic variation in this triploid clone is extremely limited due to its sterility, which poses risks upon large-scale cultivation and significantly limits further progress through breeding [6,9,[11][12][13]. In contrast, great and largely untapped genetic diversity is harboured within and among natural populations of M. sacchariflorus and M. sinensis, which have adapted to a wide range of geographical conditions [6,13].
One of the key challenges that currently impedes the wide-scale commercialization of cellulosic ethanol production resides is our inability to efficiently deconstruct plant lignocellulose into fermentable sugars. The development of feedstocks with improved biomass quality is envisioned to contribute to the economic feasibility of cellulosic biofuel technologies [5,[14][15][16]. Lignocellulosic feedstocks are composed of cellulose, hemicellulosic polysaccharides and lignin [17]. High contents of cellulose and hemicellulosic polysaccharides are desirable, as these constituents can be hydrolyzed and subsequently fermented to produce biofuels. Lignin, on the other hand, cross-links to hemicellulosic polysaccharides and forms a highly impermeable and complex matrix that shields cell wall polysaccharides from degradation, and impedes the extraction of fermentable sugars from the cell wall [18][19][20][21]. Genotypic variation in cell wall composition has been reported in M. sinensis and M. sacchariflorus, providing ample scope for improving biomass quality in these species through breeding [22,23].
Compared to annual crops, progress in breeding of perennials, such as miscanthus, is slowed-down by the need to evaluate genotype performance in multi-year field trials. Miscanthus typically matures in 3 years and selection at a premature stage, specifically during its first year of establishment, has proven unreliable [24]. Therefore, the application of marker-assisted selection could substantially increase the efficiency of breeding in miscanthus, as selections could be done at the seedling stage using marker data. Genetic maps form the basis for finding marker-trait associations, but their construction in miscanthus is complicated by the large genome size and the high levels of heterozygosity inherent to its obligate outcrossing nature [11,13]. Nonetheless, a few genetic maps of miscanthus have been published to date [25][26][27][28][29].
So far three of these genetic maps have been used for the identification of quantitative trait loci (QTL) for different traits of interest, but none of these studies focused on biomass quality for biofuel production. The randomized amplified polymorphic DNA (RAPD) marker-based map by Atienza et al. has been used for identification of QTL associated with agronomic performance and combustion quality [30][31][32][33][34]. The simple-sequence repeat (SSR) marker-based map by Swaminathan et al. was used for the identification of QTL associated with agronomic performance [35]. This map was recently extended with simple nucleotide polymorphism (SNP) markers, obtained through restriction siteassociated DNA (RAD) sequencing, and was used for the identification of QTL underlying the zebra stripe phenotype that is desirable for the use of miscanthus as an ornamental grass [29]. Currently, no markertrait associations have been reported in miscanthus for traits relating to cell wall composition or biomass quality for the production of cellulosic biofuels.
Here we report the construction of a new genetic map for M. sinensis using SNP markers obtained through a genotyping-by-sequencing (GBS) approach. The mapping population used in this study segregates for biomass quality traits, as it was derived from a cross between two parental lines with contrasting cell wall composition. The objectives of this study were (1) to detect QTL for biomass composition and quality in miscanthus regarding its use as a lignocellulose feedstock for biofuel production and (2) to align marker sequences to the sorghum reference genome to facilitate crossspecies comparisons.

Mapping population
A mapping population of 182 F1 progeny was generated by crossing two M. sinensis genotypes with contrasting cell wall composition. The male parent, hereafter referred to as P1, was a genotype (H0227) originating from the miscanthus collection of Wageningen University and Research (WUR). The female parent, hereafter referred to as P2, was derived from a cross between two genotypes from the BIOMIS mapping population (H0012 × H0163) described by Atienza et al., [25]. Both H0012 and H0163 (grandparents) were also included in the field trial and are hereafter referred to as G1-P2 (H0012) and G2-P2 (H0163), respectively. A random sample of seeds was sown in August 2011 in trays in a heated greenhouse; seedlings were subsequently potted and raised to vigorous plants by the end of the winter of 2011/2012. These were split by the end of May 2012 into four roughly equally sized clonal pieces (ramets). Three ramets of each genotype were immediately used to establish a field trial in May 2012; one spare ramet per genotype was potted to replace possible fall-outs. The trial was located at an experimental site of WUR at Wageningen (The Netherlands) and had a randomized block design with the individual ramets used as experimental units. The ramets were planted in rows with a distance between and within rows of 75 cm. The trial was surrounded by two rows of medium-sized M. sinensis plants in order to minimize possible border effects. In the second and third growth season, heading date was scored per plant. At the end of the second and third growth season (December 2013 and 2014), all plants were harvested separately, dried to constant weight using ventilated air (dm%~92%) and weighed. A random sample of each plant was subsequently taken, from which leaves and inflorescences were separated from the stem material. The stem fraction of each plant was then chopped into~2 cm chips, and air-dried at 60°C for 72 h in a forced-air oven. Stem samples (n = 186 genotypes × 3 replicates × 2 years = 1104, minus fall-outs) were ground using a hammer mill with a 1-mm screen and used for biomass quality analyses.

Biomass quality analysis
Neutral detergent fiber (NDF) and acid detergent fiber (ADF) contents of stem dry matter were determined by detergent fiber analysis using an ANKOM 2000 Fiber Analyzer (ANKOM Technology Corporation, Fairpoint, NY). Acid detergent lignin (ADL) contents were determined after 3-h hydrolysis of the ADF residue in 72% H 2 SO 4 with continuous shaking. All analyses were performed in triplicate and fiber fractions were expressed in gram per kg dry matter. Fiber fractions were used to calculate the concentrations (in g/kg dm) of cell wall (NDF), cellulose (CEL, equals ADF -ADL), hemicellulosic polysaccharides (HEM, equals NDF -ADF) and acid detergent lignin (ADL) on a dry matter basis. The residual NDF material of the replicated fiber analyses was pooled per sample and used for the determination of neutral sugar and Klason lignin (KL) content as described previously [36]. Briefly, 30 mg of NDF material was hydrolysed for 1 h in 0.3 ml 72% H 2 SO 4 at 30°C, after which the acid concentration was diluted to 4% and samples were autoclaved for 60 min at 121°C. Autoclaved samples were cooled and centrifuged, after which the supernatant was used for determination of glucose (GLU), xylose (XYL) and arabinose (ARA) contents using high performance anion exchange chromatography (HPAEC) on a Dionex system (Dionex, Sunnydale, CA) equipped with a CarboPac PA1 column and a pulsed amperometric detector. The pellet remaining after centrifugation was vacuum-filtered through a pre-weighed glass fibre filter (AP25, Fischer Scientific, Loughborough, UK). The residue was dried overnight at 103°C and weighed for the determination of KL.
Separate analyses of ground stem samples were performed for the characterization of saccharification efficiency by two different methods. The first method was used for the high-throughput, small-scale quantification of the rate of glucose release during enzymatic hydrolysis of hot-water pretreated samples, as described previously [37]. The release of glucose was expressed as the concentration in nmol of reducing sugars released per mg of biomass per hour of digestion; hereafter referred to as saccharification rate (SacR). The second method was aimed at quantifying the final yield of fermentable sugars using a highly controlled lab-scale alkaline pretreatment and enzymatic saccharification setup, as described by van der Weijde et al. [36]. The released amounts of glucose and xylose are expressed either (1) as a weight percentage of the amount of glucose and xylose present in the untreated sample as determined by neutral sugar analysis (i.e., referred to as glucose conversion (GC) or xylose conversion (XC)) or (2) as a weight percentage of the amount of cellulose and hemicellulose present in the untreated sample as determined by fiber analysis (i.e., referred to as cellulose conversion (CC) and hemicellulose conversion (HC)).
To allow high-throughput analysis of all biomass quality traits we used near-infrared spectroscopy (NIRS) technology. Multivariate prediction models that combined near-infrared (NIR) spectral data and biochemical data were developed for all traits except for SacR. Nearinfrared absorbance spectra of stem samples were obtained using a Foss DS2500 near-infrared spectrometer (Foss, Hillerød, Denmark) and processed by weighted multiplicative scatter correction and mathematical derivatization and smoothing treatments (2,6,4,1) using WinISI 4.9 statistical software (Foss, Hillerød, Denmark). Different prediction models were developed for different traits, depending on the number of samples that could be biochemically analyzed and on the availability of existing data for creating robust prediction models (containing a range of miscanthus samples from different experiments) ( Table 1). All models contained at least 140 samples from the first growing season of the mapping population. The quality of the prediction models was validated using the squared Pearson coefficient of correlation (r 2 ) between predicted and biochemical data and by evaluating for these samples the standard error of cross-validation (SECV) for each of the traits (Table 1). Subsequently, the developed prediction models were used to determine biomass composition and conversion efficiency of all 1104 stem samples (minus fall-outs).

Genotyping-by-sequencing
Genomic DNA from young leaf tissues was extracted following a CTAB based protocol [38]. DNA concentration and quality were checked using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and standardized using a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA). DNA integrity was confirmed on 1% agarose gels. Libraries were prepared for GBS using the restriction endonuclease ApeKI (five-cutter) to digest the genomic DNA for complexity reduction. Each digested DNA sample was ligated to a set of uniquely barcoded sequencing adaptor pairs, following PCR amplification with adapter-specific primers, and amplicons between 300 and 500 bp were extracted from an agarose gel and sequenced in four single lanes of Illumina HiSeq2000 using a 100 bp paired-end protocol. DNA digestion, adapter ligation, library construction, and sequencing were carried out by the Beijing Genomics Institute (BGI), China.
The de-multiplexed sequence reads obtained from BGI were filtered by removing those reads that did not start with the 5'-CWCG-3' site pattern, typically resulting from ApeKI digestion, or that contained undefined ('N') nucleotides. Reads were right-trimmed to a length of 82 nucleotides and clustered in order to count the number copies per unique read sequence. Note that this clustering was not only done for each sample individually, but also separately for the forward and reverse reads. Only unique reads that occurred at least four times were kept. Unique reads from all samples were jointly clustered using the RADSNP program (RADNPGTv1.1 package, BGI, China). Our initial approach to classify genotypes was to assign a genotypic score to the studied genotypes with a cluster size of at least five reads by applying a set of classification rules to separate clustered reads. The first classification rule was that if the genotype had a frequency of 0.8 or higher for the most abundant read in the cluster, this was considered to be present in homozygous condition. The second classification rule was applied when the two most abundant reads in a cluster if both had frequencies of at least 0.2. The genotype was then classified to be heterozygous. If for a particular cluster neither rule 2 nor 3 held true, no genotypic assignment was given. Unfortunately this approach did not result in acceptable data for map construction, because the average cluster size was too small to allow for a proper genotypic classification due to insufficient sequencing depth. Therefore we refrained from this approach and focused on segregation analyses for single reads. The number of reads for each selected sequence was in this case the basis for genotypic classification using a dominant way of scoring. Genotypes with one or more reads were considered to be either homozygous dominant or heterozygous for this short-sequence marker, whereas the ones showing no reads were supposed to be homozygous recessives. A missing value was assigned to genotype-marker combinations when both the number of reads for this marker over genotypes as well as the average number reads over all markers for a genotype was low. This was done to prevent misclassification of genotypes.

Map construction
A genetic map was constructed following the two-way pseudo test-cross strategy [39], using the dominantly scored SNP markers. To this end, suitable markers were first filtered out of all available markers (49102) based on segregation ratio, with only uniparental single-dose markers, i.e., markers that segregated in a 1:1 ratio in the population, used for further analysis. A total of 1145 markers remained and were coded according to segregation type following the coding scheme for cross pollinated populations as used in JoinMap [40]. Male simplex × female nulliplex markers were classified as lm × ll, while male nulliplex × female simplex markers

Statistical analysis and QTL mapping
General analyses of variance (ANOVA) were performed to determine the significance of genotype differences (p < 0.05) in the mapping population for cell wall composition and saccharification efficiency. Variance analyses were performed separately for both growing seasons, taking into account the randomized complete block design of the trial. Estimates of genotypic (σ g 2 ) and residual (σ e 2 ) variance were used to calculate broad sense heritability (h 2 ) estimates following h 2 = σ g 2 /(σ g 2 + σ e 2 ). To visualize associations amongst traits, a principal component analysis was performed on genotype means for all traits evaluated in both growth seasons. Origin centered, normalized scores for the first two principal components were plotted in a principal-component biplot. All statistical analyses were performed using Genstat for Windows, 18th edition software package (VSN International, Hemel Hempstead, UK).
Quantitative trait loci (QTL) analysis was performed with MapQTL 6.0 (Kyazma, Wageningen, Netherlands) using a maximum likelihood mixture model. An interval mapping approach was used with a step size of 1.0 cM. Significance of a QTL was called based on a LOD score higher than a genome-wide significance threshold based on 1000 permutations [41], which was determined to be 3.561 for the male and 3.655 for the female map. One-LOD and two-LOD support intervals were determined to show the uncertainty on the QTL position. The percentage of variance explained (PVE) by the QTL was calculated by 100 × ([residual variance with no QTL fittedresidual variance with QTL fitted]/population variance) [42].

Results and discussion
Genotypic variation for biomass quality traits Significant heritable variation was observed in the mapping population for all stem biomass quality traits determined after the second and third growth seasons, as shown by the population statistics and parental and grand-parental values summarized in Table 2. Cell wall material (NDF) comprised the largest fraction of biomass and ranged from~815 to 911 g/kg dm in the second and from~877 to 918 g/kg dm in the third growth season. The main cell wall components were CEL and HEM, with variation in the population in the second growth season ranging from~446-527 and~304-365 g/ kg dm, respectively. In the third growth season plants had on average higher CEL and lower HEM contents compared to the second growth season and ranged from 474-532 and~282-345 g/kg dm, respectively. Particularly large variation in cell wall glucose content (GLU) was also found, ranging from~35 to 50% of the cell wall fraction in the second and from~21 to 39% in the third growth season.
Variation in ADL ranged from~42 to 82 g/kg dm in the second and from~75 to 110 g/kg dm in the third growth season. ADL/cw and KL ranged from~5-9% and~12-15% of the cell wall in the second and from 8-12% and 12-16% in the third growth season, respectively. Variation in lignin content is of particular interest for improving biomass quality of miscanthus, and variation in both ADL and KL was extensive. KL values are higher than ADL values, as during the quantification of ADL, detergents are used that likely dissolve a fraction of the total lignin. However, KL values might overestimate lignin as it is more likely to be contaminated with protein [43,44]. Both methods provide different but valuable insights into biomass quality [36].
The mapping population also harbored extensive variation in conversion efficiency. Particularly for SacR, GC and XC, considerable variation was observed among genotypes. Variation in SacR ranged from~11 to 24 nmol reducing sugars per mg biomass per hour. Variation in GC ranged from~42 to 55% in the second and from 33 to 46% in the third growth season. These ranges are comparable with the ranges observed in other highly diverse sets of miscanthus genotypes [23,36,[45][46][47], indicating that variation in conversion efficiency in this population created by crossing two highly compositionally distinct parents is substantial. Conversion efficiency values in the third growth season were substantially lower than those found in the second, which is presumably associated with the increase in lignin content observed with increasing plant age ( Table 2).
Genotype performance for most of the evaluated traits was highly reproducible across replicated blocks. As a result, for most traits a high heritability (h 2 > 0.5) was observed, with the highest heritability for quality traits across years observed for lignin (ADL/cw) (h 2 = 0.62-0.72). For all traits, heritability estimates in the third growth season were reduced compared to those in the second. The lower heritabilities in the third growth season can be caused by the lower range of observed genetic variation, environmental effects, errors in biochemical analyses, and/or lower  [48][49][50]. Frequency distributions of all traits evaluated in the third growth season were reasonably uniform and showed continuous unimodal histograms (Fig. 1). For all traits, with the exception of CEL, parental and grandparental performance were contrasting and for most traits population variation extended beyond parental and grand-parental values in both directions. For KL and GLU, the performance of P1 was very near the low-end population extreme; hence genetic variation leading to concentrations lower than those observed for P1 for these traits is not expected in this population.
Principal components analysis revealed that approximately 58% of the observed genotypic variation in biomass quality resolved into two composite variables (Fig. 2). The first principal component summarized 32% of the observed genotypic variation and predominantly discriminated genotypes based on differences in the content of cellulosic and hemicellulosic polysaccharides. The second component, which summarized 26% of observed variation, discriminated genotypes mostly based on differences in lignin and conversion efficiency characters. As the angle between vectors is representative of correlations between traits, from this plot it can be deduced that the different conversion efficiency characters are positively associated to each other and negatively associated with lignin. It is also evident that SacR was more strongly correlated with the content of cellulosic polysaccharides than was cellulose conversion. These trait associations are consistent with other reports on miscanthus biomass composition and quality for biofuel production [36,45,47].

Synteny with Sorghum bicolor and coding of linkage groups
The DNA sequences of the mapped markers were aligned to the Sorghum bicolor (L.) Moench genome (version 'sbi1') from the Plant Genome Database using NCBI BLASTN [51]. Only hits with an identity score greater than 85% and an alignment length of at least 50 nucleotides were retained and used to label the miscanthus linkage groups according to which sorghum chromosome the markers in each linkage group mapped (Additional file 1: Figure S1). Linkage groups of the female map were designated by the corresponding Sorghum bicolor chromosome numbers, followed by an 'a' or 'b' , as well as linkage groups of the male map, but followed by a 'c' or 'd'. These suffixes were randomly appointed to the two homologous miscanthus linkage groups of each map that are syntenic to each sorghum chromosome, as the genome of M. sinensis consists of two sub-genomes with a high level of synteny to the sorghum genome [26,27]. In both, the male and the female map, there was one linkage group that aligned with two sorghum chromosomes; these groups were designated '4b7b' and '4d7d'. The occurrence of this phenomenon in miscanthus has been reported previously and is ascribed to an ancient chromosome fusion or translocation event between two miscanthus chromosomes syntenic to sorghum chromosomes 4 and 7. This event explains why miscanthus has a basic chromosome number of 19 and not 20 (twice the basic chromosome number of sorghum) [27,28].

QTL mapping of miscanthus biomass quality traits
QTL analysis was performed to investigate associations between genomic regions and stem composition and conversion traits. In a combined QTL analysis carried out on the male and female map simultaneously a total of 86 QTLs were found to be associated with cell wall composition and conversion efficiency characters with LOD scores ranging from 3.58 to 9.02 (Table 3).
Heterozygosity was uncovered in 58 loci of the male parent and 28 loci of the female parent, but these may be partly the same loci if the male and the female map would be combined. Twenty out of 86 QTLs were found in both growth seasons. In the combined analysis, significant QTLs were located across 21 out of the total of 33 male and female linkage groups (Fig. 3). For several traits, QTLs were observed to be present in roughly the same genomic position in presumably homeologous linkage groups in both parental maps (e.g., QTLs for ARA on groups 2c and 2d).
Out of the 86 QTLs that were observed, 9 were associated with stem cell wall, 5 with cellulose, 6 with hemicellulosic polysaccharides, 22 with lignin and 23 with neutral sugar contents (Table 3).. The QTLs on the male map for ARA were numerous; this high number was likely due to factors such as the inherent imprecision of the QTL mapping approach, marker multicollinearity and the height of the LOD score thresholds used. The large number of QTLs found to be associated with lignin content could be partly explained by the fact that three different lignin characters (ADL, ADL/cw and KL) were evaluated. Notably, QTLs associated with KL did not colocalize with QTLs for ADL or ADL/cw (Fig. 3). Two major-effect QTLs were identified for CEL/cw (CEL/cw 3 and CEL/cw 4) in linkage groups of the male parent, each respectively accounting for 29% and~17% of the observed genotypic variation during the third growth season (Table 3). These may be interesting targets for further study.
A total of 20 QTLs were found for conversion efficiency characters with LOD-scores ranging from 3.75 to 7.00, among which 7 for SacR, 4 for cellulose conversion, 2 for hemicellulose conversion and 7 for glucose conversion (Table 3).. QTLs for SacR and GC colocalized on linkage group 3c, QTLs for SacR and HC co-localized on linkage groups 6b and 6c (potentially homologous groups) and QTLs for CC and GC colocalized in linkage group 8c. However, many QTLs for the different conversion characters did not co-localize and seem to be independently controlled characters (Fig. 3). On linkage groups 1b, 1c, 3c, 4d7d, 6c, 6d and 8c QTLs for conversion efficiency characters colocalized with QTLs for lignin characters. Particularly strong evidence for co-localization of QTLs for these traits was found on linkage group 1b and 8c, where QTLs for lignin (KL, ADL and ADL/cw) and conversion characters (CC and GC) co-localized in both growth seasons. On linkage groups 4b7b, 4d7d, 6b, 6c and 8c QTLs for conversion efficiency characters co-localized with QTLs for accumulation of hemicellulosic polysaccharides. A big clustering of co-localized QTLs were observed on linkage groups 6b and 6c, possibly indicating   the presence of a master-regulator affecting cell wall biosynthesis. QTLs for the same traits co-localized in both clusters, suggesting that 6b and 6c are homologous linkage groups. Several QTLs for conversion efficiency characters did not co-localize with any of the QTLs for compositional characters evaluated in this study (e.g., SacR on 3b and 10b), suggesting that other, unidentified compositional characters are affecting conversion efficiency. One such character, for example, could be the content of hydroxycinnamic acids, such as para-coumaric or ferulic acids, which were recently identified as key factors affecting the conversion efficiency of miscanthus biomass [47].

Comparative analysis of QTL in miscanthus and sorghum
In addition to identifying QTLs for miscanthus biomass composition and conversion characters, an objective of this study was to demonstrate that by aligning the genetic map of miscanthus to the physical map of Sorghum bicolor, the exchange of information from genetic studies across species is facilitated and a wealth of information becomes available for the genetic improvement of miscanthus. For this particular objective, the heading date of the genotypes used in this study was scored in both growth seasons, as this is a trait that normally has a high heritability in miscanthus and was previously mapped in miscanthus [35]. Due to the high level of synteny between miscanthus and sorghum, QTL found in one species might have corresponding QTL in homologous regions in the other. In this study, 3 QTLs were identified for heading date, located on linkage groups 1c, 3c and 6d (Table 3). A QTL for heading date on the linkage group that aligns with Sb03 was also identified by Gifford et al., [35] on the same position at the end of the chromosome arm (position 6-9 cM) as HD2 in this study. In addition, a QTL for heading date was consistently reported in sorghum on the end of the chromosome arm of Sb06 [50,52,53], which is in accordance with HD3 found in this study. Similarly, QTLs for NDF are reported in sorghum on chromosomes Sb02, Sb03, Sb04 and Sb06 [50,54], which may correspond to QTLs for NDF found in this study on the corresponding linkage groups 2c, 3a, 3c, 4c and 6c (Table 3, Fig. 3). The QTL on chromosome Sb03 was reported to have a strong effect and explained a large fraction of the observed variation in a sorghum mapping population [50]. The strong effect of this QTL in sorghum may explain why the presumably corresponding QTL was detected on both the female and the male map in both growth seasons (NDF2 on linkage group 3a and LG linkage group, PVE percentage variance explained