Next Article in Journal
The Genetic Complexity of Prostate Cancer
Next Article in Special Issue
Transcriptome Dynamics during Black and White Sesame (Sesamum indicum L.) Seed Development and Identification of Candidate Genes Associated with Black Pigmentation
Previous Article in Journal
Role of PB1 Midbody Remnant Creating Tethered Polar Bodies during Meiosis II
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Genes
Volume 11
Issue 12
10.3390/genes11121395
genes-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide Association Study of Root System Development at Seedling Stage in Rice

by
Hongjia Zhang
1,
Mar Lar San
1,
Seong-Gyu Jang
1,
Ja-Hong Lee
1,
Na-Eun Kim
1,
Ah-Rim Lee
1,
So-Yeon Park
2,
Fang-Yuan Cao
3,
Joong-Hyoun Chin
4,* and
Soon-Wook Kwon
1,*
1
Department of Plant Bioscience, College of Natural Resources and Life Science, Pusan National University, Miryang 50463, Korea
2
National Institute of Crop Science, Rural Development Administration, Miryang 50463, Korea
3
Key Laboratory of Silkworm and Mulberry Genetic Improvement, Ministry of Agriculture, School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang 212008, China
4
Department of Integrative Biological Sciences and Industry, Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Korea
*
Authors to whom correspondence should be addressed.
Genes 2020, 11(12), 1395; https://doi.org/10.3390/genes11121395
Submission received: 7 October 2020 / Revised: 20 November 2020 / Accepted: 23 November 2020 / Published: 25 November 2020
(This article belongs to the Special Issue Genetic Research and Plant Breeding)
Browse Figures
Versions Notes

Abstract

:
Root network structure plays a crucial role in growth and development processes in rice. Longer, more branched root structures help plants to assimilate water and nutrition from soil, support robust plant growth, and improve resilience to stresses such as disease. Understanding the molecular basis of root development through screening of root-related traits in rice germplasms is critical to future rice breeding programs. This study used a small germplasm collection of 137 rice varieties chosen from the Korean rice core set (KRICE_CORE) to identify loci linked to root development. Two million high-quality single nucleotide polymorphisms (SNPs) were used as the genotype, with maximum root length (MRL) and total root weight (TRW) in seedlings used as the phenotype. Genome-wide association study (GWAS) combined with Principal Components Analysis (PCA) and Kinship matrix analysis identified four quantitative trait loci (QTLs) on chromosomes 3, 6, and 8. Two QTLs were linked to MRL and two were related to TRW. Analysis of Linkage Disequilibrium (LD) decay identified a 230 kb exploratory range for detection of candidate root-related genes. Candidates were filtered using RNA-seq data, gene annotations, and quantitative real-time PCR (qRT-PCR), and five previously characterized genes related to root development were identified, as well as four novel candidate genes. Promoter analysis of candidate genes showed that LOC_Os03g08880 and LOC_Os06g13060 contained SNPs with the potential to impact gene expression in root-related promoter motifs. Haplotype analysis of candidate genes revealed diverse haplotypes that were significantly associated with phenotypic variation. Taken together, these results indicate that LOC_Os03g08880 and LOC_Os06g13060 are strong candidate genes for root development functions. The significant haplotypes identified in this study will be beneficial in future breeding programs for root improvement.

1. Introduction

Rice (Oryza sativa L.) is one of the most significant staple food crops in the world, supporting the major energy requirements of more than half the global population [1]. Rapid societal developments have led to rice becoming a model economic crop plant, and rice yield and quality are closely associated with economic and political stability in many developing countries [2]. However, societal development is accompanied by increases in building and concomitant decreases in agricultural land availability. Moreover, climate deterioration, increasing of global temperature, and abiotic stresses, such as cold, drought, and salinity, result in the decreases in rice yield and quality [3,4,5]. Rice breeding programs to improve yield and quality are crucial to meet these challenges [6].
Rice has diverse tissues that perform different functions [7]. Root tissues play vital roles in plant growth and transport of water and mineral nutrition [8]. When plants suffer drought or other stresses, the adaptable root system supports elongation growth to obtain more water from deeper soil [9], mediated by auxin or other phytohormone responses [10]. Rice, as a monocotyledonous model plant, has root structures distinct from those of dicotyledonous models such as Arabidopsis thaliana [11]. In Arabidopsis, relatively few lateral roots are generated from the primary root. In rice, major functions are performed by the crown root instead of the primary root. The crown root generates the main lateral root, which then divides into larger and smaller lateral roots [12]. These advanced lateral roots are critical for the acquisition of moisture and nutrients and are highly adaptable to changes in the environment.
Rice growth encompasses four stages: germination, vegetative growth, reproductive growth, and maturity [13]. Root development has substantially different functions during these four stages. During the germination stage, the primary root mainly absorbs moisture and provides structural support [14]. During the vegetative stage, primary roots elongate to increase moisture absorption and root structure becomes more mature, with the crown root rather than the primary root acquiring nutrients from soil and supporting tillering [15]. During the reproductive and mature stages, seed setting rate, spikelet number, and seed weight, and therefore yield, are influenced by lateral root function [16]. However, root-shoot ratio has the opposite effects: increases in root content and root metabolism consume more photosynthetic resources, with negative effects on grain-filling and yield [16].
Research into root development has advanced considerably in the past two decades [17,18,19]. The plant hormones auxin and cytokinin were shown to have important roles in the development of root systems [20], and the first auxin root regulatory factor to be identified was Crown rootless1/Adventitious rootless1 (CRL1/ ARL1). CRL1/ARL1 was characterized as a novel transcription factor, containing an AS2/LOB domain, that responded to auxin and mediated the auxin signaling pathway to positively regulate crown root development in rice. CRL1/ARL1 mutants showed decreased numbers of lateral roots, reduced sensitivity for auxin, and reduced gravitropism [21,22]. Other regulatory factors such as OsPIN family factors OsPIN1a and OsPIN1b, which were shown to modulate auxin transport [23,24], exhibited conserved sequences. A pin1a/pin1b CRISPR/Cas9 double mutant had shorter primary roots, smaller crown roots, and reduced gravitropism. However, the phenotypes of pin1a or pin1b CRISPR/Cas9 single mutants did not differ from the wild type [25]. The CRL4/OsGNOM1 protein was shown to assist with polarity positioning of OsPIN1 [26]. OsCAND1 was found to be homologous to Arabidopsis CADN1, which participated in the development of crown root primordia via the auxin signaling pathway [27]. Overexpression of microRNA393a (miR393a) or miR393b resulted in decreased expression of auxin receptor genes OsTIR1 and OsAFB2, which disturbed the auxin signaling pathway and hampered primary crown root development [28,29]. In addition to auxin-related factors, cytokinin response factors were also found to regulate or influence root system development via cytokinin response pathways [30,31,32,33]. Expression of the metallothionein gene OsMT2b was induced by cytokinin, and overexpression led to increased adventitious root growth, larger lateral root development, and decreased cytokinin content. An RNA interference (RNAi) line showed opposite effects, indicating that OsMT2b was involved in cytokinin regulation [34]. The key regulator of adventitious root development, OsWOX11, was shown to interact with ERF3 protein, leading to expression repression of the cytokinin response factor OsRR2 and promotion of crown root growth [35].
Genome-wide association studies (GWAS) can be used for the efficient discovery of genomic loci associated with phenotypic traits [36]. The rapid advancement of next-generation sequencing (NGS) technologies has allowed single nucleotide polymorphisms (SNPs), rather than simple sequence repeat markers, to become the marker of choice for use in genotyping studies with GWAS [37]. The large amounts of data that can be accumulated using NGS allows comprehensive genome coverage and confident SNP identification. Several GWAS studies have been published that accurately detected variation in cultivars and phenotypes [38,39].
In this study, 137 rice cultivars were collected and their maximum root length and total root weight at the two-leaf stage were assessed. GWAS analysis of root traits involving two million high-quality SNPs, Principal Component Analysis (PCA), and kinship matrices, identified five previously characterized and four novel genes as candidates for root-related traits. Promoter sequence analysis indicated that two of the candidate genes contained noteworthy SNPs in root-related motifs, and haplotype analysis of these genes showed natural variation associated with rice subspecies.

