Genome-wide association study identified novel genetic loci controlling internode lengths and plant height in common wheat under different nitrogen treatments

Nitrogen is an important nutrient for crop growth and development. Plant height-related traits can be affected by nitrogen supplementation. In this study, we performed a genome-wide association study (GWAS) on plant height, spike length, length of different internodes, and lodging resistance strength at the grain-filling stage based on wheat local varieties subjected to low nitrogen and normal (CK) treatments. GWAS analysis showed that a total of 86 quantitative trait locus (QTLs) were detected, including 13 QTLs for plant height, 10 QTLs for spike length, 19 QTLs for the length of the first internode from the top of the plant, 6 QTLs for the second internode length, 11 QTLs for the third internode length, 13 QTLs for the fourth internode length, and 14 QTLs for the fifth internode length. Compared to the CK treatment, the plant height, spike length, and fourth and fifth internode lengths were significantly affected by the low nitrogen treatment. A total of 18 QTLs responding to low nitrogen level were detected, including three QTLs for the fourth internode length detected on 3A, 6A, and 6D chromosomes, eleven QTLs for the fifth internode length on 1A, 1B, 1D, 2A, 2B, 3A, 3B, 4A, 5B and 7B chromosomes, one QTL for spike length on 3A chromosome, and one QTL for plant height on 5B chromosome. These QTLs will enhance our understanding of the genetic basis of plant height responses to nitrogen deficiency and will benefit genetic reactions to nitrogen fertilization.


Introduction
Nitrogen is an important nutrient for crop growth and development. Large quantities of nitrogen fertilizer are commonly used to obtain high yields for many crops (Cormier et al. 2016;Peng et al. 2010). In addition, modern breeding programs have mostly aimed at obtaining high yields under a high nitrogen supply (Liu et al. 2021). The high input of nitrogen fertilizers has resulted in low nitrogen use efficiency and has caused a series of environmental problems, such as soil acidification and water eutrophication (Shi Abstract Nitrogen is an important nutrient for crop growth and development. Plant height-related traits can be affected by nitrogen supplementation. In this study, we performed a genome-wide association study (GWAS) on plant height, spike length, length of different internodes, and lodging resistance strength at the grain-filling stage based on wheat local varieties subjected to low nitrogen and normal (CK) treatments. GWAS analysis showed that a total of 86 quantitative trait locus (QTLs) were detected, including 13 QTLs for plant height, 10 QTLs for spike length, 19 QTLs for the length of the first internode from the top of the plant, 6 QTLs for the second internode length, 11 QTLs for the third internode length, 13 QTLs for the fourth internode length, and 14 QTLs for the fifth et al. 2010). To solve these kinds of problems and to improve nitrogen use efficiency, the breeding of new cultivars with improved productivity in low nitrogen environments provides an effective approach to increase nitrogen use efficiency (Marino et al. 2011;Xu et al. 2014).
Plant height is an important agronomic trait in wheat, as it is a key parameter affecting the lodging and grain yield of wheat. A reduction in crop height can be associated with reduced grain yield (Law et al. 1978), and the aim of breeding programs is therefore the identification of variants with reduced height without adversely affecting yield (Wurschum et al. 2015). In general, along with a lower incidence rate of lodging, an appropriate grain number per spike, and an improved harvest index, an appropriate plant height could increase grain yield and quality (Hedden 2003).
Classical genetic studies have indicated that the plant height of bread wheat is a complex trait and that genes located on almost all 21 chromosomes can control plant height (Law et al. 1978;Snape et al. 1977). Most dwarfing genes have now been well characterized. Rht-B1 and Rht-D1, which have been successfully used in wheat breeding programs worldwide, are located on the short arms of chromosomes 4B and 4D, respectively (Borner et al. 2002;Huang et al. 2003). Rht13 on chromosome 7BS has been verified to produce a marked reduction in height (Ellis et al. 2005). In addition, the photoperiod-insensitive alleles of the major photoperiod regulator Ppd-1 have also been reported to exert pleiotropic effects on plant height in wheat (Borner et al. 1993).
Previous studies showed that different nitrogen levels had significant effects on plant height in wheat (Hussain et al. 2006;Yu et al. 2020). With increasing nitrogen levels, the plant height will increase accordingly. Plant height in wheat includes the spike length plus all internode lengths above the ground. The nitrogen response of the components of plant height in wheat is probably different. Understanding which components of plant height have significant responses to different nitrogen levels would be beneficial for wheat breeders to efficiently select related traits.
Quantitative trait loci (QTL) mapping based on linkage analysis is often used to analyze the molecular basis of complex traits. The objectives of previous studies were to map QTLs for agronomic and yield traits, uptake, and utilization efficiency of nitrogen in a recombinant inbred line (RIL) population under different supplemental nitrogen environments Xu et al. 2014;Cui et al. 2016;Zhang et al. 2019;Sun et al. 2012). However, QTLs have rarely been reported from local varieties in environments with different levels of nitrogen supplementation.
With the rapid development of genotyping technologies and computational methods, genome-wide association studies (GWAS) are now becoming a powerful tool for detecting underlying complex natural variation traits in crops. Recently, important QTLs for agronomic traits have been found in crops based on GWAS methods (Lou et al. 2020;Pang et al. 2020;Sansaloni et al. 2020;Zhou et al. 2015;Hu et al. 2016). Furthermore, GWAS has more advantages than RIL population mapping, as it provides higher resolution and allows more variation to be detected without the need to construct populations (Huang and Han 2014). The wheat 90K SNP genotyping assay comprised approximately 90,000 gene-associated SNPs covering the whole genome ). Compared to SSR markers, the 90 K SNP assay is generally more abundant, stable, efficient, and costeffective (Sun et al. 2017;Gao et al. 2016). In this study, we performed a GWAS using a set of 120 winter wheat cultivars and a 90 K SNP assay to identify QTLs associated with plant height-related traits under normal nitrogen environments and low nitrogen environments, with a focus on QTLs that respond to nitrogen stress.

