Identification of candidate genes associated with cell wall digestibility and eQTL (expression quantitative trait loci) analysis in a Flint × Flint maize recombinant inbred line population

Background Cell-wall digestibility is the major target for improving the feeding value of forage maize. An understanding of the molecular basis for cell-wall digestibility is crucial towards breeding of highly digestible maize. Results 865 candidate ESTs for cell-wall digestibility were selected according to the analysis of expression profiles in 1) three sets of brown-midrib isogenic lines in the genetic background of inbreds 1332 (1332 and 1332 bm3), 5361 (5361 and 5361 bm3), and F2 (F2, F2 bm1, F2 bm2, and F2 bm3), 2) the contrasting extreme lines of FD (Flint × Dent, AS08 × AS 06), DD1 (Dent × Dent, AS11 × AS09), and DD2 (Dent × Dent, AS29 × AS30) mapping populations, and 3) two contrasting isogenic inbreds, AS20 and AS21. Out of those, 439 ESTs were assembled on our "Forage Quality Array", a small microarray specific for cell wall digestibility related experiments. Transcript profiles of 40 lines of a Flint × Flint population were monitored using the Forage Quality Array, which were contrasting for cell wall digestibility. Using t-tests (p < 0.01), the expression patterns of 102 ESTs were significantly different between high and low quality groups. Using interval mapping, eQTL (LOD ≥ 2.4) were detected for 20% (89 of 439) of the spotted ESTs. On average, these eQTL explained 39% of the transcription variation of the corresponding ESTs. Only 26% (23 of 89) ESTs detected a single eQTL. eQTL hotspots, containing greater than 5% of the total number of eQTL, were located in chromosomal bins 1.07, 1.12, 3.05, 8.03, and 9.04, respectively. Bin 3.05 was co-localized with a cell-wall digestibility related QTL cluster. Conclusion 102 candidate genes for cell-wall digestibility were validated by genetical genomics approach. Although the cDNA array highlights gene types (the tested gene and any close family members), trans-acting factors or metabolic bottlenecks seem to play the major role in controlling heritable variation of gene expression related to cell-wall digestibility, since no in silico mapped ESTs were in the same location as their own eQTL. Transcriptional variation was generally found to be oligogenic rather than monogenic inherited due to only 26% ESTs detected a single eQTL in the present study. One eQTL hotspot was co-localized with cell wall digestibility related QTL cluster on bins 3.05, implying that in this case the gene(s) underlying QTL and eQTL are identical. As the field of genetical genomics develops, it is expected to significantly improve our knowledge about complex traits, such as cell wall degradability. Comprehensive knowledge of the lignin pathway and cell wall biogenesis will allow plant breeders to choose the best genomic targets controlling these characters, for improving forage digestibility through genetic engineering or marker-assisted selection.