2. Materials and Methods

2.1. Plant Materials and Whole Genotype Collection

In total, 137 rice accessions were collected and used in this study (Supplementary Table S1): 19 tropical japonica, 62 temperate japonica, 43 indica, eight aus, three aromatic, and two admixture. Rice accessions were collected from 28 countries and were selected from 25,604 rice accessions in the Genebank Information Center of the Rural Development Administration (RDA-Genebank, Republic of Korea) [40].
The rice accessions were sequenced using an Illumina HiSeq 2500 Sequencing Systems Platform (Illumina Inc., San Diego, CA, USA). Average genome coverage was 8× and filtered reads were aligned to the rice reference genome (IRGSP 1.0). The following parameters were used to filter genotypes for GWAS: minor allele frequency (MAF) > 5%, missing data < 1%, and heterozygosis ratio < 5%, as implemented using Plink software [41]. Finally, approximately 2 million SNPs were selected from 6.5 million raw data SNPs.

2.2. Evaluation of Root System Development

To evaluate the root system, the 137 rice accessions were planted in 96-well PCR plates and cultivated in a hydroponics system, with three replicates in random arrangements. Ten seeds were used per accession for each replicate. Seed dormancy was broken by incubation at 50 °C for 4 days, after which seeds were sterilized for 10 min with 10% hypochlorite and then washed with double distilled water three times. Seeds were germinated at 28 °C for 3 days on filter paper in Petri dishes and then moved to 96-well PCR plates for growth to the seedling stage. PCR plates were incubated at 28 °C/26 °C, 12 h/12 h (day/night) at 70% humidity. Yoshida culture solution was used and was replaced daily during seedling growth [42].
Seedlings were assessed after 10 days of incubation. Maximum root length (MRL) was measured in the ten seedlings for each accession. Outlier measurements were discarded, and the average length of the remaining seedlings was determined. The average values from each of the three replicate experiments were combined to produce a final phenotype value. Next, whole seedling roots were dried at 65 °C for 5 days and used to determine total root weight (TRW). As with root length, average root weights were used to determine a phenotype value [43]. Correlation analysis was performed using Origin software to detect correlation of MRL and TRW. Bar plots showed differences according to rating levels, and box plots showed variance between diverse subgroups.

2.3. Population and Genotype Analysis

A series of genotype analyses were performed to complement the GWAS analysis and detect differences between populations. Population structure analysis was conducted using ADMIXTURE [44], with subgroups assigned according to delta K value, and results were visualized using Pophelper [45]. Plink software was used to generate a PCA matrix, and the optimized number (6) of principal components (PCs) was used as a Q-matrix for GWAS correction. PCA plot visualization was performed in R. Neighbor-joining trees (NJ-Tree) were generated using MEGA X [46], with the Newick format file output used for modified visualization in the online tool Interactive Tree of Life (iTOL) [47]. Linkage disequilibrium (LD) decay analysis to identify candidate regions was performed using PopLDdecay [48]. The average R2 values of pairwise SNP markers were calculated for all SNPs in the genome, and the candidate region was identified where average R2 decreased to half of the maximum value.

2.4. GWAS Analysis

GWAS analysis and plot visualization was performed using the GAPIT package (version 3.0) in R [49]. The Fixed and random model Circulating Probability Unification (FarmCPU) model was adopted, with PCA (Q-matrix) and kinship (K-matrix) as covariates, to correct GWAS. Due to the fact that many SNPs contain strong LD in genotype data, the thresholds decided by the total number of SNPs were too rigorous for detection of association loci [50], thus the genotype was filtered by software PLINK. Non-independence SNPs were removed and a total of 97,469 effective and independent SNPs remained, association threshold was calculated by the formula: “−log10 (1/number of independent SNPs)” [51]. Finally, the threshold was set as −log10(P) = 4.989 for identification of association loci, and SNP markers located at locus peaks were designated as lead SNPs for the detected loci. LD decay analysis identified 230 kb around the lead SNPs as candidate genomic regions for gene identification.

2.5. Bioinformatics Analysis and Candidate Gene Prediction

Candidate genomic regions were identified from LD decay and GWAS analyses. Functional gene annotations were assigned to candidate regions from the rice genome website Gramene (www.gramene.org/). RNA-seq data from the GEO database in NCBI (www.ncbi.nlm.nih.gov) were used to assess tissue-specific expression of candidate genes. GEO accessions followed GSE6893 and GSE56463, each sample was analyzed by three replicates, respectively. Data results of root, leaf, shoot apical meristem (SAM), inflorescence, seed, mature seed from dataset GSE6893, analysis process combination with results of seeds, flower buds, flowers, flag leaves, and root before and after flowering were from dataset GSE56463. Gene ontology (GO) annotations for GO enrichment analysis were obtained from agriGO (systemsbiology.cau.edu.cn/agriGOv2). Heatmaps and bubble plots were constructed using TBtools and R, respectively [52].

2.6. RNA Extraction, cDNA Synthesis, and Expression Analysis

Gene expression was assessed in three rice varieties exhibiting extreme ‘favorable’ root-related traits and three exhibiting extreme ‘unfavorable’ traits. Total RNA was extracted from seedling roots using a RNeasy Plant Mini Kit (QIAGEN) and samples were treated with RNase-Free DNase (QIAGEN) to remove genomic DNA. Complementary DNA (cDNA) was synthesized using a SuperScript III Kit (Thermo), with primers for target genes designed using Primer5 (Supplementary Table S2). SYBR Green PCR Kit (QIAGEN) reagents were used for qRT-PCR analysis. Reactions were performed in triplicate, and actin expression was used as an internal baseline control. The −ΔΔCt method was used to calculate relative gene expression [53].

2.7. Sequence and Haplotype Analysis

Sequence and haplotype analysis were performed for each gene to assess genetic and subspecies variation. For haplotype analysis, all SNP markers from the gene were used, but missing and heterozygote data were excluded. Average score and variety count were determined from phenotype data for each subspecies, and haplotypes were identified that were significantly associated with phenotype. LD analysis was conducted using HaploView 4.2 [54]. SNP markers that were favorable for haplotype analysis were also used in LD analysis for analysis of conservation level. Haplotype variation analysis was performed using PopART software [55]. The online tool Gene Structure Display Server 2.0 [56] was used for visualization of gene structure and SNP position.

3. Results

3.1. Phenotype Evaluation of Root System Development

Six rice subspecies encompassing 137 rice varieties were used for assessment of root growth and development: Tropical japonica (Trj), Temperate japonica (Tej), Indica (Ind), Aus, Aromatic (Aro), and Admixture (Adm). Maximum root length (MRL) followed a continuous distribution with a range of 29.8–104.2 mm (Figure 1A). Total root weight (TRW) was in the 7.8–36 mg−1 range and also followed a continuous distribution (Figure 1B). These results indicated that MRL and TRW were quantitative traits controlled by multiple genes. Subspecies variation of MRL was assessed. The Tej group had slightly longer roots than the Aus and Aro groups, and substantially longer roots than groups Trj and Ind (Figure 1C). TRW showed similar trends, with the exception that root weight in group Ind was similar to Tej (Figure 1D).
Five varieties with well-developed roots (‘favorable’) and five with poorly developed roots (‘unfavorable’) were examined further, and the MRL and TRW phenotypes exhibited significant differences (Figure 2B). These differences were observed during the seed germination stage as well as in the seedling stage (Figure 2A). Correlation analysis showed a slight positive correlation between MRL and TRW (Supplementary Figure S1), with a correlation coefficient of 0.206 (p < 0.001). These results suggested that some varieties would be able to absorb more nutrients from the earliest germination stage. The favorable root development phenotype would be predicted to support improved growth during the seedling stage and in subsequent developmental stages. Correlation between MRL and TRW suggested that some varieties benefited from root branch growth and that vigorous lateral root development was beneficial to the development of the root system.

