Genome-Wide Association Study Reveals Candidate Genes for Root-Related Traits in Rice

Root architecture is a determinant factor of drought resistance in rice and plays essential roles in the absorption of water and nutrients for the survival of rice plants. Dissection of the genetic basis for root structure can help to improve stress-resistance and grain yield in rice breeding. In this study, a total of 391 rice (Oryz asativa L.) accessions were used to perform a genome-wide association study (GWAS) on three root-related traits in rice, including main root length (MRL), average root length (ARL), and total root number (TRN). As a result, 13 quantitative trait loci (QTLs) (qMRL1.1, qMRL1.2, qMRL3.1, qMRL3.2, qMRL3.3, qMRL4.1, qMRL7.1, qMRL8.1, qARL1.1, qARL9.1, qTRN9.1, qTRN9.2, and qTRN11.1) significantly associated with the three traits were identified, among which three (qMRL3.2, qMRL4.1 and qMRL8.1) were overlapped with OsGNOM1, OsARF12 and qRL8.1, respectively, and ten were novel QTLs. Moreover, we also detected epistatic interactions affecting root-related traits and identified 19 related genetic interactions. These results lay a foundation for cloning the corresponding genes for rice root structure, as well as provide important genomic resources for breeding high yield rice varieties.


Introduction
Rice (Oryz asativa L.) is one of the most important crops, providing staple food for more than half of world population [1]. With the rapid growth of population and increasing scarcity of water resource, abiotic stresses have posed great threats to crop yield, making it an important task to increase and stabilize food production [2]. The foundation of the continued improvement of rice cultivars is the rich genetic diversity within domesticated populations and wild relatives [3,4]. Rice is known for its rich species diversity, Xian (indica)-Geng (japonica) subspecies differentiation, and subgroup differentiation [5][6][7][8][9]. In addition, rice has undergone decades of intensive improvement since the Green Revolution of the 1950s. Ancient human domestication and modern breeding share a common factor of both selection to improve productivity and adaptation to specific environments. Strong root structure is important for the growth, yield, and stress resistance, as well as optimization of nutrient and water absorption of rice [10,11]. Therefore, it is of great significance to understand the genetic mechanism of rice root structure [11].
Root growth is a complex process affected by multiple genes, environmental factors, and hormones [12], such as auxin (IAA) and cytokinin (CTK) [13]. Indole-3-acetic acid is the main active form of IAA in plants. Mutations in genes related to IAA metabolism, transport, or signal transduction will significantly alter the root development of rice. OsAUX1 is involved in auxin polar transport. Under hydroponic conditions, the OsAUX1 mutant exhibited longer main roots, shorter root hairs, and lower sensitivity to IAA and synthetic auxin compared with the wild type [14], while opposite phenotypes were observed in OsAUX1 overexpression lines [15]. OsYUCCA1 plays an important role in the tryptophandependent IAA biosynthesis pathway. Overexpression of OsYUCCA1 increased the IAA sociated with the phenotype. Compared with QTL mapping, GWAS can explain more extensive genetic variations [40]. For example, GWAS analysis was performed on 795 rice accessions, and 44 and 97 loci for root length and diameter were identified, respectively [41]. Zhao [42] revealed that glucosyltransferase OsIAGLU regulates rice root growth in 178 rice accessions through GWAS. GWAS has also been used to identify the related genes Nal1, OsJAZ1 and Nced of three root traits [43,44]. Although several genes related to rice root length have been cloned, more candidate genes are required to meet the needs of rice breeding. Epistasis plays an important role in complex and multi-group genetics [45]. Some methods and software packages have been developed to detect the epistatic interaction in GWAS, including iLOCi [46], BiForce Toolbox [47], BOOST [48], and 3VmrMLM [49] for genome-wide research.
In this study, 391 rice accessions were used to investigate the main root length (MRL), total root number (TNR), and average root length (ARL) of rice. Combined with rice resequencing genotypes, GWAS of rice was performed using the general linear model (GLM) and mixed linear model (MLM). The newly developed 3VmrMLM software package was used to analyze the interaction of GWAS results. SNPs related to rice root length traits were identified, and the candidate genes were analyzed to further investigate the epistatic interaction between GWAS loci. In addition, the abundant genetic heterogeneity and epistatic interaction among root traits genes in rice were revealed. The findings provide an excellent resource for genetic improvement of rice root structure.

