Lnglng M, Chunyn Qing, Ursul Frei, You Shen*, Thoms Lübberstedt,*

aDepartment of Agronomy,Iowa State University,Ames 50010,USA

bKey Laboratory of Biology and Genetic Improvement of Maize in Southwest Region,Maize Research Institute,Sichuan Agricultural University,Chengdu 611130,China

A B S T R A C T

1. Introduction

Nitrogen(N)is a major element required for plant growth.As a C4plant, maize (Zea mays L.) accumulates biomass efficiently under abundant N supply with high photosynthetic efficiency [1]. However, increasing N input is not a viable strategy for increasing grain yield, owing to high N prices and environmental pollution [2,3]. In contrast, improving nitrogen use efficiency (NUE) is a strategy for maximizing economic return while minimizing environmental impact[1,4].

In maize and other plant species, the root system plays a major role in NUE [5]. Root system architecture (RSA)traits are also closely associated with grain yield [6,7].Thus, it is desirable to understand the genetic basis of RSA response to diverse N conditions. Maize has five main types of roots: crown, seminal, primary, lateral, and brace roots [8]. Primary and seminal roots are components of the embryonic root system. Brace and crown roots are postembryonic shoot-borne roots, formed above and below the soil surface, respectively [9]. Lateral roots are derived from the pericycles of other roots. Lateral roots are important not only for plant performance but also for water and nutrient uptake, such as N uptake in maize[10,11].

Numerous quantitative trait loci (QTL) have been found to be associated with RSA or NUE in maize. Thirty QTL for RSA traits were detected using a BC4F3(Ye478×Wu312)population[7].In a genome-wide association study(GWAS)in a subset of the Ames panel [12], 263 and 4 significant markers were associated with RSA traits using a general linear model(GLM)and mixed linear model(MLM),respectively[11].In the BGEM(B indicating Iowa State University inbred lines, and GEM-DH indicating Germplasm Enhancement of Maize Double Haploid project from USDA-ARS) panel, Sanchez et al. [5] reported 35 significant SNPs associated with RSA as well as shoot traits using GLM + Q, FarmCPU, and MLM GWAS models. QTL analyses using a F2:3mapping population revealed two to six QTL associated with ear-leaf area, plant height, grain yield,ears per plant, kernel number per ear, and kernel weight under high or low nitrogen(HN,LN)levels[13].Eight QTL were identified for grain yield under LN conditions, two of them also under HN conditions. Five QTL for anthesis-silking interval or ears per plant were consistent across two LN environments[14].Only a few studies have focused on genetic dissection of RSA under diverse N conditions. For example, a total of 17 QTL were detected for RSA traits under HN and LN levels,with one major QTL explaining 43.7%of the phenotypic variation for average axial root length(AARL)[1].Owing to the difficulty in RSA trait measurement, most previous studies have measured only a few RSA traits that were easy to measure, including axial root length, axial root number,lateral root length, total root length, total root surface area,root dry weight,and shoot dry weight[1,7].There is a need to conduct further analyses in relation to a larger array of RSA traits.

Various genes controlling RSA have been cloned in maize: Rtcs (rootless concerning crown and seminal roots),Rtcl (rootless concerning crown and seminal roots like), Rth1(roothairless 1), Rth3 (roothairless 3), Rth5 (roothairless 5),Rth6 (roothairless 6), Bige1 (big embryo 1), and Rum1(rootless with undetectable meristems) [15]. Genes Rtcs and Rtcl encode LOB domain proteins, which are key elements in auxin signal transduction. Rtcs controls shootborne and seminal root initiation; Rtcl controls shoot-borne root elongation [16,17]. Rth1, Rth3, Rth5, and Rth6 affect root hair elongation by encoding SEC3 (secretory3), COBL(COBRA-like), NADPH (nicotinamide adenine dinucleotide phosphate) oxidase, and CSLD5 (cellulose synthase-like D5)proteins, respectively [18-22]. Gene Bige1 encodes a multidrug and toxin extrusion transporter that causes an increased number of crown roots [23]. The function of Rum1 is controlling seminal and lateral root initiation via auxin (IAA) signal transduction [24].