3.2. Population Structure and LD Decay Analysis

Population structure was assessed using PCA and neighbor-joining tree analyses. Cross-validation (CV) analysis indicated that K = 6 was the optimal population grouping, with the lowest error ratio compared with other K values (Figure 3A). A structure plot for the K = 6 structure was generated using Pophelper (Figure 3B). As expected, most of the rice varieties were distinguished according to their subspecies (Supplementary Figure S2). PCA (total PCs = 6) of all the genotype data was performed using Plink software. PC1 and PC2 explained most of the variation (61.86% and 25.12%) and were selected for visualization. Significant clusters were observed for rice subspecies groups Tej, Trj, Ind, and Aus, whereas the Aro and Adm groups exhibited more ambiguous separation (Figure 3C). NJ-Tree analysis of genetic distance exhibited similar results, with most varieties clearly clustered by subspecies and only Adm and Aro varieties showing dispersion among diverse clusters. LD decay was analyzed in the full panel of rice varieties. The R2 declined with increasing physical distance. The threshold value for candidate regions was determined as half the maximal R2 value, 0.26, which produced a genomic candidate region of 230 kb (Supplementary Figure S3).
Taken together, the structure, PCA, and NJ-Tree analyses divided the 137 rice accessions into six optimal subgroups. With the exception of a few varieties, the majority of accessions were grouped according to their subspecies. Candidate regions for further GWAS analysis were identified.

3.3. GWAS Analysis

GWAS analysis was performed using the GAPIT package in R. The previous results of PCA and Kinship as the Q-Matrix and K-Matrix respectively, as the covariate co-analyzed by combination with genotype and phenotype. Manhattan plots for MRL and TRW are shown in Figure 4. Four association loci, two each for MRL and TRW, were identified (Table 1). The two association loci for MRL, named qMRL6 and qMRL8, were located on chromosomes 6 and 8, and had −log10(P) of 5.135 and 5.059, respectively. The two association loci for TRW, named qTRW3 and qTRW6, were located on chromosomes 3 and 6 and had −log10(P) of 5.148 and 5.967, respectively. LD decay analysis indicated that the MRL and TRW loci on Chr6 were in close proximity, with a significant overlapping region. The distance between the lead SNPs (chr06_7017401 and chr06_6908583) was approximately 108 kb. The association loci were searched for relevant genes, and five previously characterized genes related to root development were identified. Of these, OsMADS47, OsAFB6, and PIN1a were previously characterized as positive regulators of root development [25,57,58], and AHP1 and OsGA20ox7 were found to negatively regulate root development [59,60]. Overall, this analysis identified four association loci associated with root development and identified several genes responsible for root development traits. Additional gene analysis suggested that several uncharacterized genes in the association loci might also be related to root development.

3.4. Candidate Genes and Bioinformatics Analysis

The GWAS analysis detected four association root-related loci on chromosomes 3, 6, and 8. Analysis of candidate genes from the association loci was conducted using RNA-seq data. The LD decay results focused the search for novel functional genes to a 230 kb region. For gene annotations, hypothetical proteins or non-protein coding transcripts were filtered, and then RNA-seq data were used to assess gene expression in different tissues [61,62]. Partial results are shown in Figure 5A. Genes were identified that showed significantly different expression in roots compared with other tissues and, alongside examination of gene annotations or other references studied, this allowed the identification of 44 candidate genes (Supplementary Table S3). GO enrichment analysis showed that most of the candidate genes were involved in binding and metabolic process functions (Figure 5B). Examination of published references, gene family descriptions, and homologous genes on other species was used to further refine the list of candidate genes. Finally, seven candidate genes were selected from three association loci (two loci and the overlap region).
Expression of candidate genes in rice varieties with favorable and unfavorable root traits was assessed using qRT-PCR, with the OsMADS47 gene used as a reference. Primers are shown in Supplementary Table S2. Three favorable (F1, F2, F3) and three unfavorable varieties (U1, U2, U3) that exhibited extremes of MRL and TRW were used (Figure 6). Five genes (LOC_Os08g44230, LOC_Os03g08754, LOC_Os03g08880, LOC_06g12320, and LOC_Os06g13060) were differentially expressed between the favorable and unfavorable varieties. Two candidate genes (LOC_Os08g44230 and LOC_Os06g13060) had higher expression levels in favorable varieties than in unfavorable varieties (Figure 6), and two candidate genes (LOC_Os03g08880 and LOC_06g12320) and OsMADS47 (LOC_Os03g08754) exhibited higher expression in unfavorable varieties than in favorable varieties (Figure 6). LOC_Os08g44230 encodes a zinc finger family protein of unknown function. LOC_Os03g08880 (OsPUP1) encodes a purine permease belonging to the OsPUP family. Previous research showed that OsPUP1 was expressed at high levels in roots and that another OsPUP family member, OsPUP7, was involved in growth and development in rice [63]. LOC_06g12320 encodes a transmembrane amino acid transporter of unknown function. LOC_Os06g13060 (OsDjC54) encodes a heat shock protein, and its Arabidopsis homolog, AT1G24120 (ARL1), was shown to be involved in root development [64].
The combination of RNA-seq, GO annotation, and qRT-PCR analysis, alongside previous research, facilitated the selection of four candidate root development genes from three loci.

3.5. Analysis of SNP Function in Promoter Region

Promoter involvement in expression regulation prompted the examination of SNPs in the promoter regions of candidate genes. Promoter regions were considered to be the 1500 bp immediately upstream of ATG start codon, except where intergenic regions necessitated shorter promoter sequences, as in LOC_Os03g08754, LOC_Os03g08880, and LOC_Os06g13060. Five root-related promoter motifs were identified using the online tool New PLACE: auxin response factor motifs ARFAT and ASF1MOTIFCAMV, and root development or expression motifs RSEPVGRP1, RAV1AAT, and RAV1BAT.
The root-related motifs were assessed in the promoter regions of the four candidate genes (Figure 7A), and three Sequence logo showed the structure of motifs (Figure 7B,D). Two genes (LOC_Os03g08880 and LOC_Os06g13060) harbored SNPs in root-related motifs. In LOC_Os03g08880, one SNP, A to C, was identified at promoter position 173 bp in motif 5 (ASF1MOTIFCAMV) (Figure 7B). A second SNP, C to T, was located at promoter position 309 bp in motif 1 (ARFAT) (Figure 7C). LOC_Os06g13060 contained one SNP, A to G, at promoter position 1065 bp in motif 3 (RAV1AAT) (Figure 7D). These results suggested that genes LOC_Os03g08880 and LOC_Os06g13060 were likely candidates for transcriptional regulation during root development.

3.6. Haplotype Analysis of Reported and Novel Candidate Genes

In total, five characterized and two novel candidate genes were identified by GWAS. Next, haplotype analysis was performed using genotyping data. After removal of heterozygote and missing data, SNPs from promoters and intragenic regions were used for haplotype, LD, and haplotype variation analysis. Results from two candidate genes are shown in Figure 8 and Figure 9. LOC_Os03g08880 contained four SNPs in the promoter and exon regions (Figure 8A), and the LD block showed significantly strong LD between each of the SNPs (Figure 8C). Haplotype analysis of the 136 rice varieties showed that the four SNPs divided into three haplotypes (Haps), with the maximum phenotypic TRW variation of 6.62 between Haps 2 and 3 (Figure 8B). TRW variation was significant with Hap 3 but was not significant between Haps 1 and 2. This was expected, as the major constituents of these two Haps were indica varieties (Figure 8D). LOC_Os06g13060 contained ten SNPs in the promoter, intron, and exon regions (Figure 9A). The LD block showed a similar result with LOC_Os03g08880, a strong LD level between all SNPs in this gene (Figure 9C). LOC_Os06g13060 was found in the overlapping genomic region between the MRL and TRW loci, and the ten SNPs divided into two Haps with phenotypic variations of 12.71 for MRL (Figure 9B). Hap 1, which included all the japonica and a few other varieties, showed large variation compared with Hap 2 (Figure 9D).
Four of the five characterized candidate genes exhibited significant phenotypic variation (Supplementary Figures S4–S7). With respect to LD blocks, LOC_Os08g44350 and LOC_Os03g08754 had high linkage levels between each SNP (Supplementary Figures S4C and S6C). LOC_Os08g44350 and LOC_Os08g44590 contained three and four Haps respectively, and exhibited significant phenotypic MRL variations of 28.16 and 30.04, respectively. The unfavorable Haps were comprised of indica or admixture varieties (Supplementary Figures S4B and S5B). LOC_Os03g08754 and LOC_Os06g12610 contained significant TRW phenotypic variations of 7.34 and 8.19, respectively (Supplementary Figures S6B and S7B). Large and complex variations were identified in six Haps for each of these genes, with clear grouping apparent for most of the japonica, indica, and other subspecies varieties. These results suggested that the root phenotype variation between diverse haplotypes emerged due to evolutionary processes.
In summary, haplotypes exhibiting significant root phenotype variation were identified for several candidate genes, including characterized and novel genes. These results will provide useful information towards rice breeding.