Plant Materials
A total of 391 rice accessions were obtained from 3000 Rice Genome Project (3KRGP) in 36 countries and regions and planted in paddy fields of Anhui Agricultural University (31.93 N, 117.39 E) in 2021 (Table S1) [50]. Each line was planted in three rows with eight individuals per row with spacing of 17 × 33 cm. Field management was carried out according to the local standard practices [51]. Seeds were harvested at 30 to 35 days after heading at maturity considering the inconsistent heading date, which were then dried, and seal preserved. The germination test was carried out uniformly after completion of the harvest of all the required materials.

Evaluation of Main Root Length, Total Number of Roots, and Average Root Length
The seeds of all 391 rice accessions were dried in a 42 • C oven for 7 days to break dormancy. Twenty mature and full seeds were selected, washed with distilled water, and germinated in a 32 • C incubator for 64 h. Eight seeds of each variety germinated evenly were placed in the germination box, and the germination test was carried out at 30 • C in the light incubator. On the 10th day of germination, the root length and total root number (TRN) of each material were measured [52]. The length the longest root was taken as the main root length (MRL), and the average root length (ARL) was calculated as the sum of root length/TRN. The phenotypic mean measured by each cultivar was used for subsequent analysis. The R packages of 'corrplot' and 'ggplot2' were used to draw the correlation coefficient matrix and violin graphs of phenotypic data across subgroups.