Numerous genes involved in N metabolism have been cloned and characterized in Arabidopsis, rice and other species. For example, in Arabidopsis, genes NLP7 (NIN like protein 7) and AtCIPK8 (CBL-interacting protein kinase 8) are involved in N metabolism [25,26]. Overexpression of OsENOD93-1 improved NUE in rice [27]. Gene OsGLN1;1 promoted a rapid conversion of ammonium to glutamine(GLN) at the surface cell layer under LN conditions [28]. Most of the N metabolism genes cloned in previous studies belong to the nitrogen transporter (NRT), glutamine synthetase (GS),asparagine synthetase (AS), and coordinate C metabolism pathways[29].However,only a few genes in maize have been demonstrated to affect NUE, such as GS1 (glutamine synthetase1),Dof1(DOF zinc finger protein 1),and Ms44(male sterile 44). Overexpressing the genes GS1 and Dof1 increased NUE[30-32]. Recently, gene Ms44 encoding a lipid transfer protein was reported to increase maize yield significantly under LN conditions[33].

In the present study, three GWAS models were applied to analyze RSA traits at two N levels in maize. The objectives were to 1) characterize the phenotypic variation of RSA traits of 14-day old seedlings from the BGEM panel under low and high N levels, 2) identify SNP markers associated with RSA traits under two N levels, and 3)identify candidate genes of seedling RSA traits affecting N response.

2. Materials and methods

2.1. Plant materials

GWAS was performed using the BGEM panel [5]. This population was constructed as follows: a backcross (BC1F1)population was developed using exotic maize landraces as donors and two expired Plant Variety Protection lines,PHZ51 and PHB47, as recurrent parents. PHZ51 and PHB47 represented non-Stiff-Stalk (NSS) and Stiff-Stalk (SS) heterotic groups, respectively. BC1F1plants were crossed with the inducer hybrid RWS 9×RWK-76(with an R-nj color marker)to produce haploid plants [34,35]. Haploid seed was planted in the greenhouse. Seedlings were treated with colchicine to promote genome doubling and transplanted to the field to produce DH lines, of which 226 (Table S1) were used for GWAS.

2.2. Root phenotyping

The 226 DH lines and the two recurrent parents (PHZ51 and PHB47)were grown in a growth chamber using a randomized complete block design. The chamber parameters were as follows: light/darkness, 16/8 h; temperature under light/temperature under darkness, 25/22 °C; photosynthetically active radiation,200 μmol photons m-2s-1;relative humidity,65%.In each trial(replication),eight kernels of each line and the two parents were treated with 6% sodium hypochlorite for 15 min and then rinsed three times with deionized water.The eight sterilized kernels per line were divided into two sets of four and then placed on two sheets of brown germination paper (Anchor Paper, St. Paul, MN, USA) soaked with captan solution (concentration: 2.5 g L-1). The paper sheets were rolled up and placed in vertical orientation in 2-l glass beakers with HN(15 mmol L-1NO3-)or LN(1.5 mmol L-1NO3-) Hoagland solution. The solution compositions were described by Abdel-Ghani et al. [6]. Three independent growth chamber trials were completed on February 14,March 4,and March 25,2018.

Phenotypic values of the traits described in Table S2 were recorded for 14-day old seedlings. Three seedlings from each line with healthy and consistent growth were selected to measure the traits in one replication of each N level. Before measuring, the seedlings were kept in 30% ethanol in a 4 °C cold room to prevent growth.Shoot length(SHL)and primary root length (PRL) were measured using a ruler. For the other root traits, the roots were scanned with a flatbed scanner(EPSON Expression 10,000 XL, Epson America, Inc., CA, USA)and the traits were measured in the images using ARIA(Automatic Root Image Analysis) software [36]. Shoot dry weight(SDW)and root dry weight(RDW)were then measured after oven drying at 80°C for 48 h.

2.3. Phenotypic data analysis

