Ploidy effect and genetic architecture exploration of stalk traits using DH and its corresponding haploid populations in maize

Doubled haploid (DH) lines produced via in vivo haploid induction have become indispensable in maize research and practical breeding, so it is important to understand traits characteristics in DH and its corresponding haploids which derived from each DH lines. In this study, a DH population derived from Zheng58 × Chang7-2 and a haploid population, were developed, genotyped and evaluated to investigate genetic architecture of eight stalk traits, especially rind penetrometer resistance (RPR) and in vitro dry matter digestion (IVDMD), which affecting maize stalk lodging-resistance and feeding values, respectively. Phenotypic correlation coefficients ranged from 0.38 to 0.69 between the two populations for eight stalk traits. Heritability values of all stalk traits ranged from 0.49 to 0.81 in the DH population, and 0.58 to 0.89 in the haploid population. Quantitative trait loci (QTL) mapping study showed that a total of 47 QTL for all traits accounting for genetic variations ranging from 1.6 to 36.5 % were detected in two populations. One or more QTL sharing common region for each trait were detected between two different ploidy populations. Potential candidate genes predicated from the four QTL support intervals for RPR and IVDMD were indirectly or directly involved with cellulose and lignin biosynthesis, which participated in cell wall formation. The increased expression levels of lignin and cellulose synthesis key genes in the haploid situation illustrated that dosage compensation may account for genome dosage effect in our study. The current investigation extended understanding about the genetic basis of stalk traits and correlations between DH and its haploid populations, which showed consistence and difference between them in phenotype, QTL characters, and gene expression. The higher heritabilities and partly higher QTL detection power were presented in haploid population than in DH population. All of which described above could lay a preliminary foundation for genetic architecture study with haploid population and may benefit selection in haploid-stage to reduce cost in DH breeding.


Background
Maize (Zea mays L.) is one of the important grain and feed crops in which the stalk, as one indispensable part of plant morphology, serves as the conductor of transporting water and nutrients. Stalk lodging lead to yield losses estimated to range from 5 to 20 % annually worldwide [1]. Rind penetrometer resistance (RPR), which is one of the reliable indicators of stalk strength, has been widely used to measure stalk strength and improve stalk lodging resistance [2,3]. Maize is also one of the most important annual forage crops. In vitro dry matter digestion (IVDMD) has been the most useful evaluating indicator for maize forage variety examination in many European countries [4]. Therefore, a further and better understanding of the molecular basis for RPR and IVDMD is crucial for breeding lodging-resistant and highly digestible maize [5].
The genetic analysis of quantitative traits is difficult and complex in maize, and quantitative traits are affected by key genes and interacting networks of small-effect genes. Therefore, different studies have provided different results including quantitative trait loci (QTL) number, distribution, and genetic effects for one trait [6,7]. This lack of conformity may also be explained by the many differences in parental materials, segregation-population types, ecological conditions, genetic maps, analytical methods and phenotype evaluation [8,9]. Moreover, high genome dosage levels have effect on genetic analysis [10,11].
Due to the advantages of time-saving and high genetic variance, doubled haploid (DH) technology is routinely used in modern maize breeding for production of homozygous parental lines for maize hybrid breeding and constructing DH populations for genetic research [12][13][14][15]. Although haploid populations possess the characteristics of genetic homozygosity and have one genome dosage, moderate to strong correlations have been identified between small size DH populations and their haploid version populations for some agronomic traits [16]. Moreover, haploid lines could react more sensitively to biotic and abiotic stresses and, therefore, they would effectively uncover susceptibility to diseases and environmental constraints. In A. thaliana, the utility and power of haploid genetics had been reported. Haploids can provide genetic analysis advantages that are not available in diploids, such as specifically pyramiding multiple mutant combinations, forward mutagenesis screens and swapping of nuclear and cytoplasmic genomes [17]. In yeast, haploid screens represent an ideal platform for negative selection since a certain genetic lesion set by mutagenesis will exert equal effects in all cells [18]. In this regard, the haploid lines may also be interesting in the genetic architecture exploration of maize quantitative traits.
Different segregating populations have been used in linkage analysis or genome-wide association study of RPR, and the genome set number of all these populations was two. The results suggested the genetic complexity of RPR. Flint Garcia et al. [19] first detected 35 RPR QTL in four F 2:3 populations, which accounted for more than 33 % of the total phenotype variation. Hu et al. [20] detected 9 QTL in a RIL population developed from the cross of B73 × Ce3005, which could explain 1.15-12.43 % of the phenotypic variation. Li et al. [21] narrowed the QTL interval which had the largest effect among the 7 QTL of RPR detected in two RIL populations by the method of haplotype analysis. Peiffer et al. [22] reported that 18 family-nested QTL and 141 significant GWAS associations were identified for RPR across NAM (nested association mapping) and IBM (intermated B73 × Mo17) families, while numerous weak associations were found in the NCRPIS (North Central Regional Plant Introduction Station) diversity panel for RPR. Mutations, brittle stalk (BK) genes exhibiting a lower proportion of cellulose, had dramatically weakened tissue mechanical strength than that of wild type stalks [23].
Moreover, whole plant digestibility, which can reflect the feeding value, has been extensively studied in forage maize, and several reports of QTL analyses with lowdensity markers for stalk digestibility in forage maize were published [24,25]. Maize mutants and/or genetically engineered plants have highlighted a few genes affecting maize cell wall degradability [26,27]. Reports have emerged on nucleotide diversity and the extent of linkage disequlibrium (LD) at the gene locus of lignin and cellulose synthesis [28][29][30].
It was well known that plant breeders are desired to choose lines based on minimizing negative effects of genotype agronomic value, so it was crucial to perform research on the genetic architecture of stalk traits, especially for RPR and IVDMD. In this study, we first used a DH population combined with the corresponding haploid population to identify QTL and observe candidate gene expression about stalk traits. Our objectives were to: (1) explore the genetic architecture of stalk traits; (2) evaluate consistence and difference in phenotype, QTL characters, and gene expression between two different ploidy populations in stalk traits; and (3) preliminary propose and illustrate a ploidy effect mechanism for RPR and IVDMD under one genome dosage situation with the QTL mapping method.

