Genome-wide association studies of leaf angle in maize

Compact plant-type with small leaf angle has increased canopy light interception, which is conducive to the photosynthesis of the population and higher population yield at high density planting in maize. In this study, a panel of 285 diverse maize inbred lines genotyped with 56,000 SNPs was used to investigate the genetic basis of leaf angle across 3 consecutive years using a genome-wide association study (GWAS). The leaf angle showed broad phenotypic variation and high heritability across different years. Population structure analysis subdivided the panel into four subgroups that correspond to the four major empirical germplasm origins in China, i.e., Tangsipingtou, Reid, Lancaster and P. When tested with the optimal GWAS model, we found that the Q + K model was the best in reducing false positive. In total, 96 SNPs accounting for 5.54–10.44% of phenotypic variation were significantly (P < 0.0001) associated with leaf angle across three years. According to the linkage disequilibrium decay distance, 96 SNPs were binned into 43 QTLs for leaf angle. Seven major QTLs with R2 > 8% stably detected in at least 2 years, and BLUP values were clustered in four genomic regions (bins 2.01, 2.07, 5.06, and 10.04). Seven important candidate genes, Zm00001d001961, Zm00001d006348, Zm00001d006463, Zm00001d017618, Zm00001d024919, Zm00001d025018, and Zm00001d025033 were predicted for the seven stable major QTLs, respectively. The markers identified in this study can be used for molecular breeding for leaf angle, and the candidate genes would contribute to further understanding of the genetic basis of leaf angle.