4. Discussion

4.1. Root System Development as the Crucial Trait for Rice

GWAS analysis was used in this study to detect association loci associated with maximum root length (MRL) and total root weight (TRW) traits in rice (Figure 1 and Figure 2). In rice, water capture, nutrient absorption, and metabolic processes are strongly influenced by root systems. The major root components, namely, the primary root, crown root, and lateral root, together comprise a complex root system, each part of which has diverse functions during different rice growth stages. For instance, primary roots are important for mineral nutrient acquisition and support shoot growth during the germination and young seedling growth stages [65]. The major functions of the crown and lateral roots are water and nutrient uptake and regulation of interactions between the plant and its environment [19].
Recent studies using GWAS or Quantitative Trait Loci (QTL) analysis identified several novel root-related QTLs. A GWAS analysis used 21 traits and 529 representative rice accessions and identified 143 significant associations, 25 of which were novel. The involvement of the multi-functional gene Nal1 and a novel functional gene, OsJAZ1, in the regulation of root development was demonstrated using overexpression and RNAi transgenic rice lines [9]. A study conducted by Wang and colleagues used 307 rice accessions for GWAS analysis and identified six significant QTLs associated with maximum root length, crown root number, total root length, and root tip number traits [43]. Together, these studies suggested that maximum root length and total root weight were optimal traits for identification of candidate genes for root development. Additional factors such as root responses to abiotic stresses in natural disasters were also shown to be important for rice growth. Previous GWAS analysis involving root system and abiotic stress-related phenotypic traits allowed the identification of loci and candidate genes related to root development under salt, cold, or abscisic acid stress [9]. Use of GWAS to identify useful traits will facilitate breeding programs to develop rice varieties that are more resilient to environmental changes.

4.2. Candidate Genes Identified by GWAS

In this study, GWAS analysis was performed with 137 diverse rice varieties and approximately 2 million high-quality SNPs (Figure 4). PCA matrix, kinship matrix, and LD decay analysis identified four loci and 230 kb candidate regions for root-related traits (Supplementary Figure S3). RNA-seq data, gene annotations, and previous research identified 44 candidate genes in the association loci, including five previously reported genes. Further qRT-PCR analysis showed that four of the candidate genes exhibited different expression patterns in varieties with extreme ‘favorable’ and ‘unfavorable’ root-related phenotypes (Figure 6). Ordinarily, use of structure and kinship results as covariates for GWAS analysis, with structure as the population structure for classification of subspecies, can decrease GWAS false-positives [41]. The rice collection used in this study included accessions from diverse rice subspecies (Tej, Trj, Ind, etc.), and structure analysis was therefore necessary for result correction. Kinship as the genetic relationship showed similar functions to structure.
Several candidate genes were identified in the association loci. RNA-seq data was used to assess gene expression in diverse samples, as genes were more likely to have root-related functions if significant differences in expression were observed between root and other tissues (Figure 5A). Studies of gene families or homologous genes in other species were also used as evidence for determination of promising candidate genes. Eventually, seven root-related candidate genes were identified. Five of these genes, LOC_Os08g44230, LOC_Os03g08880, LOC_Os03g08920, LOC_Os03g09170, and LOC_Os06g12320, showed significant differences in expression between root and other tissues in two RNA-seq datasets. Arabidopsis homologs of the other two genes, LOC_Os06g13060 and LOC_Os06g13090, had root-related functions in Arabidopsis. Analysis of these genes by qRT-PCR revealed that four genes, LOC_Os08g44230, LOC_Os03g08880, LOC_Os06g12320, and LOC_Os06g13060, had distinct expression patterns associated with the root traits in the diverse panel of rice accessions. These four genes were thus promising candidates for further analysis.

4.3. LD, Haplotype, and Functional SNP Analysis

Two of the four promising candidate genes, LOC_Os03g08880 and LOC_Os06g13060, contained functional SNPs in root-related promoter motifs (Figure 7). Furthermore, haplotype analysis of these two genes and five previously characterized candidate genes revealed diverse haplotypes that were associated with phenotypic variation or related to rice subspecies. Functional cis-acting elements or motifs in the promoter regions play crucial roles in signaling pathways through the activity of transcription factors that identify and bind to specific promoter motifs of downstream target genes. Research by Inukai and colleagues showed that the crl1 gene was the key regulator of lateral root development. Two auxin-response motifs (AuxREs) in the promoter region of crl1 were specifically bound by auxin-response family factors, linking regulation of root development to the auxin signaling pathway [21]. Five diverse root-related motifs were identified in the promoters of the four main candidate genes. LOC_Os03g08880 had SNPs in promoter motifs ASF1MOTIFCAMV and ARFAT, and LOC_Os06g13060 had a SNP in the RAV1AAT motif. Changes in these motifs may have altered their binding by upstream transcription factors, which would impact gene expression and lead to positive or negative changes in the regulation of root development.
In the genome region, LD was constructed by multiple SNPs. Due to the LD degraded by the crossover between genes, therefore, LD level was dependent on genome recombination rate. In a steady population, if some SNPs that located in one region never occur recombination in the genome, these SNPs will show a strong LD level between each SNP, the gene located in this region will more likely to be prominent and conserved. Popularly, LD decay analysis is aimed to find the proper LD level for detecting candidate genes, set as when the R2 decreases to half of maximum [40]. The diversity of rice accessions meant that diverse LD levels associated with genetic variation were found across the genome. LD decay was therefore assessed across the whole genome and the physical map for the appropriate regions. Also, the LD block could show the LD level between each SNP, as well the multiple SNPs make a haplotype block by strong LD level, which is contributed to identified novel haplotype and functional genes. In Yuan’s study, they identified multiple strong LD regions in candidate regions using the diverse model and diverse subspecies [66]. Two candidate genes, OsSTL1 and OsSTL2, were identified from strong LD regions in GWAS results. In the present study, as the result of LD decay analyzed, candidate genes were selected in the 230 kb region around lead SNPs, and two novel genes and four reported genes were identified. Among intragenic and promoter regions, all SNPs showed a complete strong LD level in genes LOC_Os03g08880, LOC_Os06g13060, LOC_Os08g44350, LOC_Os08g44590, and LOC_Os03g08754, which suggested that these genes contain a conserved function during diverse evolutionary or crossing processes. Haplotype analysis of most of the candidate genes showed clear grouping by rice subspecies, suggesting that the major groups emerged alongside the historical development of rice. Societal development, migration to new environments, and climate alterations may have increased the variation between subspecies and enhanced the variation in root phenotypes, but further research is needed to address this hypothesis.

5. Conclusions

