Genetic dissection of maize plant architecture using a novel nested association mapping population

Abstract

The leaf angle (LA), plant height (PH), and ear height (EH) are key plant architectural traits influencing maize (Zea mays L.) yield. However, their genetic determinants have not yet been well‐characterized. Here, we developed a maize advanced backcross‐nested association mapping population in Henan Agricultural University (HNAU‐NAM1) comprised of 1,625 BC₁F₄/BC₂F₄ lines. These were obtained by crossing a diverse set of 12 representative inbred lines with the common GEMS41 line, which were then genotyped using the MaizeSNP9.4K array. Genetic diversity and phenotypic distribution analyses showed considerable levels of genetic variation. We obtained 18–88 quantitative trait loci (QTLs) associated with LA, PH, and EH by using three complementary mapping methods, named as separate linkage mapping, joint linkage mapping, and genome‐wide association studies. Our analyses enabled the identification of ten QTL hot‐spot regions associated with the three traits, which were distributed on nine different chromosomes. We further selected 13 major QTLs that were simultaneously detected by three methods and deduced the candidate genes, of which eight were not reported before. The newly constructed HNAU‐NAM1 population in this study will further broaden our insights into understanding of genetic regulation of plant architecture, thus will help to improve maize yield and provide an invaluable resource for maize functional genomics and breeding research.

Abbreviations
• BLUP
• best linear unbiased prediction
• Chr
• chromosome
• EH
• ear height
• GWAS
• genome-wide association studies
• HNAU-NAM1
• advanced backcross-nested association mapping population in Henan Agricultural University
• JLM
• joint linkage mapping
• LA
• leaf angle
• LD
• linkage disequilibrium
• LOD
• logarithm of odds
• LRT
• likelihood ratio test
• MAGIC
• multiparent advanced generation intercross populations
• NAM
• nested association mapping
• PH
• plant height
• PVE
• phenotypic variation explained
• QTLs
• quantitative trait loci
• SLM
• separate linkage mapping
• SNPs
• single nucleotide polymorphisms.

INTRODUCTION

Maize (Zea mays L.) is one of the most important crops worldwide, representing a valuable source of bioenergy material and food for the human (Li et al., 2017). The domestication of maize occurred ∼9,000 yr ago from ancient teosinte (Z. mays subspecies parviglumis) and its continuous cultivation over the millennia led to the vast diversity in a number of important architectural traits, including leaf angle (LA), plant height (PH), and ear height (EH) (van Heerwaarden et al., 2011). These traits have been subjected to extensive improvement that significantly increased maize yield (Duvick et al., 2010; Tian et al., 2019). Importantly, the emergency resulting from a shortage of food and energy associated with the continuous growth of the world population (Gerland et al., 2014) make it an urgent priority to develop maize cultivars with enhanced architecture that may result in higher yields.

The three key components contributing to plant architecture, LA, PH, and EH, are complex quantitative traits under the control of multiple QTLs. In general, the development of a segregating population is a commonly used strategy for QTL mapping. In the past, biparental populations, such as recombinant inbred lines, near-isogenic lines, and F₂ lines, were extensively used to dissect the genetic mechanisms behind agronomic traits by linkage mapping in maize (Andorf et al., 2019). Given the fact that limited recombination events occur within such populations, the QTL mapping resolution is quite low and often makes it challenging to identify the causal genes controlling different traits. In contrast, the low levels of linkage disequilibrium (LD) in natural population comprising hundreds of accessions collected worldwide has enabled the QTL mapping to a narrow region or even single gene (Mamidi et al., 2020). However, intrinsic population structure and kinship between diverse accessions can confound association mapping results (Hoffman, 2013).

