Genome-Wide Identification of Height-Related Genes Using Three Maize Dwarfs and RNA-Seq

: Plant height is an important grain yield-associated trait in maize. To date, few genes related to plant height have been characterized in maize. To better understand the genetic mechanisms of plant height in maize, we revealed the transcriptional changes of three dwarf mutants compared to the wild type. By ethyl methane sulfonate treatment of the wild-type maize cultivar PH6WC, we obtained three dwarfs—PH6WCdwarf1 (pd1), PH6WCdwarf2 (pd2), and PH6WCdwarf3 (pd3)—and their plant heights were reduced by 42%, 38%, and 24%, respectively. RNA-Seq data suggested that 1641 differentially expressed genes (DEGs) overlapped with each other among the three dwarfs at the seedling stage. Further analysis showed that the DEGs were divided into four groups with different expression patterns. Functional analysis revealed that these DEGs were commonly enriched in 47 GO terms mainly involved in cytokinesis, hormone, and energy metabolism pathways. Among them, An1 , involved in the GA biosynthesis pathway, and mutations in An1 result in reduced plant height. EREB182 encodes ethylene-responsive element binding protein 2, which is critical for internode elongation. Microtubule-related genes Zmtub2 , Zmtub3 , Zmtub5 , Zmtub6 , and TUBG2 were commonly enriched among the three comparisons. Previous studies have shown that mutations in microtubule-associated genes cause the dwarf phenotype. However, nearly half of the common DEGs had no functional information, such as Zm00001d000107 , Zm00001d000279 , etc., implying their novel and specific functions in maize. Overall, this study identifies several potential plant height-related genes and contributes to linking genetic resources with maize breeding.


Introduction
Maize (Zea mays) is one of the most important crops worldwide.In recent decades, the increase in maize yield per unit area has been attributed to increased planting density rather than an improvement in yield potential per plant [1].However, high-density planting often leads to a tendency toward lodging, which may decrease maize yields [2].Plant height is an important grain yield-associated trait in maize and other crops.Reducing plant height is the main strategy used by breeders to prevent individual plants from lodging.During the "Green Revolution" in the late 1960s, the introduction of semidwarfing traits in breading programs produced high-yield varieties [3], of which short and strong stems confer the capacity for lodging resistance under high fertilizer input [4,5].
Plant height is a complex quantitative trait and is controlled by polygenes with minor effects [6].Unremitting efforts were made by researchers to explore the genetic mechanism of plant height.Several typical plant height-related genes have been cloned and well characterized, including na1, dwarf8(d8), d9, d3, br2, an1, and ZmGA3ox2.Among them, D8 and d9 encode DELLA proteins, which are key genes in the gibberellin GA signal transduction cascade that act as negative regulators of the GA response [7].The Anther ear1 (An1) gene product is involved in the cyclization of geranylgeranyl pyrophosphate to ent-kaurene, the first tetracyclic intermediate in the GA biosynthetic pathway [8].D3, a maize gene identified by transposon tagging, encodes a cytochrome P450-mediated step early in GA biosynthesis [9].The maize brachytic2 (br2) gene, which belongs to the MDR class of P-glycoproteins, influences polar auxin transport [10].ZmGA3ox2 encoding a GA3 β-hydroxylase was positionally cloned and was indicated to be significantly associated with plant height through association mapping [11].To date, numerous QTLs related to plant height have been identified by using linkage mapping and genome-wide association analysis, and more than 40 genes involved in various biosynthesis pathways have been cloned [12].However, identification of these genes is not sufficient for the identification of the molecular genetic mechanisms of plant height.Compared with Arabidopsis thaliana and Oryza sativa, the genetic and molecular mechanisms of maize plant height development remain poorly defined.Therefore, it is imperative that we explore plant height-related genes at the transcriptome level to provide more evidence for identifying the genetic mechanism of plant height development.
RNA sequencing (RNA-Seq) is a robust and convenient technology to uncover potential functional genes and gene regulatory mechanisms [13,14].In recent years, impressive advances have been made in high-throughput gene expression studies by RNA-Seq.For instance, Zhou et al. used RNA-Seq technologies to examine transcriptome changes between two inbred maize lines with distinct seed sizes [15].Naoumkina et al. used RNA-Seq to identify common fiber elongation-related genes in developing fibers of Li1 and Li2 mutants.The results provide new insights into the osmoregulation of short fiber mutants and the role of aquaporins in cotton fiber elongation [16].Here, to better understand the genetic mechanisms of plant height in maize, three dwarf mutants with the same genetic background were selected and defined as PH6WCdwarf1(pd1), PH6WCdwarf2(pd2), and PH6WCdwarf3(pd3).Samples were collected at the seedling stage.RNA-Seq was performed to determine plant height-related genes commonly expressed in three dwarf mutants.We found many differentially expressed genes common to the three dwarfs, and from those, the majority were unannotated genes, resulting in the identification of new plant height genes.These research results contribute to our understanding of the mechanisms underlying plant height regulation.

