Genetic Basis Underlying Tiller Angle in Rice (Oryza sativa L.) by Genome-wide Association Study

Plant materials

The 164 japonica varieties used in this study were derived from were selected from the 3K RGP (Wang et al. 2018). Detailed information regarding these varieties, including their geographical origin, is shown in Additional File 2: Table S1.

Genotyping data and SNP filtering

The 3 K RGP 4.8mio SNP dataset was downloaded from the Rice SNP-Seek Database (https://snp-seek.irri.org/) (Alexandrov et al. 2015). A total of 1,135,519 cleaned SNPs were called in the population. SNPs were filtered utilizing the software PLINK (Purcell et al. 2007) with missing rate < 20% and minor allele frequency (MAF) > 0.05.

Evaluation of tiller angle and population structure

All varieties were grown in the experimental field at density of 17 cm × 20 cm. A completely randomized block design with three replicates was performed. Six plants in the middle of every row were selected for evaluation of tiller angle at maturity stage. The mean value of three replicates was used for analysis. The field experiment was performed in December 2015 to April 2016, Sanya city, Hainan (HN) Province, China (18°15′N, 109°30′E). Principal component analysis (PCA) using PLINK (Purcell et al. 2007) was used to assess the population structure.

Genome‑wide association mapping

GWAS was performed on SNPs as described above using the factored spectrally trans- formed linear mixed models (Fast-LMM) (Lippert et al. 2011), GLM (Wang and Zhang 2021) and FarmCPU (Liu et al. 2016). Suggestive thresholds were calculated using the formula “-log10(1/effective number of independent SNPs)” as described previously (Wang et al. 2016). The SNP with the minimum P value was considered the lead SNP. Regions were considered as QTLs when SNPs exceeding the P value within a 200-kb interval of the leading SNP. Manhattan plots for the GWAS results were drawn using the R package “qqman” (https://www.r-project.org/).

Haplotype analysis and identification of candidate genes

Candidate genes were scanned within the 200 kb region centered on the lead SNP of each QTL (using the reference Nipponbare genome (http://rice.plantbiology.msu.edu/cgi-bin/gbrowse/rice/). Haplotype analysis of the candidate genes were carried out using non-synonymous SNPs in exons, while the PAY1 and TIG1 were carried out for all available SNPs in the promoter region (2 kb), non-synonymous SNPs in exons, introns, and 3′ untranslated region were selected from the Rice SNP-Seek Database (https://snp-seek.irri.org/). Only haplotypes shared by at least 5 varieties were used for multiple comparisons. Duncan’s multiple comparison tests followed by a one-way analysis of variance were completed with the R package ”agricolae” (https://myaseen208.github.io/agricolae/https://cran.r-project.org/package=agricolae).

Quantitative Reverse Transcription PCR (qRT‑PCR) Analysis

Tiller base from the varieties used in this study, were collected in triplicate at the tillering stage.

Total RNA was extracted using the Trizol reagent (Invitrogen), following the manufacturer’s instructions. cDNA was synthesized from total RNA using the Fast-Quant RT Kit (with gDNase; Tiangen). Quantitative real time PCR (qRT-PCR) was performed in a two-step reaction using SuperReal PreMix Color (SYBR Green; Tiangen) on a Roche Light Cycler 2.10 system with three technical replicates. Expression levels were normalized to OsActin1 (LOC_Os03g50885).

The sequences of the candidate genes were downloaded from the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/analyses_search_locus.shtml). Primer sequences for candidate genes were downloaded from the qPrimerDB-qPCR Primer Database (https://biodb.swu.edu.cn/qprimerdb/best-primers-ss) and designed by NCBI except for LOC_Os08g05900, which could not find the appropriate primer. Primer sequences are given in Additional File 2: Table S6.

Natural variation of rice tiller angle

To dissect the genetic structure of rice tiller angle, a total of 164 japonica varieties from 3K RGP were selected for tiller angle phenotyping at the heading stage in Hainan, China (Additional File 2: Table S1). The tiller angle observed at the maturity stage in this panel ranged from 6.3° to 35.6°, which exhibited a nearly normal distribution (Fig. 1a), indicating that tiller angle was regulated by many genes with small or moderate effects (Zhang et al. 2018; He et al. 2021). Tiller angles of about half of the panel (46.2%, 73 of 164) were between 18 and 24˚, with the median value of 19.4˚ (Fig. 1a-b). Principal component analysis (PCA) of the 164 varieties using high-confidence SNPs showed that all tested varieties are closely linked, and that there is no obvious structural pattern existed in this panel (Additional File 1: Figure S1), avoiding population structure-induced false positives. Taken together, these results indicate that the rice panel used in this study has abundant genetic diversity and is suitable for GWAS.

GWAS for tiller angle

To dissect all possible genomic loci associated with tiller angle, we performed GWAS with FastLMM, FarmCPU and GLM models using 1,135,519 high-confidence SNPs. The genome-wide suggestive thresholds P ≤ 8.81E^− 7 were calculated using the formula “-log10(1/effective number of independent SNPs)”. Totally, 18 QTLs were identified to be significantly associated with tiller angle when significant SNPs located within 200-kb range (Fig. 1). Among these QTLs, 8, 2, 8 QTLs were identified by Fast-LMM (Lippert et al. 2011), GLM (Wang and Zhang 2021) and FarmCPU (Liu et al. 2016), respectively. Notably, 2 QTLs, qTA8.3 and qTA8.4, identified in this study were co-localized with PAY1 (Zhao et al. 2015) and TIG1 (Zhang et al. 2019), respectively (Table 1), confirmed the reliability of our GWAS results. Among the remaining 16 QTLs newly identified, 2 QTLs (qTA6.1 and qTA8.1) were simultaneously identified by both GLM and FastLMM approaches, which are highly valuable for further identification of the causal genes for tiller angle trait.

PAY1, the candidate gene of qTA8.3

The QTL qTA8.3 was found to contain PAY1 gene, a single nucleotide change, G to A, at position 19468238 in exon 4 of which affected plant architecture via affecting polar auxin transport activity and altering endogenous indole-3-acetic acid distribution (Zhao et al. 2015). To understand how PAY1 sequence may affect tiller angle phenotype, we analyzed SNPs in the genomic coding region and 2 kb of upstream promoter region of PAY1 across 164 varieties, which revealed 2 major haplotypes; tiller angles of varieties containing PAY1^HapA were significantly lower than those of PAY1 ^HapB, five SNP variations including only one in the first exon were observed between these two haplotypes (Fig. 3a-b; Additional File 2: Table S1). However, none of those variations observed in these two haplotypes was the same as the functional mutation previously reported in PAY1 (Zhao et al. 2015). The SNP variations observed in both promoter region and CDS region of PAY1 indicated that both may play important role in regulating tiller angle in rice. To further investigate the causative SNP variation in PAY1 responsible for tiller angle, we used qRT-PCR to examine expression level of PAY1 between varieties harboring PAY1^HapA or PAY1^HapB, and results showed that expression levels of PAY1 are not associated with haplotypes (Fig. 3c), indicating that the only variation in the first exon, a non-synonymous SNP (T to C) identified in PAY1 could be a new functional allele that leads to spread phenotypes. In generally, japonica varieties always have a relatively smaller tiller angle and show a more compact plant architecture than indica varieties (Dong et al. 2016a). Therefore, to further understand the breeding utilization of “T” and “C” alleles in PAY1 between the japonica and indica varieties. we investigated the changes of allelic frequencies using 2584 varieties from 3K RGP (Fig. 3d and Additional File 2: Table S2). The results showed that nearly all of the indica varieties (98%) contained “C” allele, while most japonica varieties (73%) contained “T” allele, indicating that PAY1 has been selected during japonica-indica differentiation.

TIG1, the candidate gene of qTA8.4

Another QTL, qTA8.4, contained TIG1, which encodes a TCP family transcriptional activator and variations in the TIG1 promoter of indica varieties (tig1 allele) decrease expression of TIG1 in the adaxial side of tiller base and reduce cell length and tiller angle, leading to the transition from inclined tiller growth in wild rice to erect tiller growth during rice domestication (Zhang et al. 2019). To investigate the causative SNP variations in TIG1 responsible for phenotypic variations for tiller angle, we analyzed SNPs in the genomic coding region and 2 kb of upstream promoter region of TIG1 across 164 varieties, which revealed 4 major haplotypes with no SNP in CDS region (Fig. 4a; Additional File 2: Table S1). Compared with previously reported SNP variations of TIG1 in rice (Zhang et al. 2019), SNPs detected in our research were completely new. Moreover, multiple comparison tests of tiller angle showed that tiller angles of varieties with TIG1^HapA were significantly higher than those with TIG1^HapC or TIG1^HapD (Fig. 4b), indicating that SNP variations identified in promoter regions may play an important role in regulating tiller angle as a new allele. To test whether the relative expression level of different haplotypes of TIG1 were associated with tiller angle, we used qRT-PCR analysis on samples from tiller base, which revealed that expression level of TIG1 in TIG1^HapA varieties was significantly higher than that in TIG1^HapD varieties (Fig. 4c), in agreement with previously report that TIG1 positively regulates tiller angle in rice (Zhang et al. 2019).

Candidate genes analysis of qTA6.1 and qTA8.1 by qRT-PCR and haplotype analysis

To identify the candidate genes in two novel QTLs (qTA6.1 and qTA8.1) that were commonly detected by both GLM and FastLMM, we selected candidate genes in a 200-kb region (100 kb flanking the left and right side of the most significantly associated SNPs) using the Nipponbare reference genome (http://rice.plantbiology.msu.edu), and found contained 60 genes including 50 functionally annotated genes (31, 19 in qTA6.1 and qTA8,1, respectively.)and 10 retro/transposon/hypothetical proteins (Additional file 2: Table S3). Two analyses were performed further to identify candidate genes for these two QTLs: 1) expression analysis using qRT-PCR analysis in the tiller base of 22 varieties: 11 spread-tiller varieties and 11 compact-tiller varieties; 2) haplotype analyses in CDS region.

For qTA6.1, LOC_Os06g06990 was the only one showing significantly lower expression level in spread-tiller varieties than compact-tiller varieties as revealed by qRT-PCR analysis (Fig. 5a). Haplotype analysis in CDS region of all candidate genes in the same region showed that 12 genes can be divided into different haplotypes, and that none of identified SNP variations among these genes were associated with significant differences in the tiller angle (Fig. 5c and Additional File 1: Figure S2), indicating that the LOC_Os06g06990 is likely the most promising candidate gene for qTA6.1.

For qTA8.1, LOC_Os08g06090 was shown to be significantly and highly expressed in spread-tiller varieties than that in compact-tiller varieties as evidenced by qRT-PCR analysis (Fig. 5b). Additional haplotype analysis in CDS region showed that 6 genes (LOC_Os08g05840, LOC_Os08g05860, LOC_Os08g05870, LOC_Os08g05890, LOC_Os08g05900, and LOC_Os08g05910) can be divided into different haplotypes, all were associated with significant differences in the tiller angle (Fig. 5c). These results indicated that abovementioned seven genes are all potential candidate genes for qTA8.1.

Candidate genes analysis of the remaining 14 novel QTLs by integrating GWAS and transcriptomics

We further analyzed the 100 kb window near the leading SNPs within these remaining 14 QTLs. After removing 142 retro/transposon/hypothetical proteins, 297 annotated genes were identified (Additional file 2: Table S4).

To further narrow down candidate gene numbers, we referred to two publicly available transcriptomic data about tiller angle regulatory genes: one contained 4204 differentially expressed genes (DGEs) between gravistimulated shoots and control shoots without stimulation (Zhang et al. 2018); the other contained 4431 DGEs between ZnF transgenic plants and the wild type (Wu et al. 2018). By combining these DEGs and our GWAS results, we finally identified 15 overlapping candidate genes (Fig. 6a). Further qRT-PCR analysis showed that among these 15 candidate genes, expression levels of 3 genes (LOC_Os01g21590, LOC_Os03g15050 and LOC_Os04g41229) are significantly lower in spread-tiller varieties than in compact-tiller varieties (Fig. 6b and Additional file 2: Table S5), showing close association with tiller angle. These 3 candidate genes are likely to be a promising candidate gene in qTA1.2, qTA3.1 and qTA4.1, respectively, negatively regulating the tiller angle in rice.

Phenotype variations and GWAS

The erect plant architecture is favored by human beings because it improves photosynthetic efficiency and increases planting density, which eventually enhances grain yield (Wang and Li 2008, 2011). The genetic basis of rice tiller angle is complex and is associated with many loci, each with small contributions to phenotypic variance (Lu et al. 2015; Wu et al. 2019), our results confirmed these previous findings and showed rice tiller angle trait tends to show a nearly normal distribution in a segregated population (Fig. 1a). The natural variation of rice tiller angle revealed in this study, ranging from 6.3° to 35.6°(Fig. 1a-b), indicated not only that the panel selected in this study is rich in genetic diversity that suitable for GWAS analysis, but also that this panel is useful for direct use as parents in rice breeding for ideal tiller angle.

We used three models to perform GWAS in this study, all generated QTL signals, some of them were repeatedly detected in more than one model while most of them were not (Fig. 2). These results indicated that when performing GWAS, if possible, different models should be used in order to identify as much as possible the QTLs. Our results are consistent with previous report that different GWAS analysis methods possess unique features and complement each other (Liu et al. 2020). In total, three models yielded 18 QTLs for tiller angle; and the identification of two known rice tiller angle genes, PAY1 and TIG1, that co-localized with qTA8.3 and qTA8.4, respectively, indicated the reliability of our GWAS results.

Novel alleles in PAY1 and TIG1 affect tiller angle in rice

Some alleles of the genes related to yield-related QTLs were selected during domestication or breeding programs (Jin et al. 2008; Tan et al. 2008; Jiang et al. 2012). The allele mining technique has been successfully employed to find important variations at various locus for understanding mechanisms underlying domestication and for breeding (Huang et al. 2018; Guo et al. 2020; Bai et al. 2021). Present study identified novel alleles of PAY1 and TIG1 through GWAS in a japonica panel that are associated with variable tiller angle, and provided genetic resource for further understanding mechanisms underlying rice tiller angle. Clearly, the two alleles have different mechanisms regulating rice tiller angle. The SNP in coding region of PAY1 plays important role in regulating tiller angle (Fig. 3), while SNPs in the promoter region of TIG1 are associated with variation of tiller angle (Fig. 4). Therefore, learning from domestication, the promoter edited are interesting targets for genome editing to create new alleles for rice breeding as did previously (Swinnen et al. 2016; Oliva et al. 2019; Huang et al. 2020). Moreover, differences in allelic frequencies in the two sub-species (Fig. 3d) suggest that PAY1 has been subjected to selection during rice breeding, which partially explains why the Indica varieties usually have a relatively larger tiller angle and display comparatively spread-out plant architecture whereas japonica varieties always have a relatively smaller tiller angle and show a more compact plant architecture (Yu et al. 2007; Zhang et al. 2019). Similarly, the major tiller angle gene TIG1, is reported to be strongly selected in indica varieties by human beings (Zhang et al. 2019). Therefore, our study suggests the idea that rice tiller angle is selected during domestication.

Notably, the previous reported variants in both PAY1 and TIG1 were not found in 3K RGP (Zhao et al. 2015; Zhang et al. 2019), the detected SNPs, particularly those in the TIG1 promoter, in this study, thus, provide valuable information for breeding rice varieties, especially japonica rice.

Novel candidate genes

Previous studies have demonstrated that rice tiller angle is strongly related to plant gravitropic responses, which also can be effected by hormones, such as auxin and brassinosteroids (BRs) (Dong et al. 2016a; Liu et al. 2016; Waite and Dardick 2018; Zhang et al. 2018). Using expression and haplotype analysis combined with GWAS, we identified 11 candidate genes for rice tiller angle in this study (Fig. 5a-c and Fig. 6b). Among them, 4 candidate genes (LOC_Os08g05840, LOC_Os08g05870, LOC_Os08g05890 and LOC_Os01g21590) attracted specific attention based on their functional annotations. Because DNA Topoisomerase 1 (OsTOP1) functions as a role in gravitropism by linking transcriptional R-loops with auxin signaling (Shafiq et al. 2017), LOC_Os08g05840, encoding DNA topoisomerase 1, likely is involved in tiller angle regulation in rice. LOC_Os08g05840 may be regulation of tiller angle in rice. LOC_Os08g05870 and LOC_Os08g05890 encode an oxysterol-binding protein-related protein 6 and a MSP domain containing protein, respectively. both protein families have been documented to transport sterols and involved in sterol regulation (Reynolds and Hand 2009; Saravanan et al. 2009; Umate 2011), Considering that sterols also act as precursors for BR hormones, we propose that LOC_Os08g05870 and LOC_Os08g05890 are two important components of BR signaling associated with rice tiller angle. LOC_Os01g21590, encodes a homeodomain containing protein, whose homologues have been reported to be associated with tiller angle (Zhang et al. 2018; Hu et al. 2020). Considering its significant higher expression in compact-tiller varieties than in spread-tiller varieties (Fig. 6b), LOC_Os01g21590 is likely connected with tiller angle in rice. Nevertheless, these candidate genes provide useful genetic resources for expand our understanding regarding hormone mediated tiller angle regulation in rice.

In summary, this study provides important genetic information for better understanding the genetic diversity of tiller angle in rice and useful genetic resources for breeding rice varieties with ideal plant architecture.