Over the last decade, multiparental populations have been gradually applied in multiple crop genetic studies, including maize (Dell'Acqua et al., 2015; Xiao et al., 2016; Yu et al., 2008), wheat (Triticum aestivum L.) (Mackay et al., 2014), barley (Hordeum vulgare L.) (Nice et al., 2016; Schnaithmann et al., 2014), rice (Oryza sativa L.) (Meng et al., 2016), sorghum [Sorghum bicolor (L.) Moench] (Mace et al., 2013; Tao et al., 2020), soybean [Glycine max (L.) Merr.] (Song et al., 2017), oilseed (Brassica napus L.) (Hu et al., 2018), and peanut (Arachis hypogaea L.) (Gangurde et al., 2020). These strategies utilized linkage and/or association mapping methods and included different implementations, such as random-open-parent association mapping, nested association mapping (NAM), and multiparent advanced generation intercross (MAGIC) populations. In maize, a total of six NAM-related populations, ranging in size of approximately 600–5,000 lines, have so far been developed (Gage et al., 2020; Li et al., 2021; Morales et al., 2020). Using these populations, over 100 different phenotypes have been characterized, ranging from outer whole-plant agronomic traits to inner kernel ionomics (Gage et al., 2020), which represents a powerful tool to investigate the traits of interest in this crop.

Previous research on maize architectural traits allowed for the identification and mapping of hundreds of QTLs associated with LA, PH, and EH (https://www.maizegdb.org/data_center/phenotype), some of which were identified in multiple populations with different parents (Chen et al., 2019; Peiffer et al., 2014; Tian et al., 2011). However, the underlying genetic mechanisms behind plant architecture have not yet been well-characterized. Hence, further accurate dissection of the QTLs/genes regulating LA, PH, and EH will significantly contribute to a profound understanding of this issue. In the present study, we developed a maize advanced backcross-nested association mapping population in Henan Agricultural University (HNAU-NAM1) population by crossing a GEMS41 inbred line with 12 diverse inbred lines, resulting in 12 subpopulations consisting of 1,625 lines. Using this population, we analyzed the level of genetic diversity and LD decay and then identified QTLs/alleles regulating LA, PH, and EH. After integrating all the mapping results, we also inferred the associated candidate genes in some major QTLs.

MATERIALS AND METHODS Choice of founder lines and population development

Twelve diverse maize inbred lines were selected from the association mapping panel AMP540 (Yang et al., 2011), specifically CML470, P178, CIMBL29, CIMBL83, CML486, CML454, CML304, TY1, CML360, DAN598, K22, and CML496. These lines were selected based on levels of genetic diversity and traits of interest. Then the lines were crossed with the common parent GEMS41, respectively, to establish HNAU-NAM1 (Supplemental Table S1). GEMS41, an inbred line from America, served as a recurrent parent. The HNAU-NAM1 population thus consisted of 12 biparental subpopulations (including five BC₁F₄ and seven BC₂F₄), which were developed as follows: (a) a female parent GEMS41 was crossed with the 12 founders to generate ∼1,800 F₁ plants; (b) F₁ plants of each combination were crossed back with the female parent GEMS41 to produce the corresponding BC₁F₁ progenies; (c) seven subpopulations were further crossed back with the female parent GEMS41 to acquire BC₂F₁ plants; and (d) finally, the plants of each combination were subjected to three rounds of the single-seed-descent selfing strategy to obtain the corresponding BC₁F₄/BC₂F₄ lines (Supplemental Figure S1).

Field trials and phenotypic characterization

The parents and their progeny were grown in experimental fields of the Changge county (34°13′ N, 113°46′ E) in the Henan province of China in the summer season (2017 and 2019), using a randomized complete block design. Every year, two consecutive blocks were used, and 12 subpopulations were randomly assigned to 12 plots per block. In every plot, every line of corresponding subpopulation was randomly planted in a row with ∼10 seeds, which were 25 cm apart. All lines of HNAU-NAM1 population were grown under regular agricultural management practices (weeds were manually removed and pests and diseases were chemically controlled throughout the growing periods). A total of three plant architectural traits were investigated in 13 founders and all of the HNAU-NAM1 lines. To generate reliable phenotypic data, five representative plants were randomly selected from each line in every plot. The LA was then determined for four leaves above the primary ear as the angle of each plant leaf was defined by the stalk below the node of subtending leaf (Mickelsona et al., 2002). The PH and EH of the same five plants were measured using a previously described method (Zhou et al., 2016). Finally, the average of five plants was taken as the corresponding phenotype for each trait.

Core Ideas

• A novel maize HNAU-NAM1 population, composed of 1,625 lines, was developed.
• Population structure and genetic diversity in HNAU-NAM1 population were analyzed.
• Mapping methods of separate linkage mapping, joint linkage mapping, and genome-wide association studies were used to detect QTLs for phenotypic traits.
• Ten QTL hot-spot regions simultaneously associated with three traits were detected.
• Eight new candidate genes and four known genes were deduced in 13 major QTLs.

Both standard ANOVA and broad-sense heritability (H²) calculations were performed using the QTL IciMapping v4.2.53 software (Meng et al., 2015) and applying the formula H² = V_G/(V_G + V_GE /e + V_ε /[e * r]) (Knapp et al., 1985), where e and r represent the number of year and replication, respectively; and V_G, V_GE, and V_ε denote the variance components of the genotype, genotype × year, and residual, respectively. The best linear unbiased prediction (BLUP) value of each trait was obtained for each line of the HNAU-NAM1 population across 2 yr using linear mixed-effects models in the R package lme4 (Bates et al., 2015). Subsequently, the BLUP value of each line was used for downstream analyses, such as phenotypic description, the calculation of Pearson correlation coefficients, and QTL mapping for the three aforementioned traits.

Genotyping data analysis

The genomic DNA (gDNA) from all of the 13 parental founders and the HNAU-NAM1 lines was extracted from fresh leaves using the CTAB method (Murray & Thompson, 1980) and evaluated by agarose gel electrophoresis. Qualified samples were sent to the China Gold Marker Biotechnology Company (Beijing, China) for genotyping using the MaizeSNP9.4K^{Affymetrix®Axiom®} Genotyping Array, including 9,433 maize single nucleotide polymorphisms (SNPs) selected from the Axiom Maize56K SNP Array (Ganal et al., 2011) (Supplemental Figure S2). The SNPs that passed the following strict filtering criteria were used for the construction of single genetic linkage maps in each of the 12 subpopulations: (a) clear physical location on the maize B73 version 4 reference genome (Jiao et al., 2017); (b) homozygous and polymorphic between two parental founders; (c) a missing rate lower than 10%; (d) a heterozygote rate lower than 10%; (e) a segregation ratio of progenies fitting the expected 23:2:7 (corresponding to the number of homozygotes for the GEMS41 parent, heterozygotes and homozygotes for the non-GEMS41 parent, respectively) for the BC₁F₄ or 55:2:7 for the BC₂F₄ populations with a p value >.001 in a chi-squared test.

To improve the mapping resolution in the genome-wide association study (GWAS) analysis, approximately 1.03 million parental lines genotypes were downloaded from http://www.maizego.org. The genotypes of the offspring in each biparental subpopulation were imputed according to the 1.03 million SNPs based on their progenitor's genotype, which was similar to a previous study (Buckler et al., 2009). The imputation method includes two steps. First, low-density markers derived from MaizeSNP9.4K array were used to build the linkage map for each biparental subpopulation using QTL IciMapping v4.2.53 software (Meng et al., 2015), as mentioned above. Then, in accordance with the obtained linkage map, the chromosomal region of the offspring was imputed from the high-density genotype of the corresponding parents. In the crossover regions of the linkage map where the source parents could not be clearly defined, the genotype of the individual was imputed as missing for loci where the parents were polymorphic or imputed as the parental genotype for loci where the two parents contained the same genotype. After imputation, a total of 1.03 million SNP genotypes were obtained across the HNAU-NAM1 population. After filtering sites with a minor allele frequency ≥0.05, a total of 570,380 high-quality SNPs were retained for GWAS analysis (Supplemental Figure S3).

Diversity and linkage disequilibrium analyses

A merged genotyping dataset (in hapmap format) was obtained from http://www.maizego.org/Resources.html. This dataset comprises 540 maize inbred lines and the 13 parental founders of the HNAU-NAM1 population. Approximately 10,000 SNPs were randomly selected from the dataset to construct a phylogenetic tree using the Neighbor-joining method as implemented in TASSEL v5.2.65 (Bradbury et al., 2007), which was then perfected using online tool iTOL (https://itol.embl.de/).

Principal component analysis and genomic LD analyses were performed on SNPs with minor allele frequency ≥ 0.05, which were extracted from the MaizeSNP9.4K BeadChip array using the TASSEL software. The extent of LD as a function of physical distance was determined by considering pairwise LD measures within 10 Mb, separately on each chromosome, and then calculating the average r² in 100 Kb windows using an in-house developed Perl script (https://github.com/dashengzhao/HNAU-NAM1). The LD halving distance was calculated for the whole genome based on the values obtained across all HNAU-NAM1 lines.

QTL mapping by SLM, joint linkage mapping, and GWAS

To perform separate linkage mapping (SLM) analysis on the 12 HNAU-NAM1 subpopulations, separate genetic linkage maps were constructed using the QTL IciMapping v4.2.53 software (Meng et al., 2015) based on quantified, nonredundant SNP markers available from the MaizeSNP9.4K BeadChip array. The map distances were calculated based on Kosambi's mapping function (Kosambi, 1944). QTL analysis on the LA, PH, and EH traits was performed in each subpopulation using composite interval mapping with the Windows QTL Cartographer v2.5 software (Wang et al., 2012). One thousand permutations were performed for each trait within each subpopulation, to determine the threshold level for logarithm of odds (LOD) score for the QTL significance test, which resulted in a range of 2.80–3.28 (α = 0.05). For convenience, the LOD threshold score ≥3.0 was used as the global cut-off point. To avoid overestimating the number of significant QTLs, adjacent peaks within neighboring genetic regions (≤10 cM) with the same effect directions were defined as a single QTL, as previously described (Xiao et al., 2016). Each QTL region was defined by the peak position of the LOD score and its surrounding two-LOD drop regions, as calculated by WinQTLCart (Wang et al., 2012).

After combining all lines of the HNAU-NAM1 population, a joint linkage mapping (JLM) analysis was performed for further QTL identification, using a similar procedure published in previous study (Xiao et al., 2016). Briefly, for each subpopulation, the missing and heterozygous SNPs in the MaizeSNP9.4K array were filled based on the surrounding SNPs genotype and physical position information. Then, for each chromosome of HNAU-NAM1 lines, redundant SNPs were merged into 8,610 nonredundant genetic bins, where no recombination occurred. The average length of these bins was 0.23 ± 2.21 Mb (mean ± SD) (Supplemental Figure S4). Subsequently, all genetic bins from 10 chromosomes were merged into one file, and a marker-based kinship matrix for the population was calculated. Finally, a liner mixed model was fitted, and restricted maximum likelihood was used to evaluate the significance of each recombination bin as Xiao et al. study (2016), where the population mean and the intercept served as fixed effects, and the marker and polygenic effects were treated as random effects. The tested bin contained 13 additive effects corresponding to the parental alleles of the 13 founder lines of the HNAU-NAM1 population. After a permutation test consisting of 500 permutated samples, the likelihood ratio test (LRT) scores threshold was calculated to be 2.76 (Type I error rate of 0.05). Accordingly, the genomic regions with LRT ≥ 2.76 were defined as QTLs.

The GWAS analysis was conducted using a fixed and random model circulating probability unification approach (Liu et al., 2016), implemented in the memory-efficient and parallel-accelerated R package rMVP (Yin et al., 2021). The estimation of the variance components was performed using the default Brent method (Burch & Iyer, 1997). The first five principal components were included in the GWAS model to correct the hidden population structure. Associations with P value lower than the 1 × 10⁻⁵ threshold were classified as significant, in accordance with previous studies (Chen et al., 2019; Kremling et al., 2018).

QTL overlaps and candidate gene prediction

For each trait, the QTL supporting regions identified in the SLM or JLM methods, were obtained based on the physical coordinates of flanking markers. By combining the 500-kb upstream and downstream regions of the significant GWAS SNPs, the overlapping information among all the QTL regions detected by three complementary mapping methods were acquired collectively. To further investigate the genetic basis of the three key traits, the potential candidate genes located within the major QTLs (phenotypic variation explained [PVE] ≥10%) were detected by three methods simultaneously. This was done by integrating the results acquired from aforementioned three mapping methods, the protein information, and the corresponding orthologs annotations in rice and Arabidopsis, which were downloaded from the genome portal of the Department of Energy Joint Genome Institute (Nordberg et al., 2014).

RESULTS Development and genetic diversity of the HNAU-NAM1 population

The newly developed HNAU-NAM1 population is comprised of 1,625 lines, obtained from 12 corresponding subpopulations, each ranged in size from 48 to 210 lines (135.4 ± 47.6 lines; Supplemental Table S1). After removal of 158 low-quality individuals (high heterozygosity or missing rates), a total of 1,467 lines were retained in the resulting HNAU-NAM1 population (Supplemental Figure S5).

The 13 selected founders represented a relatively broad molecular diversity within the AMP504 population (Figure 1a), and also showed high levels of phenotypic variation in the plant architecture traits of interest, specifically LA, PH, and EH (Figure 1b). Principal component analysis based on a total of 6,017 SNPs showed that the HNAU-NAM1 lines closely clustered with the common parent GEMS41, when compared with the other 12 parents in both PC1 and PC2, which explained a total of 9.5% of the molecular variance (Figure 1c). In addition, the whole-genome LD decay distance in the HNAU-NAM1 population was estimated to be approximately 2.59 Mb (Figure 1d).

FIGURE 1. Genetic and phenotypic diversity of a maize advanced backcross-nested association mapping population (HNAU-NAM1) population. (a) Diversity of 13 HNAU-NAM1 founder lines in 540 diverse inbred lines using data reported in Yang et al. (2011). (b) Phenotypic distribution of leaf angle, plant height and ear height. Red dashed line stands for mean value. Colorful arrows represent the 13 founder lines according to legend in Panel a. (c) Outcome of principal component analysis on HNAU-NAM1 population and 13 founder lines by MaizeSNP9.4K BeadChip data. (d) Genomic linkage disequilibrium (LD) decay analysis in whole genome

Phenotypic variation, heritability and correlations analyses

A broad range of phenotypic variation was observed in the HNAU-NAM1 population (Figure 1b), as well as in each of the 12 subpopulations in order to investigate the genetic architecture of the LA, PH and EH traits (Supplemental Figure S6 and Table S1). Specifically, among the HNAU-NAM1 lines, LA values varied from 13.39 to 63.02°, with an average of 27.02°; PH ranged from 114.27 to 216.68 cm, with a mean value of 159.76 cm; and EH fluctuated from 40.88 to 97.26 cm, with a mean value of 62.89 cm (Table 1).

TABLE 1 ANOVA, broad-sense heritability, descriptive statistics, and correlation analysis of plant architecture traits in the HNAU-NAM1 population

	ANOVA						Descriptive statistics		Correlation analysisa
Trait	Block	Genotype	Environment	Genotype × environment	Error	Broad-sense heritability	Range	Mean ± SD	LA	PH	EH
LA	0.79 ***	52.83***	3.01***	7.98***	12.45***	0.88	13.39–63.02	27.02 ± 6.47	–	0.07*	0.21*
PH	26.68***	321.89***	35.41***	59.67***	79.91***	0.86	114.27–216.68	159.76 ± 15.79	0.07*	–	0.73**
EH	12.62***	92.24***	86.23***	30.37***	39.75***	0.78	40.88–97.26	62.89 ± 7.73	0.21*	0.73**	–

Note. ANOVA, analysis of variance; EH, ear height; LA, leaf angle; PH, plant height.

a

The values mean the Pearson correlation coefficient.

^*significant at P < .05. ^**significant atP < .01. *** significant at P < .001.

Broad sense heritability (H²) was estimated to be 0.88 for LA, 0.86 for PH, and 0.78 for EH (Table 1). The ANOVA showed that all these traits were significantly affected by genotype, environment, and the genotype × environment interaction. Among all sources of phenotypic variation, the largest part was linked to genotypes (Table 1). The correlation between PH and EH (r = 0.73, P < .01) was much stronger than those between EH and LA (r = 0.21, P < .01) or between PH and LA (r = 0.07, P < .05) (Table 1). And in each subpopulation, PH and EH showed the most significant positive correlation as well (Supplemental Figure S6).

Linkage map construction from twelve subpopulations

Filtered SNPs and the high-quality individuals from each subpopulation were used to construct linkage maps. Twelve genetic maps were generated with an average length of 1,920.73 ± 262.87 cM, ranging from 1,400.80 cM in CML470 × GEMS41 to 2,235.55 cM in CIMBL29 × GEMS41 (Supplemental Figure S7 and Supplemental Table S2). On average, each linkage map was comprised of 1,031 high-quality SNP markers, and the recombination rate along the chromosome varied from 0.67 to 1.06 cM/Mb among the 12 subpopulations of the HNAU-NAM1 population.

Genetic dissection of LA using three different methods

A total of 41 QTLs associated with LA were detected among the 12 subpopulations by SLM analysis (Figure 2; Supplemental Figure S8 and Supplemental Table S3). The corresponding absolute additive effects values ranged from 1.37 to 4.52^°. Close to half of the QTLs (46.34%) explained more than 10% of the phenotypic variation. The average length of QTLs was approximately 10.84 ± 14.19 Mb. Notably, seven QTL regions located on Chromosome (Chr) 2, 4, 4, 7, 8, 8, and 9, separately, were consistently identified in two subpopulations, and a QTL region located on Chr1 was consistently identified in three subpopulations (Supplemental Table S3). Using JLM analysis, 84 minor-effect QTLs associated with LA were identified, and the average physical distance of these QTLs was 1.45 ± 2.30 Mb (Figure 2; Supplemental Figure S9a and Supplemental Table S4). Moreover, these QTLs showed a narrower mapping region compared with that obtained using the SLM method. Furthermore, GWAS analysis identified a total of 22 SNPs distributed on nine chromosomes, which were significantly associated with LA (Figure 2; Supplemental Figure S10a and Supplemental Table S5). Of these candidate SNPs, the most robust marker-trait association was found at the distal end of Chr1 (CHR1.S_272218482: 276,987,222 bp) with a P value of 6.81 × 10⁻¹⁴. This SNP also showed the maximum additive effect, with a value of 4.42. A subsequent comparative analysis on the mapping results obtained for LA using the three different methods (i.e., SLM, JLM, and GWAS) revealed that 11–24 QTLs were detected by at least two methods, and 10–14 QTLs were identified by all three methods (Figure 2; Supplemental Table S3).

FIGURE 2. Quantitative trait loci (QTL) results for leaf angle in the maize advanced backcross-nested association mapping population . (a) The number of QTLs identified by separate linkage mapping (SLM), joint linkage mapping (JLM), and genome-wide association study (GWAS). (b) Overlaps in QTL results identified by three methods. Numbers in italics and bold represents the count of QTL identified by JLM and GWAS, respectively. (c) Overview of QTL results in 10 chromosomes. Top panel: the colored dots show the significance of genome-wide bins estimated by JLM method, in which LRT means likelihood ratio test; black triangles represent the physical positions of 30 known leaf angle-related genes. Middle panel: colored rectangles represent QTL regions identified by SLM across the 12 biparental subpopulations, and the color density of the rectangles shows the magnitude of the logarithm of the odds (LOD) values. Bottom panel: the colored dots indicate significant single nucleotide polymorphisms (SNPs) identified by GWAS

Genetic dissection of PH by three different methods

In the SLM analysis, we detected 31 QTLs in the 12 subpopulations (Supplemental Figure S8, S11, and Supplemental Table S3), including 14 QTLs (45.16%), explaining more than 10% of the phenotypic variation. The absolute values of the observed additive effects ranged from 3.61 to 8.82 cm for all of the identified QTLs. The average QTL physical length was 11.55 ± 20.21 Mb. While most of the QTLs were found in only one subpopulation, three QTL regions located on Chr3, 4, and 5 were identified in at least two subpopulations (Supplemental Table S3). The JLM analysis revealed a total of 78 QTLs associated with PH that were distributed on all chromosomes of maize (Supplemental Figure S11 and Supplemental Table S4). Among these, the JLM_qPH6-2 QTL located at Chr6: 98.75–100.30 Mb showed the highest significance, with the LRT score of 18.34. The average interval length of these QTLs was approximately 1.18 ± 1.61 Mb, which was about nine-fold shorter than those detected by the SLM method. At the same time, all of the QTLs detected by the JLM method represented small-effect alleles in 13 founder parents, with absolute values lower than 0.29 (Supplemental Figure S9b). The GWAS analysis revealed 23 SNPs scattered on nine chromosomes that were significantly associated with PH (Supplemental Figure S10b, S11, and Supplemental Table S5). Among these SNPs, the strongest correlation was SYN21611, which was located at Chr5: 1,575,005 bp, with a P value of 1.34 × 10⁻¹¹. The PH-associated QTLs detected by the three methods indicated that 10–24 QTLs were identified by at least two methods, while 9–11 QTLs were detected by all three methods (Supplemental Figure S11 and Table S3).

Genetic dissection of EH by three methods

Using SLM method, a total of 26 QTLs with an average region length of 10.42 ± 17.44 Mb and the proportion of PVE varying from 6.87 to 24.64%, were mapped on all chromosomes except for Chr9 (Supplemental Figure S8, S12, and Supplemental Table S3). Specifically, a QTL region located at Chr2: 0.04—3.55 Mb was detected in four subpopulations, and two QTL regions (Chr10: 140.70-149.26 Mb, Chr2: 39.89-61.21 Mb) were both detected in two subpopulations of HNAU-NAM1 population (Supplemental Table S3). The JLM analysis revealed 88 small-effect QTLs with an average interval length of 1.05 ± 1.50 Mb and absolute additive effect values less than 0.48 (Supplemental Figure S9c, S12). Approximately half of these loci (45.45%) were distributed on Chr2, 3, and 10 (Supplemental Table S4). Among these QTLs, the most significant was JLM_qEH2-1, with the LRT score of 42.33 and located in the physical region of 0.02–3.69 Mb on Chr2. In addition, the GWAS results showed a total of 18 SNPs that were significantly associated with EH located across all chromosomes, and with estimated absolute additive effect values ranging from 0.84 to 3.64 (Supplemental Figure S10c, S12, and Table S5). Interestingly, the most significant SNP, CHR2.S_1812565 (P = 5.01 × 10⁻²⁵), was also located at the beginning of Chr2 (at position 1,804,762 bp). When comparing the EH-related QTL regions mapped by the three aforementioned methods, 7–23 QTLs were identified by at least two methods and 6–8 QTLs were detected by all three methods (Supplemental Figure S12 and Table S3).

Analysis on overlapping genomic regions and candidate gene prediction

A total of 18–88 QTLs controlling LA, PH and EH were identified by applying the SLM, JLM and GWAS methods in the HNAU-NAM1 population, and some of these QTLs were clustered in hotspots across the genome (Figure 3; Table 2). The co-localization of these QTLs was further analyzed in order to dissect the genetic basis and relationship between LA, PH and EH. As expected, 10 QTL hot-spot regions distributed on all chromosomes except for Chr6, were observed using 10 Mb sliding windows (with 1 Mb step), of which two were both found on Chr8 (Table 2). Moreover, each of these QTL hot-spot regions harbored QTLs controlling three traits simultaneously, and QTLs for each trait were detected by at least two mapping methods.

FIGURE 3. Summary of quantitative trait loci (QTL) regions on maize genome for leaf angle (LA), plant height (PH) and ear height (EH) by three methods. SLM: separate linkage mapping; JLM: joint linkage mapping; GWAS: genome-wide association study. For GWAS, only significant single nucleotide polymorphisms (SNPs) were showed here. QTL hot-spot regions means that wherein LA, PH and EH simultaneously coincide and QTLs for each trait could be identified by at least two mapping methods. Chr, chromosome; HNAU-NAM1, maize advanced backcross-nested association mapping population

TABLE 2 QTL hot-spot regions wherein LA, PH, and EH simultaneously coincide within 10 Mb and each identified by at least two mapping methods

			LA			PH			EH
QTL region	Chr	Region (Mb)	SLM	JLM	GWAS	SLM	JLM	GWAS	SLM	JLM	GWAS
Hotspot 1	1	6.0–16.0	×	√	√	√	√	√	×	√	√
Hotspot 2	2	0–10.0	√	√	√	√	√	√	√	√	√
Hotspot 3	3	163.0–173.0	√	√	√	√	√	√	√	√	√
Hotspot 4	4	231.0–241.0	√	√	×	√	√	√	√	√	√
Hotspot 5	5	173.0–183.0	√	√	√	√	√	√	√	√	×
Hotspot 7	7	160.0–170.0	√	√	×	√	√	√	√	√	×
Hotspot 8_1	8	10.0–20.0	√	√	√	√	√	×	√	√	√
Hotspot 8_2	8	165.0–175.0	√	√	×	×	√	√	×	√	√
Hotspot 9	9	145.0–155.0	×	√	√	√	√	×	×	√	√
Hotspot 10	10	140.0–150.0	√	√	√	√	√	×	√	√	√

Note. Chr, chromosome; EH, ear height; GWAS, genome-wide association study; JLM, joint linkage mapping; LA, leaf angle; PH, plant height; QTL, quantitative trait loci; SLM, separate linkage mapping. “√” represents the QTL hotspot was detected by corresponding mapping methods, while “×” was not.

Among all QTLs obtained for three plant architecture traits, a large subset of these only exerted small effects on the phenotypic composition (Supplemental Figure S9 and Supplemental Table S3–S5). Given the potential application of these QTLs to improve future molecular assisted selection breeding of maize, and also to get in depth insights into plant architecture traits, 13 major QTL regions were focused because each had a PVE ≥ 10% and were simultaneously detected by the three mapping methods (Table 3). To further deduce the candidate genes, the positions of 87 known genes affecting LA, PH, and EH were compared with these major QTL regions (Supplemental Figure S13 and Table S6). As a result, three LA related genes, including asp1 (Cao et al., 2020), lg1 (Moreno et al., 1997), and lg2 (Walsh et al., 1998), were found to be located within corresponding QTL regions, which likely represent causal genes (Supplemental Figure S14 a–c). As for the gene associated with PH and EH, te1, whose mutant plants exhibited shorter internodes and had decreased statures (Veit et al., 1998), likely represents the causal agent for two traits (Supplemental Figure S14d and e). Moreover, three, two, and three newly candidate genes, potentially associated with LA, PH, and EH, respectively, were deduced in the remaining eight major QTL regions, based on the results of three mapping approaches and gene annotations (Table 3).

TABLE 3 Functional annotation of potential candidate genes in major QTL regions detected by three methods simultaneously

Trait	Chr	Biparental subpopulation	PVE (%)	QTL interval (Mb)	Potential candidate gene	Annotation	Annotation of best hit in Arabidopsis	Annotation of best hit in rice	References/sources
LA	1	CIMBL83/GEMS41	12.7	219.76–225.55	Zm00001d032295	DRE-binding protein 1; AP2/ERF domain	encodes a member of the DREB subfamily A-6 of ERF/AP2 transcription factor family	ethylene-responsive transcription factor, putative, expressed	https://www.maizegdb.org/
	1	CML360/GEMS41; CIMBL29/GEMS41	11.74	272.72–292.63	Zm00001d033981 (asp1)	ATP sulfurylase1	pseudouridine synthase/archaeosine transglycosylase-like family protein	bifunctional 3-phosphoadenosine 5-phosphosulfate synthetase, putative, expressed	(Cao et al., 2020)
	2	CML496/GEMS41; TY1/GEMS41	28.15	2.55–4.67	Zm00001d002005 (lg1)	liguleless1 protein; squamosa binding protein	squamosa promoter binding protein-like 8	OsSPL8 - SBP-box gene family member, expressed	(Moreno et al., 1997)
	3	CIMBL83/GEMS41	10.92	178.90–181.60	Zm00001d042777 (lg2)	liguleless2 protein; bZIP domain protein	bZIP transcription factor family protein	transcription factor, putative, expressed	(Walsh et al., 1998)
	8	CML304/GEMS41	10.76	16.43–21.38	Zm00001d008749	auxin-responsive Aux/IAA family member	indoleacetic acid-induced protein 8	OsIAA15 - Auxin-responsive Aux/IAA gene family member, expressed	https://www.maizegdb.org/
	10	CML470/GEMS41	18.12	143.57–146.90	Zm00001d026322	serine/threonine-protein kinase receptor	protein kinase superfamily protein	tyrosine protein kinase domain containing protein, putative, expressed	https://www.maizegdb.org/
PH	3	CIMBL83/GEMS41	17.58	162.51–170.03	Zm00001d042445 (te1)	protein terminal ear1	terminal EAR1-like 1	AML1, putative, expressed	(Veit et al., 1998)
	5	CML470/GEMS41	10.23	164.56–174.72	Zm00001d016708	cell wall invertase1	cell wall invertase	cell wall invertase 4	https://www.maizegdb.org/
	7	DAN598/GEMS41	15.96	136.69–149.23	Zm00001d021237	plant auxin-responsive proteins	SAUR-like auxin-responsive protein family	phytochrome-interacting factor 4, putative, expressed	https://www.maizegdb.org/
EH	2	TY1/GEMS41; CIMBL83/GEMS41; CML486/GEMS41	23.74	1.77–3.55	Zm00001d001879	ARF-transcription factor 3	auxin response factor 8	auxin response factor, putative, expressed	https://www.maizegdb.org/
	3	CIMBL83/GEMS41	10.72	162.11–170.59	Zm00001d042445 (te1)	protein terminal ear1	terminal EAR1-like 1	AML1, putative, expressed	(Veit et al., 1998)
	7	CML360/GEMS41	14.7	172.49–178.90	Zm00001d022416	early response to dehydration 15-like protein	dehydration-induced protein (ERD15)	early response to dehydration 15, putative, expressed	https://www.maizegdb.org/
	10	CML496/GEMS41; CML304/GEMS41	18.12	146.13–149.26	Zm00001d026595	early nodulin-like protein 1	early nodulin-like protein 5	plastocyanin-like domain containing protein, putative, expressed	https://www.maizegdb.org/

Note. Chr, chromosome; EH, ear height; LA, leaf angle; PH: plant height; PVE: phenotypic variation explained.

DISCUSSION

Multiparental populations are often more powerful to study the genetic basis of complex traits than biparental and natural populations. For NAM-type populations, the ideal set of founders should be selected to maximize genetic diversity (Gage et al., 2020). In this study, we constructed a novel maize HNAU-NAM1 population consisting of 1,625 lines, by crossing an inbred GEMS41 line with 12 representative inbred lines. GEMS41 line was selected as a common and recurrent parent in 12 backcross subpopulations, due to its moderate growth period (55–60 d), high combining ability, and wide environmental adaptability in China. Both the parents and the progenies in HNAU-NAM1 population exhibited high levels of genetic and phenotypic diversity in LA, PH, and EH, which provided a good prerequisite for further genetic dissection of plant architecture in maize.

Considering the special backcross population design used to obtain the BC₁F₄/BC₂F₄ lines, only a limited amount of genetic structure was observed in the HNAU-NAM1 population, which was similar to the levels observed in the ‘Halle Exotic Barley 25′ (HEB-25) NAM population but different from those estimated in the MAGIC (which lacked any significant stratification) (Dell'Acqua et al., 2015) or the random-open-parent association mapping (which showed significant structure) (Xiao et al., 2016) populations. The LD decay distance is generally regarded as an indicator of the mapping resolution, with shorter distances being beneficial for effective and narrower QTL mapping results. The average genome-wide LD decay distance calculated in the HNAU-NAM1 population was nearly 2.59 Mb, which is comparable to 1–4 Mb of the maize MAGIC population (Dell'Acqua et al., 2015) and allows for a relatively high mapping resolution.

Three complementary mapping methods, specifically SLM, JLM, and GWAS, were integrated to decipher the QTLs associated with LA, PH, and EH, as seen in many studies involving multiple-parent populations (Chen et al., 2019; Li et al., 2021; Xiao et al., 2016). During the GWAS analysis, the widely used P-value of 1 × 10⁻⁵as a threshold level to declare significant marker-trait association was adopted here, given the potential linkage disequilibrium among correlated SNPs, which could increase the proportion of false negative findings when using Bonferroni correction (Marees et al., 2018). Identification of 6–14 common QTLs for each of three traits by integration of three mapping approaches showed the strong reliability of these loci. Further, we also observed 4–62 QTLs by using these methods individually. Similar phenomena have been observed in previous maize NAM populations (Chen et al., 2019; Kump et al., 2011). The relative deficiency of overlap between distinct methods may be due to lack of coincident segregation between the causal SNPs and QTLs, different allele frequencies, varying QTL effect sizes, and intrinsic distinct QTL model assumptions for each of the mapping methods (Xiao et al., 2016). Hence, by applying three complementary methods in the HNAU-NAM1 population, we were able to comprehensively dissect the genetic basis of maize plant architectures.

While considering the mapping results for LA, PH, and EH together, we obtained a total of 10 QTL hotspot regions. These regions are likely to include QTLs that exhibit pleiotropic effects for the three aforementioned traits. Alternatively, the genes affecting each trait might be closely linked, as suggested by the significant correlation coefficients observed. In addition, several genes were reported to mediate multiple plant architecture traits (Supplemental Figure S13 and Table S6), such as na2 and acs7, which were known to participate in the regulation of LA and PH (Li et al., 2020b; Tao et al., 2004). Both te1 and pin1 could influence PH and EH (Carraro et al., 2006), and, interestingly, te1 was also found to be located in the hotspot3 region, but whether it can control the development of maize LA needs to be further investigated. In addition, some genes, like brd1 and br2, could simultaneously regulate the development of LA, PH, and EH (Li et al., 2020a; Tian et al., 2019). Hence, it is likely that some of the pleiotropic genes located in these QTL hot-spot regions play a crucial role in regulating maize plant architecture.

In retrospect, the process of maize domestication and improvement involved pyramiding favorable alleles and including some key genes and major effect QTL alleles influencing important agronomic traits (Wang et al., 2018). These genes include tb1, which inhibits the development of axillary organs and promotes the formation of female inflorescences (Doebley et al., 1997), and tga1, which mainly controls the naked kernels on the maize cob (Dorweiler et al., 1993; Wang et al., 2005). Even though a total of 18–88 QTLs associated with LA, PH, and EH have been mapped by each of the three mapping methods (Figure 2; Supplemental Figure S11, S12), but we mainly focused on the 13 large-effect QTLs simultaneously detected by all methods (Table 3). Interestingly, from these 13 QTL genomic regions, we disclosed four genes, namely Zm00001d033981 (asp1), Zm00001d002005 (lg1), Zm00001d042777 (lg2), and Zm00001d042445 (te1), which are known to be involved in the regulation of LA, PH, and EH (Cao et al., 2020; Moreno et al., 1997; Veit et al., 1998; Walsh et al., 1998). These most probably represent fundamental genes underlying plant architecture, and their location within major QTL regions strongly demonstrates the power and efficiency of our strategy to dissect the genetic basis of quantitative traits using HNAU-NAM1 population. Furthermore, we analyzed the eight remaining QTL regions and identified potential candidate genes, for example, four hormone-related genes (Zm00001d021237, Zm00001d001879, Zm00001d032295, and Zm00001d008749) were found to be located within four major QTL mapping regions. Importantly, phytohormones are known to play vital roles across the different stages of plant growth and development, including LA (Kir et al., 2015; Li et al., 2020c), PH (Carraro et al., 2006; Liu et al., 2007), and EH (Li et al., 2017; Liu et al., 2019). However, further fine-scale mapping and functional gene validation is required to confirm whether these candidate genes truly represent the causal agents underlying QTLs for different traits.

Maize germplasms exhibit remarkable genetic diversity, and both SNPs and structural variation play significant roles in phenotypic divergence (Dooner & He, 2008; Jiao et al., 2017). The availability of high-quality reference genomes, combined with accurate and complete genome annotations, provide fundamental information to study genetic and functional differences in plants (Edwards et al., 2013). Similar to the ongoing developments in the US–NAM parents sequencing project (https://nam-genomes.org/), deciphering the parental genome sequences of the HNAU-NAM1 populations in the future will greatly facilitate the cloning of the causal genes underlying the QTLs associated with plant architecture traits. This will broaden our insights into genetic and regulatory mechanisms of maize yield and ultimately allow us to achieve higher productivity in modern maize breeding programs.

DATA AVAILABILITY STATEMENT

The relevant Perl scripts, phenotypic and genotyping data of HNAU-NAM1 population, have been deposited in the Github (https://github.com/dashengzhao/HNAU-NAM1).

ACKNOWLEDGMENTS

Sheng Zhao, Xueying Li, and Junfeng Song contributed equally to this work. The authors would like to express their sincere appreciation to all the collaborators from Henan Agricultural University, who had made some efforts in the phenotyping of some traits, and from Huazhong Agricultural University and Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, who had given some advice in the data analysis. This work was supported by National Natural Science Foundation of China (31970379), and the Elite Young Scientists Program of Chinese Academy of Agricultural Sciences.

AUTHOR CONTRIBUTIONS

Sheng Zhao: Formal analysis; Methodology; Validation; Visualization; Writing-original draft; Writing-review & editing. Xueying Li: Investigation; Methodology; Resources; Junfeng Song: Investigation. Huimin Li: Investigation; Methodology. Xiaodi Zhao: Investigation. Peng Zhang: Methodology. Zhimin Li: Investigation. Zhiqiang Tian: Investigation. Meng Lv: Investigation. Ce Deng: Investigation. Tangshun Ai: Investigation. Gengshen Chen: Methodology. Hui Zhang: Methodology. Jianlin Hu: Methodology; Software. Zhijun Xu: Methodology; Software. Jiafa Chen: Methodology. Junqiang Ding: Conceptualization; Resources; Supervision; Writing-review & editing. Weibin Song: Resources; Supervision. Yuxiao Chang: Conceptualization; Writing-review & editing.

CONFLICT OF INTEREST

The authors declare no competing interests.