1 3 to achieve high yield in maize. Leaf angle, which is a major factor in determining the ideal plant-type, directly affects the distribution of light and CO 2 in the canopy and the utilization of light energy of the population and then affects the growth and development of the plant and ultimately the yield (Pendleton et al. 1968). At the high planting density, the hybrids with upright leaves always have higher grain yield than those with flat leaves (Lambert and Johnson 1978). It is one of the effective ways to increase yield by breeding ideal plant-type hybrids with suitable leaf angle in maize.
The ligular region has been proved to be the key part to regulate leaf angle (Kong et al. 2017). Biosynthesis, transport and signal transduction of plant hormones (Zhang et al. 2014b;Tian et al. 2019), and cell division and elongation in the ligular region played a crucial role in regulating the leaf angle in maize. Genes for leaf angle have been cloned mainly using mutants or homology-based cloning method in maize, for instance, rld1, dwil1, lbl1, lg1, lg2, lg3, and lg4 (https:// www. maize gdb. org/). ZmTAC1 in bin 2.07-2.08, which explained 7.30% of phenotypic variation for leaf angle and was the ortholog gene of OsTAC1 in rice, was cloned using comparative genomics (Ku et al. 2012). Only four QTLs for leaf angle have been isolated by mapbased cloning (Zhang et al. 2014b;Tian et al. 2019;Ren et al. 2020). Map-based cloning was firstly used to successfully identified ZmCLA4 on Chromosome 4 responsible for leaf angle by regulating polar auxin transport and the shape and number of cells in the ligular region (Zhang et al. 2014b). The brd1 (brassinosteroid C-6 oxidase1) gene, which underlined the QTL on Chromosome 1 (UPA1), was mapbased cloned and controlled leaf angle by altering endogenous brassinosteroid content (Tian et al. 2019). UPA2, which was a QTL for leaf angle on Chromosome 2, might repress the expression of a B3-domain transcription factor (ZmRAVL1) through a two-base sequence polymorphism. ZmRAVL1 regulated the expression of brd1 and then affecting the proliferation of cells in ligular region (Tian et al. 2019). The ZmILI1 gene underlining another QTL for leaf angle on Chromosome 2 was an ortholog of OsILI1 encoding a basic helix-loop-helix leucine zipper family transcription factor in rice (Ren et al. 2020).
Complementing the traditional linkage analysis, association mapping analysis can simultaneously evaluate the effects of multiple alleles. Moreover, association mapping analysis using diversity populations with rapid LD decline can detect the functional sequence that causes phenotypic variation in maize (Yan et al. 2011). For example, using the nesting association mapping population of maize with 4,892 RILs, Tian et al. (2011) detected 30 QTLs for leaf angle by the method of GWAS. Among them, seven SNPs in lg1, lg2, and lg4 contributed to more upright leaves. A total of 42 SNPs, which located in bins 1.02, 1.03, 1. 06, 1.11, 2.05, 3.04, 5.03, 5.04, 8.03, 8.04, and 10.07, were detected significantly associated with leaf angle by the methods of GWAS in four environments using 289 maize inbred lines (Sun et al. 2018). Among them, QTLs for leaf angle located in bins 1.11, 5.03 and 8.03 could be detected in two or more environments (Sun et al. 2018). Although some QTLs/genes associated with leaf angle have been mapped and cloned in maize, there is still a lack of deep understanding of genetic bases of leaf angle.
In the present study, three GWAS models were used to analyze leaf angle across 3 years in maize. The objectives were to (1) characterize the phenotypic variation of leaf angle from the 285 inbred lines across 3 consecutive years; (2) identify SNP markers associated with leaf angle across 3 years, and (3) identify candidate genes for leaf angle.

Plant materials
The present study involved an diverse panel of 285 inbred lines suitable for planting in Tianjin, which were selected from three sets of inbred lines, i.e., the first from the core collection established previously by Li et al. (2004), including 242 diverse accessions which were historically used in maize breeding, the second from research institutions or companies including 648 elite inbred lines, most of which are used in current maize breeding, and the last from 400 lines from USA, CIMMYT, Africa and other countries . The final maize panel included 254 lines from the three main maize production regions of China, i.e., the North China, the Huanghuaihai Valley, and the Southwest China, 24 lines from the USA, three from Europe, two from Africa, one from Australia, and one from Mexico. The pedigrees or sources of the lines are listed in Table S1.

Phenotypic evaluation
All the 285 accessions were planted in Tianjin Academy of Agricultural Sciences Innovation Base, Wuqing,Tianjin (39.43° N,116.95° E) in 2017, 2018, and 2019. Each line was grown in single rows that was 3.0 m in length and 0.6 m apart, with 24.7 cm between each plants in rows at a planting density of 67,500 plants ha −1 , using a randomized complete block design with two replications each year. Five representative plants from the middle of each plot were chosen to measure the leaf angle as follows at the 10th day after anthesis. Leaf angle was measured as the angle between the vertical and the midrib of the first three leaves above ear positon, and the means were used as the phenotypic data to analysis.

Phenotypic data analysis
The descriptive statistical analysis and Pearson correlation analysis for leaf angle across 3 years were performed using SPSS 22.0 software. Analysis of variance (ANOVA) of leaf angle across 3 years was performed using the lmer function of the lme4 package of R based on the following model: y ijk = μ + g i + e j + r k + ge ij + ε ijk , where y ij is the observation from the ijth plot, μ is the grand mean over all years, g i is the genotypic effect of the ith genotype, e j is the year effect of the jth year, r k is the replication effect of the kth replication, ge ij is the genotype × year interaction effect, and ε ij is the residual error, with genotype, year, genotype by year, and replications as random effects. The broad-sense heritability (H 2 ) of leaf angle were calculated using the following formula: H 2 (%) = σ 2 g / (σ 2 g + σ 2 ge /n + σ 2 e /nr) × 100%, where σ 2 g is the genotypic variance, σ 2 ge is the interaction variance of genotype with year, σ 2 e is the error variance, n is the number of years, and r is the number of replications (Knapp et al. 1985). To minimize environmental effects, best linear unbiased prediction (BLUP) values for leaf angle of each line across all years were calculated with the same ANOVA model. The BLUP values were also used as the phenotypic data of the final GWAS. The phenotypic distribution of leaf angle across 3 years was performed using R software.

Genotyping and SNP filtering
The 285 lines were genotyped by MaizeSNP50 Bead-Chip, which included 56,110 SNPs (http:// suppo rt. illum ina. com/ array/ array_ kits/ maize snp50_ dna_ analy sis_ kit. html). Totally 39,827 SNPs passed quality control via the software PLINK 1.90 (Purcell et al. 2007) according to the following SNP screening criteria: (1) the minor allele frequency (MAF) exceeded 0.05 and (2) the missing rate is less than 0.2. The source sequences of 39,827 SNPs were identified through BlastN search against the reference genome sequence of B73 (B73_RefGen_v4.). SNPs with ambiguous physical position or multiple blasthits were excluded from the genotyped dataset. After that, 39,074 high-quality SNPs were used in further analysis.

Population structure analysis
The model-based clustering algorithm of ADMIX-TURE v1.3 was used (Alexander et al. 2009) to investigate subpopulation structure of the 285 maize accessions. A subset of 6,527 genome-wide SNP markers filtered by PLINK v.1.90 software using the parameters -indep-pairwise 50 10 0.5, MAF ≥ 0.05 and missing data < 5% were used. A preliminary analysis was performed in multiple runs by inputting successive values of K from 1 to 10. A fivefold CV procedure was performed for each K value. The most likely K value was determined using ADMIXTURE's CV values. The stacked bar-charts were plot using barplot command in the R software. Inbred lines with membership probabilities of more than 0.5 were assigned to corresponding clusters (Liu et al. 2012a, b).
To obtain an in-depth picture of the genetic relationships among 285 inbred lines, the genetic distance between pair-lines was calculated by using 3,9074 SNPs. Modified Euclidean was defined as follows: D = 1 -identity by state (IBS) similarity (Legesse et al. 2007). IBS means the probability that alleles derived at random from two individuals at identical loci are the same. For any two lines, the probability of IBS was averaged over all non-missing loci. The cladogram was constructed based on the distance matrix described above based on the neighbor-joining (NJ) algorithm. Clusters were observed from the result of the phylogenetic tree. All calculations were carried out in TASSEL software 5.2.69 (Bradbury et al. 2007). Principal component analysis with 39,074 SNPs of the 285 individuals was carried out using the TASSEL software 5.2.69.
Association study and candidate gene identification Principal component analysis (PCA) was implemented in TASSEL 5.2.69, and the top five principal components were used to estimate the population structure matrix to control population structure. Kinship matrices were calculated using the dominance normalized IBS method in TASSEL 5.2.69 to estimate genetic relatedness among individuals. To balance false-negative and false-positive SNPs in GWAS, three models were performed in TASSEL 5.2.69: (1) Q model: GLM (General Linear Model) + PCA; (2) K model: MLM (Mixed Linear Model) + kinship; and (3) Q + K model: MLM + PCA + kinship. The quantile-quantile plots (QQ plot) obtained from each analysis were plotted using the R software, and the optimal statistical model was identified for the GWAS analysis. Because several SNPs should be in LD, the GEC software tool (Li et al. 2012) was used to calculate the effective number of independent SNPs (Ne). The recommended threshold (1/Ne) calculated by the software was used as the basis for the significant correlation between leaf angle and SNPs.
Based on the average LD decay across the whole genome, significant SNPs within the LD decay distance were binned into one QTL, and the most significant SNP were selected as the lead SNP. All the potential candidate genes within LD decay distance upstream and downstream of the lead SNP of all QTLs were identified from Gramene database (http:// ensem bl. grame ne. org/) based on the maize B73 Ref-Gen_V4 genome. The biological functions of the candidate genes were predicted from Gramene database, maize GDB database (http:// www. maize gdb. org), and Gene Ontology database. In addition, some candidate genes were also annotated by the homologous genes in rice (RAP database) or Arabidopsis thaliana (Tair database). Combined with the above gene annotation and related literatures, the most possible candidate gene for leaf angle was identified in each QTL interval.

Phenotypic variation and heritability of leaf angle
The basic statistical analysis for leaf angle of the 285 inbred lines showed that extensive and continuous variation were observed for leaf angle in different years, and the variation range of leaf angle in 2017 and 2018 was larger than that in 2019 and BLUP for the 3 years ( Fig. 1, Table 1). It indicated that the leaf angle was a typical quantitative trait and was controlled by polygenes. The variances of genotype (σ 2 g ), year (σ 2 y ) and genotype × year (σ 2 gy ) were significant at P < 0.01, and the broad-sense heritability of leaf angle was 93.94% (Table 2). Significant correlation (P < 0.01) for leaf angle among 3 years and BLUP values (r 2 = 0.79-0.97) was also observed. It indicated that leaf angle was mainly controlled by genetic factors.

3
Model-based population structure and relative kinship A total of 39,074 high-quality SNPs were mapped in the reference genome B73_RenGen_v4, and the number of SNPs on individual chromosome ranged from 2754 on Chromosome 10 to 6259 on Chromosome 1 (Fig. S1). When counting the number of SNPs per 1 Mb across the whole genome, the distribution density of SNPs at both ends of each chromosome was relatively high, while the distribution density of SNPs near the centromere was relatively low. The highest density region contained 63 SNPs/Mb, and the lowest contained 1 SNPs/Mb. The average SNP density was 20 SNPs/Mb across the whole genome.
The population of 285 accessions based on the model-based algorithm was analyzed using the even distribution of 6527 SNPs. As the ADMIXTURE software overestimated the number of subgroups for inbred lines, it was difficult to choose the "correct" K from the cross-validation error of data ( Fig. 2A). Thus, the results were compared with the known pedigree of the inbred lines for each run of different Ks. The results showed that when K = 4, the modelbased subgroups were largely consistent with known pedigrees of the inbred lines. The four groups corresponded to the four major germplasm origins in China, i.e., Tangsipingtou (TSPT), Reid, Lancaster, and P ( Fig. 2B; Table S1).
The TSPT group was the largest group with 91 inbred lines, which were mainly derived from Huangzaosi, one of the founder parents in maize breeding of China. In addition, some lines derived from the landrace of Ludahonggu and the inbred line of "Zi330" were also designed into this group. The Reid group included 56 inbred lines, which were mainly selected from Pioneer hybrids "3382" and "U8." The Lancaster group included 41 inbred lines, which were mainly derived from "Mo17" and "Va35." The P group included 32 inbred lines, which was mainly derived from the Pioneer hybrid of "P78599." Additionally, 65 inbred lines that had < 0.5 membership in each of the four subgroups and had a mixture of two or more subgroups were assigned to a mixed subgroup (Table S1).   Estimation of relative kinships showed that 94.16% of paired relative kinship ranged from 0.05 to 0.30, 0.33% of paired relative kinship equaled to 0, and 1.16% of paired relative kinship ranged from 0 to 0.05 (Fig. 3). This analysis indicated that a weak or various relative kinship existed in this collection of germplasm.
Neighbor-jointing tree and principal component analysis A neighbor-joining tree was constructed based on the modified Euclidean distance and is shown in Fig. 4. The 285 accessions were firstly clustered into three major groups according to their origins: TSPT, the mixed group of P and Reid, and Lancaster, which contained 81, 125, and 79 accessions, respectively. The TSPT group could be further clustered into two major subgroups according to their origins: the Chang7-2 cluster with 45 lines and the Zi330 cluster with 36 lines. The P + Reid group was then clustered into the Ye8112 cluster including 61 lines, the 599 cluster including 51 lines, and the Va35 cluster including 13 lines. The Lancaster group could be clustered into two subpopulations, i.e., the C416 cluster with 39 lines and the Fengke1 cluster with 40 lines (Fig. 4). These results showed moderately common with some difference compared to that evaluated by using model-based method.
The PCA results showed obvious population differences among the groups of TSPT, Reid, Lancaster and P, and a good agreement with both modelbased population structure and clustering analyses (Fig. 5). The 95% confidence intervals of Lancaster and TSPT partly overlapped, which may be resulted by the lager introgression existing between Lancaster and TSPT accessions and lower power of PCA in population structure analysis by using only two PCs. The representative lines Chang7-2, Ye8112, and 599 for TSPT, Reid, and P groups exhibited obvious differentiation.

Linkage disequilibrium analysis
The LD decay indicated that the mean of r 2 in the panel dropped rapidly at approximately 10 Kb, and the mean of r 2 at 27.5 Kb was 0.2 (Fig. S2). The average LD decay distance was 383 kb in the panel. Combined with the average LD decay distance and the physical location of peak SNP, the candidate genes at each QTL can be further inferred.

Genome-wide association analysis
The effective marker number of 39,074 high-quality SNPs was analyzed by using the GEC software. The results showed that the effective marker number Ne was 24,166, and the suggested P value was 4.14E-5, that is, -log P was 4.38. Therefore, this value was used as the threshold for significant marker-trait association in the whole genome association analysis of leaf angle.
The leaf angles in the three years and their BLUP values were analyzed by GWAS using three models (Q, K, and Q + K). QQ diagram showed that Q model and K model had high type I errors (false positive), while Q + K model was the best in reducing false positive (Fig. 6). Therefore, a mixed linear model based on Q + K model was used for leaf angle genome-wide association analysis.
Genome-wide association analysis was carried out for the leaf angles across the 3 consecutive years and their BLUP values of 285 diversity inbred lines combined with 39,074 SNPs markers covering the genome. A total of 96 SNPs significantly associated with leaf angle were detected based on the Q + K model (P < 4.14E-5), which were distributed on all chromosomes except Chr. 8 (Fig. 7, Table S2). Based on the average LD decay distance across the whole genome, 96 SNPs significantly associated with leaf angle were distributed in 43 QTLs, ranging from one QTL (Chr. 3) to nine QTLs (Chr. 2).
In the present study, the QTLs contained SNP with R 2 > 8% in at least 1 year were considered as major QTLs, and the major QTLs detected in at least 2 years and also for BLUP values were called "stable major QTLs." Thus, seven stable major QTLs located in four chromosome segments of bin 2.01, bin 2.07, bin 5.06, and bin 10.04 were detected (Fig. 7, Table 3).
The first stable major QTL was QTL11 located in bin 2.01. This locus could be detected in all 3 years and joint analysis (BLUP), explaining 7.18-9.29% of the phenotypic variation for leaf angle.
Two stable major QTLs, QTL17 and QTL18, were both detected in bin 2.07. The two QTLs could be detected in all 3 years and for their BLUP values, explaining 8.13-9.24% and 6.50-8.10% of the phenotypic variation for leaf angle, respectively. The fourth stable major QTL (QTL29) was located in bin 5.06 and could be detected in 2017, 2018, and for BLUP values, explaining 8.69-9.01% of the phenotypic variation for leaf angle.
Three stable major QTLs (QTL41, QTL42, and QTL43), including 24 SNPs significantly associated with leaf angle, were detected in bin 10.04. QTL41 was detected in 2018, 2019 and for BLUP values, explaining 6.79-8.69% of the phenotypic variation for leaf angle. Both QTL42 and QTL43 could be detected in 3 years and their BLUPs, explaining 6.79-10.44% and 6.40-10.02% of the phenotypic variation for leaf angle, respectively.

P+Reid TSPT Lancaster
Ch an g7 -2 C h a n g 7 -2 Z i3 30 Y e 8 1 1 2 Y e 8 1 1 2 5 9 9 5 9 9 C 4 1 6 C 4 1 6 F e n g k e 1 F e n g k e 1 V a 3 5 Prediction of candidate genes for leaf angle Candidate genes within the range of LD decay distance between the upstream and downstream of a given lead SNP were identified for each QTL. The most promising candidate gene was predicted for each of 43 QTLs for leaf angle (Table 3, Table S2). The candidate genes are mainly related to the imbalanced

Population structure of the 285 diverse inbred lines in China
The main heterotic groups according to the pedigree information in China varied with different stages, that is Golden Queen, Huobai, TSPT, Ludahonggu, and Lancaster in 1970s (Wu 1983); TSPT, Ludahonggu, Lancaster, and Reid in 1980s (Zeng 1990); Lancaster (Mo17 subgroup and Zi330 subgroup), Reid, TSPT, Ludahonggu, and other germplasms in the 1990s (Wang et al. 1998;Teng et al. 2004); and Reid, P, Ludahonggu, TSPT, and Lancaster in the early 2000s (Teng et al. 2004;Sun et al. 2014 In the present study, 285 inbred lines were assigned to four major groups by the model-based method and PCA, i.e., TSPT, Reid, Lancaster, and P group, which were consistent with known pedigrees of the inbred lines and previous studies Wu et al. 2014a). The TSPT group played a dominant role in Chinese maize breeding and occupied an important position in Chinese maize heterotic patterns. For example, Zhengdan958 derived from the pattern of Reid × TSPT was the most important Chinese maize hybrid with the largest annual planting area in consecutive 16 years since 2004. Ludan981   (Li and Wang 2010), most of which and their descendants were clustered into the TSPT group in the present study. These results were consisted with Xie et al. (2007Xie et al. ( , 2008, who integrated the Ludahonggu and TSPT group into the D group. The reciprocal genetic improvement of TSPT and other local germplasm, e.g., Ludahonggu and Golden queen were also founded based on the pedigree of some excellent lines, i.e., 5237 derived from Huangzaosi × Dan340 (TSPT × Ludahonggu), Jing2416 derived from 5237 × Jing24, Xun92_8 derive from Chang7-2 × 5237, and Shuang741 derived from [(3Tuan-11 × Huafeng100) × Huangzaosi] × Aijin525. In additionally, Zi330 and its descendants were also always used to improve lines from TSPT by Chinese breeders. For example, Ji853, a famous parent in TSPT, was derived from Huangzaosi × Zi330. Furthermore, several reports classified Zi330 and Ludahonggu into one group (Teng et al. 2004;Wang et al. 2008). Therefore, lines from Ludahonggu, Zi330, and local landraces were important donors to improve the TSPT germplasm. Complex genetic basis of Lancaster was confirmed by the present study. Some lines of Lancaster according to pedigree were assigned to TSPT and P + Reid by neighbor-jointing tree analysis, such as lines from the Zi330 subgroup and the Va35 subgroup. The PCA results also showed that the Lancaster germplasm generally close to other germplasm. One reason was that the Lancaster Surecrop germplasm was less than 50% in the lines of Lancaster in China (Zheng et al. 2002).
The Reid group still occupied an important position in Chinese maize heterotic patterns. The Reid germplasms, such as M14, W24, W59E, and B73, were firstly introduced to China in the 1950s. Later, using a batch of hybrids introduced from the USA, such as XL80, 3147, U8, and 3382, Chinese breeders developed several excellent inbred lines, such as Ye107, U8112, Ye478, Tie7922, and Shen5003 (Li and Wang 2010). Recently, the female parent of Denghai605, the fourth Chinese maize hybrids with planting area more than 0.85 million hectares in 2019, was also derived from Reid.  The P inbreds, such as 178, P138, 87,001, Qi319, Dan598, Dan599, Ye107, and 18-599, have been developed from the Pioneer hybrid of "78599'' or similar hybrids since the late 1980s (Xie et al. , 2008Li and Wang 2010). Due to the existence of tropical and subtropical germplasm in the P group, the resistance to biotic and abiotic stress is often strong. Therefore, the introgression of the P germplasm to Reid, SS, or NSS germplasm attracted more attention of breeders. Recent years, the success of several hybrids with extensive planting in China, e.g., Denghai605, Yufeng303, and Zhongkeyu505, was partly attributed to the P germplasm.

Comparative analysis of the genetic basis for leaf angle
The heritability of leaf angle varied, depending on genetic populations and environments in which populations are evaluated. It was found that the heritability of leaf angle was 90.5% in different years (Wassom 2013), 54.56-87.00% at different locations (Mickelson et al. 2002;Chen et al. 2015;Ding et al. 2015), and 68.00-81.15% in different years combined with locations (Ku et al. 2010;Li et al. 2015). In addition, the heritability of leaf angle was 84.00% ) under different planting densities (52,500 plants/ha, 67,500 plants/ha and 82,500 plants/ha). Except that the heritability of leaf angle in Xiema, Beibei and Hechuan in Chongqing was 54.56-58.75% , more than 68.00% was observed for the heritability of leaf angle in different environments in other reports. In this study involving 285 diverse inbred lines of maize, it was found that the heritability of leaf angle was also very high (93.94%), indicating that leaf angle was mainly controlled by genetic factors and less influenced by environments. Zhao et al. (2018b) integrated the results of 27 QTL mapping studies for leaf angle, leaf orientation, leaf length, leaf width, leaf area and leaf length/ width, and 16 meta-QTLs (mQTLs) for leaf angle were found to be distributed on 8 chromosomes except Chr. 9 and Chr. 10, suggesting the complexity of the genetic basis of leaf angle in maize. In this study, a total of 96 SNPs significantly associated with leaf angle (P < 4.14E-5) were distributed in 27 bins on the nine chromosomes except Chr. 8. We further identified 43 QTLs for leaf angles, of which 7 stable major QTLs (explaining more than 8% of the phenotypic variation for leaf angle and could be stably detected in at least 2 years and also for BLUP values) were distributed in bin 2.01, bin 2.07, bin 5.06, and bin 10.04. QTL11 located in bin 2.01 were also identified by Lu et al. (2007), Ding et al. (2015), Li et al. (2015), and Wang et al. (2017) using different linkage mapping populations, explaining 2.86%, 4.06%, 4.30%, and 8.85% of phenotypic variation, respectively. Ku et al. (2012) and Zhao et al. (2018b) integrated the previous QTLs for leaf angle and detected a mQTL in bin 2.01 by meta-analysis. This locus in bin 2.01 stably expressed in different environments and different genetic backgrounds might be an important target for further study on leaf angle.
Two stable major QTLs (QTL17 and QTL18) were detected in bin 2.07, and a total of 8 SNPs significantly associated with leaf angle within these two QTLs could be stably inherited in different years. Liu et al. (2012a) detected a major QTL explaining 10.8% of phenotypic variation for leaf angle in bin 2.02-2.07. Zhang and Huang (2021) also found a major QTL in bin 2.07, which accounted for 14.2% of phenotypic variation for leaf angle. Due to the limited number of parents, the interaction of QTL × genetic background, population size and marker density in linkages analysis, and population structure and rare alleles in association analysis, few studies identified QTLs for leaf angle in this bin.
QTL29 located in bin 5.06 was a hot spot locus for leaf angle founded by previous studies (Mickelson et al. 2002;Lu et al. 2007;Zhang et al. 2014a;Zhao et al. 2018b) and could explain 4.17-9.40% of phenotypic variation. In addition, Ku et al. (2012) integrated the previous QTL mapping results of leaf angle and found an mQTL in bin 5.06 by meta-analysis. The mQTL affected both leaf angle and leaf orientation, and four candidate genes, yabby9, SE, LIC, and yabby15, were predicted. Therefore, bin 5.06 might play an important role in leaf angle.
In previous studies, QTLs for leaf angle detected were mainly located in bin 10.07 (Yu et al. 2006;Zhang et al. 2014a;Sun et al. 2018;Wang et al. 2017), while bin 10.01-10.02 (Zhao et al. 2018b), bin 10.06 , and bin 10.08-10.09 (Yu et al. 2006) were also reported. However, the QTL in bin 10.04 was the most important locus found in the present study, which carried three stable major QTLs (QTL41, QTL42, and QTL43), and contained 24 SNPs significantly associated with leaf angle. Therefore, the present study found for the first time that bin 10.04 plays an important role in the regulation of leaf angle in maize, and there may be at least three major genes that control leaf angle within this region. It was necessary to further map-based clone the genes underline by constructing large-scale secondary segregation populations.

Analysis of candidate genes for leaf angle
In recent years, a few genes regulating leaf angle or tiller angle were cloned in rice, maize, etc. It has been clear that leaf ligular region is the key tissue to regulate leaf angle (Zhu et al. 2015;Kong et al 2017;Hu et al. 2019). Factors such as the imbalanced division and expansion of adaxial and abaxial cells, the strength of mechanical tissue, and the gravity response in the ligular region will affect the leaf angle (Zhu et al. 2015;Hu et al. 2019). The molecular regulation mechanism of leaf angle mainly includes (1) the regulatory pathway of plant hormone synthesis and signal transduction (brassinosteroid, auxin, and gibberellin have been deeply studied on the regulation of leaf angle, while other hormones such as jasmonic acid, strigolactones, and ethylene have also been found to be involved in the regulation of ligular region development) and (2) other regulatory pathways (the gravity response of the ligular region (such as LPA1, which affects the division of adaxial cells in ligular region, and LAZY1, which regulates the polar auxin transport) or the strength of mechanical tissue in ligular region (such as ILA1, which controls the formation of vascular bundles and the content of cellulose and xylan in the cell wall; OsCesA7 and OsCD1 involved in cellulose synthesis; DL, the YABBY gene family, which controls the formation of leaf blade-midribs; OsIG1 involved in the division and differentiation of vesicular cell and the formation of ligular region)) (Zhu et al. 2015;Luo et al. 2016;Hu et al. 2019).
The progress of identifying the regulatory genes for leaf angle in maize is relatively slow. Through the research on mutants, a few genes in regulating leaf angle was cloned, such as LG1, encoding squamosa promoter binds to protein (SBP); LG2, encoding a basic leucine zipper protein (bZIP); and lg3 and lg4, encoding homeo domain protein. The mutants of lg1 and lg2 (Moreno et al. 1997;Walsh et al. 1998) and the dominant mutants of lg3-O and lg4-O showed that the abnormal development of ligular region resulted in a small leaf angle (Bauer et al. 2004). Through homologous cloning, ZmTAC1, a homologous gene of TAC1 that regulated leaf angle in rice, had been proved to be able to positively regulate the leaf angle (Ku et al. 2012). In recent years, ZmCLA4 underlined qLA4-1 in bin 4.03 (Zhang et al. 2014b), brd1 underlined UPA1 in bin 1.08 and two-base sequence polymorphism underlined UPA2 in bin 2.03 (Tian et al. 2019), and ZmILI1 underlined qLA2 in bin 2.02 (Ren et al. 2020) for leaf angle have successively map-based cloned. Seven stable major QTLs for leaf angle were identified in the present study, which will be useful for the improvement of leaf angle by molecular breeding and provide a basis for the cloning of the genes.
The candidate gene of QTL11 might be Zm00001d001961, which encoded SAUR-like auxinresponsive protein. The majority of SAUR genes can induce growth in leaves, stems, and floral organs via cell elongation (Stortenbeker and Bemer 2019). In addition, the overexpression of SAUR could also influence the auxin levels, polar auxin transport, and/ or expression of auxin pathway genes (Stortenbeker and Bemer 2019). Auxin can also played role on the leaf angle via the regulation of the growth of the abaxial and abaxial cells of ligular region and auxin signal transduction (Song et al. 2009;. So Zm00001d001961 might regulate leaf angle by the division and expansion of adaxial and abaxial cells. Zm00001d006348, the predicted candidate gene for QTL17, encoded a growth-regulating factor (GRF) with two conserved domains, i.e., WRC and QLQ (Choi et al. 2004). Through the interaction between QLQ domain in the N-terminal of GRFs and the SNH domain in GIF (GRF-interaction factor) protein, a functional complex which regulated the expression of genes downstream was formed . In Arabidopsis thaliana, GRFs could promote or inhibit the growth of organs by regulating the cell size (Kim et al 2003) or the cell division (Kim and Lee 2006;Vercruyssen et al. 2015). Over-expression of BrGRF8 from Chinese cabbage  and BnGRF2 from rapeseed (Liu et al 2012b) in transgenic Arabidopsis plants increased the sizes of the leaves and other organs by regulation of cell proliferation. It was also found that genes for GRFs could regulate several traits, such as leaf size, plant height, and 1 3 flowering traits in maize (Wu et al 2014b;Zhang et al. 2008). Their transcriptional expression was affected by gibberellin (GA) ) and miR396 (Jones-Rhoades et al. 2006). Because the protein encoded by the candidate gene of Zm00001d006348 was the homologous sequence of that encoded by AT2G45480 (growth-regulating factor 9) in Arabidopsis thaliana, the function of Zm00001d006348 can be predicted by AT2G45480. AT2G45480 negatively regulated leaf growth by activating the expression of a bZIP transcription factor gene, i.e., OBP3-responsive gene (ORG3), and restricted cell proliferation in the early stage of leaf development (Omidbakhshfard et al. 2018). Therefore, Zm00001d006348 might negatively regulate the proliferation of cells in the ligular region and further affect the leaf angle.
Zm00001d006463, the candidate gene for QTL18, encoded C2C2-Dof-transcription factor with a highly conserved DNA-binding domain of Cys2/Cys2 (C2/ C2) zinc finger (Yanagisawa, 2004). The DOF proteins were found involved in many biological processes, such as tissue differentiation, seed development, and regulation of metabolism (Noguero et al. 2013). It is worth noting that the DOF transcription factor AtDOF5.1 could regulate leaf adaxial-abaxial polarity by controlling the expression of Revoluta (REV), a Class III homeodomain-leucine zipper (HD-ZIPIII) transcriptional regulatory protein . Therefore, Zm00001d006463 might also regulate the leaf adaxial-abaxial polarity in maize.
Zm00001d017618, a possible candidate gene for QTL29, encoded a B3-domain transcription factor that might regulate the expression of genes in the BR synthesis pathway and the signal transduction pathway (Tian et al. 2019;Je et al. 2010). Recently, Tian et al. (2019) found that the B3-domain transcription factor gene Zm00001d002562 (ZmRAVL1) on Chrom. 2 activated the expression of brd1 downstream, increased endogenous brassinosteroid content, promoted proliferation of cells in ligular region, and enlarged leaf angle. Therefore, Zm00001d017618 was predicted to regulate the leaf angle via brassinosteroid content in the ligular region.
Zm00001d024919 encoding an adenylate kinase was the possible gene for QTL41. Two near isogenic lines with different leaf angle were analyzed by the method of comparative proteomic analysis, and adenylate kinase complicated in glycometabolism might be associated with leaf angle formation and the physical and mechanical properties of midribs .
The candidate gene for QTL42 might be Zm00001d025018 encoding alpha expansin. Expansins responsible for cell-wall relaxation are required for cell elongation (Kong et al. 2017). Furthermore, increased expansin expression at early stages of organ development induced cell proliferation and consequently increased growth (Wyrzykowska et al. 2002). Both the size of the cells determined by the cell elongation and the number of cells determined by the cell division influence the size of the auricle and ultimately leaf angle (Kong et al. 2017).
Zm00001d025033 encoding a TCP transcription factor was the possible candidate gene of QTL43. TCP transcription factors with TCP conserved domain formed a bHLH structure had been found to regulate many aspects of plant development, such as branching in maize (Doebley et al. 1995;Bai et al. 2012). TB1 (Teosinte Branched 1) can inhibit the growth and development of lateral branches, and tb1 mutant gene leaded to increased number of lateral branches in maize (Doebley et al. 1995). BAD1, which is another TCP transcription factor gene and could promote cell proliferation in a lateral organ, i.e., the ligular region, may influence inflorescence architecture by regulating the angle of lateral branch emergence (Bai et al. 2012). Zm00001d025036 encoding indole-3-acetate beta-glucosyltransferase also located in this QTL. The glycosylation of auxin might influence the phenotype of plant via the regulation of auxin levels (Jackson et al. 2002).
Code availability Not applicable.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.
Competing interests The authors declare no competing interests.