Analysis of variance (ANOVA) of the phenotypic data was based on the model yij=μ+Ri+Gj+Eij,where yijrepresents the observation from the ijth plot, μ is the overall mean, Riis the effect of the ith replication,Gjis the effect of the jth genotype,and Eijis the experimental error. The ANOVA table, expected mean squares, and least-square means were obtained using the function PROC GLM of SAS (Statistical Analysis System,version 9.3, SAS Institute, Cary, NC, USA). Type 3 sums of squares were used to account for missing data. The statisticsand H2(heritability)were calculated by entry means as follows[11]:

Here,MSG,MSE,and rep denote mean square for genotype,mean square error, and the number of replications (3) in the experiment. Eleven traits studied previously and having heritabilities ＞0.3 (Table S3) were used for GWAS. These were shoot length (SHL), total root length (TRL), lateral root length (LRL), network area (NWA), primary root length (PRL),median (MED), total number of roots (TNR), root dry weight(RDW), total plant biomass (TPB), surface area (SUA), and shoot dry weight(SDW).

SPSS (Statistical Product and Service Solutions, version 21.0, IBM, Armonk, NY) was used to compare means of phenotypes for the two N levels for each trait based on ttests, calculate Pearson correlation coefficients, and evaluate phenotypic trait correlations.

2.4. Genotypic data analysis

BGEM lines were genotyped using genotyping by sequencing(GBS)resulting in 955,690 SNPs(Genomic Diversity Laboratory,Cornell University). After filtering and corresponding, 62,077 SNPs were yielded by Sanchez et al.[5].In this study,a total of 61,634 markers remained for GWAS after removal of SNPs with ＞20% missing data, ＞20%heterozygosity, or minor allele frequency ＜0.05.

2.5. Population structure, linkage disequilibrium, and association study

Based on 61,634 SNPs, principal component analysis (PCA)implemented in GAPIT (Genome Association and Prediction Integrated Tool) [37] was used to estimate the number of subgroups in the BGEM population. TASSEL 4.0 [38] was used to calculate linkage disequilibrium (LD) among all 61,634 markers, and a threshold of r2= 0.2 was considered to determine the genetic distance within which LD decay occurred in the panel. HaploView software (http://www.broad.mit.edu/mpg/haploview) was used to calculate the LD decay between the SNPs in target regions.

To balance false negative and false positive SNPs in GWAS,three models were used: GLM (General Linear Model) + PCA,FarmCPU (Fixed and random model Circulating Probability Unification), and MLM (Mixed Linear Model). PCA was calculated by GAPIT and used as a covariate.The models GLM+PCA,FarmCPU,and MLM were conducted using TASSEL 4.0[38],the FarmCPU R package [39], and the GAPIT R package [37],respectively. Models FarmCPU and GAPIT were both run in the R studio(version 3.4.1)environment.

The simpleM program [40,41] was used to account for multiple testing in R studio and to calculate the effective number of independent tests (Meff_G). The significance threshold value was calculated using a Bonferroni correction:P-value=α/Meff_G(α=0.05).The calculation process of Meff_Gwas as follows: first, a correlation matrix for all 61,634 SNPs was constructed and corresponding eigenvalues for each SNP locus were calculated.A composite LD(CLD)correlation was then calculated directly from SNP genotypes [37], and once the SNP matrix was created, the Meff_Gwas calculated. Here,Meff_Gwas 23,760, so that the P-value was 0.05/23,760(2.10×10-6).

The identification of candidate genes with highest priority was based on the following principles:1)a gene harboring the significant SNP(s) was identified as a candidate gene; 2) if a significant SNP was located in an intergenic region, the LD region of the SNP was scanned for all the gene models.The r2between the significant SNP and each of the SNPs in the above gene models was then calculated. If the r2≥0.8, the corresponding gene model was considered a candidate gene[42].

3. Results

3.1. Phenotypes of the BGEM-DH lines