Plant Materials and Phenotype Evaluation
The three dwarf maize mutants, pd1, pd2, and pd3, were derived from the same background (PH6WC) by ethyl methane sulfonate (EMS) treatment.All experimental materials, including the three dwarfs and PH6WC, were grown in a field on the farm of Shanxi Academy of Agricultural Sciences in Yuci, Jinzhong, China (112.75 • E, 37.68 • N), with a distance of 30 cm between rows and 25 cm between individuals.Field management followed normal agricultural practices.The middle ten plants in each row were used for data collection.Plant and ear height were measured as the height to the tassel top and the height to the base of the ear node, respectively.Leaf area was calculated by the formula leaf length × leaf width × 0.75.Ears were harvested at physiological maturity for phenotyping, and ear length and ear width were measured using standard procedures.The internode number was calculated from the number of internodes between the top brace root internode and the top of the plant, including the tassel [17].All phenotypic data were determined as the mean of measurements from ten individual plants.

RNA-Seq Library Construction and Sequencing
Young stems from the wild-type and mutant plants at the V3 stage were harvested for RNA preparation.Two biological replicates were used per genotype in this study.Total RNA was extracted using an RNAiso Plus Kit (TaKaRa, Otsu, Japan).An Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) was applied to detect the purity and integrity of the total RNA.mRNA samples were purified using an oligo (dT)-coupled magnetic bead system.Enriched mRNA was fragmented and then transcribed into cDNA using random primers with M-MuLV Reverse Transcriptase.Blunt-ended cDNA was generated by end repair and phosphorylation, and adenine and adaptors were added to generate ligation products.The ligation products were purified and amplified by PCR to create the final library.Illumina-based RNA sequencing was performed on an Illumina HiSeqTM 4000 platform (Personalbio, Shanghai, China).

GO and KEGG Enrichment Analysis of DEGs
Gene ontology (GO) enrichment was conducted using the R software (v.4.3.0)package topGO.All the genes with a given GO category in the B73 genome were considered the background for calculating the significance [21].GO terms with a corrected p-value ≤ 0.05 were considered significantly enriched among the DEGs.Significant pathway enrichment analysis was performed using KOBAS 2.0, and a false discovery rate (FDR)-adjusted p-value ≤ 0.05 indicated that the DEGs were significantly enriched in the pathway.

Validation of DEGs by Quantitative Real-Time PCR
Nine genes were randomly selected for qRT-PCR analysis to verify the RNA-Seq results.Gene-specific primers were designed using Primer Premier 6.0, and qRT-PCR was performed on a StepOne Plus real-time PCR machine (Applied Biosystems, Foster City, CA, USA).Each reaction contained 2.0 µL diluted cDNA, 10 µL SYBR Green qPCR Master Mix and 10 mM each of forward and reverse gene primers in a final volume of 20 µL.The reaction was performed with three technical and three biological replicates.The thermal cycling conditions were as follows: 3 min at 95 • C, followed by 45 cycles for 3 s at 95 • C, 30 s at the annealing temperature, and 30 s at 60 • C. The relative expression levels of the candidate genes were determined with the 2 −△△CT method [22].