Genotyping and Population Structure Analysis
The resequencing data of 391 rice accessions have been published in NCBI, and the variation site information can be obtained by SNP-Seek database (https://snp-seek.irri. org/index.zul, accessed on 15 January 2021) [50,53]. We used PLINK (version 1.9) [54], MAF (0.05), and GENO (0.2) to further filter SNPs and screened 1,821,713 SNPs with minor allele frequency (MAF) > 5% and missing data rate (MDR) < 0.2 (Supplementary Figure S1). The distance matrix was constructed using VCF2Dis (https://www.opensourceagenda. com/projects/vcf2dis, accessed on 20 February 2021) based on filtered files, and the phylogenetic tree was constructed using iTol (https://itol.embl.de/, accessed on 15 April 2021). GCTA [55] (version 1.93) was used for principal component analysis (PCA) to estimate the number of subgroups in the GWAS panel. In order to avoid the influence of linkage SNP in the process of population structure analysis, SNP loci were filtered according to linkage disequilibrium (LD), and unlinked SNP was retained to obtain 45,311 SNPs. The population structure analysis was carried out using Admixture (https://speciationgenomics.github.io/ ADMIXTURE/, accessed on 18 May 2021) software, and the genetic structure of the whole population was predicted. In order to construct the kinship matrix required for the mixed linear model as the covariance matrix, this study used the software Tassel (version 5.0) [56] to analyze the kinship. The population structure diagram, PCA diagram, and kinship diagram were drawn by R.

Genome-Wide Association Study
A total of 1821713, 1184751, and 882400 SNPs with MAF > 5% and MDR < 0.2 were screened for association analysis of the whole population, Xian rice population, and Geng rice population. MRL, TRN, and ARL of the whole population, Xian, and Geng rice were analyzed by GWAS using R software package rMVP [57]. Two models, GLM (only considering population structure) and MLM (considering population structure and relative kinship), were used to analyze the association between phenotype and genotype datasets. The suggested p value of dominance was 1.0 × 10 −5 to control the genetic false positive error rate of the population. The Manhattan and quantile graphs of GWAS results were created using 'rMVP' package of R. The LD attenuation of different populations, different chromosomes, and different genome regions in the same chromosome was different. Hence, for a certain site, the resolution is determined by the local LD attenuation level. The average LD decay distance of accession resources was analyzed by software PopLDdecay (version 3.4), and the screening standard was the physical distance when R 2 was reduced to half of the maximum. The LD decay distance of the whole population, Xian population, and Geng population was 350 kb, 125 kb and 346 kb, respectively (Supplementary Figure S2), which is consistent with previous findings in rice LD [34]. Therefore, the associated signal may appear in the upstream 350 kb or downstream 350 kb of the pathogenic gene. Hence, we chose a 700 kb distance as the recognition overlap marker-trait associated signal range.

Selection and Analysis of Candidate Genes
In order to identify candidate genes related to MRL, TRL, and ARL, the rice genome annotation project (http://rice.plantbiology.msu.edu, accessed on 15 January 2021) was used to search candidate genes in the 700 kb genomic region of selected important SNPs. Among all candidate genes, four types of genes, including expression protein, hypothetical protein, retrotransposon, and transposon, were excluded. ANNOVAR software was used to further select the polymorphism significantly related to GWAS trait variation and predicted to induce amino acid exchange or change splicing connection as candidate genes [58].

Haplotype Analysis
SNP data of candidate genes were extracted based on genotype of SNP with MAF > 0.05. This dataset only contained double-allele SNPs. Further haplotype analysis excluded heterozygotes and missing alleles. Haplogroups consisting of fewer than 10 accessions were deleted. For the genes found in QTLs, haplotype analysis only used the non-synonymous SNPs in the coding region of these genes for haplotype analysis of R, and the student t test was used to determine whether this locus could cause changes in rice root traits. R was used to visualize the results.

Epistatic Interaction Analysis
We selected 10,801, 4227, and 6056 SNPs from GWAS (−Log 10 (p) > 3.0) of MRL, TRN, and ARL under GLM for epistasis interaction analysis. Epistasis method of 3 V multi-locus random SNP effect mixed linear model (3VmrMLM) [49] was used to detect epistatic interaction of main effect QTLs.

Population Structure and Root Traits
A phylogenetic tree and principal component analysis were used to analyze the genetic structure of SNPs in 391 rice accessions. The phylogenetic tree showed that 391 accessions were clustered into four groups (Figure 1a). Similar results were obtained by principal component analysis, and the first two principal components can explain most of the genetic variations between accessions (Figure 1b). For the first two components, most accessions were clustered into four groups. Then, the 391 rice accessions were divided into four groups, namely Xian (n = 213), Geng (n = 147), Aus (n = 21), and admix population (n = 10) ( Table 1 and Supplementary Table S1).

Genome-Wide Association Study of MRL, ARL, and TRN Traits
The general linear model (GLM) and mixed linear model (MLM) were used to conduct GWAS on three root traits (MRL, ARL, and TRN). Considering the attenuation distance of linkage disequilibrium (LD) in rice, adjacent SNPs with spans less than 700 kb were defined as single QTLs, and the SNP with the lowest p value was used as the lead SNP to reduce the redundancy of association signals of different traits. The results showed  As shown in Table 1, the mean MRL of 391 rice accessions was 10.11 cm, among which the mean MRL of Aus and Xian rice was 15.14 cm and 9.71 cm, respectively. There was no significant difference in mean MRL between Xian rice and Geng rice ( Figure 1c and Table S2). For ARL, the average value of Xian rice was 5.88 cm, whereas that of Aus was 8.68 cm, and that of Geng rice and admix was 6.39 cm and 6.88 cm, respectively ( Figure 1d). The average TRN of Aus was 3.83, which was significantly different from that of other three clusters (Figure 1e). The variation of MRL, ARL, and TRN in rice was extensive with a normal distribution (Supplementary Figure S3), indicating that the root length traits are controlled by multiple genes, which provides a basis for further understanding the genetic structure of roots.

Genome-Wide Association Study of MRL, ARL, and TRN Traits
The general linear model (GLM) and mixed linear model (MLM) were used to conduct GWAS on three root traits (MRL, ARL, and TRN). Considering the attenuation distance of linkage disequilibrium (LD) in rice, adjacent SNPs with spans less than 700 kb were defined as single QTLs, and the SNP with the lowest p value was used as the lead SNP to reduce the redundancy of association signals of different traits. The results showed that 22, 11, and 12 QTLs detected by GLM were significantly correlated with MRL, ARL, and TRN, respectively. A total of 9, 2, and 3 QTLs were detected for MRL, ARL, and TRN using the MLM approach (−log 10 (p) ≥ 5), respectively (Supplementary Tables S3 and S4). The GLM of MRL contains 88.9% of QTLs of MLM, and all QTLs of MLM of ARL and TRN are located in GLM. GLM had a high false positive rate, although MLM could effectively reduce the false positive rate. To effectively reduce the false positives, we selected the QTLs co-located by MLM and GLM as candidate regions (Table 2). In general, eight QTL regions were determined to be significantly correlated with MRL, including qMRL1.  Table S5). Five MRL-related QTLs (qMRL1.1, qMRL1.2, qMRL3.1, qMRL3.2 and qMRL3.3) and one TRN-related QTL (TRN11.1) were detected in the Xian rice population, and no QTL co-located with ARL was detected (Supplementary Figure S4 and Table S5). In Geng rice population, no QTL was detected to be co-located with that of the whole population (Supplementary Figure S5 and Table S5). There was no significant difference between Xian and Geng in three phenotypic traits (Figure 1c-e). In our study, the number of Xian rice was greater than that of japonica rice. Through GWAS of Xian and Geng rice subgroups, it can be concluded that Xian rice provides the main genetic variation for rice root length.  [42]) associated with MRL, among which five were novel QTLs (qMRL1.1, qMRL1.2, qMRL3.1, qMRL3.3, and qMRL7.1). Further annotation with the AN-NOVAR software screened the genes that can predict the induction of amino acid exchange or change the splicing point. In qMRL1.1, qMRL1.2, qMRL3.1, qMRL3.3, and qMRL7.1, there were respectively two, three, one, fifteen, and ten genes causing non-synonymous mutations (Supplementary Table S6). Interestingly, the lead SNPs of qMRL1.1 and qARL1.1 were 7 kb apart, sharing the same candidate genes (Supplementary Figure S6). The QTL region contained two candidate genes annotated as starch synthase (LOC_Os01g52250) and pathogen-related protein (LOC_Os01g53090). We then focused on LOC_Os01g52250 and performed a haplotype analysis. Twelve non-synonymous SNP loci were found in the region (30.034-30.042 Mb), and haplotypes with fewer than 10 accessions were deleted, which formed five haplotypes (Figure 3a-c). Multiple comparisons revealed that HapB had significantly longer MRL and ARL than other four haplotypes (Figure 3g,f). HapA was the main haplotype of Xian rice, with nearly exclusive distribution in Xian rice, which might be related to the significant reduction of MRL and ARL. The major difference between HapA and HapB was in two non-synonymous mutations. HapB was mainly distributed in Au's rice and Xian rice, which was related to the significant increase in MRL and ARL in Aus rice. HapC was mainly distributed in Geng rice, and the major difference from HapB was a non-synonymous mutation (Figure 3e). These results explained the response of different subpopulations to MRL and ARL variations. Similarly, we also focused on the gene LOC_Os03g22830 in the qMRL3.1 region that can cause amino acid changes, which was annotated as a zinc finger protein containing the C3HC4 domain. LOC_Os03g22830 had three haplotypes. The accessions carrying HapC showed longer MRL than those carrying HapA or HapB (Supplementary Figure S7). In the qMRL3.3 region, it is worth noting that the reported lipoxygenase gene OsLOX5 was proved to be the first key enzyme in JA biosynthesis pathway. OsLOX5 had five major haplotypes, with HapA, HapD, and HapE mainly existing in Aus, Geng and Xian rice, respectively (Supplementary Figure S8). Multiple comparisons revealed that HapA had significantly higher MRL than HapD and HapE, indicating that OsLOX5 is a causal gene in MRL. Similarly, OS-ETR4 was a cloned gene in the qMRL7.1 region, which mainly had four haplotypes. The accessions carrying HapC and HapD had significantly longer MRL than those carrying HapA and HapB (Supplementary Figure S9). These results indicated that OS-ETR4 is also involved in the regulation of MRL in rice.