For most traits, high variation was observed in the panel.Under HN conditions, the traits TRL and LRL of the BGEM-DH panel showed the highest standard deviations (SDs) of 94.9 cm and 93.3 cm, respectively (Table S3). SDs of TRL and LRL were 80.4 cm and 78.8 cm under LN conditions (Table S3). The increased percentage of phenotypic values ranged from-50.0% to 18.0% compared to that under the two N levels(Table S3). Heritability (H2) estimates of the 24 seedling traits ranged from 0.39 to 0.85 and 0.20 to 0.83 under HN and LN conditions, respectively (Table S3). For the 11 traits used for GWAS, the phenotypic values of these traits all followed normal distributions (Figs. S1, S2). Among the 24 traits, 19 showed significant phenotypic differences (P ＜ 0.05) between the two N levels (Table S3), suggesting that the two N treatments for the BGEM-DH panel were effective.

At each N level,most of the 11 traits displayed significant correlations (P ＜0.05) with one another, and the Pearson correlation coefficients(r)ranged from 0.10 to 1.00 at HN and 0.04 to 1.00 at LN levels (Figs.S1,S2).For example,traits TRL and NWA showed significant correlations under both N conditions, and the relationship between them was extremely close (r = 0.99) (Figs. S1, S2). For most of the 11 seedling traits, a significant correlation (P ＜0.05) between phenotypic values at HN and LN levels was observed (Table S4). The correlation coefficients ranged from 0.02 to 0.79(Table S4).

3.2. Population structure

The 226 BGEM-DH lines were divided into two major subpopulations using PCA implemented in GAPIT (Fig. S3). This division is consistent with the construction of the BGEM population.One subgroup with PHB47 background(stiff stalk)contained 112 (49.6%) lines, and the other subgroup with PHZ51 background (non-stiff stalk) contained 80 (35.4%) lines(Fig.S3).Thirty-four(15.0%)lines,BGEM-0007-S,BGEM-0017-S,BGEM-0052-S, BGEM-0053-S, BGEM-0078-S, BGEM-0092-S,BGEM-0094-S, BGEM-0111-S, BGEM-0112-S, BGEM-0113-S,BGEM-0114-S, BGEM-0115-S, BGEM-0116-S, BGEM-0117-S,BGEM-0118-S, BGEM-0165-S, BGEM-0166-S, BGEM-0171-S,BGEM-0175-S, BGEM-0189-S, BGEM-0220-S, BGEM-0266-S,BGEM-0269-S, BGEM-0005-N, BGEM-0085-N, BGEM-0107-N,BGEM-0121-N, BGEM-0129-N, BGEM-0132-N, BGEM-0215-N,BGEM-0227-N, BGEM-0232-N, BGEM-0248-N, and BGEM-0260-N were mis-grouped (Fig. S3). Here, N and S represent the PHZ51 (NSS) and PHB47 (SS) backgrounds. This misclassification was consistent with findings of a previous study[43]and was most likely due to the contribution of exotic germplasm of BGEM-DH lines.

3.3. Linkage disequilibrium decay

All 61,634 SNPs, evenly distributed over all 10 chromosomes(Fig. S4-A) were used to estimate LD decay using TASSEL 4.0.According to the LD threshold of r2= 0.2, LD decay in the BGEM-DH panel among all 10 chromosomes reached or exceeded 2 Mb(Fig.S4-B).The development of the population from two recurrent parents with large contributions was speculated to be the main reason for the slow LD decay.

3.4. Genome-wide association study

Using the MLM, FarmCPU, and GLM models, 33 and 51 SNPs were found to be significantly associated with 11 seedling traits under HN and LN conditions, respectively (Tables 1, 2).Among them, LD regions of 9 HN-SNPs and 22 LN-SNPs overlapped with the intervals of QTL that were identified in previous studies as controlling root system development(RSD) and N response (Table 3). The QTL that reside in the respective intervals bnlg1025/umc2029 (223,622,305/233,892,356 bp), bnlg278/phi048 (194,870,616/204,606,322 bp),bnlg278/phi048 (194,870,616/204,606,322 bp), umc2043/umc1061 (134,791,672/139,123,323 bp), and umc2043/umc1061(134,791,672/139,123,323 bp) and control RSD and N response were located in the LD regions of five significant SNPs,S1_223834059, S5_204522218, S5_205927835, S10_134650981,and S10_138694384 detected under the LN treatment (Table 3). According to the significant SNPs, 43 and 68 candidate genes were obtained for the HN and LN treatments, respectively (B73 genome, RefGen_v2, Tables 4 and S5). The annotations of the candidate genes were shown in Tables 4 and S5.