Plant Materials and Experimental design
A total of 120 bread wheat cultivars collected from the Yellow and Huai Valley Regions of China, including Henan, Hebei, Shaanxi, and Shandong provinces, were selected as a collection of local varieties (Table S1). These cultivars were planted and harvested in 2017-2018 and 2018-2019. In each cropping season, all the surveyed cultivars were planted in experimental fields of Shandong Agricultural University, Tai'an city (116°36′ E, 36°57′ N). Between the two cropping season, the statistic of temperature and rainfall was similar between 2018 and 2019, according to official data (http:// sd. cma. gov. cn). Two Page 3 of 15 146 Vol.: (0123456789) treatments were applied: low nitrogen (LN) and a normal fertilized control (CK). Hereafter, "2017LN", "2017CK", "2018LN", and "2018CK" represent the treatment trials. Each block had four rows with a length of 1.5 m spaced 25 cm apart.
To analyze the soil nitrogen content, samples were selected at a depth of 0-40, 40-60 cm using a fivepoint sampling method in each block during the seedling stage. The nitrogen content of the soil samples was measured using the alkali-hydrolyzable nitrogen method (Fan et al. 2005). The soil nitrogen contents were the average of the samples. In the depth of 0-40 cm, the soil nitrate nitrogen content in the low nitrogen environments was 90 mg/kg, while that in the CK environments was 135 mg/kg. Significant differences (p < 0.01) existed between the LN and CK blocks. In the depth of 40-60 cm, there was no different soil nitrate nitrogen content (72.5 mg/kg) between LN and CK blocks. In the LN and CK blocks, urea were applied about 300 kg/hm 2 at elongation stage of wheat, and no other fertilizer was applied at the crop growth stages. Crop management was performed according to the local cultivation practices.

