Integrating a genome-wide association study with a large-scale transcriptome analysis to predict genetic regions influencing the glycaemic index and texture in rice

Abstract

Reliably generating rice varieties with low glycaemic index (GI) is an important nutritional intervention given the high rates of Type II diabetes incidences in Asia where rice is staple diet. We integrated a genome-wide association study (GWAS) with a transcriptome-wide association study (TWAS) to determine the genetic basis of the GI in rice. GWAS utilized 305 re-sequenced diverse indica panel comprising ~2.4 million single nucleotide polymorphisms (SNPs) enriched in genic regions. A novel association signal was detected at a synonymous SNP in exon 2 of LOC_Os05g03600 for intermediate-to-high GI phenotypic variation. Another major hotspot region was predicted for contributing intermediate-to-high GI variation, involves 26 genes on chromosome 6 (GI6.1). These set of genes included GBSSI, two hydrolase genes, genes involved in signalling and chromatin modification. The TWAS and methylome sequencing data revealed cis-acting functionally relevant genetic variants with differential methylation patterns in the hot spot GI6.1 region, narrowing the target to 13 genes. Conversely, the promoter region of GBSSI and its alternative splicing allele (G allele of Wx^a ) explained the intermediate-to-high GI variation. A SNP (C˃T) at exon-10 was also highlighted in the preceding analyses to influence final viscosity (FV), which is independent of amylose content/GI. The low GI line with GC haplotype confirmed soft texture, while other two low GI lines with GT haplotype were characterized as hard and cohesive. The low GI lines were further confirmed through clinical in vivo studies. Gene regulatory network analysis highlighted the role of the non-starch polysaccharide pathway in lowering GI.

Keywords: final viscosity; gene-regulatory network analysis; genome-wide association study; glycaemic index; transcriptome-wide association study; whole-genome bisulfite sequencing.

Figure 1