Characterization of Three Dwarf Mutants in Maize
Three dwarf mutants-pd1, pd2, and pd3-showed 42%, 38%, and 24% reductions in plant height, respectively, compared with wild-type PH6WC, especially pd1, which had a height that was less than half that of the wild-type plant (201.8 cm) (Figure 1).The dwarf phenotype of these mutants was attributed to both shortened internodes and decreased internode numbers.In addition to plant height, other traits of the mutants were also reduced significantly, including ear height, ear length, ear width, and leaf area (Table 1).Among the three mutants, pd1 displayed the most severe developmental abnormalities, such as dark green, shorter, broader and shriveled leaves; delayed maturity; and long and thickened husks.The pd2 and pd3 mutants with normally developed tassels were fertile under open-pollination conditions.In contrast, it was difficult to self-pollinate the pd1 plant due to its long and thickened husks, which seriously delayed the silk emergence of the dwarf mutant.
ties, such as dark green, shorter, broader and shriveled leaves; delayed maturity; and lo and thickened husks.The pd2 and pd3 mutants with normally developed tassels w fertile under open-pollination conditions.In contrast, it was difficult to self-pollinate t pd1 plant due to its long and thickened husks, which seriously delayed the silk emergen of the dwarf mutant.

Illumina Sequencing and Sequence Read Analysis
A total of 348,456,350 paired-end reads were produced for PH6WC, pd1, pd2, a pd3, which generated 49 Gb of data.These reads had a Q20 percentage of over 96% (p centage of sequences with a sequencing error rate lower than 1%).After removing lo quality reads, adapters, ambiguous reads, and duplicates, 346,420,920 sequences w successfully obtained from the 8 libraries, with 40-46 million reads per library.Of the high-quality reads, 96.04% were uniquely aligned to the maize B73 reference geno

Illumina Sequencing and Sequence Read Analysis
A total of 348,456,350 paired-end reads were produced for PH6WC, pd1, pd2, and pd3, which generated 49 Gb of data.These reads had a Q20 percentage of over 96% (percentage of sequences with a sequencing error rate lower than 1%).After removing low-quality reads, adapters, ambiguous reads, and duplicates, 346,420,920 sequences were successfully obtained from the 8 libraries, with 40-46 million reads per library.Of these high-quality reads, 96.04% were uniquely aligned to the maize B73 reference genome (Zea_mays.AGPv3.25).The percentage of Illumina sequencing reads mapped to exon regions was 98.21%.Principal component analysis revealed that the eight samples could be clustered into four groups: WT, pd1, pd2, and pd3 (Figure 2A).The correlations of each pair of biological replicates were evaluated using Pearson's correlation tests.The correlation coefficient (r) of biological replicates ranged from 0.93 to 0.97, which demonstrated a high level of reproducibility of the RNA expression patterns (Figure 2B).
(Zea_mays.AGPv3.25).The percentage of Illumina sequencing reads mapped to exon regions was 98.21%.Principal component analysis revealed that the eight samples could be clustered into four groups: WT, pd1, pd2, and pd3 (Figure 2A).The correlations of each pair of biological replicates were evaluated using Pearson's correlation tests.The correlation coefficient (r) of biological replicates ranged from 0.93 to 0.97,which demonstrated a high level of reproducibility of the RNA expression patterns (Figure 2B).