Under HN conditions, 2, 11, and 20 significant SNP-trait associations(P threshold=2.10×10-6)were detected using the MLM,FarmCPU,and GLM+PCA models,respectively(Table 1).One SNP,S9_2483543 located on chromosome 9,was found to be associated with traits SDW (P = 8.01 × 10-7, SNP effect =-0.02), TPB (P = 1.79 × 10-7, SNP effect = -0.02), and SDW (P =1.57 × 10-6, SNP effect = -1.30) by MLM (Fig. 1-A, B), FarmCPU(Fig.1-C,D)and GLM+PCA(Fig.1-E,F),respectively(Table 1).This SNP is within gene model GRMZM5G833563(2,479,762-2,484,723 bp), which encodes a DNA polymerase lambda (POLL) (Table 4). SNP S3_10516162 that was identified as controlling trait TPB, was detected by both FarmCPU and GLM + PCA (FarmCPU: P = 2.39 × 10-7, SNP effect = -0.01;GLM + PCA: P = 9.45 × 10-7, SNP effect = -1.15) (Table 1). This SNP is within gene model GRMZM5G802971, located between 10,515,400 and 10,518,765 bp on chromosome 3 and annotated as encoding a hypothetical protein(Table 4).Under FarmCPU,SNP S10_146939958 was significantly associated with LRL(P =1.81 × 10-9, SNP effect = 34.0), NWA (P = 3.86 × 10-9, SNP effect=0.20),and TRL(P=2.01×10-6,SNP effect=24.7)(Table 1). This SNP is within a QTL region (umc2122/umc1344)located on chromosome 10, which was previously reported to control root system development (Table 3). SNP S7_29655312 was associated with PRL (P = 6.89 × 10-7, SNP effect = -10.04) and SUA (P = 1.35 × 10-6, SNP effect = -6.53)using GLM + PCA (Table 1), placed within unknown gene model GRMZM5G821968 (29,652,133-29,656,462 bp) (Table 4).Another SNP(S7_96016906)was associated with both PRL(P =1.35 × 10-7, SNP effect = -8.68) and SUA (P = 5.16 × 10-8, SNP effect=-5.46)based on the GLM+PCA model(Table 1),located within unknown gene GRMZM2G035451(96,014,968-96,017,273 bp) (Table 4).

Table 1-Significant SNPs detected by GWAS using three models under HN treatment.

Under LN conditions, 24 and 27 SNPs were detected (P threshold = 2.10 × 10-6) as significantly associated with seedling traits using the FarmCPU and GLM + PCA models,whereas no significant SNP was identified using the MLM model (Table 2). Three significant SNPs, S3_15473839,S5_168222967, and S2_5836079 were identified as controlling multiple traits using FarmCPU or the GLM+PCA model(Table 2). SNP S3_15473839, located in the intergenic region of chromosome 3, was associated with three traits, namely LRL(P = 9.75 × 10-10, SNP effect = 35.2), NWA (P = 9.64 × 10-8, SNP effect=0.18),and TRL(P=5.81×10-8,SNP effect=32.1)(Table 2). Three seedling traits, LRL (P = 9.53 × 10-8, SNP effect =-29.0), TNR (P = 7.37 × 10-9, SNP effect = -1.66), and TRL (P =4.96×10-8,SNP effect=-31.3)were all significantly associated with SNP S5_168222967 (Table 2). The candidate gene for this SNP is GRMZM2G059851 (168,222,641-168,226,086 bp), which encodes a HSF-transcription factor 6 (hsftf6)(Table 4).Marker S2_5836079 was associated with five seedling root traits,NWA(P = 3.66 × 10-7, SNP effect = -11.69), TNR (P = 3.69 × 10-7, SNP effect = -9.72), LRL (P = 7.73 × 10-7, SNP effect = -17.86), TRL(P = 8.81 × 10-7, SNP effect = -18.17), and MED (P = 1.94 × 10-6,SNP effect=-5.83)using the GLM+PCA model(Table 2).This marker is located within cQTL2_1Y (TIDP3756/gpm762a),which was identified as controlling yield-related traits under LN by a previous study(Table 3).

