Combined Linkage and Association Mapping Reveals QTL and Candidate Genes for Plant and Ear Height in Maize

Plant height (PH) and ear height (EH) are two very important agronomic traits related to the population density and lodging in maize. In order to better understand of the genetic basis of nature variation in PH and EH, two bi-parental populations and one genome-wide association study (GWAS) population were used to map quantitative trait loci (QTL) for both traits. Phenotypic data analysis revealed a wide normal distribution and high heritability for PH and EH in the three populations, which indicated that maize height is a highly polygenic trait. A total of 21 QTL for PH and EH in three common genomic regions (bin 1.05, 5.04/05, and 6.04/05) were identified by QTL mapping in the two bi-parental populations under multiple environments. Additionally, 41 single nucleotide polymorphisms (SNPs) were identified for PH and EH by GWAS, of which 29 SNPs were located in 19 unique candidate gene regions. Most of the candidate genes were related to plant growth and development. One QTL on Chromosome 1 was further verified in a near-isogenic line (NIL) population, and GWAS identified a C2H2 zinc finger family protein that maybe the candidate gene for this QTL. These results revealed that nature variation of PH and EH are strongly controlled by multiple genes with low effect and facilitated a better understanding of the underlying mechanism of height in maize.


INTRODUCTION
Maize (Zea mays L.) is one of the most cultivated crops worldwide, owing to its versatility and wide adaptability, and serves as food, animal feed, and raw material for various industrial products. Plant density is considered an important factor to increase yield (Cardwell, 1982;Zhang et al., 2014). Plant height (PH) and ear height (EH) are two important traits that directly affect plant density (Cai et al., 2012). A better understanding of the genetic architecture of PH and EH could help breeders to develop varieties with high density tolerance and, consequently, increased yield.
The mapping of quantitative trait loci (QTL) based on bi-parental populations is an effective method for detecting the genetic variation of in crops. Molecular markers have been widely used in QTL mapping for various traits in maize, including maturity (Buckler et al., 2009), disease resistance (Zwonitzer et al., 2010;Chen et al., 2012;Ding et al., 2012), abiotic stress (Agrama and Moussa, 1996), plant morphology (Xu et al., 2009;Ding et al., 2015;Chen et al., 2016), and plant color (Chen et al., 2014). More than 280 QTL have been identified for PH and EH in maize (Cai et al., 2012); Sibov et al. (2003) detected four QTL for PH and five QTL for EH that explained 24.76 and 20.91% of the phenotypic variance, respectively, using an F 2:3 population with 400 families; Tang et al. (2006) detected six QTL for PH, of which four QTL were located in the same genomic region of the QTL affecting average internode length using a recombinant inbred line (RIL) population and suggested that average internode length was the main contributor to PH in maize; and Wang et al. (2006) detected 127 QTL for PH based on the maize integrated QTL map and confirmed 40 of them using meta-analysis. Additionally, Peiffer et al. (2014) showed that PH in maize has a highly polygenic genetic architecture using nestedassociation mapping population. Therefore, multiple models of genetic architecture, differing in polygenicity and effect sizes, can plausibly explain population variation for height in maize. All the previous studies provided valuable information to determine the genetic basis of height in maize, but only few QTL were validated in different populations. Therefore, the underlying mechanism of plant development remains unclear.
The available single nucleotide polymorphisms (SNPs) (Elshire et al., 2011;Chia et al., 2012) have been widely used in genome-wide association study (GWAS) for discovering loci controlling various traits in maize, including leaf architecture (Tian et al., 2011), disease and insect resistance (Zila et al., 2014;Samayoa et al., 2015;Li et al., 2016;Mahuku et al., 2016), and other important agronomic traits (Yang et al., 2014). Compared with traditional QTL mapping, GWAS has higher resolution and allows the identification of genes for multiple traits in a single population (Breseghello and Sorrells, 2006;Yu and Buckler, 2006;Yan et al., 2011). However, population substructure and genotyping errors can result in high false positive rates (Andersen et al., 2005;. Therefore, the traditional QTL mapping powerfully compares pairs of alleles with low resolution, whereas association analysis provides a highresolution evaluation of numerous alleles with uneven statistical power (Wilson et al., 2004). In this study, we combined linkage mapping and GWAS to identify QTL for PH and EH and exploit the strengths and weaknesses of both approaches for discovering causal loci across the genome.
Previous studies have cloned several genes that strongly influence PH in maize, such as brachytic2 on Chromosome (Chr.) 1, which influences polar auxin transport (Multani et al., 2003;Xing et al., 2015); ZmGA3ox2 on Chr. 3 and dwaft3 on Chr. 9 related to gibberellin synthesis (Winkler and Helentjaris, 1995;Teng et al., 2013); and dwarf8 on Chr. 1 and dwarf9 on Chr. 5 involved in gibberellin signal transduction pathways (Thornsberry et al., 2001;Lawit et al., 2010). In maize, over 40 genes have been reported for PH that are related to different biosynthesis pathways such as hormone synthesis, transport, and signaling (Wang and Li, 2008;Peiffer et al., 2014). However, the underlying mechanism of plant development remains unclear. The identification of additional loci/genes related to PH in elite inbred maize lines might help to better understand the mechanism of height development and optimize selection in breeding programs.
Here, we evaluated two bi-parental populations and one GWAS population for PH and EH in different environments and used both linkage mapping and GWAS to identify and verify QTL and candidate genes associated with the natural variation of height in maize. Our objectives were to: (1) identify QTL that significantly affect PH and EH in the two bi-parental populations; (2) discover top SNPs and their candidate genes associated with PH and EH by GWAS; (3) verify QTL in a near-isogenic line (NIL) populations; and (4) discuss the function of candidate genes.

Germplasm and Phenotype Evaluation
An F 2:3 bi-parental population, consisting of 225 families were developed from a cross between Zheng58 and Chang7-2, which are the parental lines of Zhengdan958, one of the most popular maize hybrids in China with an annual acreage of 4.0 × 10 6 ha, accounting for 20% of total maize production (Du et al., 2006). Parental line Chang7-2 has higher PH and EH than another parental line Zheng58. The population and their parents were evaluated in Zhengzhou (34 • 52 ′ N 113 • 37 ′ E), Jiyuan (35 • 34 ′ N 112 • 5 ′ E), and Xichang (27 • 32 ′ N 102 • 10 ′ E) in 2007 (07Zhengzhou, 07Jiyuan, 07Xichang) and Zhengzhou and Jiyuan in 20008 (08Zhengzhou, 08Jiyuan). Another bi-parental recombinant inbred line (RIL) population consisting of 250 lines was produced from a cross between BT-1 and N6. BT-1 is an inbred line reformed by tropical Asia material with relatively high PH and EH and high resistance to various diseases, whereas N6 is a Tangsipingtou inbred line with relatively low PH and EH and low resistance to diseases. The population and their parents were evaluated in Zhengzhou in 2007 and 2008 (07Zhengzhou, 08Zhengzhou), Xuchang (34 • 2 ′ N 113 • 81 ′ E) in 2009 (09Xuchang), and Zhengzhou in spring and summer of 2009 (09Zhengzhou1, 09Zhengzhou2).
NIL population was developed from a cross between the low-PH line Zheng58, and the donor line Chang7-2 through four cycles of advanced backcrosses and one generation of selfcrossing using marker assisted selection for the target genomic region. Two NILs, one homozygous for the QTL of the donor parent (+) and another homozygous for the QTL of the recurrent parent (−) were evaluated in Zhengzhou and Wenxian (34 The GWAS population, consisting of 258 maize inbred lines including the parental inbred lines of the two bi-parental populations, from the heterotic populations Tangsipingtou and Reid and some tropical lines from International Maize and Wheat Improve Center (CIMMYT), was evaluated in Zhengzhou, Wenxian and Hainan (18 • 21 ′ N 109 • 10 ′ E) in 2015 (15Zhengzhou, 15Wenxian, and 15Hainan, respectively).
All the populations were arranged in a randomized complete block design with three replications for each location. Sixteen plants were planted in 4 m row plots with 0.67 m row spacing. Field management was performed according to the standard agronomic practices in each location. PH was measured from the soil level to the tip of main inflorescence, and EH was measured from the soil level to the node attachment of the primary ear. The final value of each plot was the mean value of each plant inside the plot.

Phenotypic Data Analysis
Analysis of variance (ANOVA) of phenotype data was performed using the SPSS software (www.spss.com). The components of variance were estimated using a complete random effects model and broad-sense heritability was calculated as defined by Knapp et al. (1985). Best linear unbiased predictions (BLUPs) of the combined PH and EH were calculated by using a mixed linear model (lmer) in R software (R Core Team, 2015), which used replication, environments which combination location and years, and entries as a random effect. The BLUP value of each line was used for GWAS analysis. Histogram of phenotypic data was generated by using hist function in R software (R Core Team, 2015).

Linkage Mapping
The two bi-parental populations were genotyped by simple sequence repeat (SSR) markers. The linkage map of RIL population contained 207 polymorphic SSR markers had a total length of 1820.8 cM and an average distance of 11.7 cM between markers . The F 2:3 population was genotyped with 180 SSR markers, and the constructed map had a total length of 1987.7 cM and an average distance of 11.0 cM between markers .
Composite interval mapping (CIM) for QTL mapping was performed by Windows QTL Cartographer 2.5 (Zeng, 1994). The threshold value was set using 1000 random permutations (Churchill and Doerge, 1994). The proportion of phenotypic variation explained by a single QTL was determined by the square of the partial correlation coefficient (R 2 ).

GWAS and Candidate Gene Annotation
A total of 955,650 SNPs were identified in the GWAS population using the genotyping-by-sequencing (GBS) method (Elshire et al., 2011;Glaubitz et al., 2014) which was conducted in Cornell University and the SNP flanking sequence and position information were available in "panzea" website (http://cbsusrv04. tc.cornell.edu/users/panzea/download.aspx?filegroupid=4). The SNPs was filtered with a missing value greater than 0.25 and a minor allele frequency (MAF) less than 0.05, resulting in 224,152 SNPs for future analysis. The principal component analysis (PCA), Kinship matrix and linkage disequilibrium (LD) between each pair of SNPs were conducted using TASSEL 5.0 software (Bradbury et al., 2007). The population structure was determined using an admixture ancestry model with correlated allele frequency in software STRUCTURE v2.3.3 (Pritchard et al., 2000). STRUCTURE was run four replicates for K (number of subpopulations) = 1-8, and two replicates for K = 9-12 with a run-length of 100,000 repetitions of Markov chain Monte Carlo following a burn-in period of 10,000 iterations. GWAS was conducted using a mixed linear model (MLM) that included BLUPs, markers, Kinship and PCA in TASSEL (Bradbury et al., 2007).
Candidate gene information were obtained from the MaizeGDB (http://www.maizegdb.org/) genome browser based on physical position of significant SNPs. Phytozome database (http://phytozome.jgi.doe.gov/pz/portal.html) were used to define relevant pathways and annotate possible functions of candidate genes.

Meta-Analysis
All QTL detected in the bi-parental populations and SNPs identified in the GWAS population were analyzed based on their physical position using R (R Core Team, 2015). The R code is provided in File S1.

Phenotype Analysis
Descriptive statistics for PH and EH in the three mapping populations are presented in Table 1. There were a wide variations were observed in each population, for example, the PH of combined environments in RIL population ranged from 105 to 273 cm, whereas in the GWAS population from 106 to 293 cm. The EH in the F 2:3 population ranged from 38.8 to 123 cm, and from 21.7 to 147 cm in GWAS population. The mean of PH and EH in the F 2:3 families were 183.6 and 73.7 cm, respectively. The frequency of phenotypic value in all three populations for PH and EH followed an approximately normal distribution ( Figure S1). The genotypic variance (σ 2 g ) and the genotype-byenvironment variance (σ 2 ge ) of PH and EH were significant in all three populations. Heritability for PH was 0.92 in the F 2:3 population, 0.97 in the RIL population, and 0.91 in the GWAS population. Heritability for EH is similar to that for PH ( Table 1). The high repeatability and heritability indicated that much of the phenotypic variance was genetically controlled in the populations and suitable for QTL mapping.

GWAS for PH and EH
Population structure analysis showed that the LnP(D) value was stable between replications when K ≤ 3 and became unstable when at K > 3 (Figure 1A), indicating that the appropriate K-value was 3. At K = 3, 98.2% of the GWAS population was divided into Subpopulations 1, 2, and 3, containing 165, 68, and 60 lines, respectively, whereas the rest 2.8% was grouped into a mixed subpopulation ( Figure 1B). Subgroup 1 was the most diverse and included tropical lines from CIMMYT; Subgroup 2 included lines from the Tangsipingtou heterotic group and Subgroup 3 included lines from the Reid heterotic group. PCA showed the lines source could separate very well by PC1 and PC2 (Figure S4A), and the first three PCs explained most of the variance (Figure S4B). PCA was in agreement with the results generated by STURCTURE, clearly showing the existence of three subgroups in the GWAS population. Hence, the first three PCs were included in the mixed model for GWAS. Single marker-based GWAS was performed using MLM incorporating both the population structure (first three PCs) and K into the model. A total of 21 significant SNPs for PH and 20 for EH was identified with p < 1.0 × 10 −4 (Figure 2), which explained 6.2-9.3% and 6.2-9.54% of the phenotypic variation for PH and EH, respectively. The most significant SNP was S6_124299082 for PH located on Chr. 6 which explained 9.3% of the phenotypic variation in the GWAS population, whereas the most significant SNP was S3_1976523 for EH located on Chr. 3 which explained 8.8% of the phenotypic variation ( Table 3). The quantile-quantile (QQ) plots showed that population structure was well controlled by PCA and K ( Figure S5). Based on the physical position of the significant SNPs in the B73 reference genome, 29 SNPs ware located in 19 unique candidate genes, whereas 12 SNPs in intergenic regions (Table 3). Candidate gene annotation showed that most of the candidate genes function belonged to the kinase family transport and signaling.

Meta-QTL Analysis for PH and EH
In order to identify the consistent loci for PH and EH, metaanalysis was conducted for both QTL mapping and GWAS results. Six common loci were identified on Chr. 1, Chr. 2, Chr. 3, Chr. 5, Chr. 6, and Chr. 9 (Figure 3). The QTL for PH in bin 1.05 was detected in both the two bi-parental populations, but no SNPs with p < 10 × 10 −4 were identified in the GWAS population; the QTL for EH in bin 6.04 was detected in both the bi-parental populations, and one significant SNP was identified in the GWAS population; the QTL for EH and PH in bin 9.04 was detected in both the bi-parental populations, and one significant SNP was identified in the GWAS population. From meta-QTL analysis, six consistent loci were detected in more than two populations, suggesting that may have stable effects at different genetic backgrounds and environmental conditions.

QTL Verification
The qPH1-1 for PH in bin 1.05/06 share the same chromosome region of QTL qEH1-1 for EH, suggesting that they might be the same QTL, controlling both PH and EH. Furthermore, qPH1-1 had the highest LOD value in the F 2:3 population. In order to fine mapping this QTL, NIL populations were developed by using marker assistance selection with flanking markers. The phenotype evaluation of NIL populations showed that positive homozygous allele (homozygous from donor parents) could increase PH by 18 cm and EH by 11 cm compared with the negative homozygous allele (Figures 4A,B), verifying QTL mapping results, and indicating that QTL in bin 1.05 could significantly increase PH and EH.
The Manhattan plot of GWAS ( Figure 4C) showed that two markers with very low p-value, S1_92787084 with p = 7.4 × 10 −4 for PH and p = 7.6 × 10 −3 for EH and S1_92245779 with p = 2.4×10 −3 for PH and p = 1.8×10 −3 for EH, represented two candidate genes, GRMZM2G161459 and GRMZM2G114667, respectively. GRMZM2G161459 encodes the peptide transporter PTR2, which transports a wide spectrum of dipeptides and  is involved in multiple pathways (Chiang et al., 2004;Xia et al., 2008), whereas GRMZM2G114667 encodes the C 2 H 2 zinc finger family protein which is related to plant development and is required for the transcription of a number of genes coding proteins with stress-protective functions (Görner et al., 1998;Iuchi, 2001). Additionally, GRMZM2G114667 is highly expressed in the internode of the plant stem and may relate to the development of internodes. Compare with GRMZM2G161459, GRMZM2G114667 has higher possibility to be the candidate gene of qPH1-1, because of it is highly expressed in the internode, has a similar high effect on both PH and EH, and is highly related to plant development.

Heredity of Height in Maize
Height is one of the most important agronomic traits and the most heritable trait in maize. However, the underlying genetic FIGURE 3 | Meta-quantitative trait locus (QTL) analysis for QTL detected by two bi-parental populations and significant SNPs detected by genome-wide association study (GWAS) in this study. Colored lines represent different QTL or SNPs detected in this study.
mechanisms of natural variation remain unclear. In the present study, the normal distribution and high heritability of PH and EH variation were detected in both the bi-parental and GWAS populations, suggesting that height in maize is a highly polygenic trait. The large number of markers with small effects that detected in GWAS confirmed that height in modern maize inbred lines is a polygenic trait controlled by multiple genes with low effect. These results were in agreement with previous studies (Peiffer et al., 2013(Peiffer et al., , 2014.

Relationship between PH and EH
PH and EH are significantly correlated in maize breeding populations, thus some QTL for these traits are mapped in the same genomic regions due to pleiotropic effects (Tang et al., 2006). In the present study, QTL for both PH and EH were located in the same genomic region such as qPH1-1 for PH and qEH1-1 for EH (bin 1.05/06) in the F 2:3 population as well as qPH3-1 for EH and qEH3-2 for EH (bin 3.04/05) in the RIL population. However, some QTL were only mapped for one of the two traits such as qEH6-2 for EH in the RIL population. These results suggested that the genetic mechanism of PH and EH may be similar, but not identical.

Genome Wide Association Study
GWAS is a powerful approach for exploring the molecular basis of phenotypic variations in plants, and the population is a foundation for GWAS. There were several GWAS populations have been reported for PH and EH research (Weng et al., 2011;Yang et al., 2011;Peiffer et al., 2014). The different population have different features, some populations used global diversity of lines which focus on global variation analysis (Yang et al., 2011;Peiffer et al., 2014), and others may only use region germplasm which focus on region variation analysis. For example, Weng et al. (2011) used a GWAS population contained Chinese elite lines identified a major QTL for height on chromosome 5, which played an important role in Chinese maize breeding programs. In the present study, the GWAS population mainly contained Chinese elite inbred lines and some tropical inbred lines from CIMMYT, it represented most of the diversity of Chinese and tropical germplasm which were the important germplasm in Chinese breeding programs. The Chinese and tropical germplasm were used for GWAS analysis, in order to reveal the genetic basis of height nature variation in these important germplasm in breeding programs. The further analysis on diversity comparison between present GWAS population and previous reported population was profitable to understand the different significant markers detected in different populations and may facilitate a better understanding of the underlying mechanism of height in maize. GWAS has been widely used for functional gene discovery that has yielded a large number of associations between markers and various complex traits (Visscher et al., 2012;Li et al., 2016). However, some markers with high significant p-value may not be genuine whereas others with moderate or even low significant p-values can be considered genuine based on the biological trait (Panagiotou and Ioannidis, 2012). The selection of the most suitable threshold for GWAS using statistical approaches, such as the Bonferroni correction (Bland and Altman, 1995) or the false discovery rate (Benjamini and Hochberg, 1995), is important for the identification of genuine markers for target traits. A fixed cut-off value based on the number of markers and traits, such as p < 10 −5 or p < 10 −4 , provides flexibility and is commonly used in GWAS. In the present study, no SNPs passed the Bonferroni correction threshold, probably because height is a complex traits in modern maize inbred lines controlled by multiple genes with low effect (Peiffer et al., 2014). Therefore, to identify the top SNP associated with PH and EH, we choose a fixed cut-off threshold of p < 10 −4 .

Comparison with QTL Identified in Previous Studies
Previous studies have shown that common QTL can be detected in different environments and populations (Beavis et al., 1991;Chen et al., 2012). In the present study, we detected nine QTL for PH and EH in the F 2:3 families population and 12 QTL for PH and EH in the RIL populations. Of these, three "stable QTL" (detected in both bi-parental populations) were identified for PH and EH on Chr. 1, Chr. 5, and Chr. 6 and also reported in previous studies (Cai et al., 2012;Park et al., 2013;Sa et al., 2014;Wei et al., 2015). These results indicated that those three loci were stable QTL for PH and EH in modern maize inbred lines.

Candidate Genes for Height in Maize
The QTL for EH on Chr. 6 (bin 6.04/05) could explain 4% of the phenotypic variation in the F 2:3 population and 10% of the phenotypic variation in the RIL population. However, no QTL was detected in this genomic region for PH in these two populations. The associated SNP S6_113668398 had the lowest p-value (p = 1.95 × 10 −5 ) and explained 7.3% of the phenotypic variation in the GWAS population for EH. Candidate gene analysis showed that S6_113668398 was located in AC202439.3_FG006, which encodes a urea transmembrane transporter. The homologous genes in Arabidopsis (AtDUR3), and rice (LOC_Os10g42960.1) also encode transporters that play an important role in the uptake and use of urea (Kojima et al., 2007;. Urea is the main source of nitrogen in modern crop production, because of its low cost and easy absorption by the plants, and provides 50% for the total weight of nitrogen fertilizer worldwide (Kojima et al., 2007;. The efficient use of the nitrogen by the plant increases its height (Lafitte and Edmeades, 1994). Therefore, this gene could be considered an important candidate gene for maize height, and further functional verification might help to better understand the underlying genetic mechanism.
The QTL for PH and EH on Chr. 3 could explain 15 and 19% of the phenotypic variance for PH and EH, respectively, in the RIL population. An associated SNP with small pvalue was detected for PH located in the annotated gene GRMZM2G322186. GRMZM2G322186 encodes a pyridoxal phosphate-dependent transferase and also has multiple catalytic functions in various organisms (Percudani and Peracchi, 2003). The protein level of GRMZM2G322186 was highly different between two parental lines BT-1 and N6 of the RIL population, and the N6 showed similar protein sequence with reference genome B73 (Figure S6). Expression analysis showed that GRMZM2G322186 has a high expression level in all the developmental stages of maize, indicating its importance in plant development (Winter et al., 2007;Sekhon et al., 2011). These results suggested that GRMZM2G322186 was the most possible candidate gene for the QTL of PH on Chr. 3.

CONCLUSIONS
In this study, two bi-parental populations and one GWAS population were used to identify and map QTL for PH and EH. A total of 21 QTL for PH and EH was detected by QTL mapping and 41 SNPs were identified by GWAS. A QTL on Chr. 1 was verified in NILs and indicated that a C 2 H 2 zinc finger family protein might be the candidate gene for PH. One candidate gene on Chr. 6 (AC202439.3_FG006), encoding a urea transmembrane transporter, was considered as the candidate gene for EH. The candidate genes for other stable QTL were also discussed. These results confirmed that nature variation of maize height is strongly controlled by multiple genes with low effect and the QTL and candidate genes identified in this study could help to better understand the genetic basis of PH and EH in maize.

AUTHOR CONTRIBUTIONS
JW, JC, JD, and XZ designed, led, and coordinated the overall study. XL, ZZ, YW, BZ, RW, JM, SW, and XZ perform the field experiment. ZZ, XL, and JC carried out the analysis. JC, ZZ, XL, and JW wrote the manuscript.

FUNDING
The work was supported by National Natural Science Foundation of China (Grant No. 31271307).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00833 Figure S1 | Phenotypic distribution of plant height (PH) and ear height (EH) in the F 2:3 population, the recombinant inbred line (RIL) population, and the genome-wide association study (GWAS) population. The density line is shown in black.   File S1 | R code for meta-quantitative trait locus (QTL) analysis.