Heterosis in Early Maize Ear Inflorescence Development: A Genome-Wide Transcription Analysis for Two Maize Inbred Lines and Their Hybrid

Heterosis, or hybrid vigor, contributes to superior agronomic performance of hybrids compared to their inbred parents. Despite its importance, little is known about the genetic and molecular basis of heterosis. Early maize ear inflorescences formation affects grain yield, and are thus an excellent model for molecular mechanisms involved in heterosis. To determine the parental contributions and their regulation during maize ear-development-genesis, we analyzed genome-wide digital gene expression profiles in two maize elite inbred lines (B73 and Mo17) and their F1 hybrid using deep sequencing technology. Our analysis revealed 17,128 genes expressed in these three genotypes and 22,789 genes expressed collectively in the present study. Approximately 38% of the genes were differentially expressed in early maize ear inflorescences from heterotic cross, including many transcription factor genes and some presence/absence variations (PAVs) genes, and exhibited multiple modes of gene action. These different genes showing differential expression patterns were mainly enriched in five cellular component categories (organelle, cell, cell part, organelle part and macromolecular complex), five molecular function categories (structural molecule activity, binding, transporter activity, nucleic acid binding transcription factor activity and catalytic activity), and eight biological process categories (cellular process, metabolic process, biological regulation, regulation of biological process, establishment of localization, cellular component organization or biogenesis, response to stimulus and localization). Additionally, a significant number of genes were expressed in only one inbred line or absent in both inbred lines. Comparison of the differences of modes of gene action between previous studies and the present study revealed only a small number of different genes had the same modes of gene action in both maize seedlings and ear inflorescences. This might be an indication that in different tissues or developmental stages, different global expression patterns prevail, which might nevertheless be related to heterosis. Our results support the hypotheses that multiple molecular mechanisms (dominance and overdominance modes) contribute to heterosis.


Introduction
Heterosis, or hybrid vigor, refers to the phenomenon in which progeny of two inbred lines (hybrids) exhibit enhanced agronomic performance such as biomass production, growth rate, fertility, and disease resistance relative to both parents [1]. Heterosis has been extensively used in agriculture, especially in the breeding of maize and rice. For example, it is estimated that approximately 95% of the United States maize acreage and 55% of rice acreage in China are planted with hybrids. Furthermore, hybrid maize technology for large-scale production has a yield advantage of 15% over the elite inbred varieties. Though the concept of heterosis has been introduced over 100 years ago and different genetic models considered [1][2][3], the genetic and molecular basis of heterosis remains elusive.
Three classical genetic hypotheses to explain heterosis have been proposed: the dominance, overdominance, and epistasis hypotheses. The dominance hypothesis states that deleterious recessive alleles cause inbreeding depression. A cross of two inbred parents will benefit from complementation of these deleterious alleles and will display a superior phenotype [2][3][4]. The over-dominance hypothesis refers to allelic interactions at one or multiple heterozygous genes resulting in superior trait expression compared to the better parent [1,5]. Both hypotheses may be insufficient to explain the molecular mechanism for heterosis [6]. Epistasis hypothesis refers to interactions of alleles at different loci from two parents in F 1 hybrids, leading to heterosis [7,8].
Maize immature ear inflorescences show heterosis in ear architectural traits [40], and significant positive correlations between grain yield and ear architectural traits, such as ear length, kernel row number, number of kernels per row, kernel number density, and cob diameter have been reported [40], which suggests that each of these components contributes to greater yields. A thorough knowledge of the genes affecting the various components and their interactions will facilitate our understanding of the molecular basis of heterosis of grain yield. In this study, we applied a highly effective approach of high throughout deep sequencing to identify genes, which are highly expressed in maize elite inbred lines of B73 and Mo17, and their F 1 hybrid (B73 × Mo17) at an early stage of ear inflorescences development. We provide first molecular evidence that regulatory mechanisms underlying the phenomenon of heterosis are very early active in maize ear development. Furthermore, we also found that differentially expressed genes between hybrids and their parents can be involved in certain regulatory networks, which suggested that complicated gene networks might be underlying heterosis. Results of the present study might help promote further understanding of mechanisms underlying heterosis.

Statistics and Analysis of Library Sequencing
Here, we sequenced three ear digital gene expression (DGE) libraries from two inbred parents (B73 and Mo17) and their F 1 hybrid (B73 × Mo17) using massively parallel sequencing on the Illumina platform at BGI-Shenzhen, China (Tables 1 and S1). A total of approximately 4.2 million raw tags per library with 259,890 distinct tag sequences were identified. After data-processing steps (see Materials and Methods), the total number of filtered, high-quality clean tags was almost the same in three libraries. Furthermore, the F 1 library had the highest number of distinct tags (259,282), followed by the B73 (242,184), and Mo17 (239,963) libraries (Table 1 and S1). These distinct tags and their genomic frequency as well as the raw data were deposited in NCBI Sequence Read Archive (SRA) database with the accession number (PRJNA248701). Copy numbers of most of the distinct tags (over 77%) ranged from 1-5. However, a small number of distinct tags (less than 3.29%) with a frequency higher than 100 make up over 62% of all clean tags in all three libraries ( Figure S1).
When sequencing depths reach 1 million total tags, the number of novel distinct tags discovered dropped dramatically in all three libraries ( Figure S2A). From that point, increasing sequence depth results in a slow and stable accumulation of new distinct tags indicating that sequencing has reached saturation. Moreover, as shown in Figure 2B, when the total tag number in B73 reached 1 million, the increase of identified genes started to level out, and stabilized when the number of tags reached 3 million. Next, the level of gene expression was determined by calculating the number of unambiguous clean tags for each gene and then normalized to the number of transcripts per million tags (TPM). The Mo17 and F 1 data showed a similar trend ( Figure S2). This suggests that only few more distinct genes would be identified when the total clean tag number reached a certain value.

Mapping Tags to the Maize Reference Genome
We used SOAP2 software [41] to map all distinct tags to the maize reference genome (B73 RefGen_v2) [42]. Mapping results showed that 51.85%, 55.49% and 58.96%, respectively, of distinct clean tags mapped to the reference database (sense or anti-sense), and 46.02%, 49.02%, and 52.29%, respectively, of the distinct clean tags mapped unambiguously to the reference genes (Table 1 and S1). Out these, we identified 20,784-21,372 genes expressed in the three genotype comparison (Table S2), our analysis revealed 17,128 genes expressed in all samples and 22,789 genes expressed collectively in the present study (Table S3). In total, 32.84%, 29.19%, and 25.43% of all distinct clean tags for F 1 , B73, and Mo17 data sets, respectively, did not map to the reference maize genome sequence or associated transcripts (Tables 1 and S1). These non-mapped tags most likely represent regions where the maize reference sequence is incomplete [43] or there are differential mRNA processing events for most maize genes, such as alternative splicing [44]. Only 0.03% of non-mapped tags matched maize chloroplast or mitochondrial genome sequences (Table S1). Because Solexa sequencing can distinguish transcripts originating from both DNA strands, we found evidence for bidirectional transcription in 14,012-14,420 of all detectable overlapping genes and 1031-1119 antisense-stand specific transcripts based on the strand-specific nature of the sequenced tags (Table S2). By comparison, the ratio of sense to antisense strand of the transcripts was approximately 1.3:1 for all libraries. As summarized in Figure S3 and Table S3, most expressed genes (approx. 19,965) are represented in fewer than one hundred copies and only a small proportion of genes are highly expressed (Table S3). We map the clean tags, which cannot be mapped to mRNA, mitochondria and chloroplasts, to the whole genome, providing start positions that could be uniquely mapped by those tags (Table S4).

Different Gene Expression Analysis
A total of 32,973 significantly changed tags entities were detected among three genotypes (see Materials and Methods) ( Figure S4A,C,E and Table S5). Then, the processed DGE data were used to determine different gene expression between inbred parents or between parental lines and their F 1 hybrid and we identified 8621 out of 22,789 genes that were differentially expressed among B73, Mo17, and their hybrid, representing 38% of the ear digital gene expression ( Figure 1A). Among them, 3401 (3344) and 2226 (2189) genes were significantly up-regulated and down-regulated, respectively, in the F 1 compared to B73/Mo17 ( Figure 1B, Figure S4B,D and Table S6). The comparison between parental line libraries also revealed significant variation in expression. A total of 3883 genes, including 1896 up-regulated and 1987 down-regulated genes, were identified in B73 compared to Mo17 ( Figure 1B, Figure S4F and Table S6). We further investigated the mode of gene action for these different genes (Table S6). 8.9% (767 of 8621) exhibited an expression pattern that was not distinguishable from additivity, while the other 91.1% (7854 of 8621) genes showed non-additive expression patterns (Figure 2A). The non-additive

Number of DEGs
differentially expressed genes from the cross were further classified into four distinct classes: high-parent dominance (1984), low-parent dominance (2559), over-dominance (1085), and under-dominance (1963) ( Figure 2B and Table 2). A sample of 30 differentially expressed genes was randomly selected for validation by qRT-PCR. The trends in the expression of these genes detected by DGE were consistent (29 genes) or partially (1 gene) consistent with those determined in qRT-PCR analyses ( Figure 3). These findings are consistent with a recent study in maize [31] or rice [16] that supports the involvement of multiple modes (dominance and overdominance) of gene action in association with heterosis.

Functional Enrichment Analysis for Differentially Expressed Genes Using Gene Ontology (GO)
To gain a better understanding of the functional roles of different genes between inbred parents or between parental lines and their F 1 hybrid, we looked for gene enrichment regarding the Gene Ontology (GO) cellular component, molecular function and biological process categories [44,45]. Functional-annotations from the maize genome (http://maizesequence.org) were used for functional classification of the 8621 different genes and performed by the official web-based tools for searching and browsing the Gene Ontology database (AmiGO) (http://www.geneontology.org/) [45] (Figure 2 and Table 3). We found that different genes showing differential expression patterns were mainly enriched in five cellular component categories (organelle, cell, cell part, organelle part and macromolecular complex), five molecular function categories (structural molecule activity, binding, transporter activity, nucleic acid binding transcription factor activity and catalytic activity), and eight biological process categories (cellular process, metabolic process, biological regulation, regulation of biological process, establishment of localization, cellular component organization or biogenesis, response to stimulus and localization) ( Table 3). We further classified different genes in more detail based on Gene Ontology (GO) by AmiGO [45] (Tables S7-S9). These different genes in GO annotation analysis of cellular component showed significant enrichment in the following components: non-membrane-bound organelle (6.40 × 10 −4 ), intracellular non-membrane-bound organelle (6.40 × 10 −4 ), nucleus (2.56 × 10 −3 ), ribosome (3.34 × 10 −2 ), and macromolecular complex (4.27 × 10 −2 ). In the GO biological process enrichment analysis, there are five significant GO terms (FDR corrected p-value < 0.05). Eighty percent were macromolecular synthesis-related, such as F 1 cellular protein metabolic process, cellular macromolecule metabolic process, protein metabolic process, and macromolecule metabolic process (Table 4).

The Expression Patterns of Biological Macromolecular Synthesis-Related Genes in the B73 × Mo17 Hybrid
To investigate pathways in which different genes were involved and enriched, pathway analysis was performed using the KOBAS 2.0 web tool [46]. 2935 out of 5627 different genes were involved in 122 pathways between F 1 and B73, 2877 out of 5533 different genes were involved in 121 pathways between F 1 and Mo17 (Table S10), Among them, the majority were present in ribosome, spliceosome and various metabolic pathways (such as arginine and proline metabolism, purine metabolism, glycolysis/gluconeogenesis, and biosynthesis of alkaloids derived from histidine and purine, etc.) (Table S10). q-value estimates [47] revealed only four pathways showing significance (q < 0.05): ribosome (q = 1.58 × 10 −3 ), arginine and proline metabolism (q = 1.58 × 10 −3 ), spliceosome (q = 1.63 × 10 −2 ) and pyruvate metabolism (q = 3.80 × 10 −2 ) (Tables 5 and S11), and the first two pathways showed extreme significance (q < 0.01). Of the top 10 differentially expressed genes enriched in pathways between parental lines and their F 1 hybrid, it should be noted that the purine metabolism pathway was the third largest pathway, following the ribosome pathway, although it does not reach the significance level (Tables 5 and S11).   [47]; the top 10 pathways with q-value are listed; * pathway with q-value < 0.05 is considered as significant; ** pathway with q-value < 0.01 is considered as extreme significant.
In the four significant pathways, 188 different genes (4.68%) were detected in the ribosome pathway, 46 different genes (1.11%) in the arginine and proline metabolism pathway, 194 (4.82%) in the spliceosome pathway and 49 (1.69%) in the pyruvate metabolism pathway. Approximately 50% of different genes in every pathway were differentially expressed between inbred and hybrid ( Table 5). In the two significant metabolic pathways, almost all (202 out of 204) different genes showed non-additive expression patterns (HPD, UDO, LPD, or UDO). Seventeen of the 46 different genes in the arginine and proline metabolism pathway were up-regulated (with expression pattern of HPD or ODO) in the F 1 hybrid, and 26 of 46 were down-regulated (with expression pattern of LPD or UDO). As compared with parental lines, the transcriptional level of one different gene in the F 1 hybrid (GRMZM2G062142) encoding ornithine carbamoyltransferase (OTC) was up-regulated by 7-and 11.6-fold, respectively. The gene (GRMZM2G169458) encoding fatty aldehyde dehydrogenase 1 and another gene (GRMZM2G366392) encoding S-adenosylmethionine decarboxylase proenzyme (SAMDC) showed 4.7-/3.4-fold and 2.1-/2.9-fold higher expression in F 1 hybrid than the parental lines, respectively. Simultaneously, the expression levels of another different gene in the F 1 hybrid (GRMZM2G035042) encoding IMP dehydrogenase/GMP reductase were down-regulated by 8.5-and 9.3-fold, respectively, comparing with parental lines, another gene (GRMZM2G010406) encoding argininosuccinate synthase showed significant down-regulation expression in F 1 hybrid. Moreover, 50 different genes involved in the pyruvate metabolism included genes encoding glyoxylatereductase (GRMZM2G166899), phosphoenolpyruvate carboxylase1 (PEPCase 1, GRMZM2G083841), oxidoreductase (GRMZM2G118770), and Pyruvate, orthophosphate dikinase 1 (GRMZM2G097457). (Table S8).

Resolving Transcription Factors (TFs) among Differentially Expressed Genes
A primary objective was to identify genes that encode TFs and to determine their modes of gene action. To test this, we retrieved putative orthologs of maize genes in our differently expressed data based on information from the Ensembl Compara gene trees [50] at maizesequence.org and gramene.org [51]. We then queried known Arabidopsis TFs in the database of Arabidopsis Transcription Factors (http://datf.cbi.pku.edu.cn/) and identified 435 maize different genes with sequence similarities to Arabidopsis TFs between parental lines and their F 1 hybrid libraries (Table  S12). We further interrogated these different genes using gene ontologies, InterPro domains, and known maize annotations. Of the 435 putative TFs, 11 exhibited additive gene action, and the majority of these TFs (n = 424) detected in this study exhibited nonadditive modes of gene expression. Most of these genes exhibited low-parent dominance (n = 138), high-parent dominance (n = 112), and underdominance (n = 112). However, overdominance (n = 42) and other gene action (n = 19) were also observed.
We also surveyed the differential expression of TFs across a wide range of transcript abundance in hybrid here and a mutant in RAMOSA3 (RA3) gene reported in a previous study [52], the latter of which showed an increased branching phenotype resulting from a loss of determinacy of basal spikelet pair meristems. A total of 39 differentially expressed putative TFs in our dataset were identified, and exhibited all nonadditive modes of gene expression (Table 6). Moreover, these TFs were also differentially expressed over a wide range of abundances in ra3 mutants [52]. These TFs includes several kinds of members of TF families associated with functions in development and meristem maintenance or identity (NAC, YABBY, GRAS and TCP), while others have roles in hormone-mediated or stress-mediated signaling by auxin (AUX/IAA), brassinosteroids (BES), or ethylene and stress (AP2/ERF). Among the 42 differentially expressed TFs, eight were characterized as AP2/ERF family proteins (Table 6). Therefore, these TFs possibly not only contribute to heterosis, but also provide insight into genetic control of branching.

A Significant Number of Genes Were Expressed in Only One Inbred Line or Absent in both Inbred Lines
A total of 5660 (24.8%) genes with no detectable expression in one inbred line or both two inbred lines were identified (Table 7). In other word, there are a significant number of genes that displayed presence-absence expression patterns. Among these genes, 46.4% (2624 of 5660) genes exhibited an expression pattern that was present in B73, and absent in Mo17, while some other 34.8% (1971 out of 5660) genes exhibited a similar expression pattern that present in Mo17, and absent in B73 (Tables 7  and S13). Moreover, these genes expressed in only one inbred line also displayed presence-absence expression patterns in their hybrid. Specifically, regarding the genes present in B73 and absent in Mo17, the ratio of genes present to genes absent in their hybrid was approximately 0.56:1, while the ratio of genes present to genes absent in their hybrid was approximately 1.31:1 for other genes present in Mo17 and absent in B73. In combination, these results suggested that the majority of genes that were expressed in only one inbred line exhibited parental effects on gene expression levels, and presence-absence expression patterns of some genes may be related with presence/absence variations in maize genes [53] (See the following analysis). Surprisingly, it was found that 18.8% (1065 of 5660) genes were not expressed in the two inbred parents. However, these genes were expressed in their hybrid (Tables 7 and S13). Additionally, nine genes expressed in only one inbred line were also found in previous study by Stupar et al. [14], and these results in two studies were consistent. For instance, GRMZM2G037255 (corresponding accession #CF629797) and GRMZM2G152258 (accession #BM073080) absent in both inbred lines in our study were not detected in both B73 and Mo17 by PCR in the study of Stupar et al. [14] (Table S13). However, the other 106 genes were expressed in only one inbred in Stupar's study [14], and were expressed in both B73 and Mo17 in present study.

Analysis of Presence/Absence Variation (PAVs) Genes Action in Maize Hybrids
Presence/absence variations (PAVs) have been described in maize genes [53], and most of the PAVs reflected true differences in gene content between the B73 and Mo17 genomes in recent research. Because of the availability of a complete list of PAV genes identified by Lai et al. [54], the B73 × Mo17 cross was first examined regarding PAV genes between the two parental lines and its relation to different gene expression among the three genotypes. As shown in Table S14, there were 104 PAV genes between the two parental lines, and only 37 PAV genes (35.6%) were mapped by identified tags in three libraries. Interestingly, most of these mapped genes products were hypothetical proteins, and the other 55.8% (58 of 104) of genes were unknown proteins except for nine genes including terpene synthase and ZCN20 (Table S14). Of these mapped genes, 45.9% (17 of 37) genes were expressed in only B73. Six genes expressed in only one inbred line were also not expressed in their hybrid. What was more puzzling is there were three genes which were not expressed in the two inbred parents, nevertheless, expressed in their hybrid. Interestingly, among the 37 PAV genes, the expression levels (TPM) of 30% genes in F 1 hybrid was the same as that of B73 or Mo17, and 50% showed higher or lower expression in F 1 hybrid than both the parental lines (Table S15). However, only five PAV genes were identified with significantly differential expression with HPD or LPD (Table S14).

Discussion
In this study, we assayed genome-wide patterns of gene expression of the maize ear at an early flower differentiation stage among two maize elite inbred lines (B73 and Mo17) and their F 1 hybrid (B73 × Mo17) using Solexa's digital gene expression (DGE) system, a tag-based novel high-throughput transcriptome deep sequencing method. Given the nature of the DGE system, we have pooled biological replicates from three varieties for each group to make representative samples for deep sequencing analysis. We obtained a sequencing depth of approximately 4.2 million tags per library (Table 1) and found 22,789 genes expressed collectively except for putative new transcripts found in the study. Levels of some genes not expressed in the present study were responsive to abiotic stress, e.g., aquaporin PIP1-6 gene (GRMZM2G136032), MADS-box transcription factor 14 gene (GRMZM2G137510), beta-fructofuranosidase 1 precursor gene (GRMZM2G394450), etc., were induced by heavy metal Pb treatment in one recent study [55]. Interestingly, we found evidence for bidirectional transcription in all datasets. By comparison of all libraries, the ratio of sense to antisense strands of the transcripts was approximately 1.3:1, which suggested that not only a high number of antisense expressions, but also the transcriptional regulation in the young ear development acted most strongly on the sense strand. A similar observation was also reported in a recent study [56]. Furthermore, approximately 20% of genes with no detectable expression in one inbred line were identified. These genes may exhibit parental effects on gene expression levels.
Based on our digital gene expression analysis, approximately 37.8% of genes exhibited differential expression between every two genotypes in the B73 × Mo17 hybrid cross. QPCR validation was performed both on the same pooled material that was used for deep sequencing and on independent RNA extractions from each sample, and almost all confirmed the direction of change detected by DGE analysis (Figure 3). These different genes exhibited additive and non-additive expression patterns ( Figure 2). A small fraction of DGs (8.9%, 767 of 8621) exhibited a mode of gene action that was indistinguishable from additivity, which is similar to recent studies in maize [14,23,25,31]. Several studies reported more nonadditively expressed genes, including many with F 1 expression levels outside the parental range [11][12][13]23,27]. Among those nonadditively expressed genes, the proportion of genes with clear over and under-dominant gene action were 12.6% and 22.8%, respectively, which is similar to results from prior studies [12,15,16,31,57]. Additionally, we further compared the differences of modes of gene action between Swanson-Wagner's study [31] and the present study. A total of 8621 different gene BLAST searches were performed using 1367 ESTs as queries. Four hundred and twenty-one different genes were matched by 547 ESTs with high-scoring segment pairs (Table S15), but only a small number of these genes (75/421) had the same mode of gene action between the two studies. Thus, different global expression patterns in different tissues or developmental stages might prevail. Our results support that multiple molecular mechanisms (dominance and overdominance modes) contribute to heterosis, which is consistent with previous reports [21,31].
No consensus has yet been reached about the genetic basis of heterosis [58]. However, some mechanisms were supported by the observations that sequence polymorphism in promoter alleles between inbred lines preferentially occurred in those differentially transcribed genes [16]. When two alleles are exposed to a common trans-acting factor, cis-elements in hybrids might differentially interact with gene regulators, resulting in allele-specific gene expression [14,16,59]. This undoubtedly is one of the causes of gene-expression changes in hybrids. There was evidence that phenotypes in hybrids resulted from the dosage effect of such regulatory genes [60]. In this study, we found 424 putative TFs, exhibiting differential expression in the hybrid compared with either parent, in agreement with recent studies [16]. Remarkably, about 9.2% of these TFs were also differentially expressed in ra3 mutants ( Table 6) and many of the differentially expressed genes that could be mapped onto metabolic pathways were associated with primary carbohydrate biosynthesis and degradation, respiration, and energy production as well as redox and nitrogen cycling processes [52].
In conclusion, the expression of TFs in maize hybrids might be important for allele-specific gene expression in heterosis. Another important finding in a recent study is that many SNPs, indel polymorphisms (IDPs) and PAV genes identified between the B73 and Mo17 genome [54] were consistent with the occurrence of insertion/deletion (indel) variants in 5' regions between the alleles of genes that are differentially expressed in different rice strains [16]. We also analyzed the expression of all 104 PAV genes between the two parental lines by DGE data. Of these genes, 67 PAV genes (64.4%) did not express in both B73 and Mo17, possibly because these genes were nonfunctional. For another 37 PAV genes, most (33/37) of their gene products were hypothetical proteins, unknown proteins, or no significant BLAST hits were obtained by using an e-value cutoff of 1 × 10 −5 , and almost 50% of these genes were expressed in both inbred lines and their hybrid. Specially, eight PAV genes were expressed only in one inbred line and their hybrid, which possibly contributed to heterosis, because the phenomenon conformed to the previous assumption that inbred lines with large differences in gene content could complement one another [54]. Moreover, the expression levels of approximately 30% (9 of 37) of these genes in the F 1 hybrid was the same as that of B73 or Mo17, and 50% showed higher or lower expression in F 1 hybrid than both the parental lines. Interestingly, only five PAV genes were identified with significantly differential expression in the hybrid. It suggests that the expression of only a part of PAV genes was consistent with the complementation hypothesis. In conclusion, it is unlikely that heterosis is the result of any single mechanism [58,61].

Plant Growth and RNA Isolation
The hybrid corn combination, B73 × Mo17 (F 1 and its parental lines, B73 and Mo17), was originally obtained from Thomas Lübberstedt (Iowa State University, Ames, IA, USA). The materials were offered for high-throughput sequencing and quantitative real time PCR (qRT-PCR) analysis. The inbred lines were cultivated at the experimental station of Sichuan Agricultural University, Chengdu, for seed propagation of the inbreds, and for production of B73 × Mo17 hybrid seed. Kernels of the combination were grown in environmentally controlled growth chambers that provided 15 h of light (25 °C) and 9 h of dark (20 °C) as described previously. Light intensity was ≈650-800 μmol·m −2 ·s −1 . Ears of 10 random healthy plants at early flower differentiation stages were manually collected as a pool for each genotype as described previously [62]. In brief, ears at stage 3 (three stamen primordia and one pistil primordium can be observed) were collected according to the plant features (number of visible leaves, leaf age index, number of unfolded and folded leaves) combined with microscopic observation. At stage 3, the number of visible leaves, leaf age index, and the number of unfolded and folded leaves are 19, 65%, 13.7 and 5, respectively. Morphological observations of ears at stage 3 in two inbred lines and their F 1 were almost the same as that in another inbred line X178 reported previously [62]. After separately grinding meristems in liquid nitrogen, RNA were extracted for constructing three digital gene expression libraries, and quantitative real-time PCR validation from ≈5 g of frozen tissue by using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions.

Digital Gene Expression Library Preparation and Sequencing
Tag library preparation for the three genotypes (B73, Mo17 and their F 1 ) was performed in parallel by using the Illumina gene expression sample preparation kit as described previously [63]. An extract of 6 μg of total RNA was obtained and treated with Oligo (dT) magnetic bead adsorption to purify mRNA. Oligo (dT) was then used as a primer to synthesize the first-and second-strand cDNA. The 5' ends of tags were generated by two endonucleases NlaIII or DpnII. The bead-bound cDNA was subsequently digested with restriction enzyme NlaIII, which recognizes and cuts off the CATG sites. The fragments apart from the 3' cDNA fragments connected to Oligo (dT) beads were washed away and the Illumina adaptor 1 was ligated to the sticky 5' end of the digested bead-bound cDNA fragments. The junction of Illumina adaptor 1 and CATG site is the recognition site of MmeI, which is a type of Endonuclease with separate recognition and digestion sites. It cuts 17 bp downstream of the CATG site, producing tags with adaptor 1. After removing 3' fragments with magnetic beads by precipitation, Illumina adaptor 2 was ligated to the 3' ends of tags, acquiring tags with different adaptors at both ends to form a tag library. After 15 cycles of linear PCR amplification, 95 bp fragments were purified by 6% TBE PAGE Gel electrophoresis. After denaturation, the single-chain molecules were fixed onto the Illumina Sequencing flowcell. Each molecule grows into a single-molecule cluster sequencing template through in situ amplification. The four types of nucleotides were labeled by four colors, and added to perform sequencing by synthesis (SBS) [64]. Each tunnel will generate millions of raw reads with sequencing length of 35 bp.

Quantitative RT-PCR and Gene Expression Analysis
In order to verify a sample of genes that exhibited statistically significant differential expression in the analysis of DGE data, we used quantitative real time PCR analysis. The RNA samples used for the qRT-PCR assays were the same as for the DGE experiments from 10 biological replicates. First, 1 μg of RNA was treated with RNase-free DNase (Promega, Madison, WI, USA), and cDNA was synthesized with PrimeScript RT reagent kit (TaKaRa, Tokyo, Japan). Then, qRT-PCR of 30 differentially expressed genes (Table S16), which were involved in arginine and proline metabolism, pyruvate metabolism or other genes, were performed using the SYBR PremixExTaq™ protocol (TaKaRa, Tokyo, Japan) on an Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). For each sample, measurements were performed in triplicate, and the average cycle thresholds (C t ) were used to determine fold-change. 18S rRNA was employed as an endogenous control. The results were calculated using the 2 −∆∆Ct method [65].

Analysis and Mapping of Digital Gene Expression Tags
Raw sequencing image data were transformed by base calling into sequences. These raw data reads were stored in FASTQ format, and their analysis conducted as described by Qin et al. [63]. In brief, prior to mapping to the reference database, all sequences were filtered to trim the 3' adaptor sequence, filter empty tags (reads with only 3' adaptor sequences but no tags) and low-quality tags containing Ns, and remove tags which are too long or too short. A virtual library containing all possible CATG + 17 base-length sequences of the maize genome database (AGPv2, release 5b.60) [43] was utilized. All clean tags were mapped to the reference sequences and a mismatch of only 1 bp was considered. Clean tags that were mapped to the maize genome reference sequences from multiple genes were filtered. The remaining clean tags were designed as unambiguous clean tags. The expression level of each gene was estimated by the frequency of clean tags and then normalized to TPM (number of transcripts per million clean tags) [66], which is a standard method and extensively used in DGE analysis [67]. KOG functional classification, Gene Ontology (GO), pathway annotation and enrichment analyses were based on the NCBI COG (http://www.ncbi.nlm.nih.gov/COG) [68], Gene Ontology Database (http://www.geneontology.org/) [69] and KEGG pathway (http://www.genome.jp/kegg/) [70], respectively. When we investigate pathways in which DGs were involved and enriched, q-value was used for aided identification according to the previous description [47].

Identification of Differentially Expressed Genes
To examine differential expression across samples (B73, Mo17 and their hybrid), the number of raw clean tags in each library was normalized to TPM to obtain normalized gene expression levels. Detection of different tags across samples were performed as previously described [71]. The false discovery rate (5%) is controlled by the Benjamini and Hochberg's procedure [72]. After multiple testing between pairwise comparisons, we use "FDR ≤ 0.001 and the absolute value of log 2 Ratio ≥ 1" as the threshold to judge the significance of gene expression difference. More stringent criteria with smaller FDR and bigger fold-change values were used to identify different genes.
In the present study, the same strategy was performed in a linear-in-genotype contrast when F 1 genotype was compared to the two parental lines as described by Zhang et al. [16]. The genes with "FDR ≤ 0.001 and the absolute value of log 2 Ratio ≥ 1" were regarded as non-additivity, when the genes with "FDR > 0.001 and the absolute value of log 2 Ratio < 1" were regarded as not statistically significantly different from additivity. To classify the genes further, the high-parent dominant genes and low-parent dominant genes were identified from the non-additive group based on the criterion, that the F 1 genotype was significantly different from one parent and not significantly different from another parent. From the non-additive group, the expression of genes was identified as over-or under-dominant, when expression in the F 1 genotype was significantly higher or lower than in both inbred parents, respectively.

Conclusions
Our analysis revealed 17,128 genes expressed in these three genotypes and 22,789 genes expressed collectively in the present study. Approximately 38% of the genes were differentially expressed in early maize ear inflorescences from heterotic cross, including many transcription factor genes and some presence/absence variations (PAVs) genes, and exhibited multiple modes of gene action. Additionally, a significant number of genes were expressed in only one inbred line or absent in both inbred lines. Comparison of the differences of modes of gene action between previous studies and the present study revealed only small number of different genes had the same modes of gene action in both maize seedlings and ear inflorescences, it might of be an indication that, in different tissues or developmental stages, different global expression patterns might prevail, which might nevertheless be related to heterosis. Our results support the hypothesis that multiple molecular mechanisms (dominance and overdominance modes) contribute to heterosis.