Among all the significant SNPs detected under the two N conditions, only one SNP, S10_66115383, located in the intergenic region of chromosome 10, was found to be associated with PRL under both N treatments (HN: P =8.92× 10-8,SNP effect= 20.52;LN: P = 1.62 × 10-9, SNP effect=29.14)(Table 1,Table 2).

4. Discussion

4.1. Phenotypic data analysis

Root system traits are associated with yield traits [6,7].However, RSA traits have not attracted extensive attention during selective breeding in maize, owing to the difficulty in their measurement. To overcome this issue, a paper roll method was developed [46]. It was used in this study toculture 14-day old seedlings.Phenotypes of root system traits were collected using a root scanner and a high-throughput image-analysis tool, ARIA. Currently, there are several image analysis software packages: DART [47], WinRhizo (Pro, 2004),EzRhizo [48],and ARIA [36].Compared with DART, WinRhizo,and EzRhizo, ARIA has several advantages: 1) the algorithms of ARIA are fast and efficient;2)It provides a large number of traits by a high-throughput image-analysis package; 3) this tool is freely available and openly accessed[36].

Under HN, line BGEM-162-S showed the highest values for TRL, LRL, MED, TNR, and NWA, whereas BGEM-109-N displayed the lowest values for TRL, LRL, TNR, and NWA.Similar phenomenon was observed under LN. Interestingly,all above traits showed positive correlations(P ＜0.05)with oneanother (Figs. S1, S2). Moreover, most of the 11 traits showed significant correlations (P ＜ 0.05) with each other (Figs. S1, S2). All these findings suggest that there are positive interactions between seedling traits and that multiple traits may be controlled by the same genes. In this BGEM-DH population, most of the traits showed 2- to 3-fold differences between the minimum and maximum values under each N condition (Table S3).However, in a study of maize root architecture by Pace et al. [11], the differences between the minimum and maximum values of 22 seedling traits for Ames panel were of almost 10- to 30-fold. The main reason is probably that the population used in this study was a BC1F1-DH panel. The background of this panel contains about 75% background from recurrent parents PHZ51 or PHB47, whereas the Ames panel harbors a wider range of genetic variation [12].

Table 3-Significant SNPs overlapping with QTL for root system development and nitrogen response identified in previous studies.

Heritability estimates of seedling traits ranged from 0.20 to 0.85 across both N levels (Table S3), higher than those reported by Pace et al. [11] and Sanchez et al. [5].Specifically, 13 of the 24 seedling traits showed heritabilities exceeding 0.50 under both N conditions (Table S3),whereas no trait exceeded heritability values of 0.50 in the two earlier studies. Environmental variation influences trait heritability estimates. In the present study, the seedlings were cultured in Hoagland solution with LN or HN.However, Hoagland solution was replaced with deionized water in the two previous studies. The different culture conditions probably account for the difference of the H2estimates between the previous studies and our study.Ribaut et al. [14] reported higher heritability estimates of grain yield-related traits under LN and HN conditions, in agreement with our findings.

As in previous studies, a heritability threshold of 0.3 of seedling traits was set for conducting GWAS analyses.Thus, 11 traits, SHL, TRL, LRL, NWA, PRL, MED, TNR,RDW, TPB, SUA, and SDW were selected for GWAS in this study.

4.2. Linkage disequilibrium decay