Correlation among in vitro glycemic index (GI), eating (amylose content; AC) and cooking (final viscosity; FV) qualities of 305 rice germplasm lines. The scatter plots showed (a) the weak correlation of AC with GI (Pearson's r ² = 0.2809), (b) practically no meaningful relationship between GI and FV and (c) correlation between AC and FV with r ² = 0.30. Genotypes that had low GI values (GI ≤ 55) were shown as red dots whereas lines with GI value between 56 and 60 (intermediate GI) were randomly selected and plotted as blue dots. All three data sets (GI, AC and FV) were transformed into a distribution with normally distributed residuals that had mean 0 and unit variance prior to genetic analysis without regard to their raw data distribution. (d and e) represent the plots for in vivo GI evaluation of two landraces GQ02522 and GQ02497 identified to be low GI based on cohort studies of at least 10 human subjects with reference to control (glucose solution, blue line). (f) Overlaid sensory and texture profiling of three low GI lines using 14 textural attributes representing the characteristic variations in the radar graph. Fourteen attributes in case of the three accessions were evaluated in comparison to ‘Dinorado’ a popular rice variety in the Philippines (depicted by grey areas, wideness of which corresponds to the differences between the profiles generated in the two sessions of sensory test).

Figure 2

Genome‐wide association study of glycemic index and transcription‐wide association study (TWAS) highlighting the expression level of GBSS I. (a) Manhattan plot shows significant association signals in chromosomes 5 and 6 indicated by the red line at –log(p) = 7.7 that represents the Bonferroni‐corrected P‐value significance threshold. The blue line at –log(p) = 4.8 is the suggestive line represented by the association P‐value that corresponds to the q‐value<0.5 (false discovery rate). A C→T synonymous SNP (1 525 361 bp) in exon 2 of LOC_Os05g03600 (encodes a signalling receptor kinase) is a significant SNP in chromosome 5. In chromosome 6, we found a 230 kb region (1.60–1.83 Mb) that we refer to in this manuscript as GI6.1. (b) Zoomed‐in GI6.1 region, where multiple significant SNPs have implicated within a total of 26 genes that have been validated in the gene‐set analysis, distributed in 19 different linkage disequilibrium (LD) blocks. (c) Manhattan plot of the transcription‐wide association study that highlighted the expression level of GBSS I. The red and blue lines were as described in (a) with the exception of the suggestive line being a little bit higher at 4.9 due mainly to the adjustment made for false positives. (d) A plot of TWAS with –log10 P‐values in the similar GI6.1 region, based on the expression QTL (eQTL) analysis, narrowing down to 13 candidates genes in the GI6.1 region.

Figure 3

Summary for genomic variation present within the gene and 1 kb promoter region of granule‐bound starch synthase I ( GBSS I). (a) The –log plot of the association P‐values for the glycemic index (GI), shown in the top graph. A red bar indicates that the negative effect of the allele on the trait, while a black bar indicates a positive effect. Bar thickness reflects the relative effect size of the respective allele, as beta value. The middle one is the –log plot of the association P‐values for final viscosity (FV). The third track is the –log plot of the P‐values of the cis‐eQTL analysis done on the expression of GBSS I using TWAS. The bottom track that shows SNP IDs with their alleles, the placement of the SNPs on the chromosome segment, and the haplotype blocks calculated using Gabriel's algorithm implemented in Haploview 4.2. The two linkage blocks within the gene were formed from the linkage break at the A→C SNP (located at chromosomal position 1 768 006 bp) in GBSS I. This SNP is situated two base pairs away from the 3′ end of exon 6. (b and c) explained the haplotypes lying within GBSS I and contributing influence the GI and FV values, respectively. The replacement of T with G allele of the splice junction SNP of intron1 predisposes with converting high to intermediate GI, whereas a conversion of T to C allele in the exon 10 SNP correlated to shifting FV from high to intermediate.

Figure 4

Presence of CpG methylation and their correlation with the GI and amylose content (AC) observed in ten resequenced germplasm lines. (a) Representation of the correlations between degree of CpG methylation present in the promoter region of genes underlying GI6.1 hotspot region with the GI and AC; evaluated in three main categories, fully‐unmethylated, partially (10%–90% methylation) and fully methylated (>90% methylated region) categories. Correlation coefficients ranged from −1 (green) to +1 (red). Horizontal red arrow depicts the candidates showing the significant correlations with the respective trait. (b) Graphs showing the degree of the CpG methylation (fully and partially methylation) patterns in the genic region of the key candidates LOC_Os06g04200 linked with AC, (c) LOC_Os06g04070, (d) LOC_Os06g04230 and (e) LOC_Os06g04240, linked with and GI values (c–e).

Figure 5

Weighted gene co‐expression analysis among Intermediate and low GI lines. (a) The heat map of differentially expressed genes between the intermediate versus low GI lines (Heatmap scale ranges from min −3 to max +3 normalized expression values which are shown as blue to red colour gradient representing the low and high expression, respectively). (b) Gene dendrogram with two clustered modules blue, turquoise in intermediate versus low; (c, d) Gene co‐expression subnetwork of Blue and Turquoise module, respectively in intermediate versus low. Nodes in each network represent the gene and edges as the interaction, the variation of node size show the different degree of connectivity and bordered wall nodes represents the hub. Nodes are colour coded based on MapMan functional categories.

Figure 6

Schematic overview summarizing the predicted regions influencing rice glycaemic index and texture. Overlaying genome‐wide association study with transcriptome‐wide association study results (shown at the bottom left) narrowed the candidate genes down to 13 in the GI6.1 region, and further validated using expression QTL (eQTL) analysis. Four important loci influencing GI trait within GI6.1 region show interesting methylation patterns. The diagnostic haplotypes identified in LOC_Os06g04070, LOC_Os06g04200, LOC_Os06g04230 and LOC_Os06g04240 are critical to lowering the glycaemic index.