Deep sequencing transcriptional fingerprinting of rice kernels for dissecting grain quality traits

Abstract

Background: Rice represents one the most important foods all over the world. In Europe, Italy is the first rice producer and Italian production is driven by tradition and quality. All main rice grain quality traits, like cooking properties, texture, gelatinization temperature, chalkiness and yield, are related to the content and composition of starch and seed-storage proteins in the endosperm and to grain shape. In addition, a number of nutraceutical compounds and allergens are known to have a significant effect on grain quality determination. To investigate the genetic bases underlying the qualitative differences that characterize traditional Italian rice cultivars, a comparative RNA-Seq-based transcriptomic analysis of developing caryopsis was conducted at 14 days after flowering on six popular Italian varieties (Carnaroli, Arborio, Balilla, Vialone Nano, Gigante Vercelli and Volano) phenotypically differing for qualitative grain-related traits.

Results: Co-regulation analyses of differentially expressed genes showing the same expression patterns in the six genotypes highlighted clusters of loci up or down-regulated in specific varieties, with respect to the others. Among them, we detected loci involved in cell wall biosynthesis, protein metabolism and redox homeostasis, classes of genes affecting in chalkiness determination. Moreover, loci encoding for seed-storage proteins, allergens or involved in the biosynthesis of specific nutraceutical compounds were also present and specifically regulated in the different clusters. A wider investigation of all the DEGs detected in pair-wise comparisons revealed transcriptional variation, among the six genotypes, for quality-related loci involved in starch biosynthesis (e.g. GBSSI, starch synthases and AGPase), genes encoding for transcription factors, additional seed storage proteins, allergens or belonging to additional nutraceutical compounds biosynthetic pathways and loci affecting grain size. Putative functional SNPs associated to amylose content in starch, gelatinization temperature and grain size were also identified.

Conclusions: The present work represents a more extended phenotypic characterization of a set of rice accessions that present a wider genetic variability than described nowadays in literature. The results provide the first transcriptional picture for several of the grain quality differences observed among the Italian rice varieties analyzed and reveal that each variety is characterized by the over-expression of a peculiar set of loci affecting grain appearance and quality. A list of candidates and SNPs affecting specific grain properties has been identified offering a starting point for further works aimed to characterize genes and molecular markers for breeding programs.

Fig. 1

Differentially modulated genes from pair-wise comparisons. a Mean expression versus log₂ Fold Change (FC) plots (MA-plots) representing the DEGs in the pair-wise comparisons, considering the expression in the second cv indicated with respect to the first one. Normalised expression mean values are plotted versus log₂FC and DEGs are visualized as coloured dots. b For each comparison the numbers of over-expressed and under-expressed genes in the second cv with respect to the first one indicated and the total DEGs are reported

Fig. 2

General representation of the clusters obtained by the co-regulation analysis. Numbers indicating each cluster are reported. Node colours corresponds to the mean expression values, calculated on the base of the number of reads mapping at each locus on rice reference genome, for all the varieties considering the three biological replicates and a scale from white, representing low expression, to purple, corresponding to high transcription

Fig. 3

Summary of the number of DEGs detected for each cluster. Bar plots represent the expression pattern of a representative gene of the cluster in all the six genotypes. Values in bar plots correspond to the mean expression values, corresponding to the number of reads mapping at each locus on rice reference genome (raw read counts as obtained from DESeq-normalized matrix of expression data), of the three biological replicates (y-axis) for each cv (ARB-Arborio, BAL-Balilla, CAR-Carnaroli, GV-Gigante Vercelli, VN-Vialone Nano and VOL-Volano; x-axis). The main enriched GO terms and functional classes identified for each cluster are reported (the numbers of DEGs belonging to each class are in brackets)

Fig. 4

Bar-plots of some tightly co-expressed loci in cluster 2. a, b and c Different sub-clusters of tightly co-regulated loci belonging to cluster 2. Mean expression values, corresponding to the number of reads mapping at each locus on rice reference genome (raw read counts as obtained from DESeq-normalized matrix of expression data), of the three biological replicates (y-axis) for each cv (ARB – Arborio, BAL – Balilla, CAR – Carnaroli, GV – Gigante Vercelli, VN – Vialone Nano and VOL – Volano; x-axis) are represented

Fig. 5

Expression of genes belonging to the Wx cluster. a Representation of the Wx cluster and Bar-plots indicating the mean expression values, corresponding to the number of reads mapping at each locus on rice reference genome (raw read counts as obtained from DESeq-normalized matrix of expression data), of the three biological replicates (y-axis) for each cv (ARB–Arborio, BAL–Balilla, CAR–Carnaroli, GV–Gigante Vercelli, VN–Vialone Nano and VOL–Volano; x-axis). b Sashimi plots of RNA-Seq reads aligned at LOC_Os06g04200, on Nipponbare Reference Genome, in the six different varieties. Only one biological replicate for each genotype is represented. The alignments were performed using the IGV software. The coverage for each alignment track is plotted as a bar graph. Coverage ranges for each variety are reported in the bars on the left. Arcs represent splice junctions. Junction depth is indicated by the number in the arc and the thickness of the arc. Graphs were built setting a Min Junction Coverage corresponding to the 1 % of the mean exons coverage. The alternative isoforms of the gene annotated in the Rice Genome browser MSU database are reported in the lower bar. Thick and intermediate boxes represent exons and 5’ and 3’-UTR regions, respectively. Lines between boxes indicate introns. The mapping coverage for each variety at the level of each exon and intron is represented by the coloured parts

Fig. 6

IGV outputs of the alignments of the RNA-Seq reads for each cv at the LOC_Os06g12450 and LOC_Os11g14220 loci on Nipponbare Reference Genome (left part). Only one biological replicate for each genotype is represented. The structures of the genes annotated in the Rice Genome browser MSU database are reported in the lower bar. Thick and intermediate boxes represent exons and 5’ and 3’-UTR regions, respectively. Lines between boxes indicate introns. The mapping coverage for each variety at the level of each exon and intron is represented by the upper grey parts. Forwards reads are in red, reverse reads are in blue. Enlargement of the sections indicated by black squares are represented on the right part of the Figure. Putative functional SNPs are indicated as coloured bars: orange for G to A and red for G to T and C to T

Fig. 7

qRT-PCR at 6, 10 and 14 DAF for the 14 selected DEGs in all the six varieties. Values are expressed as log₂FC of expression at 10 and 14 DAF in comparison to 6 DAF, respectively. FCs are expressed as mean of the three biological replicates. Bars indicate standard errors across for the three biological replicates considering the three technical replicates for each of them