Phenotyping
For each cultivar surveyed, five representative tillers from different individual plants were investigated. Plant height and its component traits ( Fig. 1), including plant height (PH), spike length (SL), length of the first internode from the top (FIRILT), length of the second internode from the top (SECILT), length of the third internode from the top (THILT), length of the fourth internode from the top (FORILT), length of the fifth internode from the top (FIFILT), length of the sixth internode from the top (SIXILT), and the peduncle length (PL), were investigated (Fig. 1). The lodging resistance strength (LRS) was measured in the field during the grain filling stage 15 days after flowering. The strength was recorded when the stem was pushed down 20 cm above the ground using a YYD-1 plant culm-strength meter (Hangzhou Wanshen Detection Technology Co., Ltd., HangZhou, China).
In the analysis of variance (ANOVA), the formula trait = σ g 2 + σ e 2 + σ ge 2 + σ er 2 was used to calculate σ g 2 , σ e 2 , σ ge 2 and σ er 2 . The formula h 2 = σ g 2 / (σ g 2 + σ e 2 + σ ge 2 /n + σ er 2 /nr) was used to calculate broad-sense heritability (h 2 ), where σ g 2 , σ e 2 , σ ge 2 , and σ er 2 represent the variance components genotype, environment, genotype × environment and error, respectively, n and r are the numbers of environments and replications, respectively (Lou et al. 2020;Yang et al. 2020). The calculation of pairwise correlation coefficients for these traits was conducted using the "Performance Analytics" package in R software, and the boxplot of traits was visualized using the "ggplot2" package under multiple comparisons methods and was utilized to compare the same trait in the different environments using the "LSD.test" (Fisher's least significant difference test) method. The formula T delt_aver = 100 * (T CK -T LN )/T LN was used to calculate the trait reaction response to LN; T CK represents the trait in the CK environment, and T LN represents the trait in the LN environment.
DNA extraction and genotyping DNA was extracted from multiple young leaf tissues from different plant of each cultivar using the SDS method (Cui et al. 2011). DNA quality was checked by electrophoresis on 1% agarose gels, and DNA concentration was determined with a Nanophotometer (NanoDrop Technologies Co., Ltd., Wilmington, USA). Genotyping was performed by Beijing EMTD Technology & Investment Co., Ltd. (Beijing, China; http:// www. emtd. com. cn) using the 90 K wheat genotyping assay ) following the manufacturer's protocol.
The genotypic clusters for each SNP were determined using Genome Studio version 2011.1 software, filtering out SNPs with a call rate > 0.35. After the filtering out of SNPs with a minor allele frequency of < 5% and SNPs with > 10% missing data, 120 cultivars with 20,689 polymorphism SNPs were retained for further analysis. This processing was implemented in PLINK software (Purcell et al. 2007). The physical maps of all SNPs were obtained from the website (http:// wheat genom ics. plant path. ksu. edu/ wpdb/ assay? from= singl emess age) and are aligned to the IWGSC reference version 1.0 genome (International Wheat Genome Sequencing et al. 2018). Among the 20,689 SNPs, 8,632 SNPs with known physical positions were remained.

Linkage disequilibrium and population structure analysis
Linkage disequilibrium (LD) among markers was calculated for the A, B, and D genomes in PLINK software. The window size for linkage disequilibrium calculation was set based on the number of SNPs located in each genome. Pairwise linkage disequilibrium was measured using squared allele frequency correlations and was assessed by calculating the squared allele frequency correlation (r 2 ) for pairs of SNP loci (Pang et al. 2020). Principal component analysis (PCA) was conducted using GAPIT with the default parameters to calculate principal components (PCs) (Lipka et al. 2012).
Genome-wide association study GWAS was implemented by GAPIT packages in R software using the mixed linear model (MLM, PCA + K), which took the population structure and relative kinship into account ). The kinship matrix (K) was automatically calculated using the VanRaden method (VanRaden 2008). The first four principal components of the SNP data were included in the GWAS model. To combine the GWAS results from all environments, a uniform suggestive genome-wide significance threshold (P value = 1.0e−3) was used. Quantile-quantile (Q-Q) plots and Manhattan plots were visualized with the R package "CMplot" (https:// github. com/ YinLi Lin/ CMplot). The R 2 value was used to evaluate the magnitude of the location's effects.

LD analysis
A total of 120 cultivars and 8632 SNPs with known physical positions were retained for further analysis. The number of SNPs within 1 MB across each chromosome is shown in Figure S1. Moreover, the number of SNPs across each chromosome is shown in Table S1. These SNPs were distributed on all 21 wheat chromosomes, with 3260, 4310, and 1062 SNPs in the A, B, and D subgenomes, respectively.
LD was estimated by calculating the r 2 among all possible pairs of markers for each of the 21 chromosomes. The obtained r 2 values were then plotted against the physical distance for each of the three genomes separately ( Figure S2) and across the whole genome. When the cutoff threshold for r 2 was defined as 0.1, the LD decay distance for the entire genome was greater than 15 MB.

Phenotypic trait analysis among subpopulations under different conditions
We analyzed each phenotypic trait among subpopulations under different conditions. The performance of all traits in the LN and CK environments is displayed in Fig. 2. From Fig. 2, we found that the same traits of subgroups were different. By comparing the same traits under the CK and LN environments within a subgroup (Table S2), we found that the PH and PH components were both greater in the CK treatment than in the LN treatment, except for the SECILT and THILT of the C subgroup. In the ANOVA results, a significant P value of σe 2 existed among the different environments (Table S3), which illustrated that the environment affected the traits. According to the comparison of all traits in the CK and LN environments, the PH, SL, FORILT, FIFILT, and LRS showed significant differences (P < 0.001). The results indicated that the decrease in PH in the LN treatment was mainly caused by decreases in the length of FORILT, FIFILT, and SL in the LN treatment. Using the ANOVA results and formula of broad sense heritability, we calculated h 2 for each trait.
The average data based on each trait were used to calculate Pearson's correlation coefficients. PH was positively correlated with SECITL, THILT, and FORILT ( Figure S3); SL was negatively correlated with FORILT and FIFILT ( Figure S3). PL showed a strong positive correlation with FIRILT ( Figure S3), and the LRS showed no correlation with other traits ( Figure S3), which indicated that LRS had no significant relationship with the length of internodes ( Figure  S3).

GWAS of plant height-related traits under low nitrogen and normal environments
To eliminate the influence of the kinship and structure of the population, we performed GWAS for all traits with the MLM method under CK and LN environments, and significant SNPs (P < 1e−3), repetitive SNPs in multiple environments (all, E > 3; CK, E = 2; LN, E = 2) and SNPs located within the LD decay distance (15 M) were defined as repetitive QTLs. Manhattan plots and QQ plots for all traits are shown in Fig. 2 and Figure S4. All detected SNPs and QTLs in the CK and LN environments are displayed in Table S4. GWAS analysis showed that a total of 86 QTLs were detected in the LN and CK environments.

GWAS of plant height-related traits in response to nitrogen stress
According to the comparison of the same traits in the CK and LN environments, different nitrogen contents in the soil caused the difference in FORILT, FIVILT, SL, PH, and LRS. T delt_aver was measured as the trait that responded to nitrogen stress. Based on GWAS and T delt_aver , QTLs that responded to nitrogen stress were detected ( Fig. 3; Table 1). A total of 18 QTLs that responded to nitrogen stress were detected.

QTLs for plant height
Classical genetic studies indicated that plant height was a complex trait controlled by both Mendelian genes and quantitative genes. Previous genetic studies have identified numerous major QTLs for PH among 21 wheat chromosomes (Yu et al. 2020). Most dwarfing genes have been characterized, and closely linked molecular markers have been developed for marker-assisted selection in wheat breeding programs. Among the dwarfing genes, Rht-B1 and Rht-D1, which have been successfully utilized in wheat breeding programs worldwide, are located on the short arms of chromosomes 4B and 4D, respectively (Wurschum et al. 2015). In additional to Rht-B1 and Rht-D1 genes, other genes associated with PH had been reported in prervious studies such as Rht 8, Rht 13, Rht 22, Rht 24 et al. (Zhang et al. 2017;Wang et al. 2014;Li et al. 2013;Herter et al. 2018). What's more, other QTLs were also identified. Griffiths et al. identify 16 Meta-QTL using four DH population, which distributed in all chromosome except 3D, 4A, 5D, 7A ,7B and 7D chromosomes (Griffiths et al. 2010). Li et al. found a QTL located in 6A chromosome, explained 8-10.4% phenotypic variation (Li et al. 2013). A total of 33 QTL were identified across eight environments, and of these QTL, 17 stable QTL were located on chromosomes 1B, 2D, 3A, 4B, 4D, 5A, 6A, 6D, 7A, and 7B, explaining more than 95% of the observed variation of PH (Guan et al. 2018). Pang conducted a large-scale genome-wide association study using a panel of 768 wheat cultivars that were genotyped with 327 609 SNPs and detected 85 QTLs associated with PH, which distributed in 19 chromosomes except 1D and 3D chromosomes (Pang et al. 2020).
In this study, many QTLs were also detected for PH. The GA-responsive dwarfing gene Rht-1b (Wilhelm et al. 2013;Li et al. 2013) on the short arm of chromosome 4B was not detected, possibly because the materials used in this study are mainly widely used cultivars that maybe possess the Rht-B1b genes. QPH.sdau-7A.1, represented by SNP BS00078359_51, was located to a position similar to Rht 22 (Li et al. 2013), and it can be detected in 2017CK and 2018CK envirnoments. QPH.sdau-7B.1 and QPH.sdau-7B.2 might be similar to Rht 13 on chromosome 7BL . The FIRILT, SECILT, and THILT internodes may have an association with maturity in reproductive growth. The QTLs QTHILT.sdau-4B, QSECILT.sdau-4B, and QFIRILT. sdau-4B, which could be represented by SNP Bob-White_c5694_1201 located on 4B, were detected for PH, FIRILT, SECILT, and THILT, which could be the same loci reported in a previous study on maturity (Zou et al. 2017).

QTLs for the development of different internodes
The development period of wheat PH is dynamic. Genetic control of PH in bread wheat is complex, and most chromosomes harbor factors that can affect it. Stem elongation occurs in an ordered sequence. When the lower internode elongates to half of its final length, the one above it begins to elongate, then the next, and so on.
In the present study, QTLs were analyzed for PH and the lengths of each internode. When the lowermost internode of wheat (FIFILT) was elongating, the QTLs that were associated with FIFILT were probably activated. As FIFILT reached half of its final length, FORILT began to elongate, and the QTLs that were associated with FORILT were activated. This phenomenon continued with the development of other internodes. In the present study, the lengths of the different internodes were correlated with each other. This means that a certain QTL might not only control a certain length of a particular internode but might also control other internodes. In other words, some QTLs were temporarily activated in the development of one internode, while some QTLs were Page 9 of 15 146 Vol.: (0123456789) consistently activated in the development of several internodes.

QTLs responded to low nitrogen level
Nitrogen is one of the most important inputs, and it is commonly the most deficient nutrient in soils (Humbert et al. 2013). PH in wheat is a complex trait; its components include SL and each internode length. As shown in Fig. 5, THILT, FORILT, and FIFILT accounted for approximately 30% of the PH length. Combined with the development period of wheat (White and Edwards 2008), the growth of these three internodes during the period of vegetative growth leads to reproductive growth. FORILT and FIFILT showed a significant difference between the CK and LN environments in our study. Combined with the different nitrogen contents in the soil, we could infer that in the period of vegetative growth to reproductive growth, internode lengths were affected by the nitrogen content to some degree. Moreover, the PH in the LN environment was shorter than that in the CK environment, and the LRS in the LN environment was stronger than that in the CK environment. It could be presumed that under a degree of low nitrogen level, the PH was reduced, and the stem became stronger. As lodging has occurred frequently in recent years, an appropriate reduction in nitrogen could increase the lodging resistance. In this study, three QTLs for FORILT in response to nitrogen stress were detected on 3A, 6A, and 6D chromosomes. Among them, the nearby area of 6A loci has been reported to LN stress amd N use efficiency and it maybe a same QTL. A total of 11 QTLs for FIFILT were detected on 1A, 1B, 1D, 2A, 2B, 3A, 3B, 4A, 5B and 7B chromosomes and among them the loci on 4A was similar to qPhdv-4A ). One QTL for SL was detected on 3A; one QTL for PH was detected on 5B; and two QTLs were detected for LRS. QTLs were detected on chromosome 3A for FORILT, FIFILT, while in the nearby area QYd-3A has been reported to increase yeild in LN, probably is a similar locus (Cui et al. 2014). Qdelt_aver_LRS.sdau-5B, which was detected for LRS, had the lowest P-value and was located at 520,552,412 bp on chromosome 5B. TraesCS5B01G334500 (UDP-glycosyltransferase) was found near this location (upstream within 5 MB) and is considered a key gene for LRS in response to low nitrogen level. UDP-glycosyltransferase, a key enzyme that may be involved in starch synthesis, is directly limited by the carbohydrate supply and has been shown to be impacted by N supply (Jiang et al. 2004).

Conclusion
In this study, we performed GWAS for plant height with spike length, the length of the first internode to the five internodes from the top during harvest, and the lodging resistance strength at the seedling stage based on local varieties in a low nitrogen and CK environments. GWAS analysis showed that a number of QTLs were detected, including 13 QTLs for PH, 10 QTLs for SL, 19 QTLs for FIRILT, 6 QTLs for SECILT, 11 QTLs for THILT, 13 QTLs for FORILT, and 14 QTLs for FIFILT. Compared to the normal environment, the of PH, FORILT, FIFILT, and SL were significantly shorter and the LRS was significantly stronger under the low nitrogen environment. GWAS was performed to detect QTLs that responded to nitrogen stress. In total, three QTLs for FORILT in response to nitrogen stress were detected on 3A, 6A, and 6D, and 11 QTLs for FIFILT were detected on 1A, 1B, 1D, 2A, 2B, 3A, 3B, 4A, 5B and 7B. In addition, one QTL for SL was detected on 3A, one QTL for PH was detected on 5B, and there were two QTLs for LRS. These QTLs may greatly deepen the understanding of the genetic basis of trait responses to nitrogen stress. The present study provided an opportunity for the detection of QTLs induced by nitrogen fertilization as well as direct evidence to dissect traits that responded to nitrogen stress induced by environmental factors.
Author's contribution Piyi Xing performed the experiments and prepared the manuscript. Xia Zhang and Dandan Li performed partial experiments. Yinguang Bao and Honggang Wang performed partial experiments and revised the manuscript. Xingfeng Li designed the experiment and prepared the manuscript. All authors reviewed and approved the manuscript.

Ethical standards
The manuscript has not been submitted to other journals for simultaneous consideration. The submitted article is original and has not been published elsewhere in any form or language. This study was not be split up into several parts to increase the number of submissions and submitted to various journals or one journal over time. Results are presented clearly, honestly, and without fabrication, falsification, or inappropriate data manipulation (including image-based manipulation). No data, text, or theories by others are presented as if they were the author's own. Authors have permission for the use of software questionnaires surveys and scales. References articles are cited appropriate and relevant literature in support of the claims made. Authors avoid untrue statements about an entity or descriptions of their behavior or actions that could potentially be seen as personal attacks or allegations about that person. Research has no threat to public health or national security. The author group, the Corresponding Author, and the order of authors are all correct. *Authors respect all of the above guide-