LD decay differs among maize populations. For example,the LD decay was below 1000 bp for maize landraces [49],approximately 2000 bp for diverse maize inbred lines [50],and even larger than 100,000 bp for commercial elite inbred lines [51]. In the present study, the LD decay exceeded 2 Mb for the BGEM panel, roughly 20-fold that of commercial elite inbred lines. The BGEM panel is a BC1F1-derived population of DH lines that was constructed by backcrossing donor accessions (exotic maize landraces from GEM project) with two recurrent parents PHZ51 and PHB47 (See “Plant materials” section). The low LD decay in the BGEM population was probably due to the high proportion of recurrent-parent genome. In our previous study [52], when only the donor alleles of the BGEM population for a given region, Comt,were considered, the LD decay was more rapid (below 100 bp), supporting our speculation. To identify the reason for low LD decay, four significant SNPs (S1_9992325,S4_197073985, S9_2483543, S9_154381179) were randomlyselected. First, SNPs located in 2 Mb up- and downstream of these four SNPs were respectively collected. Then, the LD decay between the SNPs was calculated in the first step for four markers, respectively. The results revealed 45, 33, 27,and 74 blocks within the LD regions of S1_9992325,S4_197073985, S9_2483543, and S9_154381179, respectively.Notably, several large (＞50 kb) blocks were identified (Figs.2, S5-S7) as block 24 (72 kb; including S9_4219970,S9_4289806, S9_4291242, S9_4292438, and S9_4292681) (Fig.2-A), block 4 (162 kb; including S9_1427784, S9_1443236,S9_1447762, S9_1487251, S9_1587620, and S9_1590534) (Fig.2-B), block 13 (478 kb; including S9_2634177, S9_2711824,S9_2825523, S9_2972263, S9_3010175, S9_3067052,S9_3082788, and S9_3112757) (Fig. 2-C), and block 2 (76 kb;including S9_999394, S9_1002574, S9_1028270, and S9_1075796) (Fig. 2-D) for LD region of S9_2483543. These findings suggest that the presence of some large blocks from the recombinant panel (BGEM-DH population) is a key factor causing low LD decay.

Table 4-Annotations of candidate genes harboring the significant SNPs based on B73(RefGen_v2)genome.

4.3. Genome-wide association study and candidate genes involved in RSA under two N conditions

In this study, three models, GLM + PCA, FarmCPU, and MLM were used for GWAS to balance false positives and false negatives identified by GWAS. As described in a previous report [53], the GLM model tends to result in increased numbers of false positives, whereas the MLM model may result in too many false negatives. FarmCPU is based on iterative algorithms and addresses the confounding problem of testing markers and using them as covariates [39]. Among the three GWAS models, MLM was most stringent and GLM was least stringent. Across 11 traits under HN conditions, 2,11,and 20 significant SNP-trait associations were detected by MLM,FarmCPU,and GLM+PCA,respectively.Under LN,0,24,and 27 SNP-trait associations were detected. These findings are consistent with those of previous studies [5,53]. Only two significant SNPs were co-detected by multiple models. The SNP S3_10516162 controlling TPB was identified by both FarmCPU and GLM + PCA. S9_2483543 was found by all three models. The candidate genes for these two SNPs were GRMZM5G802971 (S3_10516162) and GRMZM5G833563(S9_2483543), which encode a hypothetical protein and a DNA polymerase lambda(POLL),respectively.

Fig.1-Significant SNPs detected by GWAS using three models under HN conditions.HN represents high nitrogen.(A),(C),and(E) represent Manhattan plots of S9_2483543 identified by MLM,FarmCPU,and GLM+PCA,respectively.(B),(D),and (F)represent QQ plots of S9_2483543 identified by MLM,FarmCPU,and GLM+PCA,respectively.SDW,shoot dry weight;TPB,total plant biomass.