Identification of Candidate Genes for ARL
In the two loci (qARL1.1 and qARL9.1) significantly associated with ARL, the candidate region in chromosome 9 was mapped from 6862 to 7642 kb (Figure 4a-c). This region contained eight candidate genes for nonsynonymous mutations in the exon region (three reported genes, two annotated as plant proteins with unknown functional domains, and three annotated as enzymes) (Supplementary Table S6). We then mainly focused on the known gene OsMYBc, which can bind to the AAANATNC sequence of OsHKT4 promoter to regulate its expression, thereby affecting the salt tolerance of rice. There were four haplotypes (Figure 4a). The accessions carrying HapA had a significantly longer ARL than the accessions carrying HapD (Figure 4e). The difference between HapD and HapB was in four non-synonymous mutations. HapD was mainly distributed in Xian rice, and HapB was mainly distributed in Geng rice (Figure 4d). The ARL of Xian rice carrying HapD was significantly lower than that of Geng rice carrying HapB (Figure 4e). Based on these data, it could be speculated that OsMYBc is a functional gene at this locus and is involved in the regulation of rice ARL.   in four non-synonymous mutations. HapD was mainly distributed in Xian rice, and HapB was mainly distributed in Geng rice (Figure 4d). The ARL of Xian rice carrying HapD was significantly lower than that of Geng rice carrying HapB (Figure 4e). Based on these data, it could be speculated that OsMYBc is a functional gene at this locus and is involved in the regulation of rice ARL.  In each haplotype network, two adjacent haplotypes are separated by mutation changes. Based on ARL (e) of OsMYBc haplotype, differences between the haplotypes were statistically analyzed using Tukey's test.