In this study, maximum root length and total root weight were assessed in a panel of 137 diverse rice accessions. Four association loci associated with the root-related traits were identified by GWAS analysis. Five previously characterized genes and four novel candidate genes were identified, and SNPs in root-related promoter motifs were assessed. Haplotype analysis of the five characterized genes and two novel genes revealed diverse haplotypes associated with root phenotypic variation. This study supported beneficial information for improvement of future breeding.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4425/11/12/1395/s1, Figure S1: Correlation analysis of MRL and TRW, Figure S2: Information of population structure, Figure S3: LD decay analysis, Figure S4: Haplotype analysis of LOC_Os08g44350, Figure S5: Haplotype analysis of LOC_Os08g44590, Figure S6: Haplotype analysis of LOC_Os03g08754, Figure S7: Haplotype analysis of LOC_Os06g12610, Table S1: List of whole panel for this study, Table S2: List of primers for qRT-PCR, Table S3: List of tentative candidate genes identified by this study.

Author Contributions

J.-H.C., S.-W.K., and H.Z. experiment design and hypothesis development; H.Z. and F.-Y.C. data collection; S.-G.J., M.L.S., and A.-R.L. data analysis; J.-H.L., N.-E.K., and S.-Y.P. material development and validation; H.Z. writing; J.-H.C. and S.-W.K. supervision and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We thank Yu-Ting Zeng for doing the contribution of GWAS analysis process; Shen Hu for doing the contribution of haplotype analysis. The research grant from the Rural Development Administration, Republic of Korea (PJ01480501).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, Q.; Wing, R. Genome studies and molecular genetics: Understanding the functional genome based on the rice model. Curr. Opin. Plant Biol. 2013, 16, 129–132. [Google Scholar] [CrossRef] [PubMed]
  2. Muthayya, S.; Sugimoto, J.D.; Montgomery, S.; Maberly, G.F. An overview of global rice production, supply, trade, and consumption. Ann. N. Y. Acad. Sci. 2014, 1324, 7–14. [Google Scholar] [CrossRef] [PubMed]
  3. Sreenivasulu, N.; Butardo, V.M., Jr.; Misra, G.; Cuevas, R.P.; Anacleto, R.; Kavi Kishor, P.B. Designing climate-resilient rice with ideal grain quality suited for high-temperature stress. J. Exp. Bot. 2015, 66, 1737–1748. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Ma, Y.; Dai, X.; Xu, Y.; Luo, W.; Zheng, X.; Zeng, D.; Pan, Y.; Lin, X.; Liu, H.; Zhang, D. COLD1 confers chilling tolerance in rice. Cell 2015, 160, 1209–1221. [Google Scholar] [CrossRef] [Green Version]
  5. Ghomi, K.; Rabiei, B.; Sabouri, H.; Sabouri, A. Mapping QTLs for traits related to salinity tolerance at seedling stage of rice (Oryza sativa L.): An agrigenomics study of an Iranian rice population. OMICS 2013, 17, 242–251. [Google Scholar] [CrossRef]
  6. Quero, G.; Gutiérrez, L.; Monteverde, E.; Blanco, P.; Pérez de Vida, F.; Rosas, J.; Fernández, S.; Garaycochea, S.; McCouch, S.; Berberian, N. Genome-Wide Association Study Using Historical Breeding Populations Discovers Genomic Regions Involved in High-Quality Rice. Plant Genome 2018, 11, 1–12. [Google Scholar] [CrossRef] [Green Version]
  7. Xiao, Y.; Liu, D.; Zhang, G.; Gao, S.; Liu, L.; Xu, F.; Che, R.; Wang, Y.; Tong, H.; Chu, C. Big Grain3, encoding a purine permease, regulates grain size via modulating cytokinin transport in rice. J. Integr. Plant Biol. 2019, 61, 581–597. [Google Scholar] [CrossRef] [Green Version]
  8. Jiang, W.; Zhou, S.; Zhang, Q.; Song, H.; Zhou, D.-X.; Zhao, Y. Transcriptional regulatory network of WOX11 is involved in the control of crown root development, cytokinin signals, and redox in rice. J. Exp. Bot. 2017, 68, 2787–2798. [Google Scholar] [CrossRef]
  9. Li, X.; Guo, Z.; Lv, Y.; Cen, X.; Ding, X.; Wu, H.; Li, X.; Huang, J.; Xiong, L. Genetic control of the root system in rice under normal and drought stress conditions by genome-wide association study. PLoS Genet. 2017, 13, e1006889. [Google Scholar] [CrossRef] [Green Version]
  10. Okushima, Y.; Overvoorde, P.J.; Arima, K.; Alonso, J.M.; Chan, A.; Chang, C.; Ecker, J.R.; Hughes, B.; Lui, A.; Nguyen, D. Functional genomic analysis of the AUXIN RESPONSE FACTOR gene family members in Arabidopsis thaliana: Unique and overlapping functions of ARF7 and ARF19. Plant Cell 2005, 17, 444–463. [Google Scholar] [CrossRef] [Green Version]
  11. Coudert, Y.; Le, V.A.T.; Adam, H.; Bès, M.; Vignols, F.; Jouannic, S.; Guiderdoni, E.; Gantet, P. Identification of CROWN ROOTLESS 1-regulated genes in rice reveals specific and conserved elements of postembryonic root formation. New Phytol. 2015, 206, 243–254. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Li, J.; Zhu, S.; Song, X.; Shen, Y.; Chen, H.; Yu, J.; Yi, K.; Liu, Y.; Karplus, V.J.; Wu, P. A rice glutamate receptor–like gene is critical for the division and survival of individual cells in the root apical meristem. Plant Cell 2006, 18, 340–349. [Google Scholar] [CrossRef] [Green Version]
  13. Das, S.; Parida, S.K.; Agarwal, P.; Tyagi, A.K. Transcription factor OsNF-YB9 regulates reproductive growth and development in rice. Planta 2019, 250, 1849–1865. [Google Scholar] [CrossRef] [PubMed]
  14. Kikui, S.; Sasaki, T.; Maekawa, M.; Miyao, A.; Hirochika, H.; Matsumoto, H.; Yamamoto, Y. Physiological and genetic analyses of aluminium tolerance in rice, focusing on root growth during germination. J. Inorg. Biochem. 2005, 99, 1837–1844. [Google Scholar] [CrossRef] [PubMed]
  15. Zhang, Y.; Wang, Z.; Li, L.; Zhou, Q.; Xiao, Y.; Wei, X.; Zhou, M. Short-term complete submergence of rice at the tillering stage increases yield. PLoS ONE 2015, 10, e0127982. [Google Scholar] [CrossRef] [PubMed]
  16. Huang, M.; Chen, J.; Cao, F.; Jiang, L.; Zou, Y. Root morphology was improved in a late-stage vigor super rice cultivar. PLoS ONE 2015, 10, e0142977. [Google Scholar] [CrossRef]
  17. Sazuka, T.; Kamiya, N.; Nishimura, T.; Ohmae, K.; Sato, Y.; Imamura, K.; Nagato, Y.; Koshiba, T.; Nagamura, Y.; Ashikari, M. A rice tryptophan deficient dwarf mutant, tdd1, contains a reduced level of indole acetic acid and develops abnormal flowers and organless embryos. Plant J. 2009, 60, 227–241. [Google Scholar] [CrossRef]
  18. Behringer, C.; Schwechheimer, C. B-GATA transcription factors–insights into their structure, regulation, and role in plant development. Front. Plant Sci. 2015, 6, 90. [Google Scholar] [CrossRef] [Green Version]
  19. Zhao, H.; Ma, T.; Wang, X.; Deng, Y.; Ma, H.; Zhang, R.; Zhao, J. OsAUX 1 controls lateral root initiation in rice (Oryza sativa L.). Plant Cell Environ. 2015, 38, 2208–2222. [Google Scholar] [CrossRef]
  20. Müller, B.; Sheen, J. Cytokinin and auxin interaction in root stem-cell specification during early embryogenesis. Nature 2008, 453, 1094–1097. [Google Scholar] [CrossRef] [Green Version]
  21. Inukai, Y.; Sakamoto, T.; Ueguchi-Tanaka, M.; Shibata, Y.; Gomi, K.; Umemura, I.; Hasegawa, Y.; Ashikari, M.; Kitano, H.; Matsuoka, M. Crown rootless1, which is essential for crown root formation in rice, is a target of an AUXIN RESPONSE FACTOR in auxin signaling. Plant Cell 2005, 17, 1387–1396. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Liu, H.; Wang, S.; Yu, X.; Yu, J.; He, X.; Zhang, S.; Shou, H.; Wu, P. ARL1, a LOB-domain protein required for adventitious root formation in rice. Plant J. 2005, 43, 47–56. [Google Scholar] [CrossRef] [PubMed]
  23. Xu, M.; Zhu, L.; Shou, H.; Wu, P. A PIN1 family gene, OsPIN1, involved in auxin-dependent adventitious root emergence and tillering in rice. Plant Cell Physiol. 2005, 46, 1674–1681. [Google Scholar] [CrossRef] [PubMed]
  24. Sun, H.; Tao, J.; Bi, Y.; Hou, M.; Lou, J.; Chen, X.; Zhang, X.; Luo, L.; Xie, X.; Yoneyama, K. OsPIN1b is involved in rice seminal root elongation by regulating root apical meristem activity in response to low nitrogen and phosphate. Sci. Rep. 2018, 8, 1–11. [Google Scholar] [CrossRef]
  25. Li, Y.; Zhu, J.; Wu, L.; Shao, Y.; Wu, Y.; Mao, C. Functional divergence of PIN1 paralogous genes in rice. Plant Cell Physiol. 2019, 60, 2720–2732. [Google Scholar] [CrossRef]
  26. Kitomi, Y.; Ogawa, A.; Kitano, H.; Inukai, Y. CRL4 regulates crown root formation through auxin transport in rice. Plant Root 2008, 2, 19–28. [Google Scholar] [CrossRef] [Green Version]
  27. Wang, X.-F.; He, F.-F.; Ma, X.-X.; Mao, C.-Z.; Hodgman, C.; Lu, C.-G.; Wu, P. OsCAND1 is required for crown root emergence in rice. Mol. Plant 2011, 4, 289–299. [Google Scholar] [CrossRef] [Green Version]
  28. Bian, H.; Xie, Y.; Guo, F.; Han, N.; Ma, S.; Zeng, Z.; Wang, J.; Yang, Y.; Zhu, M. Distinctive expression patterns and roles of the miRNA393/TIR1 homolog module in regulating flag leaf inclination and primary and crown root growth in rice (Oryza sativa). New Phytol. 2012, 196, 149–161. [Google Scholar] [CrossRef]
  29. Xia, K.; Wang, R.; Ou, X.; Fang, Z.; Tian, C.; Duan, J.; Wang, Y.; Zhang, M. OsTIR1 and OsAFB2 downregulation via OsmiR393 overexpression leads to more tillers, early flowering and less tolerance to salt and drought in rice. PLoS ONE 2012, 7, e30039. [Google Scholar] [CrossRef]
  30. Debi, B.R.; Taketa, S.; Ichii, M. Cytoinin inhibits lateral root initiation but stimulates lateral root elongation in rice (Oryza sativa). J. Plant Physiol. 2005, 162, 507–515. [Google Scholar] [CrossRef]
  31. Cheng, X.; Jiang, H.; Zhang, J.; Qian, Y.; Zhu, S.; Cheng, B. Overexpression of type-A rice response regulators, OsRR3 and OsRR5, results in lower sensitivity to cytokinins. Genet. Mol. Res. 2010, 9, 59. [Google Scholar] [CrossRef] [PubMed]
  32. Kitomi, Y.; Ito, H.; Hobo, T.; Aya, K.; Kitano, H.; Inukai, Y. The auxin responsive AP2/ERF transcription factor CROWN ROOTLESS5 is involved in crown root initiation in rice through the induction of OsRR1, a type-A response regulator of cytokinin signaling. Plant J. 2011, 67, 472–484. [Google Scholar] [CrossRef] [PubMed]
  33. Zhao, Y.; Hu, Y.; Dai, M.; Huang, L.; Zhou, D.-X. The WUSCHEL-related homeobox gene WOX11 is required to activate shoot-borne crown root development in rice. Plant Cell 2009, 21, 736–748. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Yuan, J.; Chen, D.; Ren, Y.; Zhang, X.; Zhao, J. Characteristic and expression analysis of a metallothionein gene, OsMT2b, down-regulated by cytokinin suggests functions in root development and seed embryo germination of rice. Plant Physiol. 2008, 146, 1637–1650. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Zhao, Y.; Cheng, S.; Song, Y.; Huang, Y.; Zhou, S.; Liu, X.; Zhou, D.-X. The interaction between rice ERF3 and WOX11 promotes crown root development by regulating gene expression involved in cytokinin signaling. Plant Cell 2015, 27, 2469–2483. [Google Scholar] [CrossRef] [PubMed]
  36. Ma, X.; Feng, F.; Zhang, Y.; Elesawi, I.E.; Xu, K.; Li, T.; Mei, H.; Liu, H.; Gao, N.; Chen, C. A novel rice grain size gene OsSNB was identified by genome-wide association study in natural population. PLoS Genet. 2019, 15, e1008191. [Google Scholar] [CrossRef]
  37. Petersdorf, E.W.; O’hUigin, C. The MHC in the era of next-generation sequencing: Implications for bridging structure with function. Hum. Immunol. 2019, 80, 67–78. [Google Scholar] [CrossRef]
  38. Yano, K.; Morinaka, Y.; Wang, F.; Huang, P.; Takehara, S.; Hirai, T.; Ito, A.; Koketsu, E.; Kawamura, M.; Kotake, K. GWAS with principal component analysis identifies a gene comprehensively controlling rice architecture. Proc. Natl. Acad. Sci. USA 2019, 116, 21262–21267. [Google Scholar] [CrossRef]
  39. Yano, K.; Yamamoto, E.; Aya, K.; Takeuchi, H.; Lo, P.-c.; Hu, L.; Yamasaki, M.; Yoshida, S.; Kitano, H.; Hirano, K. Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice. Nat. Genet. 2016, 48, 927. [Google Scholar] [CrossRef]
  40. Kim, T.-S.; He, Q.; Kim, K.-W.; Yoon, M.-Y.; Ra, W.-H.; Li, F.P.; Tong, W.; Yu, J.; Oo, W.H.; Choi, B. Genome-wide resequencing of KRICE_CORE reveals their potential for future breeding, as well as functional and evolutionary studies in the post-genomic era. BMC Genom. 2016, 17, 408. [Google Scholar] [CrossRef] [Green Version]
  41. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.; Bender, D.; Maller, J.; Sklar, P.; De Bakker, P.I.; Daly, M.J. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 2007, 81, 559–575. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Yoshida, S. Routine Procedure for Growing Rice Plants in Culture Solution; Laboratory Manual for Physiological Studies of Rice, International Rice Research Institute: Los Baños, Philippine, 1976; pp. 61–66. [Google Scholar]
  43. Wang, F.; Longkumer, T.; Catausan, S.C.; Calumpang, C.L.F.; Tarun, J.A.; Cattin-Ortola, J.; Ishizaki, T.; Pariasca Tanaka, J.; Rose, T.; Wissuwa, M. Genome-wide association and gene validation studies for early root vigour to improve direct seeding of rice. Plant Cell Environ. 2018, 41, 2731–2743. [Google Scholar] [CrossRef] [PubMed]
  44. Alexander, D.H.; Novembre, J.; Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009, 19, 1655–1664. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Francis, R.M. pophelper: An R package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 2017, 17, 27–32. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Kumar, S.; Stecher, G.; Li, M.; Knyaz, C.; Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 2018, 35, 1547–1549. [Google Scholar] [CrossRef]
  47. Letunic, I.; Bork, P. Interactive Tree of Life v2: Online annotation and display of phylogenetic trees made easy. Nucleic Acids Res. 2011, 39, W475–W478. [Google Scholar] [CrossRef]
  48. Zhang, C.; Dong, S.-S.; Xu, J.-Y.; He, W.-M.; Yang, T.-L. PopLDdecay: A fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics 2019, 35, 1786–1788. [Google Scholar] [CrossRef]
  49. Liu, X.; Huang, M.; Fan, B.; Buckler, E.S.; Zhang, Z. Iterative usage of fixed and random effect models for powerful and efficient genome-wide association studies. PLoS Genet. 2016, 12, e1005767. [Google Scholar] [CrossRef]
  50. Zhao, Y.; Zhang, H.; Xu, J.; Jiang, C.; Yin, Z.; Xiong, H.; Xie, J.; Wang, X.; Zhu, X.; Li, Y.; et al. Loci and natural alleles underlying robust roots and adaptive domestication of upland ecotype rice in aerobic conditions. PLoS Genet. 2018, 14, e1007521. [Google Scholar] [CrossRef]
  51. Yang, W.; Guo, Z.; Huang, C.; Duan, L.; Chen, G.; Jiang, N.; Fang, W.; Feng, H.; Xie, W.; Lian, X.; et al. Combining high-throughput phenotyping and genome-wide association studies to reveal natural genetic variation in rice. Nat. Commun. 2014, 5, 5087. [Google Scholar] [CrossRef]
  52. Chen, C.; Chen, H.; Zhang, Y.; Thomas, H.R.; Frank, M.H.; He, Y.; Xia, R. TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol. Plant 2020, 13, 1194–1202. [Google Scholar] [CrossRef] [PubMed]
  53. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  54. Barrett, J.C.; Fry, B.; Maller, J.; Daly, M.J. Haploview: Analysis and visualization of LD and haplotype maps. Bioinformatics 2005, 21, 263–265. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Leigh, J.W.; Bryant, D. POPART: Full-feature software for haplotype network construction. Methods Ecol. Evol. 2015, 6, 1110–1116. [Google Scholar] [CrossRef]
  56. Guo, A.-Y.; Zhu, Q.-H.; Chen, X.; Luo, J.-C. GSDS: A gene structure display server. Yi Chuan 2007, 29, 1023–1026. [Google Scholar] [CrossRef] [PubMed]
  57. Duan, K.; Li, L.; Hu, P.; Xu, S.P.; Xu, Z.H.; Xue, H.W. A brassinolide-suppressed rice MADS-box transcription factor, OsMDP1, has a negative regulatory role in BR signaling. Plant J. 2006, 47, 519–531. [Google Scholar] [CrossRef] [PubMed]
  58. He, Q.; Yang, L.; Hu, W.; Zhang, J.; Xing, Y. Overexpression of an auxin receptor OsAFB6 significantly enhanced grain yield by increasing cytokinin and decreasing auxin concentrations in rice panicle. Sci. Rep. 2018, 8, 1–11. [Google Scholar] [CrossRef] [PubMed]
  59. Sun, L.; Zhang, Q.; Wu, J.; Zhang, L.; Jiao, X.; Zhang, S.; Zhang, Z.; Sun, D.; Lu, T.; Sun, Y. Two rice authentic histidine phosphotransfer proteins, OsAHP1 and OsAHP2, mediate cytokinin signaling and stress responses in rice. Plant Physiol. 2014, 165, 335–345. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Liu, X.; Cai, W.-J.; Yin, X.; Yang, D.; Dong, T.; Feng, Y.-Q.; Wu, Y. Two SLENDER AND CRINKLY LEAF dioxygenases play an essential role in rice shoot development. J. Exp. Bot. 2020, 71, 1387–1401. [Google Scholar] [CrossRef] [PubMed]
  61. Jain, M.; Nijhawan, A.; Arora, R.; Agarwal, P.; Ray, S.; Sharma, P.; Kapoor, S.; Tyagi, A.K.; Khurana, J.P. F-box proteins in rice. Genome-wide analysis, classification, temporal and spatial gene expression during panicle and seed development, and regulation by light and abiotic stress. Plant Physiol. 2007, 143, 1467–1483. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Wang, H.; Niu, Q.W.; Wu, H.W.; Liu, J.; Ye, J.; Yu, N.; Chua, N.H. Analysis of non-coding transcriptome in rice and maize uncovers roles of conserved lnc RNA s associated with agriculture traits. Plant J. 2015, 84, 404–416. [Google Scholar] [CrossRef] [PubMed]
  63. Qi, Z.; Xiong, L. Characterization of a Purine Permease Family Gene Os PUP 7 Involved in Growth and Development Control in Rice. J. Integr. Plant Biol. 2013, 55, 1119–1135. [Google Scholar] [CrossRef] [PubMed]
  64. Cui, H.; Hao, Y.; Kong, D. SCARECROW has a SHORT-ROOT-independent role in modulating the sugar response1. Plant. Physiol. 2012, 158, 1769–1778. [Google Scholar] [CrossRef] [Green Version]
  65. Qin, H.; Wang, J.; Chen, X.; Wang, F.; Peng, P.; Zhou, Y.; Miao, Y.; Zhang, Y.; Gao, Y.; Qi, Y. Rice Os DOF 15 contributes to ethylene-inhibited primary root elongation under salt stress. New Phytol. 2019, 223, 798–813. [Google Scholar] [CrossRef]
  66. Yuan, J.; Wang, X.; Zhao, Y.; Khan, N.U.; Zhao, Z.; Zhang, Y.; Wen, X.; Tang, F.; Wang, F.; Li, Z. Genetic basis and identification of candidate genes for salt tolerance in rice by GWAS. Sci. Rep. 2020, 10, 1–9. [Google Scholar] [CrossRef] [PubMed]