The LD regions of 9(27.27%)SNPs under HN and 22(43.14%)SNPs under LN overlapped with the intervals of QTL for RSD and N response that were reported previously (Table 3). In particular, the known genes that affect root system development or NUE were located in the LD regions of four SNPs(S1_9992325, S9_154381179, S4_197073985, and S9_151726472)(Table 3). The distance between SNP S1_9992325 and Rtcs is 0.83 Mb and that between S9_151726472 and Rtcl is 0.75 Mb(Fig.3-A,B;Table 3).In addition,S9_154381179 is 1.90 Mb from Rtcl, and S4_197073985 is 2.0 Mb from Ms44 (Fig. 3-C, D; Table 3). SNPs S1_9992325 and S9_154381179 were located within GRMZM2G062841 and GRMZM2G092776 (annotation: putative DUF26 domain receptor-like protein kinase), respectively(Table 4). Only one SNP (S10_66115383) was detected under both N conditions, suggesting that this SNP controlled PRL independent of N level. The candidate gene for SNP S10_66115383 was GRMZM2G043749, which encodes a brain protein 44 (Table S5). The other SNPs that affected seedling trait development only at specific N levels was speculated to regulate N metabolism.

According to the significant SNPs detected in this study, a total of 43 and 68 genes were identified under HN and LN conditions, respectively (Tables 4, S5). Several candidate genes were associated with seedling development, seed development, root system development, or N metabolism.Among them, GRMZM2G139811 (S2_37861383), associated with PRL is annotated as a C2 and GRAM domaincontaining protein (At5g50170; Table 4). Remarkably, several GRAM domain-containing proteins, Abr1 (ABA-responsive protein 1), Abr2 (ABA-responsive protein 2), and Abr3(ABA-responsive protein 3) play important roles in ABAmediated inhibition of seed germination in Arabidopsis [54].The SDW-associated gene, GRMZM2G314898, (S3_8300032)encodes a HXXXD-type acyl-transferase family protein(Table 4). Its ortholog in Arabidopsis, AT2G39980, was induced to be expressed in roots by ACC (1-aminocyclopropane-1-carboxylic acid) treatment, and shown to control root cell elongation [55]. Gene model GRMZM2G054050(S3_24354497), which was associated with trait SHL encodes a multicopper oxidase Lpr2 (low phosphate root 2) (Table 4),whose homologous gene, Lpr1 was reported [56] to control the length of primary roots under Pi (Potassium)-deficient conditions. Lpr1 and Lpr2 act as key elements in Pi sensing at root tips [57]. GRMZM2G173682 (S10_148533344), associated with TNR, encodes Atg4b (autophagy 4b) (Table 4). A previous study [58] using a T-DNA insertion double mutant of Atg4a and Atg4b (atg4a4b-1) demonstrated that Atg4a and Atg4b contribute to root system development under N stress conditions. A TNR-associated gene, GRMZM2G470914(S3_5849337), is annotated as a receptor-like serine/threonine-protein kinase Ale2 (Table 4). Tanaka et al. [59]reported that Ale2 controls shoot development by specifying epidermis cells in Arabidopsis. Another PRL-associated gene, GRMZM2G462325 (S7_19083078), is annotated as Ost1(oligosaccharide transferase1) (Table 4), reported [60] to regulate seed development by participating in ABA signaling. GRMZM2G416184 (S7_116386029), associated with PRL,encodes a BTB/POZ domain-containing protein Npy1 (Table 4), playing an essential role in root gravitropic responses[61]. GRMZM2G064302 (S9_20906793), which was associated with TNR, is annotated as Eno1 (enolase1). Two T-DNA insertion Eno1 mutants were previously reported to reduce numbers of root hairs and distort trichomes in Arabidopsis[62].

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.11.004.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Fig.3-Manhattan plots of four significant SNPs linked to known genes.(A)Rtcs located within the LD region of S1_9992325;(B)Rtcl located within the LD region of S9_151726472;(C)Rtcl located within the LD region of S9_154381179;(D)Ms44 located within the LD region of S4_197073985.Rtcs,and Rtcl are associated with root development,and Ms44 is associated with nitrogen response.TNR,total number of roots;SHL,shoot length;TPB, total plant biomass.

Acknowledgments

The authors thank the China Scholarship Council (CSC) for Langlang Ma's funding. The authors also thank USDA's National Institute of Food and Agriculture (IOW04314,IOW01018), Hatch Multistate Project NC-007, as well as the R.F.Baker Center for Plant Breeding at Iowa State University,for supporting this work.