Background
Lignin content is well known as a major factor affecting forage quality in maize. However, correlations between lignin content and forage quality can be variable according to the genetic background [1]. Moreover, breeding for a higher digestibility of maize involves also other, so far unknown mechanisms [1]. Correlations between maize whole plant digestibility and cell-wall (or stover) digestibility ranged from 0.60 to 0.96, with average values close to 0.80, whereas correlations between whole-plant digestibility and grain or ear content were close to 0.4 [2]. Therefore, cell-wall digestibility is the major target for improving the feeding value of forage maize. An understanding of the molecular basis for cell-wall digestibility is crucial towards breeding of highly digestible maize.
An important first step to elucidate the mechanisms underlying cell-wall degradability is to identify causative genome regions. Five major QTL clusters involved in cellwall digestibility, located on chromosomal bins [3] 1.03, 3.05/06, 6.06, 8.05, and 9.02, were identified by several QTL analyses, but the genes underlying these QTL are not yet known [1]. A second, more recent approach is transcriptome analysis to simultaneously measure the expression of thousands of genes. In comparison with normal maize genotypes, brown-midrib (bm) mutants show a significantly reduced lignin content, altered lignin composition, and/or a significantly higher cell-wall digestibility [4]. Molecular mechanisms underlying cell wall digestibility in maize have been studied in three sets of maize brown-midrib isogenic lines in the genetic background of inbreds 1332 (1332 and 1332 bm3), 5361 (5361 and 5361 bm3), and F2 (F2, F2 bm1, F2 bm2, and F2 bm3) [5]. 53 ESTs were differentially expressed in all three isogenic bm3 comparisons, whereas 32 ESTs were consistently differentially expressed in different bm isogenic lines in F2 background. Moreover, gene expression studies can be conducted on phenotypically extreme lines from mapping populations. Replicate pools of extreme lines can be profiled independently, so that differences in gene expression will be specific to the differing pools. This strategy was recently used to identify candidate genes for drought response QTL in rice [6]. Thus, by synthesis of expression profiling data from bm mutants and extreme lines of a mapping population segregating for cell wall digestibility, it should be possible to identify candidate genes related to cell-wall degradability, and to construct a microarray enriched for candidate genes underlying cell wall digestibility. Transcriptome analysis using such microarrays would provide a fingerprint of cell-wall metabolism in maize.
Genetic and gene expression approaches have been joined in the concept of "genetical genomics" [7], which aims to detect eQTL (expression quantitative trait loci) control-ling gene expression differences. Often, eQTL map to the genetic position of the respective gene itself, indicating that cis changes (within the gene) are responsible for the different levels of expression. In contrast, genes revealing (trans) eQTL at positions different from the genetic position of the respective gene are thought to be regulated by, e.g., trans-acting factors controlling their expression levels [8]. Detection of the master regulators, affecting expression levels of groups of genes, is a major feature of eQTL studies [9]. In plants, this strategy has been successfully applied to 76 maize lines in a F3 population [10] and 91 poplar lines in a backcross population [11], respectively. Altogether, combining expression profiling with genetic analysis could enrich our understanding of regulatory networks underlying cell-wall digestibility and assist plant breeders to choose the most relevant genomic targets for improvement of silage maize digestibility.
The objectives of our study were to 1) select candidate ESTs for cell-wall digestibility to establish a "Forage Quality Array", 2) identify ESTs differentially expressed between low and high digestible lines in a Flint × Flint mapping population, and 3) detect eQTL using the "genetical genomics" approach.

Selection of candidate ESTs in association with cell wall digestibility
In order to identify genes in association with cell-wall digestibility in maize, three sources of genetic material were used, including 1) three sets of brown-midrib isogenic lines in the genetic background of inbreds 1332 (1332 and 1332 bm3), 5361 (5361 and 5361 bm3), and F2 (F2, F2 bm1, F2 bm2, and F2 bm3) [5], 2) the contrasting extreme lines of FD, DD1, and DD2 DH mapping populations, and 3) two isogenic inbreds, AS20 and AS21, significantly differing in cell wall digestibility. Two complementary approaches, SSH (suppression subtractive hybridization) and microarray-based expression profiling, were used to isolate and identify candidate genes in all comparisons. About 70% of the ESTs isolated by SSH were absent on the unigene microarray [5].
The number of at least two-fold induced ESTs in stems from SSH and microarray analyses, respectively, were 246 and 1417 in FD-pop (Figure 1), 122 and 317 in DD1-pop, 71 and 1805 in DD2-pop, as well as 225 and 630 in AS20 vs. AS21 (Data not shown). In FD-pop, 34 ESTs were jointly differentially expressed both in SSH and microarray experiments, as well as 8, 33, and 31 ESTs in DD1-pop, DD2-pop, and AS20 vs. AS21, respectively. The numbers of genes simultaneously differentially expressed in pairs of mapping populations were 145 (between FD-pop and DD1-pop), 374 (between FD-pop and DD2-pop), and 135 (between DD1-pop and DD2-pop), respectively. 58 ESTs were consistently differentially expressed across the three mapping populations.
In total, 5460 distinct ESTs differentially expressed in one or multi-comparisons were identified. In order to verify the most interesting of those genes at low costs, a small microarray specific for cell wall digestibility (Forage Qual-ity Array) was designed. The first selection criterion (Figure 2) covered jointly differentially expressed ESTs both in  SSH and microarray experiments for each comparison  (287 ESTs). Secondly, 460 ESTs were selected with a more than five-fold change in respective comparisons. Thirdly, 91 EST homologues of lignin-related genes, PAL, 4CL, C3H, CCoAOMT, CCR, COMT, and CAD, were spotted.