Genes 11 01395 g001 550
Figure 1. Evaluation and statistical analysis of root system development in 137 rice varieties. (A) Evaluation of maximum root length (MRL). (B) Evaluation of total root weight (TRW). (C) Statistical analysis of maximum root length (MRL) according to rice subspecies. (D) Statistical analysis of total root weight (TRW) according to subspecies. Tej: temperate japonica, Trj: tropical japonica, Ind: indica, Aro: aromatic, Adm: admixture. Box plot: solid lines within boxes represent median values.
Figure 1. Evaluation and statistical analysis of root system development in 137 rice varieties. (A) Evaluation of maximum root length (MRL). (B) Evaluation of total root weight (TRW). (C) Statistical analysis of maximum root length (MRL) according to rice subspecies. (D) Statistical analysis of total root weight (TRW) according to subspecies. Tej: temperate japonica, Trj: tropical japonica, Ind: indica, Aro: aromatic, Adm: admixture. Box plot: solid lines within boxes represent median values.
Genes 11 01395 g001
Genes 11 01395 g002 550
Figure 2. Phenotypes of rice varieties with extreme root phenotypes. (A) Root development at germination stage. Upper and lower rows represent ‘favorable’ and ‘unfavorable’ root development, respectively. Scale bar = 1 cm. (B) Root development at the seedling stage. Left and right seedlings are favorable and unfavorable varieties, respectively. Scale bar = 1.5 cm.
Figure 2. Phenotypes of rice varieties with extreme root phenotypes. (A) Root development at germination stage. Upper and lower rows represent ‘favorable’ and ‘unfavorable’ root development, respectively. Scale bar = 1 cm. (B) Root development at the seedling stage. Left and right seedlings are favorable and unfavorable varieties, respectively. Scale bar = 1.5 cm.
Genes 11 01395 g002
Genes 11 01395 g003 550
Figure 3. Population structure analysis. (A) Cross-validation (CV) error of diverse groups (K). The dotted transverse line represents the lowest level. (B) Structure analysis outcome (K = 6). (C) Principal Component Analysis (PCA) (PC1 and PC2). Pink, orange, blue, green, yellow, and red represent the Tej, Trj, Ind, Aus, Aro, and Adm rice subspecies, respectively. (D) NJ-Tree of the rice population. Pink, orange, blue, green, yellow, and red represent the Tej, Trj, Ind, Aus, Aro, and Adm rice subspecies, respectively. Grey dots represent bootstrap information.
Figure 3. Population structure analysis. (A) Cross-validation (CV) error of diverse groups (K). The dotted transverse line represents the lowest level. (B) Structure analysis outcome (K = 6). (C) Principal Component Analysis (PCA) (PC1 and PC2). Pink, orange, blue, green, yellow, and red represent the Tej, Trj, Ind, Aus, Aro, and Adm rice subspecies, respectively. (D) NJ-Tree of the rice population. Pink, orange, blue, green, yellow, and red represent the Tej, Trj, Ind, Aus, Aro, and Adm rice subspecies, respectively. Grey dots represent bootstrap information.
Genes 11 01395 g003
Genes 11 01395 g004 550
Figure 4. Manhattan plots for maximum root length (MRL) and total root weight (TRW) in 137 rice varieties. (A) GWAS of MRL using the FarmCPU model. (B) GWAS of TRW using the FarmCPU model. X axis indicates physically mapped chromosomes. Y axis indicates significance as calculated by −log10 (P).
Figure 4. Manhattan plots for maximum root length (MRL) and total root weight (TRW) in 137 rice varieties. (A) GWAS of MRL using the FarmCPU model. (B) GWAS of TRW using the FarmCPU model. X axis indicates physically mapped chromosomes. Y axis indicates significance as calculated by −log10 (P).
Genes 11 01395 g004
Genes 11 01395 g005 550
Figure 5. Analysis of candidate root-related genes. (A) Gene expression in different tissues using RNA-seq data. Red and green represent the high and low expression, respectively. Seed-L represents mature seed. (B) GO enrichment analysis of 44 candidate genes.
Figure 5. Analysis of candidate root-related genes. (A) Gene expression in different tissues using RNA-seq data. Red and green represent the high and low expression, respectively. Seed-L represents mature seed. (B) GO enrichment analysis of 44 candidate genes.
Genes 11 01395 g005
Genes 11 01395 g006 550
Figure 6. Analysis of candidate genes in varieties with ‘favorable’ and ‘unfavorable’ root development with qRT-PCR. F1, F2, and F3 represent favorable varieties, and U1, U2, and U3 represent unfavorable varieties. Three biological and three technical experimental replicates were performed. Red and green represent high and low expression, respectively.
Figure 6. Analysis of candidate genes in varieties with ‘favorable’ and ‘unfavorable’ root development with qRT-PCR. F1, F2, and F3 represent favorable varieties, and U1, U2, and U3 represent unfavorable varieties. Three biological and three technical experimental replicates were performed. Red and green represent high and low expression, respectively.
Genes 11 01395 g006
Genes 11 01395 g007 550
Figure 7. Promoter analysis of four candidate genes. (A) Structure of candidate genes. Structure description followed the legend, the shorter promoter depends on distance of intergenic. Motif 1: ARFAT, Motif 2: RSEPVGRP1, Motif 3: RAV1AAT, Motif 4: RAV1BAT, Motif 5: ASF1MOTIFCAMV. (B–D) Sequence logo of root-related motifs. The height of each base corresponds to its frequency of occurrence in the population.
Figure 7. Promoter analysis of four candidate genes. (A) Structure of candidate genes. Structure description followed the legend, the shorter promoter depends on distance of intergenic. Motif 1: ARFAT, Motif 2: RSEPVGRP1, Motif 3: RAV1AAT, Motif 4: RAV1BAT, Motif 5: ASF1MOTIFCAMV. (B–D) Sequence logo of root-related motifs. The height of each base corresponds to its frequency of occurrence in the population.
Genes 11 01395 g007
Genes 11 01395 g008 550
Figure 8. Haplotype analysis of LOC_Os03g08880. (A) Schematic representation of gene structure and SNPs positions in LOC_Os03g08880. Green, gray, and blue blocks represent promoter, untranslated region (UTR), and exon region, red vertical bars represent SNPs. (B) Results of haplotype analysis. Hap: Haplotype. Letters: a, b represent different significance level at *** p < 0.001 (Duncan’s test). (C) Linkage disequilibrium (LD) analysis of SNPs in LOC_Os03g08880. D’ was used to indicate LD level, with LD blocks defined using the Confidence Intervals function in the analysis software. Red blocks indicate complete LD between each SNP. (D) Haplotype variation analysis. Colors indicate rice subspecies as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between two Haps.
Figure 8. Haplotype analysis of LOC_Os03g08880. (A) Schematic representation of gene structure and SNPs positions in LOC_Os03g08880. Green, gray, and blue blocks represent promoter, untranslated region (UTR), and exon region, red vertical bars represent SNPs. (B) Results of haplotype analysis. Hap: Haplotype. Letters: a, b represent different significance level at *** p < 0.001 (Duncan’s test). (C) Linkage disequilibrium (LD) analysis of SNPs in LOC_Os03g08880. D’ was used to indicate LD level, with LD blocks defined using the Confidence Intervals function in the analysis software. Red blocks indicate complete LD between each SNP. (D) Haplotype variation analysis. Colors indicate rice subspecies as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between two Haps.
Genes 11 01395 g008
Genes 11 01395 g009 550
Figure 9. Haplotype analysis of LOC_Os06g13060. (A) Schematic representation of gene structure and SNP positions of LOC_Os06g13060. Green, gray, and blue blocks represent promoter, untranslated regions, and exon regions. Solid lines indicate intron regions. Red vertical bars represent SNPs. (B) Results of haplotype analysis. Hap: haplotype. Letters: a, b represent different significance level at *** p < 0.001 (Duncan’s test). (C) LD analysis of SNPs in LOC_Os06g13060. D’ was used to indicate LD level, with LD blocks defined using the Confidence Intervals function in the analysis software. Red blocks indicate complete LD between each SNP. (D) Haplotype variation analysis. Colors indicate rice subspecies as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between two Haps.
Figure 9. Haplotype analysis of LOC_Os06g13060. (A) Schematic representation of gene structure and SNP positions of LOC_Os06g13060. Green, gray, and blue blocks represent promoter, untranslated regions, and exon regions. Solid lines indicate intron regions. Red vertical bars represent SNPs. (B) Results of haplotype analysis. Hap: haplotype. Letters: a, b represent different significance level at *** p < 0.001 (Duncan’s test). (C) LD analysis of SNPs in LOC_Os06g13060. D’ was used to indicate LD level, with LD blocks defined using the Confidence Intervals function in the analysis software. Red blocks indicate complete LD between each SNP. (D) Haplotype variation analysis. Colors indicate rice subspecies as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between two Haps.
Genes 11 01395 g009
Table
Table 1. Summary of genome-wide association loci for maximum root length (MRL) and total root weight (TRW).
Table 1. Summary of genome-wide association loci for maximum root length (MRL) and total root weight (TRW).
TraitsQTLsLead SNPsChrPosition−log10(P)Reported Genes
MRLqMRL6chr06_7017401670174015.135OsPIN1a
qMRL8chr08_278688398278688395.059OsAHP1, OsGA20ox7
TRWqTRW3chr03_4592524345925245.148OsMADS47, OsAFB6
qTRW6chr06_6908583669085835.967OsPIN1a
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, H.; San, M.L.; Jang, S.-G.; Lee, J.-H.; Kim, N.-E.; Lee, A.-R.; Park, S.-Y.; Cao, F.-Y.; Chin, J.-H.; Kwon, S.-W. Genome-Wide Association Study of Root System Development at Seedling Stage in Rice. Genes 2020, 11, 1395. https://doi.org/10.3390/genes11121395

AMA Style

Zhang H, San ML, Jang S-G, Lee J-H, Kim N-E, Lee A-R, Park S-Y, Cao F-Y, Chin J-H, Kwon S-W. Genome-Wide Association Study of Root System Development at Seedling Stage in Rice. Genes. 2020; 11(12):1395. https://doi.org/10.3390/genes11121395

Chicago/Turabian Style

Zhang, Hongjia, Mar Lar San, Seong-Gyu Jang, Ja-Hong Lee, Na-Eun Kim, Ah-Rim Lee, So-Yeon Park, Fang-Yuan Cao, Joong-Hyoun Chin, and Soon-Wook Kwon. 2020. "Genome-Wide Association Study of Root System Development at Seedling Stage in Rice" Genes 11, no. 12: 1395. https://doi.org/10.3390/genes11121395

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop