Genome-Wide Association Study (GWAS) Reveals The Genetic Basis of Brace Root Angle and Diameter in Maize

Background:Brace roots are an important part of the maize root system. Among brace root traits, brace root angle (BRA) and brace root diameter (BRD) are important components that affect plant growth and development. However, there are no reports on the genetic basis of maize BRA and BRD. Results:Here, a genome-wide association study (GWAS) was conducted using 508 associated populations with extensive natural variation. The broad heritability of BRA and BRD reached 0.91 and 0.82, respectively. The analysis of different subgroups showed that there were signi�cant differences in BRA traits in different subgroups, whereas there was no signi�cant difference in BRD. Evaluation of phenotypic diversity in three different environments showed that BRA and BRD exhibit a wide range of natural variability. In the GWAS, the BRA and BRD were combined with 55,8629 single nucleotide polymorphisms, and four candidate genes were found for BRA within the threshold of P < 1.78×10 -6 that were signi�cantly related to BRD). These genes may (1) participate in maize brace root cell wall synthesis through cell transport (GRMZM2G479243); (2) involve hormone signaling pathways in the horizontal expansion of brace root cells (GRMZM2G101928 and GRMZM2G174736); or (3) involve the PLETHORA (PLT1/2) gene (GRMZM2G151934) to promote stem cells and transport expanded cells to affect the growth of root meristems. Conclusions:These results provide theoretical information for understanding the genetic basis of brace root development. Further research on candidate genes will help clarify the molecular pathways regulating BRA and BRD in maize.


Background
The maize root system, as the main organ supporting the above-ground part of the plant, can seriously affect the lodging resistance of maize. Maize root lodging resistance will greatly improve the yield of grain or dry matter and the quality of plant harvest [1]. The maize root system is composed of main radicle and lateral radicle. It participates in the absorption of water, N, and other mineral nutrients [2], synthesizes physiologically active substances, and improves the lodging resistance of roots. Sex [3][4][5][6][7][8][9]under drought conditions, the maize root system has a partial impact on maintaining grain yield and providing anchorage in the soil [10,11]. The postembryonic root system is composed of continuous bud nodes and lateral roots. The bud roots or node roots formed in the ground are called crown roots, and those formed on the ground are called brace roots [10]. The brace roots of maize appear from continuous stem nodes, and their contribution to the development of maize plants depends on their related traits. These traits consist of the tier number of brace roots (TN), the number of brace roots (RN), the radius of brace roots (RBR) [12], brace root fresh weight, and brace root dry weight [13], among others. Maize brace roots can improve lodging resistance and can also have a signi cant impact on grain yield by absorbing nutrients and water under limited soil water conditions [4,14,15]. Root lodging is related to brace root traits, such as root volume [16], number of upper internodes [17], root diameter [18], and number of internodes [17]. In particular, the strength of plant anchoring is determined by the number, diameter, and root angle of brace roots [19,20]. This in turn affects the lodging resistance of the plant.
Studies have shown that plants with higher phosphorus acquisition e ciency exhibit a higher root-to-stem weight ratio [21,22], reduced root diameter [23], longer and denser root hairs [11], and increased lateral roots [24]. In addition, brace roots are different to crown roots, and the molecular mechanisms affecting brace root development are still poorly understood [25]. The maize AP2/ERF transcription factor decreased the number of brace roots by inserting a mutant RAP2.7-Mu into ZmRAP2.7, thereby revealing the function of ZmRAP2.7 in brace root development. Multi-drug and toxin extrusion (MATE) transporters may lead to the development of brace roots [26].
Genome-wide association analysis, also known as linkage disequilibrium mapping (LD mapping) or association mapping, is based on linkage disequilibrium, and the materials used must contain a wide range of natural mutations. The analysis method of nding the target trait and the gene locus that controls the trait has been widely used in genetic mechanism analysis of complex quantitative traits [27][28][29][30]. At present, the genetic basis of maize brace root grip radius, brace root layer number, and quantitative traits of brace roots [12] has been analyzed through GWAS. In addition, aranzana successfully found key genes that affect the owering period and control disease resistance in maize using GWAS analysis [31]. Cui analyzed the genetic structure of husk length, Husk layer number, husk width, and husk thickness using an inbred line composed of 508 related populations with extensive natural variation [12,32]. However, the genetic basis of BRA and BRD has not yet been studied.
The brace root angle and diameter used in this study are quantitative traits and are regulated by multiple genes. Under three environments, 508 maize inbred lines with 558,629 SNPs genotypes were subjected to GWAS analysis. The purpose of this study was to analyze BRA and the phenotypic diversity and genetic basis of BRD. This study also identi ed a series of candidate genes related to BRA and BRD, which provided a bene cial resource for further functional research.

BRA and BRD diversity and broad sense heritability
Diversity and broad heritability of BRA and BRD The phenotype data of 508 inbred lines were measured independently in Shenyang City, Liaoning Province, in 2016 and 2017 and Siping City, Jilin Province, in 2016, and the phenotypic values of three environmental eld trials were used as random effects to calculate the BLUP values (Supplementary Table S1). In the three environments, the average heritability of BRA and BRD was about 0.91% and 0.82%, respectively (Table 1). The correlation between the angle and diameter of brace roots and the other 16 agronomic traits is shown in (Fig. 2 Population structure, BRA and BRD of four major subgroups and different environments From the violin diagram of BRA of the four major subgroups ( Fig. 3-a), it can be seen that brace root angle traits were signi cantly different in different subgroups, and the mean brace root angles in SS and TST subgroups were signi cantly lower than those in the average MIXED and NSS subgroups, indicating that population structure had a certain in uence on this trait. As seen in the violin diagram of BRD of the four major subgroups ( Fig. 3-b), distributions of the MIXED, NSS, and SS subgroups were relatively similar, while the distributions of TST subgroups were special and widely distributed. In addition, there were no signi cant differences in BRD among the four subgroups, indicating that the population structure had no signi cant effect on this trait. There were signi cant differences in the angle and diameter of brace roots in different environments (Fig. <link rid=" g3">3</link>-c and 3-d).

Genome-wide association analysis
To minimize the effect of environmental variation, phenotypic BLUP values across three environments were used for association studies. GWAS was performed using a mixed linear model (MLM) and both kinship relation-ship (K matrix) and population structure (Q matrix) were taken into account to avoid spurious associations. In total, we identi ed 2 and 3 SNPs signi cantly associated with BRA and BRD, respectively ( Fig. 4; Fig. 5; Table 2). The percentage of phenotypic variation explained by the identi ed SNPs (R 2 ) for BRA and BRD were 1.9% and 14.68%, respectively ( Table 2). Candidate gene expression analysis is conducted through the analysis of the FPKM value of the candidate gene on the Maize GDB website, and the statistical analysis results were developed in the form of tissue-speci c expression heat maps. Four candidate genes related to BRA and BRD were identi ed through MLM (Table 2). Among them, four candidate genes were divided into eight functional types: anabolic pathways, protein binding, hormone signaling pathways, transport regulation, DNA binding, cell transport, structural proteins, and unknown functions (Supplementary Table S2). Among them, four candidate genes were signi cantly related to BRA and BRD (Table 2). To determine whether genes represented by signi cant SNPs were speci cally expressed in brace root angle and diameter tissues, the published RNA-Seq data set was used from tissue-speci c expression heat maps of candidate genes compiled in 11 different organs/tissues, including brace root tissues (Fig. 6). Among them, there were two candidate genes. Compared with other tissues, GRMZM2G479243 and GRMZM2G174736 had a lower expression trend in brace roots, while GRMZM2G151934 and GRMZM2G101928 had a higher expression trend in brace roots compared with other tissues. The high and low expression levels of these candidate genes in brace roots further indicate their relevance to brace roots.

Discussion
Novel revealing the genetic basis of BRA and BRD The heritability of maize brace root-related traits has been reported; the heritability of brace root layer (TN), brace root number (RN), and brace root grip radius (RBR) ranged between 0.69 and 0.80 [40]. In this study, BRA and BRD showed extensive phenotypic variation and followed a normal distribution, and showed a very signi cant positive correlation (0.27; Fig. 1), indicating that brace root traits are quantitative traits. Under multi-gene regulation, there may be co-regulated gene loci. Genetic analysis showed that BRA and BRD had high heritability (h 2 =0.91 and 0.82, respectively; Table 1), indicating that the quantitative trait is suitable for GWAS [41]. In the three environments, the genotypic effects were signi cant for BRA and BRD, demonstrating that these genes were involved in the control of BRA and BRD ( Table 1). The environment × genotype interaction also had a signi cant impact on BRA and BRD, which indicated that the degree of transmission of BRA and BRD in maize was different from one location to another (Table 1). In other words, the genetic effects of BRA and BRD involved in inheritance were different depending on the environment. Therefore, different lines can be selected for speci c environments to improve BRA and BRD in breeding projects.
Analyzing correlation between maize BRA, BRD and other agronomic traits As a component of the root system of maize, brace roots can not only play a supporting role but also have a certain impact on yield traits. As shown in (Fig. 2), BRA has a signi cant positive correlation with ear length (EL) and ear diameter (ED), with correlation coe cients of 0.16 and 0.13, respectively. BRD had a signi cant positive correlation with grain width (KW), with a correlation coe cient of 0.16. As an important nutrient, N seriously affects the growth and development of maize plants. Studies have shown that nitrogen uptake by plants affects the number, diameter, angle, volume, and dry weight of brace roots [42]. The angle of maize canopy roots and brace roots was positively correlated with rooting depth. As the rooting depth increases, the plant's ability to absorb nitrogen increases and the more nitrogen it absorbs [43]. Nitrogen can effectively provide the nutrients required by maize, and higher nitrogen utilization can increase the dry matter quality of maize plants and grain yield and provide energy for the growth and development of maize. The indicators of maize growth and development include the number of ear rows, the number of grains, ear length, and ear thickness. Therefore, by affecting the angle of brace roots, the ear length and ear thickness can be further affected in maize [44,45]. The diameter of brace roots is highly correlated with brace root fresh weight (r=0.730) and brace root dry weight (r=0.729) [13]. The increase in fresh weight and dry weight of brace roots directly affects the nitrogen uptake by maize brace roots [42], and nitrogen uptake will increase grain yield by affecting the increase in sucrose and grain starch content [46]. Grain width is an important indicator of grain yield traits [47], as explained by the signi cant positive correlation between brace root diameter and grain width (0.16).

Changes in the angle and diameter distribution of brace roots among the four major subgroups
Maize originated in the tropics and was then domesticated in subtropical and temperate regions. Therefore, the morphological structure of maize is strongly in uenced by population structure [48]. To study the in uence of population structure on BRA and BRD, the phenotypic variation of maize BRA and BRD between different subgroups was compared (Fig. <link rid=" g3">3</link>-a and 3-b). As shown in (Fig. 3-a), there were signi cant differences in brace root angles among the four subgroups. The two subgroups of NSS and SS had a similar distribution of brace root angles. The average and scale of brace root angles of the two subgroups of MIXED and TST increased, and the two subgroups of NSS and SS were comparable to the two subgroups of MIXED and TST. There were signi cant differences in brace root angles, indicating that the brace root angles of maize inbred lines of tropical and subtropical origin had a wider distribution range. This may be due to the external environment. To prevent the excessive evaporation of water from maize plants due to the high temperature in tropical regions, natural selection prefers maize brace roots with a wide range to extract more water. In contrast, due to the relatively low temperature in temperate regions, transpiration is relatively weak, and generally, the angle of brace roots has a small distribution range. Therefore, maize from temperate regions tends to have inbred lines with more concentrated angles. As shown in (Fig. 3-b), there was no signi cant difference among the four subgroups of brace root diameter. During the evolution of maize, the trait genes that control the diameter of the brace roots were distributed in each subgroup after a high degree of recombination, and during the long domestication process, the trait genes were less affected by genetic relationships. In addition, the angle and Role of genes speci cally expressed in brace roots causing morphological changes classi cation and pathways of candidate genes involved in BRA and BRD Brace roots can not only improve lodging resistance but also have a signi cant impact on grain yield by absorbing nutrients and water under limited soil moisture conditions [4,14,15,49]. However, the genetic basis of BRA and BRD is still unclear. In this study, the MLM method was used to identify four genes that were signi cantly related to BRA and BRD (Fig. 6). Functional annotations indicated that these candidate genes were mainly placed in several functional groups, such as DNA binding, hormone signaling pathways, and cell transport.
The plant cell wall is a dynamic structure, and its change is a response to the external environment. Two Leu-rich-like receptor kinases (FEI1 and FEI2) have been identi ed in Arabidopsis thaliana with a high degree of interaction with plant cell wall synthesis. The FEI system can process the turgor signal transmitted by the SOS5 system and stimulate cellulose synthase to synthesize the cell wall. The gene mutations of FEI1 and FEI2 lead to the break of the signal transmission chains, and plants sense the changes in the external environment and its biological effects on cellulose. Synthesis also plays a role in inhibition [50]. Debarati Basu (2016) compared the growth indicators of FEI mutants in different environments and found that compared with wild-type Arabidopsis thaliana, the hypocotyls of FEI mutants had shorter radicles and a larger diameter [51]. Fluorescent protein staining showed that FEI had a higher expression in root tips. FEI1 was highly expressed in the speci c expansion of Arabidopsis root cells. GRMZM2G479243 is highly homologous to the FEI coding gene AT1G31420.2, which is predicted to play a role in the synthesis of brace root cell walls in maize.
The morphogenesis of plant organs depends on the intercellular ow of plant hormones, and directional signal transduction is determined by the polar subcellular localization of the auxin transporter synthesized by the PINFORMED (auxin export carrier PIN) gene [52]. The auxin synthesized from the brace roots is transported to the roots via vascular tubes. After auxin transport reaches the root tip, it is redistributed to the elongation zone of the root system. At this time, IAA enters the cell through the AUX1 carrier in the form of auxin, then enters the cell through the PIN carrier, and nally exports in the form of ions [53]. During this process, auxin is unevenly distributed on both sides of the elongation zone under the in uence of gravity, and a geotropic reaction occurs. The phosphorylation kinetics of the PIN protein are affected by protein phosphokinase PP2A and PINOID, and the ROTUNDA3 (RON3) protein, as a regulator of PP2A phosphokinase, can play a role in regulating PP2A activity [54]. RON3 is a unique higher plant-speci c gene. The map location technology identi ed 18 genes, including AT4G24500, that were signi cantly associated with RON3, which encodes a proline-rich protein and phenotypes the RON3-1 mutant. The identi cation showed that the number of lateral roots at the seedling stage was signi cantly reduced, and geotropic growth, which is affected by gravity, was signi cantly weakened [55]. Furthermore, Zhan (2012) identi ed AT4G24500 as the SIC gene, which encodes a hydroxyproline-rich protein [56]. The SIC gene is a functional gene necessary for the plant to maintain normal growth and development. The phenotype of SIC mutants was identi ed. It showed a variety of developmental defects, including reduced plant height, delayed owering, increased leaf edge serrations, and abnormal roots. In maize, GRMZM2G174736 is a homologous gene of AT4G24500. Therefore, GRMZM2G174736 may play a role in the polar transport of auxin and affect the morphogenesis of maize roots.
Polar transport of the phytohormone auxin is the key to regulating plant growth and development. Polar transport of auxin in Arabidopsis roots requires the action of the MFS transporter, which is zinc-induced promoter-like 1 (ZIFL1). Reverse genetics showed that the ZIFL1.1 transporter regulates various root growth hormone-related processes and may indirectly regulate growth through cytokines used in the process of auxin transport; it may act by regulating the abundance of PIN2 in the plasma membrane of the root tip [57]. Auxin regulates various unrelated processes by directing cell division and expansion, such as embryo, root, and vascular bundle formation, post-embryonic organogenesis, and geotropism[58]. The chemical growth hypothesis better describes the basis for the movement of auxin cells, in which the proton gradient is mainly generated by the plasma membrane H + -ATPases between the neutral cytoplasm and the outside of acidic cells to drive the uptake and e ux of auxin [59]. This hypothesis presumes that there are local auxin in ux and e ux carriers in the plasma membrane; thus, the asymmetric positioning of the coupling between adjacent cells provides directionality for cytokine ow[60]. In addition, auxin can also control the growth of Arabidopsis roots by regulating the cellular response to GA[61]. Therefore, GRMZM2G101928, as a homologous gene of AT5G13750.1, may affect the horizontal expansion of maize brace root cells and then affect morphogenesis.
Root growth and development is controlled by root tip meristems, and the size of root tip meristems affects the size of roots. In this study, GWAS was used to identify a candidate gene related to BRD traits and DNA binding, GRMZM2G151934, whose DA1 homologous protein 2 (DA1-related protein 2, DAR2) can regulate the size of root apical meristems. Relevant studies have shown that cytokinin and auxin antagonism can affect cell proliferation and differentiation, thereby regulating the size of roots by affecting the abundance of SHORT HYPOCOTYL2 (SHY2/IAA3). SHY2 affects the distribution of auxin in root meristems by inhibiting the auxin-induced expression of the PIN-FORMED (PIN) auxin transport gene. The PLETHORA (PLT1/2) gene affects the growth of root meristems by promoting stem cells and transporting and expanding cells [62]. In addition, DA1-related protein 2 (DAR2) acts downstream of cytokinin and SHY2 but upstream of PLT1/2 to affect the size of the root meristem [63]. Therefore, GRMZM2G151934, as a homologous gene of AT4G24500.1, may play a role in the growth of brace roots and affect the morphogenesis of maize roots.

Conclusions
This study revealed the genetic structure of maize BRA and BRD and the mechanism of controlling natural variation through GWAS. The 558,629 SNPs markers from 508 inbred line genotype populations around the world showed extensive variation. Studies have shown that BRA and BRD traits have a high level of inheritance. In addition, the four candidate genes of BRA and BRD may be involved in maize brace root cell wall synthesis, lateral expansion of brace root cells, cell transport, DNA binding, and the hormone signaling pathway. Candidate genes provide valuable resources for further dissecting the molecular network that regulates BRA and BRD in maize. The identi cation of SNPs will help promote the marker-assisted selection of BRA and BRD in maize breeding programs.

Plant material
The genome-wide association population consisted of 508 inbred lines with extensive polymorphisms. Among them, there were 60 inbred lines from the American Maize Germplasm Innovation Project, 223 inbred lines from the International Center for Maize Improvement (CIMMYT), and 225 inbred lines from Chinese germplasm resources. All materials has preserved College of Biological Science and Technology, Shenyang Agricultural University. Most of the inbred lines from CIMMYT were derived from tropical or subtropical regions, and most of the inbred lines from the United States and China were of temperate origin.
The genome-wide related population was clearly divided into four subgroups: the SS subgroup had 27 inbred lines; the NSS subgroup had 70 inbred lines; the TST subgroup had 196 inbred lines; and the remaining 215 inbred lines were divided into MIXED subgroups [33,34]. SNPs data from maize SNP50Beadchip and RNA-Seq developed by Illumina provided the genotypes for all inbred lines, with a total of 558,629 SNPs markers. Each SNP marker underwent quality testing, and unquali ed SNPs were recorded and deleted [33].

Field trial
The 508 inbred lines of the genome-wide association group were grown in three environments in China: Shenyang City, Liaoning Province (SY) (123°25E, 41°48N) in 2016 and 2017, and Siping, Ji lin Province, in Northeast China in 2016. All lines were planted in a randomized complete block design and repeated twice. Each inbred line was planted in a single row, with a length of 3 meters and a width of 0.67 meters, with a 0.4-meter aisle in the middle.

Method
The angle and diameter of the root system were investigated 30 days after planting. At this time, the brace roots were in good condition without cracking or wilting, and the measured brace root traits were more accurate, improving the accuracy of the measurement and facilitating subsequent correlation analysis. Measurement selection was based on the rst level of the root system from top to bottom. Two plants with similar growth were selected for measurement in each row. BRD was measured in the middle of the brace roots using an electronic vernier caliper (mm). BRA was measured between the occurrence point of brace roots and the main stem using a protractor (degrees).

Statistical analyses
The mixed linear model (MLM) was used to calculate the analysis of variance of the brace root traits of maize in related populations: yijk = µ + e l +rk(l)+ +(fe)il + εlik, where µ represents the total average phenotype, and is the rst "i genotype effect" of the "family," e l represents the "l" environmental effect, (fe)il is the interaction effect between the "i" family and the "l" environment rk, (l) is the repeated in uence within the environment, and εlik is a random error. The "PROCMIXED" program in SAS software was used to analyze the phenotypic variation of maize brace root traits. The Best Linear Unbiased Prediction (BLUP) of maize brace root traits also used a mixed linear model, and the average value was added to the estimated value to obtain the nal BLUP value.
The generalized heritability calculation formula is as follows: h 2 = σg 2 /(σg 2 + σge 2 /e + σε 2 /re) [35], where σg 2 represents genetic variance, σge 2 represents the genotype by environment interaction variance, σε 2 represents a residual error, and e and r represent the number of environments and the number of repetitions in each environment, respectively.

Candidate gene prediction and function annotation
To select independent and signi cant candidate genes, genes with R 2 <0.2 between signi cant SNPs sites on the same chromosome were considered to not be in the same segment. The speci c physical location of SNPs refers to the v2 version of the maize B73 genome sequence. In the range of 50 kb upstream and downstream of the signi cantly related SNPs sites, genes that have been annotated or analyzed according to functional domains, such as genes in this segment or their homologous genes in Arabidopsis participating in root morphogenesis, were searched. Then, the gene was predicted as a candidate gene.
The B73 v2 versiB website. Functional annotations were performed on the selected maize BRA and BRD candidate genes based on the orthologous gene information of Arabidopsis and maize[36].

Genome-wide association map and phenotypic variation contribution of important loci
The original number and the obtained BLUP value were saved in tab delimiter format, and 558,629 SNPs markers (MAF≥0.05) were used for GWAS analysis from two genotyping platforms (RNA-seq and SNPs chip). The MLM was used to perform an association analysis on BRA and BRD[8] and to combine the in uence of population structure (K) and kinship (Q) to remove false-positive results [37]. The method by Hao (2015) was used to count the number of independent SNPs with R 2 <0.2 in all LDs [38]. According to the Bonferroni correction threshold MLM model of the same association group published in the article, P <1/n (n=558,629) [36,39]; therefore, the MLM model threshold of this research topic was P <1.79×10 −6 . The contribution rate of the SNPs phenotype was calculated using the R language "anova ()" function. After considering the population structure, the phenotypic contribution rate (R 2 ) of each signi cant SNP site was calculated using the following linear model: The phenotypic contribution rate of all maize BRA-and BRD-related SNPs was calculated using the following linear model: In the model, Y and X represent phenotypic variables and SNPs genotype variables, respectively; P is a subpopulation variable, α is a SNP effect, β is a subpopulation effect, and ε is a random effect.

Annotation of candidate genes
The SNPs with the most signi cance within the same LD block d (r 2 <0.2) was selected to represent the locus. The physical locations of the SNPs were recorded according to the B73 RefGen_v2 (www.maizesequence.org). The corresponding genes were annotated by performing BLASTP search through NCBI website, and the candidate genes were assigned into different biological processes on the basis of the literatures describing the function of their homologs in other species or the knowledge in conserved domain database (CDD).

Heat-map of candidate genes
Raw datasets of RNA-Seq from different maize tissues were downloaded from NCBI's Sequence Read Archive (SRA) database. RNA-Seq reads were aligned to the maize B73 reference genome v2 (www.maizesequence.org) using the TopHat2 pipeline with the built-in Bowtie mapping program. The unique mapped reads were counted by htseq-count (HTSeq). To normalize the RNA-Seq data across the eleven samples from different maize tissues, we used the scaling normalization method provided in the edgeR package, based on a trimmed mean of M-values algorithm to compute the scaling factors according to the library size of each sample. After edgeR normalization, the RPKM values were averaged from replicates and used in the further analysis. The values used in the Fig. 6 are the log10 transformed ratio of normalized RPKM count in husk relative to other tissues. The values greater than +2 or less than -1 are adjusted to 2 or -1, respectively.

Additional les
Additional le 1: Agricultural University. These plant materials don't include any wild species at risk of extinction. No speci c permits are required for sample collection in this study. We comply with relevant institutional, national, and international guidelines and legislation for plant study. Correlation between brace root angle, diameter, and 16 other agronomic traits.*Signi cant when P ≤0.05; **Very signi cant when P ≤0.01