FD-pop
Using hierarchical cluster analysis based on the EST tree, distinct expression patterns between high and low quality groups were found for most ESTs included on the Forage Quality Array (Figure 4). Using t-tests with a p value cutoff of 0.01, 102 ESTs were significantly different between the two groups of high and low quality lines, including 39 down-and 63 up-regulated ESTs in the high quality group, respectively (  sified into "chaperone activity", "catalytic activity", and "binding", respectively, whereas 1, 11, 2, 1, and 4 ESTs of up-regulated ESTs were classified into "motor activity", "catalytic activity", "structural molecular activity", "transporter activity" and "binding", respectively. According to the current annotation in the TIGR Maize Gene Index [12], four ESTs encoded enzymes involved in lignin biosynthesis, CD970581 (Cinnamyl-alcohol dehydrogenase), CD973094 (Phenylalanine ammonia-lyase), AW120445 (Class III peroxidase 67 precursor) and CF243853 (Caffeoyl CoA 3-O-methyltransferase). In addition, AI622068 (Cellulose synthase-2) is a homologue encoding a cellulose biosynthesis related enzyme, and BM333894 (YABBY-like transcription factor) and BU093700 (transcription factor and jumonji family protein) are putative transcription factors.

Genetical genomics analysis
40 extreme lines in the FF population with respect to dNDF were selected to represent two contrasting high and low quality groups ( Assignment of the 439 ESTs spotted on the "Forage quality array" to functional classes by GO vocabularies Figure 3 Assignment of the 439 ESTs spotted on the "Forage quality array" to functional classes by GO vocabularies. Each EST was assigned with GO vocabularies according to the TIGR Maize Gene Index [12]. Unclassified ESTs (66%) were not taken into account.
The expression profiling of 439 ESTs in 2 × 20 recombinant inbred lines contrasting for dNDF   High quality Low quality   Maize genetic maps have been divided into 100 segments (=bins) of approximately 20 cM length [3]. In total, 271 eQTL were detected in 24 bins covering all chromosomes except for chromosome 7 (Figure 5b). The percentage of eQTL per bin was plotted against chromosomal bins (Figure 5). eQTL hotspots, containing more than 5% of the total number of eQTL, were detected on bins 1.07, 1.12, 3.05, 8.03, and 9.04. Out of those, bin 3.05 co-localized with one of five major QTL clusters for cell-wall digestibility ( Figure 6), which were identified on bins 1.03, 3.05/ 06, 6.06, 8.05, and 9.02 [1].
Out of 89 ESTs resulting in eQTL, for 14 ESTs bin information from the Maize GDB [13] and map positions from Maize Genome Mapping project [14] (Figure 6) are known. For another 10 ESTs, only chromosomal but not bin assignment were provided from Maize Genome Mapping project. In silico mapping information indicate that the ESTs were located in various linkage groups, and generally not in the genetic location of the respective eQTL. AW438124 (nucleosome/chromatin assembly factor C), BM080577 (nodulin MtN3 family protein), and AI881378 (maize clone contig234 mRNA sequence) were co-localized with a cell-wall digestibility related QTL cluster on bin 3.05.

Power of eQTL analysis
An essential first step in QTL analysis is to assess how many samples must be collected in order to achieve sufficient power to detect a hypothetical effect. eQTL mapping is different from QTL mapping of complex inherited traits, since we expect (and are mainly interested in) major eQTL due to 1) eQTL in cis (then most of the variation should be explained by this one eQTL), and 2) major regulatory eQTL (one or few) in trans. Thus, we expect a simpler inheritance of expression patterns of single ESTs with only one or few eQTL involved, each explaining a comparatively large percentage of the phenotypic variation, in contrast to "true" quantitative characters such as grain yield. Major eQTL should be detectable even with this very limited number of RILs. In addition, this requirement is particularly salient in eQTL analysis due to the high expenses of gene-expression analyses. In plants, both 76 F3 families in maize [10] and 91 BC1 families in poplar [11] were used in respective eQTL analysis. In this study, 40 RIL lines were used to detect eQTL in the FF population. Because of an increased heritability when using homozygous and homogeneous RIL lines as compared to segregating F2:3 families, an eQTL explaining 28% of the variation in RNA abundance among the F2 population  [9]. In addition, at least 73% of the experiments with 40 RIL lines should detect an eQTL explaining 18% of the phenotypic variance [9]. Since 79% (214 of 271) identified eQTL explained more than 28% of the transcription variation of the corresponding ESTs and minimum eQTL explaining 24.1% transcription variation, 40 extreme lines should be sufficient to detect major eQTL. However, the explained phenotypic variance by individual eQTL is very likely overestimated due to the small population size [15].

Genetic architecture of transcription
Genetical genomics permits quantitative assessment of the proportion of gene loci (ESTs) displaying co-segregation with the respective eQTL (cis-eQTL), relative to those that produce unlinked trans-eQTL [13]. In maize ear leaf, 34% (6481 of 18,805) differentially expressed genes produced cis-eQTL with LOD >3.0 [10]. In contrast, no in silico mapped ESTs were in the same location as their own eQTL in this study ( Figure 6). Accordingly, no key genes involved in lignin biosynthesis co-localized with cis-eQTL in poplar [11]. Therefore, trans-acting factors could play a major role in regulating heritable variation of gene expression in cell wall digestibility.
Only 26% (23 of 89) ESTs detected a single eQTL in the present study. For the remaining 74%, two or more eQTL were found ( Figure 5). A similar distribution of eQTL was found in yeast [16]. Only 3% transcripts were consistent with single locus inheritance, 18% suggested control by two loci, and >50% required at least five loci under an additive model. Thus, transcriptional variation was gener- LG1 LG3 LG2 LG4 LG5 LG6 LG7 LG8 LG9 LG10 Bin   [13] and the IDP mapping project [27], according to map bins [3]. The x axis represents ten linkage groups of the Flint × Flint genetic map arranged end-to-end. The middle line represents LOD score threshold 2.4. The "× " indicates five major QTL clusters for cell wall digestibility on bins 1.03, 3.05/06, 6.06, 8.05, and 9.02 [1]. The dashed rectangle shows some eQTLs co-localizing with one major QTL cluster for cell wall digestibility on bin 3.05.  ally found to be oligogenic rather than monogenic inherited.
A common feature of eQTL studies is the detection of "hotspots" or hubs of trans-acting eQTL: chromosomal regions that affect the expression of a much larger number of genes than expected by chance [9]. eQTL hotspots, containing more than 5% of the total number of eQTL, were found on bins 1.07, 1.12, 3.05, 8.03, and 9.04 ( Figure 5). The strong clustering in hubs of eQTL reflects highly correlated expression levels of many gene transcripts in association with cell-wall digestibility. The five hotspots could contain important transcript factors for cell wall digestibility. However, no obvious transcription factors were found among those genes currently mapped to these bin regions [13]. Similarly, eQTL identified in yeast frequently resulted from genetic variants that altered activity of a signaling or metabolic pathway without directly altering gene transcription [17]. 4CL2, F5H, and CCR1 genes involved in lignin biosynthesis map to bin 1.07 [1]. These genes are candidates for the eQTL found in the same genome region due to potential metabolic feedback mechanisms acting through the same enzyme(s) in lignin biosynthesis. However, due to the limited understanding of genetic and physiological control of gene expression so far, any conclusions with regard to hotspots for gene regulation should be interpreted with caution [9]. Kirst et al. [11] found that two eQTL hubs for lignin-related genes co-localized with growth QTL on linkage groups 4 and 9 in poplar. One eQTL hotspot was co-localized with cell wall digestibility related QTL cluster [1] on bins 3.05 (Figure 6), implying that in this case the gene(s) underlying QTL and eQTL are identical.