Identification of Candidate Genes for TRN
GWAS identified three TRN-related QTLs (qTRN9.1, qTRN9.2 and qTRN11.1), and screened eight (qTRN9.1), one (qTRN9.2), and 17 (qTRN11.1) non-synonymous mutation genes in 350 Kb upstream and downstream of the lead SNP in the three QTLs, respectively (Supplementary Table S6). According to the known genes controlling rice root traits, we focused on LOC_Os09g07900 in the qRN9.1 region, LOC_Os09g15400 in the qTRN9.2 region, and LOC_Os11g17080 in the qTRN11.1 region. A haplotype analysis of the known gene OsMPK15 (LOC_Os11g17080) showed that OsMPK15 was located at 179 kb of the lead SNP of qTRN11.1 (Figure 5b,c). OsMPK15 had four major haplotypes, and HapA was nearly exclusively distributed in Geng rice. The difference between HapC and HapA was in 35 non-synonymous mutations (Figure 5a,d). Multiple comparisons demonstrated that HapC had a significantly lower TRN than other three haplotypes (Figure 5e). These results suggested that OsMPK15 may be an important gene regulating TRN.
screened eight (qTRN9.1), one (qTRN9.2), and 17 (qTRN11.1) non-synonymous mutation genes in 350 Kb upstream and downstream of the lead SNP in the three QTLs, respectively (Supplementary Table S6). According to the known genes controlling rice root traits, we focused on LOC_Os09g07900 in the qRN9.1 region, LOC_Os09g15400 in the qTRN9.2 region, and LOC_Os11g17080 in the qTRN11.1 region. A haplotype analysis of the known gene OsMPK15 (LOC_Os11g17080) showed that OsMPK15 was located at 179 kb of the lead SNP of qTRN11.1 (Figure 5b,c). OsMPK15 had four major haplotypes, and HapA was nearly exclusively distributed in Geng rice. The difference between HapC and HapA was in 35 non-synonymous mutations (Figure 5a,d). Multiple comparisons demonstrated that HapC had a significantly lower TRN than other three haplotypes (Figure 5e). These results suggested that OsMPK15 may be an important gene regulating TRN.

GWAS of Rice Root Traits Reveals Abundant Epistatic Interactions
We used the recently developed 3VmrMLM software package to conduct epistatic interaction analysis on the GWAS (−log 10 (p) > 3) results of GLM. As a result, 19 significant loci were detected in MRL (n = 9), ARL (n = 6) and TRN (n = 4), among which 13 were associated with GWAS results whereas six were not ( Figure 6). The loci in the qMRL3.2 region and the position of 5,013,178 bp on chromosome 1 were found to have the greatest contribution rate to MRL (10.35%) (Tables 3 and S3). When only single loci were considered, 8, 2, and 3 QTLs were associated with MRL, ARL, and TRN, but when interactions were considered, more association loci emerged. In GWAS analysis, interaction loci are often overlooked, resulting in the absence of favorable genes. In the GWAS analysis process, when only one single locus was considered, neither of the two interacting genes had a significant effect, but the interaction between showed a significant effect. These problems could be solved by epistasis. These results show that epistatic interaction is an important part of GWAS and may help to understand the heritability of complex traits in GWAS. separated by mutation changes. Based on TRN (e) of OsMPK15 haplotype, differences between the haplotypes were statistically analyzed using Tukey's test.

GWAS of Rice Root Traits Reveals Abundant Epistatic Interactions
We used the recently developed 3VmrMLM software package to conduct epistatic interaction analysis on the GWAS (−log10(p) > 3) results of GLM. As a result, 19 significant loci were detected in MRL (n = 9), ARL (n = 6) and TRN (n = 4), among which 13 were associated with GWAS results whereas six were not ( Figure 6). The loci in the qMRL3.2 region and the position of 5,013,178 bp on chromosome 1 were found to have the greatest contribution rate to MRL (10.35%) (Tables 3 and S3). When only single loci were considered, 8,2, and 3 QTLs were associated with MRL,ARL, and TRN, but when interactions were considered, more association loci emerged. In GWAS analysis, interaction loci are often overlooked, resulting in the absence of favorable genes. In the GWAS analysis process, when only one single locus was considered, neither of the two interacting genes had a significant effect, but the interaction between showed a significant effect. These problems could be solved by epistasis. These results show that epistatic interaction is an important part of GWAS and may help to understand the heritability of complex traits in GWAS.  and TRN (c). The red arrow represents the QTLs co-localized by GLM, MLM, and epistatic interactions. The −log10(p) value is reported on the left y-axis, which was obtained from the whole genome scanning of all markers in 3VmrMLM. The LOD score is reported on the right y-axis, which was obtained from the significant likelihood ratio test and suggested the dominant locus in the second step of 3VmrMLM, with the threshold LOD = 3.0 (black dotted line). If LOD score ≥ 20, the LOD scores obtained are transformed as LOD′ = 20 + (LOD − 20)/100. These LOD scores, along with their known genes, are shown in points with straight lines. All the main-effect genes with grey and blue characters are identified in one and at least two datasets. All the genes with underline are around suggested loci.  and TRN (c). The red arrow represents the QTLs co-localized by GLM, MLM, and epistatic interactions. The −log 10 (p) value is reported on the left y-axis, which was obtained from the whole genome scanning of all markers in 3VmrMLM. The LOD score is reported on the right y-axis, which was obtained from the significant likelihood ratio test and suggested the dominant locus in the second step of 3VmrMLM, with the threshold LOD = 3.0 (black dotted line). If LOD score ≥ 20, the LOD scores obtained are transformed as LOD = 20 + (LOD − 20)/100. These LOD scores, along with their known genes, are shown in points with straight lines. All the main-effect genes with grey and blue characters are identified in one and at least two datasets. All the genes with underline are around suggested loci.

Discussion
Rice root is a very important and complex system controlling rice yield and adaptation to a complex environment [61]. However, it is difficult to observe root growth since it occurs underground and to choose the ideal plants, which may be solved by indirect selection using markers closely related to genes controlling the root traits [62,63]. In this study, MRL, ARL, and TRN of 391 rice accessions were measured, and statistical epistasis and association mapping were used to excavate the excellent alleles regulating rice root system.
GWAS is a genetic survey of the entire genome to detect variations associated with traits in natural populations and is a powerful method for profiling complex traits. Its efficiency is largely determined by marker density, population size, and statistical methods [38,39]. Previously, population sizes ranging from 200 to 3000 have been used for GWAS in rice [64][65][66]. Here, we performed GWAS on 1821713 high-confidence SNPs and MRL, ARL, and TRN of 391 rice accessions using GLM and MLM from rMVP package (Supplementary Figure S1 and Table S1). Although the natural population of 391 accessions was not large enough, there were significant phenotypic variations in root traits, as the MRL ranged from 2.64 to 24.81 cm, the ARL from 1.99 to 12.38 cm, and the TRN from 2.36 to 6.86 (Table 1 and Supplementary Figure S3). These significant phenotypic variations may be associated with high genetic diversity.
Rice roots are mainly composed of seed roots, adventitious roots, and crown roots, which are of great significance to the fixation, nutrition, and water absorption of rice plants.
Here, we detected 13 QTLs related to rice root traits (MRL (n = 8), ARL (n = 2), TRN (n = 3)) ( Table 2). In reference to the previously reported QTLs and genes, three known root traitrelated QTLs and genes were mapped in this study. To be specific, a gene OsGNOM1 affecting the adventitious root formation in rice was mapped in qMRL3.2; a gene OsARF12 mediating root elongation in rice was mapped in qMRL4.1; and a qRL8.1 site that has been proved to affect the length of main roots in rice was mapped in qMRL8.1. These results confirm the reliability of GWAS in this study and can be used for further candidate gene analysis. In addition, 10 previously undescribed QTLs were identified as novel QTLs ( Table 2).
Rice root length and root growth angle are also involved in the resistance ability of rice to drought [67], which are key traits that directly affect rice yield. In this study, 12 SNPs were found in the coding region (qMRL1.1 and qARL1.1) of LOC_Os01g52250 in Xian rice (Supplementary Figure S7 and Table S5) and the whole population ( Figure 2 and Table 2), and their variations led to differences in MRL and ARL between rice accessions. In the qMRL3.3 region, the lipoxygenase gene OsLOX5 is the first key enzyme in the JA biosynthesis pathway [68]. JA is also a signal substance in plants and plays an important role in plant development. It has been reported that rice development is controlled by JA. Therefore, it can be speculated that OsLOX5 also affects rice MRL through the JA pathway. We found an ethylene receptor gene OS-ETR4 associated with MRL in the qMRL7.1 region. The positive regulation of IAA on OS-ETR4 mRNA level is mediated by an acetylene dependent pathway [69]. It has been reported that OsRTH1 may be a homologue of RTE1 in rice, which regulates ethylene response in rice [70]. OsRTH1 overexpression in rice could significantly inhibit the ethylene-induced growth and development changes and affects the development of adventitious roots. Therefore, we speculate that OS-ETR4 is directly involved in MRL regulation by reducing ethylene sensitivity. We also identified the cloned gene OsMYB (qARL9.1) associated with rice ARL. Recent studies have reported that OsHKT4 mutant plants are highly sensitive to salt stress. Compared with the wild type, the mutant exhibited retarded growth and reduced fresh weight and stem length. OsMYBc can regulate the expression of OsHKT4 and affect the salt tolerance of rice, thereby affecting its development [71,72]. Therefore, according to our results, we speculated that OsMYBc regulates salt tolerance of rice by controlling ARL. These results provide a new theoretical basis for research on the regulation mechanism of rice root length in the future. It has been reported that the slender rod gene SLR1 encodes a DELLA protein, which is a negative regulator of GA signaling [73]. It can also integrate and amplify SA-and JA-mediated defense signals, control root length and root number of rice, and play multiple roles in regulating rice growth and innate immune response. Interestingly, we found a mitogen-activated protein kinase OsMPK15 associated with TRN in our candidate genes. LOC_Os11g17080 (OsMPK15) is a cloned gene on chromosome 11 and annotated as a mitogen-activated protein kinase. Knockout of OsMPK15 can lead to the constitutive expression of disease persistence-related (PR) genes, increase the accumulation of reactive oxygen species triggered by chitin, enhance the accumulation of SA and JA, and significantly promote the resistance of rice [74]. We speculate that OsMPK15 also regulates rice TRN by controlling SA and JA. Our results provide a reliable resource for exploring genes related to rice root length traits. In Oryza sativa, there are ancient and definitive differences between the two major subspecies, Xian and Geng, but the breeding history suggests a finer genetic architecture. In this study, there were significant differences in genetic diversity between 213 Xian and 147 Geng rice accessions. A total of 5 QTLs were detected in Xian and the total subpopulation, whereas no QTLs were detected in Geng. During the domestication process, most homozygous individuals of rice were composed of homozygous individuals, and infiltration from other populations was limited, resulting in genome diversity within and among rice populations.
Rice root is a complex system affected by genetic factors, environment, and population stratification. In GWAS, SNPs are usually detected to understand the association between the whole genome and the traits of interest. Unit point analysis is usually used to estimate the effect of single SNPs. However, it can only identify SNPs with relatively strong effects, while missing some important SNPs with minor effects [75,76]. We used the VmrMLM method to detect the epistasis of SNP loci exceeding the GLM threshold p = 1.0 × 10 −3 . We detected a total of 19 pairs of genes with significant epistasis, among which six were not detected in the two models of GWAS (Table 3). These genes may be new genes affecting rice roots. For a pair of genes, each individual gene may have weak or even no association with the trait, but their interaction is strongly associated with the trait. Overall, we detected more dominant loci related to rice roots by epistasis, providing abundant genetic resources for rice breeding.
In general, we obtained a set of loci significantly associated with rice root traits by GWAS and epistasis analysis of 391 rice accessions. The candidate genes were further screened by haplotype block structure analysis of SNPs that are significantly associated with agronomic traits and functional annotation of each gene. As expected, although a large linkage disequilibrium contains many SNPs in candidate regions, our results demonstrate that the number of candidate genes can be significantly reduced by combining haplotype block structure with gene function annotation. Further epistatic analysis improved the accuracy of genome prediction. Generally, our results should be helpful for future gene functional analysis and provide valuable information for rice gene cloning research.