Results
Performance of parental lines, F1 generation and DH and haploid populations derived from each DH line Performance of parents and derived DH and haploid populations across five environments was presented in Table 1. RPR, water content (WC), acid detergent fiber (ADF), neutral detergent fiber (NDF), and cellulose(Cel) of the male parent Chang7-2 (C7-2) showed significantly higher values than those of the female parent Zheng58 (Z58) in both DH and haploid populations. In contrast, for IVDMD and WSC (water soluble carbohydrate), Z58 had a higher value than the male parent C7-2 in both populations. There was no significant difference in lignin (Lig) content between two parents in the DH and haploid populations. RPR and IVDMD showed a normal distribution in both two ploidy populations (Fig. 1). For all traits investigated in this study, coefficients of variation (CV) in the DH and haploid population ranged from 7.56 to 49.48 % and from 8.28 to 35.28 %, respectively. The genotypic variance (σ G 2 ) was significant at P < 0.01 in both the DH and haploid populations (

Inter-population and intra-population phenotypic correlation
The phenotypic correlation coefficients of all stalk traits between the DH and haploid populations ranged from 0.38 to 0.69 (Fig. 2). Coefficients of phenotypic correlation among different traits in DH population showed similar patterns to those in haploid population. In both populations, ADF, NDF and Cel showed high positive correlation among themselves, significantly positively correlated with RPR but negatively correlated with IVDMD, Lig and WSC. RPR negatively correlated with IVDMD but with different correlation coefficients in DH and haploid populations, respectively (Table 3).

Constructing a linkage map and the characteristics of markers
A total of 190 DH lines were used for genotyping with MaizeSNP3K chip, which was carried out on the Illumina Golden-Gate SNP genotyping platform [31] and 2956 high-quality SNPs were detected. The missing rate for these SNPs ranged from 0 to 20.00 % (average 1.50 %), the heterozygosity ranged from 0 to 14.21 % (average 2.06 %). A total of 4.74 % (9/190) of the DH lines with SNP heterozygosity ≥ 10 % were excluded in further analysis. Minor allele frequency (MAF) for these SNPs ranged from 0 to 0.50 (average 0.42) (Additional file 1: Table S2). Of these high-quality SNPs, 1318 SNPs were polymorphic between the two parental lines, and the marker distribution frequency for the two parents ranged from 30 to 65 % (Additional file 1: Figure S3). After quality control, 1137 SNPs were left and used to construct a linkage map using the Joinmap4.0 instructions [32]. The total length of the linkage map was 1426.83 cM with an average interval of 1.26 cM (Additional file 1: Table S3).

QTL characteristics in the DH and haploid populations
Across five environments, the number and position of QTL detected in the DH and haploid populations was shown in Fig. 3. For each trait evaluated in this study, one or more QTL were identified in one region or even shared the same support intervals with the distance of less than 20 cM between the DH and haploid populations.
In the haploid population, four QTL for RPR were detected on chromosomes 1, 2, 3 and 5, two of which were identified on chromosomes 1 and 5 using the DH population (Table 4 and Fig. 3). The position of QTL identified in the haploid population on chromosome 1 was close to that detected in the DH population, which accounted for   For IVDMD, three QTL were identified in each population, which explained 8.60-18.50 % of total genetic variation in the DH population and 6.80-18.60 % in haploid population These QTL were detected on chromosomes 1, 2, and 8 in the DH population and on chromosomes 5, 6 and 8 in the haploid population. Two QTL detected on chromosome 8 were tightly linked, which explained the 16.00 and 18.60 % of IVDMD genetic variation, respectively, and both were contributed by the IVDMD-higher parent Z58 in the DH and haploid populations.
In the DH population, QTL of RPR, IVDMD, ADF, NDF and WSC shared the same region ranging from 39.91 cM to 59.43 cM on chromosome 1. The QTL intervals for RPR, ADF, NDF and Cel detected in the

Candidate gene identification for RPR and IVDMD in the DH and haploid populations
With a relatively high mapping resolution, some QTL representing the small genomic regions and the linear B73 genome can be used for searching candidate genes related to RPR and IVDMD. Based on the available annotation of the B73 reference sequence Version 5b.60 (http://ftp.maizesequence.org/release-5b/filtered-set/), we applied the MapMan BIN classification [33] and maizeGDB website (http://www.maizegdb.org/) to search for candidate genes. In this study, four QTL which equally assigned for RPR and IVDMD in the two populations could account for more than 15.00 % of the genetic variation. The number of genes located in the four QTL of RPR and IVDMD were different from each other based on gene screening using qteller3 (http://qteller.com/qteller3/index.php) (Additional file 1: Table S5-S8). Nineteen genes have previously been demonstrated to be associated with cell wall formation mainly involved with cellulose and lignin synthesis (Table 5) [34][35][36][37][38][39][40][41][42][43]. Candidate genes participating same bioprocess were predicated between the DH and haploid populations for RPR and IVDMD, which were consistent with the results proposed by previous studies [20,21,44]. Moreover, although some evidence illustrated  Transcriptional expression analyses of key genes involved in lignin and cellulose synthesis for haploid and diploid version of parental lines To determine whether key genes involved in lignin biosynthesis were ploidy-modulated at a transcriptional level, relative expression levels of five genes, PAL, COMT, ccoAOMT, CCR and CAD, were analyzed (Fig. 4). In parental line Z58, transcript levels of the five genes increased 1.57-5.30folds in haploid plant relative to diploid plants.
Particularly, COMT had the largest change fromZ58 in diploids to Z58 in haploids, while with no significant change from C7-2 in diploids to C7-2 in haploids. In the other parental line C7-2, higher expression levels of haploids than diploids were also observed across five genes, however, the changes of expression level (1.61-2.12-fold) from diploids to haploids were lower than what observed in Z58.
We also examined expression of cellulose synthesis genes encoding glycosyltransferase which detected in both populations for RPR and IVDMD (Fig. 4), and CesA11 and CesA12 were consistently co-expressed at all developmental stages of and were predominantly associated with the deposition of the secondary cell wall in maize stems even after the anthesis stage [53]. For C7-2, the two genes, CesA11 and CesA12, were up-regulated 2.08-fold and 3.04-fold, respectively, in the haploid version relative to the diploid version. For Z58, the expression of CesA12 increased 1.48 fold in haploids in comparison to diploids, however, haploids had lower expression level for CesA11.

Discussion
Performance and heritability of stalk traits in DH and haploid populations The previous genetic investigations on RPR and IVDMD in diploid populations revealed that the traits were likely polygenic in maize and were affected by several mechanisms and complicated by confounding factors.
In this study, two parent lines presented consistent trends on RPR, as well as other traits, between the DH and haploid populations. This result showed each parent contributed coherent negative or positive allele effects even under different genome dosages and additive effect may played important role for partial phenotypic variation.
Heritability estimation depended on genetic background of the material, population types surveyed, interaction with environments and experimental design [54]. In this study, the heritability of RPR was 0.72 and 0.78 estimated in DH and haploid population respectively, Fig. 4 Expression levels of lignin and cellulose synthesis genes in FIAG rind of haploid and diploid parental lines at milky stage. Quantitative RT-PCR analysis for lignin and cellulose synthesis genes was shown in the first five pictures and the last two pictures, respectively (ACTIN as an internal control). Four bars in each picture presented Z58 haploid, Z58 diploid, C7-2 haploid and C7-2 diploid from left to right. Different lowercase indicated that statistical significant difference (P < 0.05). Error bar ± SD which were all in agreement with a high heritability of RPR reported in previous studies [19][20][21]. However, Peiffer et al. [22] reported that the heritability of RPR estimated from 26 RIL families of maize nested association mapping (NAM) population ranged from 0.08 to 0.34 (average 0.21), which may be due to a wide range in the flowering time among NAM families. Likewise, a high heritability was obtained for IVDMD in both DH (0.81) and haploid (0.89) populations, which were different from a moderate heritability, ranging from 0.55 to 0.68, reported in previous studies [55][56][57][58]. Moderate or high heritability values were obtained for other stalk traits as well. In conclusion, the high heritability values of stalk traits evaluated in this study could provide solid basis for QTL mapping analysis.
For most of traits evaluated in this study, the heritability estimated from the haploid population was higher than that from the DH population. Our results were consistent with heritability of DH and its haploid populations in maize reported by Geiger et al. [16], although different traits were studied. The higher heritability estimated from the haploid population can be explained by the relatively smaller σ G × E 2 and σ e 2 in haploid population than in DH population, which can be explained by that haploid lines mainly reacted more sensitively than DH lines to biotic and abiotic stress and therefore effectively uncover susceptibility to diseases and outer constraints, which had been proposed by Chase et al. [59] and Geiger et al. [16]. In addition, all traits evaluated in this study were measured with high precision and then had a solid genetic basis in two ploidy populations, fundamentally suggesting that the haploid population as well as the DH population could be used in QTL analysis.

Phenotypic correlations and QTL co-localization for the same trait between the DH and haploid populations
Geiger et al. [16] reported moderate to high correlations between the DH and haploid lines from three material sets (KWS, SWS, and MON) for early vigor, silking, plant height, and stover weight per plant. We also observed a significant (P < 0.0001) moderate to strong positive correlation (r = 0.38-0.69) between the DH and haploid populations for all stalk traits (Fig. 2). This could suggest that moderate to strong correlations can occur independently of material background and trait restrictions. This high correlation between haploids and corresponding DH lines may provide reference information for maize breeders to select desirable lines at haploid stage, which could reduce breeding costs. However, the genetic mechanism on the connection between the DH and haploid populations has not yet been studied and therefore, is still unclear. In this study, through QTL mapping studies conducted in DH and its haploid population, we intended to understand this issue in term of genetic architecture. We first identified common QTL regions between the DH and haploid populations for each stalk trait, which could be considered as the genetic reason for the phenotypic correlation. Other QTL located on different chromosomes or having larger distance (>20 cM) may be partially caused by the change of genome dosage and explained by the different population size. Ming et al. [60] reported that many QTL for sugar content detected in sugarcane autopolyploids were not consistent with known candidate genes and suggested that other approaches will be necessary to isolate the genetic determinants of high sugar content of vegetative tissues. Until now, QTL detection in haploid population has not been reported.

Phenotypic correlations and QTL co-localizations among different traits
In an attempt to further understand the genetic architecture of RPR and IVDMD in maize, genomic regions for RPR, IVDMD and other stalk component traits were compared and phenotypic correlations between RPR, IVDMD and other stalk component were evaluated. Forty-seven QTL were identified in the DH and haploid populations (Additional file 1: Table S4, Fig. 3 and Table 4). The incidence of QTL clusters in similar genomic regions reflected trait associations [61].
Two studies have proposed that genes associated with the biosynthesis of cell wall components were considered as candidate genes for RPR [19,20]. We also observed the positive correlations and QTL co-location of RPR with ADF, NDF and Cel, which were consistent with previous studies. RPR was negatively correlated with WC in a high-oil RIL population [20]. The same correlation trend of RPR with WC and WSC were observed in the DH and haploid populations, except that RPR had no correlation with WC in the haploid population. In addition, Hu et al. [20] reported that the internode diameter, fresh weight of internode and dry weight of internode were also significantly positively correlated with RPR, and the difference in planting years, densities and maize varieties led to different stalk RPRs [62].
IVDMD showed the opposite correlation direction as the correlations of RPR with ADF, NDF, Cel and WSC, and had the same correlation direction as the correlations of RPR with WC. Therefore, WC may be one of the improved elements for practical breeding for stalk lodging resistance and forage maize. Several QTL associated with IVDMD and other stalk components were located in the same bins as identified in our studied [57,58,63]. Lig was positively correlated with IVDMD and was not correlated with RPR, which was not in agreement with previous studies [20]. This may be due to the no-forage background materials used in this study. No reports were available on QTL both for RPR and IVDMD. Only in the DH population evaluated in this study, we first detected one QTL cluster for RPR and IVDMD at bin 1.10 and bin 1.07, respectively. However, we found more than one QTL of RPR or IVDMD sharing common regions or flanking markers with the QTL of other stalk components, which suggested close linkage or pleiotropy as the explanation for the correlations and some common genes had effects on RPR and IVDMD. The QTL clusters could be deployed for improving RPR and IVDMD in maize through marker-assisted selection.

Compare QTL identified in this study with those identified in previous studies in diploid populations
We have identified additive QTL for RPR on chromosomes 1, 2, 3 and 5. Flint-Garcia et al. [19] detected one QTL region on chromosome 3 contained overlapping support intervals across four F2:3 maize populations. Also, in less than four populations, other QTL were detected at bins 1.07-1.09, 2.02, 2.06-2.07, 3.04-3.08, and 5.02. Similarly, our mapping study of the DH and haploid populations identified five RPR QTL located near bins 1.07, 1.10, 2.02, 3.09 and 5.01. Hu et al. [20] investigated RPR in a RIL population derived from a high-oil population and reported that RPR QTL were detected on all chromosomes except for chromosome 5 and the QTL located in bin 3.06 was the most important one and it accounted for 12 % of the phenotypic variation. Li et al. [21] identified seven RPR-associated QTL in two RIL populations. Among these QTL, the largest-effect QTL accounted for 18.9 % of the phenotypic variation was located at bin 3.06, and other QTL for RPR were observed at bins 2.10, 3.08, 9.03-9.04, 4.06, 6.05, and 6.07, explaining 4.40-13.80 % of the phenotypic variation. In the present study, the QTL location on chromosome 3 were only detected in the haploid populations and accounted for 10.30 % of the genetic variation with highly detected frequency in 1000 runs cross-validation (Additional file 1: Figure S4). Moreover, it is worth noting that RPR-associated QTL, which were observed at bin 2.02 in the haploid populations and explained more than 15.00 % of the contribution to genetic variation, were located in the same region as QTL detected by Flint-Garcia et al. [19]. The QTL detected at bin 5.05 in the DH population could account for the highest percentage of RPR genetic variation (up to 16.90 %) and were not located in the QTL cluster with other traits, and this QTL has not been proposed by previous studies. Since these two newly discovered QTL were also detected with high frequencies in the 1000 crossvalidation, this confirmed our conclusion that QTL at bins 2.02 and 5.05 likely carried major candidate genes for RPR (Additional file 1: Figure S4).
Six QTL for IVDMD in total were detected in DH and its haploid population in this study. Two QTL detected in the DH and haploid populations were located in adjacent bins 8.04 and 8.05 with a genetic distance of less than 3 cM. These two QTL also showed high detection frequencies in cross-validation (Additional file 1: Figure S4). Similarly, Wei et al. [57] reported that a IVDMD QTL located at bin 8.06-8.07 were detected in Pop2 combined analysis, which was adjacent to QTL for IVDMD on chromosome 8 detected in this study. Other QTL for IVDMD identified in the DH and haploid populations were distributed on chromosomes 1, 2, 5 and 6. IVDMD QTL located at bin 1.07 can explain 18.50 % of the genetic variation. Previous reports showed that QTL on chromosome 1 had a great effect on stalk digestibility [57,63]. The QTL located at bins 5.02-5.03 and 5.03-5.06 were detected in Xuchang and Luoyang Pop2, respectively, by Wei et al. [57]. Wang et al. [58] suggested that IVDMD QTL explained more than 10 % of the genetic variation in both F3 and F4 generations were mapped on the same genomic position on chromosome 6, which were the same as QTL detected in maize recombinant inbred line progeny of F288 × F271 [64]. One IVDMD QTL detected in our study was also on chromosome 6. These QTL described above were closely linked under high-density SNP markers and deserve further investigation for finding candidate genes underlying IVDMD in a no-forage genetic background.
The role of genome dosage changes on gene expression of lignin and cellulose synthesis in inbred and haploids of two parental lines Most candidate genes were involved in lignin and cellulose synthesis which affect the stalk cell wall structure. Lignin was a phenolic polymer that imparted mechanical strength of the plant secondary cell wall, and therefore, was considered to confer stalk rot resistance and involve in plant evolution [65]. Particularly, genes participating in lignin synthesis were identified only in haploid population in our study. Therefore, based on the gene function annotations for RPR and IVDMD QTL detected in the DH and haploid populations, we analyzed the key gene expressions of lignin and cellulose synthesis and genome dosage regulation. The expression levels and phenotypes showed several interesting results, suggesting a partial explanation for ploidy effect mechanisms in the haploid condition.
In the one dosage genome, compared to the inbred, CesA11 and CesA112 gene expressions were up-regulated except for the CesA11 and CesA12 gene in Z58, which was consistent with the decreased Cel content in haploid Z58 and higher content in C7-2 haploids (Fig. 4). Unlike other gene expressions in lignin synthesis, the COMT gene showed significantly lower expression levels in C7-2 haploids than that in Z58 haploids. All these results illustrated the existence of genetic variation in morphological responses to ploidy changes, which was consistent with the results which reported by Riddle et al. [66].
All lignin and synthesis gene expression showed dramatically inversed expression levels, except for CesA11 in Z58, as genome dosage number increased (Fig. 4), which was consistent with higher Lig contents for two parents and increased Cel contents for parent line C7-2 in the haploid condition and can be explained by dosage compensation [67]. These findings implied that the major effect of genome dosage changes were not simply proportional to copy number of regulatory genes and may not be directly related to the phenotype, but rather result from various regulatory components which resulted in compensation for gene expression in the one dosage genome. Previous studies had also reported that many target loci exhibited dosage compensation, such as Adh expression with reduced alleles number [68]. Quantitative traits trended to be regulated by genes exhibiting dosage effects and most were transcription factors [69,70]. Fortyseven dosage-dependent modifiers operated in macromolecular complexes in a regulatory network [71]. Many other mechanisms, such as siRNA, DNA methylation and so on, may be concerned with expression changes in dosage-sensitive region genes, which can also sharpen the evolution of copy-number varied regions [72][73][74].

Conclusions
Using DH and its corresponding haploid populations, this analyse revealed important genomic regions associated with eight stalk-related traits. These QTL explained a large extent of phenotypic variance. One or more QTL sharing common region for each stalk-related trait were detected between this two different ploidy populations. The heritabilities in haploid population were higher than DH population in all stalk traits except WC. This study identified candidate genes involving in lignin and cellulose synthesis for RPR and IVDMD, which were the two most important stalk traits. The expression levels of most of these candidate genes were significant higher in haploid parents than that in corresponding diploid parents. Haploid population may be used as one of platforms providing information on the genetic basis for stalk-related traits, and the genetic connection between DH and haploids for traits can facilitate the selection of materials at haploid-stage which can boost the practical maize breeding.
The dosage compensation mechanism and dosagesensitivity genes may be further examined by analyzing the genetic architecture of certain traits by comparing QTL mapping results between different ploidy populations, which may also have important implications for understanding gene regulatory networks and genome evolution.

Materials and population construction
The parental lines of two populations belong to two distinct maize germplasm groups. The maternal inbred line Zheng58 (Z58) is dent corn of the Reid heterotic group, and the paternal inbred line Chang7-2 (C7-2) is flint corn. The single cross (ZD958) of these two parental lines is one of the most widely planted commercial varieties in China.
The DH population was constructed by production of haploids with an in vivo haploid induction procedure and followed by chromosome doubling (Additional file 1: Figure S1).
In total, the DH population consisted of 190 DH lines. Accordingly, the haploid population was constructed by crossing each DH line, used as source germplasm, with a haploid inducer, used as pollinator, and we finally produced 170 haploid lines. To obtain enough DH and haploid seeds for evaluating phenotypes in field trial, in the experimental station of the China Agricultural University (Yacheng in Hainan Province), we planted two replications of 190 DH lines with 21 plants in each plot that again pollinated with inducer line CAU5 at the anthesis stage [75]. The haploid seeds were manually picked up from the harvested ears based on the R1-nj color markers [76].

Field experiments
In 2013 and 2014, the two populations, their parental lines (inbred and haploid) and F 1 generation were sown in two experimental stations of the China Agricultural University (Shangzhuang in Beijing and Quzhou in Hebei Province), and Shijiazhuang experimental station of Hebei Academy of Agriculture and Forestry Sciences. At each location, the two populations were adjacently planted to reduce the influence of field heterogeneity. We used a randomized complete block design with two blocks for both populations. Within each block, each line was assigned to a single row plot (Additional file 1: Figure S2).

Trait evaluations
The fourth internode above ground (FIAG) was measured for nine traits described as follows. Three randomly selected plants in each row (i.e., plot) were chosen for trait evaluation. RPR was measured with an electronic penetrometer (AWOS-SL04, Aiwoshi Company, Hebei, China) at the milky stage. The measurement of RPR and water content (WC) of the internode, as well as measurement of IVDMD, acid detergent fiber (ADF), neutral detergent fiber (NDF), lignin (Lig), cellulose (Cel) and water soluble carbohydrate (WSC) followed a standard procedure described by Hu et al. [20].

Phenotypic data analysis
In order to evaluate phenotypic traits, all DH lines and haploid lines, were planted in 3 locations in 2013 and 2014, which were treated totally as 5 independent macro environments (not including Quzhou 2014 due to large proportion of missing data points,). Then a linear model was used to perform an analysis of variance (ANOVA), genotypic value estimation (BLUE) and variance components estimation for the DH and haploid populations for each trait: where μ is the grand mean, G is genotypic effect, E is the environment effect, G × E is the genotype-byenvironment interaction, R(E) is the effect of block within environment, and ε is the random error. All statistical analyses were performed using SAS version 9.3 [77]. The genotypic effect G was estimated when it was treated as a fixed effect using LSMEANs with PROC GLM. The genotypic value estimates (BLUEs) were used to calculate correlation coefficients by PROC CORR with Pearson's method. The variance components were estimated using PROC VARCOMP with the method of restricted maximum likelihood (REML). The estimates of genotypic variance (σ G 2 ), genotype-by-environment interaction (σ G × E 2 ) and random error (σ e 2 ) were used to estimate heritability based on the formula [78]: where l is the number of macro environments and r is the number of blocks within each environment, which equals 5 and 3, respectively in the study.

Genotyping and genetic map construction
Genomic DNA was extracted from young leaves of DH and the parental lines using a CTAB method [79], and was then purified. The purified DNA was genotyped with the maizeSNP3K chip (3,072SNPs), which is a subset of the Illumina MaizeSNP50 BeadChip [80]. SNP genotyping was performed on the Illumina Golden-Gate SNP genotyping platform at the National Maize Improvement Center of China of the China Agricultural University. Checking the quantity of each SNP was carried out manually as described by Yang et al. [81]. The qualified SNPs were reserved for further screening and creation of linkage maps. After genotyped, 1228 polymorphic markers between the two parental lines were determined. Subsequently, heterozygosity of each line, the missing rate, minor allele frequency (MAF) and heterozygosity of each SNP were calculated. The DH lines with heterozygosity ≤ 0.1 and the SNPs having polymorphisms between two parents with a missing rate ≤ 0.2 and MAF ≥ 0.05 were selected to construct a genetic linkage map with software package Joinmap4.0. SNPs and lines that shared 100 % similarities were deleted according to the Joinmap4.0 instructions [32]. The linkage groups were composed at a minimum Logarithm of odds (LOD) of 7, and the Kosambi mapping function and regression-mapping algorithm were used for calculating map distances.

QTL analysis
Since the segregating populations used in this study were DH and haploid populations, an additive genetic model was chosen for QTL analysis, using the BLUEs across environments as phenotypic data. Composite interval mapping (CIM) with a regression approach [82] in combination with the use of cofactors was employed [83]. A two-step procedure was utilized for QTL detection as described by Hu et al. [84]. A threshold of LOD =3.0 was used to determine QTL for all traits based on 2000 permutations [85]. 1-LOD support intervals were calculated from the significant peak to a certain position on both sides along the chromosome at which the LOD score had a 1.0 unit decrease [86]. To facilitate comparisons among linked QTL, two QTL were designated as overlapping when they were separated by less than 20 cM [87]. The total proportion of genotypic variation (p G ) explained by the detected QTL was calculated by the formula p G = R adj 2 /h B 2 [78]. Fivefold cross-validation was used to assess the reliability of QTL mapping results with 1000 runs [88].
Physical positions of QTL were obtained from the physical positions of their flanking markers. All genes between the flanking markers of each QTL were extracted from the filtered gene set of the maize genome sequence (http://ftp.maizesequence.org/release-5b/filtered-set/) by qteller3 (http://qteller.com/qteller3/index.php), assuming a linear relationship between recombination and physical distances within this interval.

RNA extraction and RT-PCR
Stalk rind used in gene expression analysis was collected from plant FIAG at the milky stage. Total RNA was extracted using Ultrapure RNA Reagent (Cat#CW0581, CWbio.Co.Ltd, Beijing, China) according to the manufacturer's instructions. Total RNA was digested with RNase-Free DNase I to remove any DNA (Cat# CW2090, CWbio.Co.Ltd, Beijing, China). cDNA was synthesized using HiFi-MMLVcDNA (Cat#CW0744, CWbio.Co.Ltd, Beijing, China). Two microliters of diluted cDNA was used for quantitative RT-PCR analysis using the primer pairs listed in Additional file 1: Table S1. An ABI7500 machine was used to conduct RT-PCR analysis [89]. Actin was used as an internal control to calculate the relative