Molecular mechanisms underlying cell-wall digestibility
Lignin restricts the degradation of structural polysaccharides by hydrolytic enzymes, thereby limiting the bioconversion of forages into animal products [18]. Out of 102 significantly differentially expressed ESTs between high and low quality groups (Table 2), four ESTs encode enzymes involved in lignin biosynthesis, CD970581 (Cinnamyl-alcohol dehydrogenase), CD973094 (Phenylalanine ammonia-lyase), BE051493 (4-coumarate-CoA ligase-like protein), and CF243853 (Caffeoyl CoA 3-Omethyltransferase). In agreement with the expectation of low levels of lignin, CD970581 and CD973094 were down-regulated in lines with high forage quality. However, BE051493 and CF243853 were up-regulated, which might be due to the complexity of cell wall biosynthesis.
The associations between one lignin characteristic (e.g., lignin composition) and cell wall degradability can be confounded and correlated to concurrent changes in other lignin properties (e.g., lignin cross-linking) that influence cell wall degradability [18]. Since reduced cellulose synthesis invoked lignification and defense responses in Ara-bidopsis [19], repressed ESTs encoding cellulose synthesis enzymes, such as AI622068 (Cellulose synthase-2) ( Table  2), could activate lignin synthesis in low quality lines. In addition, four regulatory genes, AI734769 (Class III HD-Zip protein 4), BM080754(Translational initiation factor eIF-4A), BM334752 (Zinc finger protein family-like), and BM333894 (YABBY-like transcription factor), were downregulated in lines with high forage quality, whereas one regulatory gene, BU093700 (transcription factor and jumonji family protein) was up-regulated (Table 2). These regulatory genes might control cell wall digestibility spatially and temporally among the two groups of lines in the FF population. Mapping of these genes will reveal, whether they colocalize with eQTL (clusters) identified in this study. Furthermore, there are a number of genes, where a connection to digestibility is unclear, some of which might be false positives.
Using cellular UV-microspectrophotometry, Shi et al. [5] demonstrated that lignin content represented by the absorbance values within the cells of sclerenchyma fiber tissue was decreasing from the periphery to the center. Lignin content was lowest in S2 (secondary cell wall) of parenchyma cells, followed by S2, CML (compound middle lamella) and CC (cell corner) of sclerenchyma fibers. Each cell type in the stem likely expresses a unique transcriptome. Transcriptome analysis using complete stems provides average gene expression levels integrated over all cell types. In order to further understand the molecular basis of cell wall degradability, new methods, for instance, laser-capture microdissection [20], are needed for efficient isolation of CC of sclerenchyma fiber from different time points to further validate 439 candidate genes.

Application of genetical genomics to plant breeding
The use of functional genomics is contributing to many aspects of QTL analysis and cloning [21]. To date, most plant QTL have been cloned using a positional cloning approach following identification in experimental crosses. Transcriptional profiling can quickly provide a list of differentially expressed genes between contrasting QTL genotypes. Subsequently, those genes functionally related to the target trait and mapping to the QTL region can be selected as candidates. Unfortunately, the number of QTL cloned to date in plants is too small to test the validity of this approach. If the difference in gene expression level between alleles is too low (approximately twofold), the candidate genes cannot be identified using standard microarray-based transcriptome analysis [21]. Less than 15% of ESTs were induced more than two-fold by microarray analyses in FD-pop (12%, 1,417 of 11,827) ( Figure  1), DD1-pop (3%, 317 of 11,827) and DD2-pop (15%, 1,805 of 11,827), respectively. A genetical genomics approach might be better suited for detecting differences in gene expression. 23% of the ESTs on the Forage Quality Array (102 of 439) were significantly different between low and high digestible (dNDF) RILs. However, this higher proportion of significantly differentially expressed genes identified by the genetical genomics approach is most likely due to successful selection of candidates by the criteria outlined above (Figure 2).
Since most of the important agronomic traits are complex inherited traits, transcript abundance may act as an intermediate phenotype between genomic DNA sequence variation and complex traits. Integration of genome-wide expression profiling with linkage analysis is a new approach to identifying genes underlying complex traits. Hubner et al. [22] has demonstrated large-scale identification of positional candidates and regulatory pathways for previously mapped physiological QTL by genetical genomics approach in the rat. Cis-acting eQTL are good candidates for physiological QTL because they show linespecific differences in gene expression that are under the control of DNA sequence variants in or close to the gene itself [23]. Similarly, eQTL identified in an intercross of inbred C57BL/6J and DBA/2J mice could accurately identify transcripts in which expression was regulated by cisacting polymorphisms [24]. Cis-acting eQTL might be detected by additional EST mapping or indirectly through the ongoing maize genome sequencing project, since so far only 24 of 89 ESTs resulting in eQTL were in silico mapped. Trans-acting eQTL represent loci that influence expression of genes or transcripts remote from the eQTL itself. Coincidental mapping of trans-acting eQTLs for multiple transcripts to the same chromosomal location, as observed on bin 3.05 ( Figure 6), may represent a shared regulatory transcriptional control mechanism by a single gene at the eQTL. The locations of trans-acting eQTL in relation to physiological QTL, together with the locations of cis-acting eQTL, may point to genes and regulatory pathways underlying individual cell wall digestibility related QTL.
However, the scope of genetic analysis of gene expression also presents enormous technical and analytical challenges. Extraordinarily large number of comparisons were involved in a genome-wide linkage scan for several hundreds or thousands of transcripts, thus technical and experimental design issues need to be addressed to handle the large data sets that are being generated, and new statistical tools are still being evaluated [25]. As the field of genetical genomics develops, it is expected to significantly improve our knowledge about complex traits, such as cell wall degradability [8]. Therefore, in the future, it is conceivable that QTL cloning will increasingly rely on candidate gene information. Comprehensive knowledge of the lignin pathway and cell wall biogenesis will allow plant breeders to choose the best genomic targets controlling these characters, for improving forage digestibility through genetic engineering or marker-assisted selection.