Global Transcriptome Changes in the Three Dwarf Mutants
To identify significantly differentially expressed genes (DEGs) in the three dwarf mutants-pd1, pd2, and pd3-relative to the wild type (PH6WC), comparisons were performed based on FPKMs (fold change ≥ 2 and p value ≤ 0.05).A total of 4055 (1751 upregulated and 2304 downregulated), 5273 (2130 upregulated and 3143 downregulated), and 8515 (3914 upregulated and 4601 downregulated) DEGs were detected for pd1 vs. WT pd2 vs. WT, and pd3 vs. WT, respectively.Then, all of the DEGs from the three pairs were pooled, redundancies were removed, and a total of 11,917 DEGs were identified.Among them, 1641 DEGs were commonly differentially expressed across the three comparisons (Figure 3).Based on the expression patterns, the 1641 common DEGs were divided into nine clusters (cluster 1-cluster 9), and these clusters were then classified into four groups (G1-G4).Genes in G1 (cluster 1) were consistently upregulated or downregulated in the three dwarf mutants, and hsp22 (GRMZM2G007729), encoding a low-molecular-weight heat shock protein precursor, was downregulated in this group.In G2 (clusters 2-3), most of the genes were downregulated in pd3 and upregulated in pd1 and pd2.The genes brittle stalk 2 precursor (bk2, Zm00001d047276), beta-carotene hydroxylase (HYD1 Zm00001d002589), and proline-rich cell wall protein 1 precursor (prcw1, Zm00001d024027) were in this group.Clusters 4-7 formed G3, in which all genes were downregulated in the three dwarf mutants.Hormone-related protein-encoding genes were identified in this group, such as genes encoding GA receptor GID1L2 and brassinosteroid insensitive 1-associated receptor kinase 1 precursor.Most of the genes in G4 were inconsistently changed in the three comparisons.The expression patterns were opposite those of the G2 genes.The auxin response factor 4-encoding gene ARF4 (Zm00001d001945) and ent-kaurene biosynthetic gene an1 (Zm00001d032961) were in this group (Figure 3).

Global Transcriptome Changes in the Three Dwarf Mutants
To identify significantly differentially expressed genes (DEGs) in the three dwarf mutants-pd1, pd2, and pd3-relative to the wild type (PH6WC), comparisons were performed based on FPKMs (fold change ≥ 2 and p value ≤ 0.05).A total of 4055 (1751 upregulated and 2304 downregulated), 5273 (2130 upregulated and 3143 downregulated), and 8515 (3914 upregulated and 4601 downregulated) DEGs were detected for pd1 vs. WT, pd2 vs. WT, and pd3 vs. WT, respectively.Then, all of the DEGs from the three pairs were pooled, redundancies were removed, and a total of 11,917 DEGs were identified.Among them, 1641 DEGs were commonly differentially expressed across the three comparisons (Figure 3).Based on the expression patterns, the 1641 common DEGs were divided into nine clusters (cluster 1-cluster 9), and these clusters were then classified into four groups (G1-G4).Genes in G1 (cluster 1) were consistently upregulated or downregulated in the three dwarf mutants, and hsp22 (GRMZM2G007729), encoding a low-molecular-weight heat shock protein precursor, was downregulated in this group.In G2 (clusters 2-3), most of the genes were downregulated in pd3 and upregulated in pd1 and pd2.The genes brittle stalk 2 precursor (bk2, Zm00001d047276), beta-carotene hydroxylase (HYD1, Zm00001d002589), and proline-rich cell wall protein 1 precursor (prcw1, Zm00001d024027) were in this group.Clusters 4-7 formed G3, in which all genes were downregulated in the three dwarf mutants.Hormone-related protein-encoding genes were identified in this group, such as genes encoding GA receptor GID1L2 and brassinosteroid insensitive 1associated receptor kinase 1 precursor.Most of the genes in G4 were inconsistently changed in the three comparisons.The expression patterns were opposite those of the G2 genes.The auxin response factor 4-encoding gene ARF4 (Zm00001d001945) and ent-kaurene biosynthetic gene an1 (Zm00001d032961) were in this group (Figure 3).

GO Enrichment Analysis of the Common DEGs
GO enrichment analysis was performed, and 47 commonly enriched GO terms were identified across the three comparisons.In the biological process category, the three enriched GO terms were cell division related, including cell cycle (GO:0007049), regulation of mitotic cell cycle (GO:0007346), and cell cycle checkpoint (GO:0000075).Three GO terms were photosynthesis related, including light harvesting (GO:0009765), light harvesting in photosystem I (GO:0009768), and light reaction (GO:0019684).For molecular function, three GO terms were commonly enriched, among which two were related to hydrolase activity (GO:0004553, GO:0016798) and another was related to chlorophyll binding.For cellular components, diverse GO terms were commonly enriched.The majority of them were involved in plastids, including plastids (GO:0009536), plastid stroma (GO:0009532), chloroplaststroma (GO:0009570), and chloroplast thylakoids (GO:0009534).GO terms relevant to the photosystem were also commonly enriched (Table 2).

GO Enrichment Analysis of the Common DEGs
GO enrichment analysis was performed, and 47 commonly enriched GO terms were identified across the three comparisons.In the biological process category, the three enriched GO terms were cell division related, including cell cycle (GO:0007049), regulation of mitotic cell cycle (GO:0007346), and cell cycle checkpoint (GO:0000075).Three GO terms were photosynthesis related, including light harvesting (GO:0009765), light harvesting in photosystem I (GO:0009768), and light reaction (GO:0019684).For molecular function, three GO terms were commonly enriched, among which two were related to hydrolase activity (GO:0004553, GO:0016798) and another was related to chlorophyll binding.For cellular components, diverse GO terms were commonly enriched.The majority of them were involved in plastids, including plastids (GO:0009536), plastid stroma (GO:0009532), chloroplaststroma (GO:0009570), and chloroplast thylakoids (GO:0009534).GO terms relevant to the photosystem were also commonly enriched (Table 2).

Validation of RNA-Seq Data by qRT-PCR
To validate the accuracy of the RNA-Seq analysis, nine DEGs were further analyzed by qRT-PCR: Zm00001d040423, Zm00001d027511, Zm00001d034897, Zm00001d009399, Zm00001d043400, Zm00001d050400, Zm00001d053513, Zm00001d015889 and Zm00001d010740.The results showed that the expression patterns of these nine genes were similar to those measured by RNA-Seq (Figure 4).The RNA-Seq results showed high relevance to the qRT-PCR results (Pearson's r = 0.77, 0.85, and 0.86), indicating the reliability of the RNA-Seq data.

Discussion
Plant height, a key component of plant architecture due to its relatedness to lodging resistance, has been demonstrated to be controlled by multiple qualitative and quantitative genes.However, few plant height genes have been well characterized in maize.Further study is required to uncover the molecular genetic basis of plant height development.Herein, three dwarf maize mutants induced by EMS and derived from a common wildtype plant, PH6WC, were selected for RNA-Seq at the seedling stage.A total of 8515, 5273, and 4055 DEGs were detected for pd1/WT, pd2/WT, and pd3/WT, respectively.Compared to WT, the largest number of DEGs was identified for pd1 among the three dwarfs, which is consistent with the most severe dwarf phenotype of pd1.A total of 1641 DEGs were commonly differentially expressed across the three comparisons, which are most likely involved in the regulation of plant height.The roles of these DEGs differed and included cytokinesis, cell wall characteristics, hormone metabolism, and energy metabolism, which are all considered to be involved in plant height development.In addition, nearly half of the common DEGs had no functional information.These DEGs were also important candidate genes for plant height.

Experimental Design for the Identification of Plant Height-Related Genes
Dwarf mutants are effectively used in the identification of plant height-related genes in many crops.To be er understand the genetic regulation of plant height in maize, we compared the transcriptomes of three dwarf mutants of maize and their common wild type, PH6WC.All three dwarfs-pd1, pd2, and pd3-were EMS-induced mutants with the same genetic background.Due to their significantly different plant heights, the phenotypes of pd1, pd2, and pd3 are most likely caused by different types of genes.In addition to plant height, other traits changed in the three dwarf plants, such as dark green leaves and delayed maturity.These altered phenotypes of the three dwarf lines show in-

Discussion
Plant height, a key component of plant architecture due to its relatedness to lodging resistance, has been demonstrated to be controlled by multiple qualitative and quantitative genes.However, few plant height genes have been well characterized in maize.Further study is required to uncover the molecular genetic basis of plant height development.Herein, three dwarf maize mutants induced by EMS and derived from a common wild-type plant, PH6WC, were selected for RNA-Seq at the seedling stage.A total of 8515, 5273, and 4055 DEGs were detected for pd1/WT, pd2/WT, and pd3/WT, respectively.Compared to WT, the largest number of DEGs was identified for pd1 among the three dwarfs, which is consistent with the most severe dwarf phenotype of pd1.A total of 1641 DEGs were commonly differentially expressed across the three comparisons, which are most likely involved in the regulation of plant height.The roles of these DEGs differed and included cytokinesis, cell wall characteristics, hormone metabolism, and energy metabolism, which are all considered to be involved in plant height development.In addition, nearly half of the common DEGs had no functional information.These DEGs were also important candidate genes for plant height.

Experimental Design for the Identification of Plant Height-Related Genes
Dwarf mutants are effectively used in the identification of plant height-related genes in many crops.To better understand the genetic regulation of plant height in maize, we compared the transcriptomes of three dwarf mutants of maize and their common wild type, PH6WC.All three dwarfs-pd1, pd2, and pd3-were EMS-induced mutants with the same genetic background.Due to their significantly different plant heights, the phenotypes of pd1, pd2, and pd3 are most likely caused by different types of genes.In addition to plant height, other traits changed in the three dwarf plants, such as dark green leaves and delayed maturity.These altered phenotypes of the three dwarf lines show interruptions to different parts of a biological metabolic process, but all cause a shorter stature.Therefore, many DEGs in the three dwarfs may be a consequence of chain reactions in response to adverse effects of the causative mutation and are not necessarily directly related to plant height.To reduce the noise in the data, common DEGs across the pd1/WT, pd2/WT, and pd3/WT comparisons at the seedling stage can be selected to identify transcripts directly related to plant height regardless of the far downstream effects of the causative mutation.

Cytokinesis-Related Genes Involved in Plant Height
Maize plant height was mainly determined by stem elongation driven by cell division and expansion within the internodes.Here, three cytokinesis-related GO terms, namely, cell cycle (GO:0007049), regulation of mitotic cell cycle (GO:0007346), and cell cycle checkpoint (GO:0000075), were commonly enriched among the three comparisons.Microtubulerelated DEGs, such as Zmtub2 (Zm00001d033850), Zmtub3 (Zm00001d013159), Zmtub5 (Zm00001d006651), Zmtub6 (Zm00001d021891), and TUBG2 (Zm00001d037582), were also enriched.Microtubules play an essential role in plant cell division, which is dependent on their organization and dynamics and regulated by a suite of diverse microtubuleassociated proteins [23].Previous studies have shown that mutations related to microtubuleassociated proteins cause dwarf phenotypes.For instance, one point mutation in an αtubulin gene created the tid1-1 mutant, showing dwarfism and twisted growth in rice [24].A single base pair change in the α-tubulin 1 gene (EtTUA1) in tef (Eragrostis tef ) resulted in semidwarfism and lodging tolerance [25].In addition, numerous cyclin-related DEGs (e.g., Zm00001d043164 and Zm00001d010656) were enriched in these GO terms.Cyclin genes play key roles in regulating the cell cycle during plant development.Moreover, genes involved in cell wall-related processes might also function in cell division and expansion during stem elongation.In the 1461 genes commonly differentially expressed in the three dwarfs, some cell wall-related genes were identified, such as ZmBK2 (Zm00001d047276) and three genes encoding probable xyloglucan glycosyltransferase (Zm00001d012666, Zm00001d020392, and Zm00001d030103).ZmBK2 is involved in cell wall biogenesis and affects the mechanical strength of maize stalks [26].Altering the xyloglucan side chain structure in Arabidopsis resulted in mutants that were dwarfed with curled rosette leaves and short petioles [27].

DEGs Related to Energy Metabolism
Photosynthesis is the fundamental process for biomass accumulation in plants.Sugars produced by photosynthesis serve as building blocks for the growth of plants.Moreover, sugars also act as signal molecules interconnected with plant hormones to finely regulate the main biological processes throughout the whole plant life cycle [28].In this study, unexpectedly, numerous genes related to photosynthesis (e.g., photosystems and chlorophyll-binding proteins) were observed among the 1641 common DEGs.These genes all fell into the same group (G3), in which all genes were downregulated in the three dwarfs.Moreover, GO categories related to photosynthesis were significantly overrepresented among the commonly enriched GO terms.Previous reports showed that defects in photosynthesis-related genes may cause the dwarf phenotype in plants.For instance, the OsPS1-Fencoded photosystem 1-F subunit mutant showed striking phenotypic changes, including a reduction in plant height [29].In addition, genes involved in height reduction have various effects on photosynthesis [30].GA also plays a positive role in enhancing photosynthetic activity.Overexpression of CcGA20ox1 causes an overall upregulation of genes encoding photosystem-related proteins and chlorophyll-binding proteins [31].Accordingly, we speculated that overrepresentation of DEGs related to photosynthesis was most likely caused by the altered GA levels in the three dwarfs.Nevertheless, their functions in plant height development need further identification.

DEGs Related to Hormone Metabolism Pathways
Among the 1461 common DEGs, some candidate genes were annotated and related to GA biosynthesis and signaling pathways.For instance, Anther ear1 (Zm00001d032961) encodes an enzyme involved in the synthesis of ent-kaurene, the first step in the GA biosynthesis pathway.Mutations in An1 result in a GA-responsive phenotype, including reduced plant height, delayed maturity, and the development of perfect flowers on normally pistillate ears [9].Four cytochrome P450-encoding genes (Zm00001d032042, Zm00001d024946, Zm00001d038511, and Zm00001d029290) were found, and previous studies indicated that cytochrome P450 could comodulate plant height by influencing cell elongation [32].GID1L2 (Zm00001d050837), encoding a GA receptor protein, plays a key role in the GA signal transduction pathway [33].In addition to GA, the expression of other phytohormone-related genes was also commonly altered in the three dwarfs.EREB182 (Zm00001d027929) encodes ethylene-responsive element binding protein 2. Ethylene is critical for growth and development.For example, ethylene induces the expression of SNORKEL1 and SNORKEL2 and then triggers remarkable internode elongation via GA in rice [34].Zm00001d037297, encoding a brassinosteroid insensitive 1-associated receptor kinase 1, was significantly downregulated in the three dwarfs.Defects in BR biosynthesis in maize cause a pleiotropic effect on plant development, including dwarf phenotypes [32].

Conclusions
We identified 1641 differentially expressed genes (DEGs) overlapped with each other among the three dwarfs at the seedling stage.These DEGs were divided into four groups with different expression patterns.Functional analysis revealed that these DEGs were commonly enriched in 47 GO terms mainly involved in cytokinesis, hormone, and energy metabolism pathways.Nearly half of the common DEGs had no functional information, implying their novel and specific functions in maize.Overall, this study identifies several potential plant height-related genes and contributes to linking genetic resources with maize breeding.

Figure 2 .
Figure 2. Evaluation of RNA-Seq data quality: (A) principal component analysis of the 8 sample datasets; (B) Correlation test of 8 sample datasets.

Figure 2 .
Figure 2. Evaluation of RNA-Seq data quality: (A) principal component analysis of the 8 sample datasets; (B) Correlation test of 8 sample datasets.

Figure 3 .
Figure 3. DEGs identified in the three dwarfs: (A) Venn diagram of the DEGs; (B) heat maps of 1641 DEGs common across the three transcriptome comparisons.

Figure 3 .
Figure 3. DEGs identified in the three dwarfs: (A) Venn diagram of the DEGs; (B) heat maps of 1641 DEGs common across the three transcriptome comparisons.

Figure 4 .
Figure 4. qRT-PCR validation of nine differentially expressed genes in the three dwarfs identified by RNA-Seq.The results for each gene are based on three biological and three technical replicates.The error bars indicate standard errors.

Figure 4 .
Figure 4. qRT-PCR validation of nine differentially expressed genes in the three dwarfs identified by RNA-Seq.The results for each gene are based on three biological and three technical replicates.The error bars indicate standard errors.

Table 1 .
Comparisons of phenotypic traits between the three dwarfs and the wild-type PH6WC.

Table 1 .
Comparisons of phenotypic traits between the three dwarfs and the wild-type PH6WC.

Table 2 .
Commonly enriched GO terms across the three comparisons.

Table 2 .
Commonly enriched GO terms across the three comparisons.