Conclusions
Our identified genetic regions associated with root development traits and the newly identified loci may contribute to marker-assisted QTL pyramid / gene introgression breeding programs in the future to stabilize and improve rice yield. It also provides a rich source of information on the natural genetic variation in the evolution, domestication, and breeding of indica and japonica rice associated with seed root length and other adaptive traits. Future research can be focused on verifying the role of these candidate genes and their functional variations. Genetic transformation and DNA insertion mutation screening will be used to determine whether these genes affect rice root structure.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cimb44100301/s1, Figure S1: Distribution of single nucleotide polymorphisms (SNPs) and nucleotide diversity across the rice Nipponbare genome in the rice association panel; Table S1: Information of 391 rice accessions; Figure S2: LD decay distance estimated for 391 rice accessions; Table S2: MRL ARL and TRN mean comparison results; Figure S3: Histograms showing the data distribution of three traits of 391 rice accessions; Table S3: All significant loci of root traits detected by two models; Figure S4: The genome-wide association plots of MRL, ARL, and TRN in Xian rice population were plotted; Table S4: Significant SNPs above the threshold line in MLM and GLM; Figure S5: The genome-wide association plots of MRL, ARL, and TRN in Geng rice population were plotted using general linear model and mixed model methods; Table S5: The significant SNP loci in Xian and Geng rice populations; Figure S6: Venn diagram of the proportion of QTL controls for single and multiple traits. Table S6: All candidate gene annotations for non-synonymous mutations in exons of 10 QTLs; Figure S7: Haplotype analysis of LOC_Os03g22830; Figure S8: Haplotype analysis of LOC_Os03g49380; Figure S9: Haplotype analysis of LOC_Os07g15540.
Author Contributions: L.L. and H.Y. conducted experiments and collected data; C.Z., J.X. and N.W. data collation and statistical analysis; Z.L., J.X. and Z.Z. graphic finishing; J.X. wrote the paper; Y.S. designed the experiment, provided intellectual guidance, and wrote and reviewed the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.