Conclusion
102 candidate genes for cell-wall digestibility were validated by genetical genomics approach. Although the cDNA array highlights gene types (the tested gene and any close family members), trans-acting factors or metabolic bottlenecks seem to play the major role in controlling heritable variation of gene expression related to cell-wall digestibility, since no in silico mapped ESTs were in the same location as their own eQTL. Transcriptional variation was generally found to be oligogenic rather than monogenic inherited due to only 26% ESTs detected a single eQTL in the present study. One eQTL hotspot was colocalized with cell wall digestibility related QTL cluster on bins 3.05, implying that in this case the gene(s) underlying QTL and eQTL are identical. As the field of genetical genomics develops, it is expected to significantly improve our knowledge about complex traits, such as cell wall degradability. Comprehensive knowledge of the lignin pathway and cell wall biogenesis will allow plant breeders to choose the best genomic targets controlling these characters, for improving forage digestibility through genetic engineering or marker-assisted selection.

Methods
Plant materials AS20, AS21, recombinant inbred lines (RILs), and doubled haploid (DH) lines were obtained from KWS Saat AG, including twenty highly digestible and twenty low digestible lines each in four mapping populations: FD (Flint × Dent DH population, AS08 × AS 06), DD1 (Dent × Dent DH population, AS11 × AS09), DD2 (Dent × Dent DH population, AS29 × AS30), and FF (Flint × Flint RIL population, AS18 × AS07) RIL mapping populations, respectively. These lines can be obtained from KWS Saat AG for non-commercial research purposes. All plants were grown and maintained in the greenhouse under a 12 h photoperiod at 23°C and 50% relative humidity. Stems were harvested 5 weeks after germination, since the highest level of COMT expression was observed in stems 5 and 7 weeks after germination [5]. For biological replication, two independent sets of three plants per set were harvested for AS20, AS21, and every line in FD, DD1, and DD2 populations, respectively, and one set of ten plants per line in the FF population. Since the maize stalk is small and covered by leaf sheath at early vegetative stages, the stem fraction investigated in this study is a mixture of stalk and leaf sheath [26]. Total RNA was extracted from the stems of each set of each line using TRIzol reagent (Invitrogen, Carlsbad, CA, USA).

Construction and hybridisation of Forage Quality Array
The macroarrays containing SSH clones were derived from the studies on three sets of maize brown-midrib isogenic lines in the genetic background of inbreds 1332 (1332 and 1332 bm3), 5361 (5361 and 5361 bm3), and F2 (F2, F2 bm1, F2 bm2, and F2 bm3) [5]. In total, 2688 clones were spotted in duplicate on each macroarray. These clones were randomly picked from five SSH libraries. In the FD, DD1, and DD2 populations, the total RNA representing high and low forage quality DH lines was pooled from the total RNA of those twenty lines in each of the populations with the highest or lowest digestibility (DNDF), respectively. The pooled total RNA and the total RNA of AS20 and AS21 were labeled with 32 P and hybridized to different arrays using the Strip-EZ RT kit (Ambion, Austin, TX, USA). Hybridization signals were detected by the Storm 860 imaging system (Amersham) with a resolution of 50 μm. After scanning, labeled cDNA probes were completely stripped from the arrays using the Strip-EZ system (Ambion). Each array was hybridized four times with four replicates per line, including two biological replications and two labeling replications in each biological replica.
Maize unigene-microarrays were provided by the laboratory of Prof. Schnable (Iowa University, USA) and contained 11,827 maize ESTs [27]. Poly (A)+ RNA was isolated from the same total RNA for probe preparation in the SSH approach via Dynabeads ® Oligo(dT) 25 [29]. All calculated gene expression ratios were log 2transformed and averaged over four replicates in each comparison. Differentially expressed ESTs at the 95% confidence level were determined using intensity-dependent Z-scores (with Z = 1.96) as implemented in MIDAS and the union of all genes identified in each comparison was considered significant in this experiment. Each EST was annotated and assigned with GO vocabularies according to the TIGR Maize Gene Index [12]. The mapped ESTs contain bin information from the Maize GDB [13] and the IDP mapping project [27].
Among 865 candidate ESTs for cell wall digestibility, 208 ESTs were ordered from the Iowa Schnable Lab [30] and 506 ESTs from the Arizona BAC/EST resource center [31]. From the stab cultures of these ESTs, plasmid minipreps were conducted by use of the R.E.A.L. Prep 96 Plasmid Kit (Qiagene AG, Germany). Using the plasmids of each EST as template, two independent re-amplification (100 μl) reactions were performed and pooled to reduce the effects of variation in PCR efficiency. Pooled PCR products were concentrated from 200 μl to approximately 25 μl using MultiScreen-PCR plates (Millipore, Billerica, MA, USA). Due to poor bacteria recovery and unspecific PCR amplification, high quality PCR products were derived from only 439 ESTs and spotted on Nexterion ® Slide A+ glass slides (SCHOTT Jenaer Glas GmbH, Germany) in a triple pattern using a QarrayMini spotter (Genetix GmbH, Germany). The total RNA of each line in the FF population was labeled and hybridized in the same way as the total RNA of pooled lines in the other three populations. Using loop design [32], each Forage Quality Array was hybridized with two contrasting lines (Table 1) labeled with different dyes. In total, 40 Forage Quality Arrays were used and each line was employed with dye-swap replication.

Data analysis of the "genetical genomics" experiment
The raw images of Forage Quality Arrays obtained from a GeneTAC™ UC4 Microarray Scanner (Genomics Solutions Ltd., USA) were imported into ArrayVision 8.0 for spot detection and quantification of hybridization signals. In ArrayVision 8.0, Lowess normalization was implemented to remove dye effects. The transcript level of each EST was represented by the average over six measures, including dye swap replication and three in-slide replications within either Cy3 or Cy5 replications.
The dataset of all 40 lines in the FF population were imported in Multiexperiment Viewer (MeV) [29]. "Hierarchical cluster analysis" [33] was conducted to discover similar expression patterns across 439 ESTs, and "between subject t-tests (p < 0.05)" [34] were utilized to identify candidate genes differentially expressed between high and low digestible lines.
A complete linkage map was constructed by KWS Saat AG with 156 loci using 270 lines, including forty extreme lines, in FF population. Low quality lines have all unfavorable QTL/alleles for dNDF, and high quality all favorable QTL/alleles for dNDF. The same mapping population was used for both QTL (Krützfeldt et al., in preparation) and eQTL mapping. Five dNDF QTLs were identified, and the positions of the QTL were in good agreement with previous QTL studies [1]. Identification of QTL for gene expression traits was performed based on the 40 extreme lines using interval mapping implemented in MapQTL [35]. An empirical LOD score threshold of 2.4 was determined by permutation testing [36,37] and adopted from Kirst et al. [11].
Original microarray expression data presented in this manuscript are available through ArrayExpress [38] with accession number E-MEXP-253.

Authors' contributions
CS carried out the most of experiments, performed the data analysis, and drafted the manuscript. AU carried out the hybridization of Forage Quality Arrays. MO conceived of the study and participated in its coordination. ML participated in the design of the study. GW conceived of the study and assisted CS with the experiments. TL conceived of the study, participated in its design, and drafted the manuscript. All authors read and approved the final manuscript.