Next Article in Journal
Neuroplasticity and Multilevel System of Connections Determine the Integrative Role of Nucleus Accumbens in the Brain Reward System
Next Article in Special Issue
The Carboxyl-Terminus of TRANSPARENT TESTA GLABRA1 Is Critical for Its Functions in Arabidopsis
Previous Article in Journal
Optimised Heterologous Expression and Functional Analysis of the Yersinia pestis F1-Capsular Antigen Regulator Caf1R
Previous Article in Special Issue
3,4-Dehydro-L-proline Induces Programmed Cell Death in the Roots of Brachypodium distachyon
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 22
Issue 18
10.3390/ijms22189787
ijms-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

RNA-Seq-Based Profiling of pl Mutant Reveals Transcriptional Regulation of Anthocyanin Biosynthesis in Rice (Oryza sativa L.)

by
Ruonan Xu
1,†,
Ronghui Pan
1,†,
Yuchan Zhang
1,
Yanlei Feng
2,3,
Ujjal Kumar Nath
4,
Yinbo Gan
1,
Chunhai Shi
1,* and
Delara Akhter
1,5,*
1
Department of Agronomy, Zhejiang University, Hangzhou 310058, China
2
School of Life Sciences, Westlake University, Hangzhou 310024, China
3
Westlake Institute for Advanced Study, Hangzhou 310024, China
4
Department of Genetics and Plant Breeding, Bangladesh Agricultural University, Mymensingh 2202, Bangladesh
5
Department of Genetics and Plant Breeding, Sylhet Agricultural University, Sylhet 3100, Bangladesh
*
Authors to whom correspondence should be addressed.
†
These authors have contributed equally to this work.
Int. J. Mol. Sci. 2021, 22(18), 9787; https://doi.org/10.3390/ijms22189787
Submission received: 29 July 2021 / Revised: 6 September 2021 / Accepted: 8 September 2021 / Published: 10 September 2021
(This article belongs to the Special Issue Plant Cell and Organism Development 2.0)
Browse Figures
Versions Notes

Abstract

:
Purple-colored leaves in plants attain much interest for their important biological functions and could be a potential source of phenotypic marker in selecting individuals in breeding. The transcriptional profiling helps to precisely identify mechanisms of leaf pigmentation in crop plants. In this study, two genetically unlike rice genotypes, the mutant purple leaf (pl) and wild (WT) were selected for RNA-sequencing and identifying the differentially expressed genes (DEGs) that are regulating purple leaf color. In total, 609 DEGs were identified, of which 513 and 96 genes were up- and down-regulated, respectively. The identified DEGs are categorized into metabolic process, carboxylic acid biosynthesis, phenylpropanoids, and phenylpropanoid biosynthesis process enrichment by GO analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) confirmed their association with phenylpropanoid synthesis, flavonoid synthesis, and phenylalanine metabolism. To explore molecular mechanism of purple leaf color, a set of anthocyanin biosynthetic and regulatory gene expression patterns were checked by qPCR. We found that OsPAL (Os02g0626100, Os02g0626400, Os04g0518400, Os05g0427400 and Os02g0627100), OsF3H (Os03g0122300), OsC4HL (Os05g0320700), and Os4CL5 (Os08g0448000) are associated with anthocyanin biosynthesis, and they were up-regulated in pl leaves. Two members of regulatory MYB genes (OsMYB55; Os05g0553400 and Os08g0428200), two bHLH genes (Os01g0196300 and Os04g0300600), and two WD40 genes (Os11g0132700 and Os11g0610700) also showed up-regulation in pl mutant. These genes might have significant and vital roles in pl leaf coloration and could provide reference materials for further experimentation to confirm the molecular mechanisms of anthocyanin biosynthesis in rice.

1. Introduction

The leaf color of the plant is considered to be one of the important morphological traits that might be the potential source of a reliable phenotypic marker in the breeding program [1]. Morphological traits provide a significant contribution to genetic advancement in rice breeding; thus, breeders are looking for the distinct phenotypes that could be used as a marker for effective selection [2,3]. Mechanisms of purple leaf color in plants are quite complex; factors such as the type of pigment as well as their concentration and distribution in leaf cells contribute to leaf colors. Three types of pigment are found in the leaves of the higher plants: (i) chlorophyll (Chl) (mainly Chla and Chlb), (ii) carotenoids (lutein and carotenoid), and (iii) flavonoids (anthocyanins) [4]. Anthocyanins are water-soluble colors and extensively disseminated among the higher plants. They contribute to physiological functions of plants and act as attractants of pollinators and seed dispersers. Recent research has confirmed that anthocyanins have protective functions against different biotic and abiotic stresses [5,6,7,8,9]. It produces coloration of the leaves and leaf sheaths in rice [10]. However, to date, detailed research on the genes involved in anthocyanin biosynthesis in rice has not been carried out.
The anthocyanin biosynthesis pathway starts from phenylalanine under synchronized regulation of structural genes in the cytoplasm, which encode some regulatory transcription factors and enzymes [11]. Glutathione S-transferases export the anthocyanins from cytoplasm to the vacuole, and they are stored there permanently [12]. The anthocyanin biosynthesis process involves three consecutive stages: firstly, phenylalanine is transformed into transcinnamic acid and 4CL (4-coumarate-CoA ligase), respectively [4,8,13]; thereafter, CHS (chalcone synthase), CHI (chalcone isomerase), and F3H (flavanone-3-hydroxylase) are catalyzed a series of enzymatic reactions and formed dihydroflavonol from one molecule of coumarate-CoA and three molecules of malonyl-CoA; finally, anthocyanidins are produced from dihydroflavonols via catalyzing of DFR (dihydroflavonol 4-reduatse) and ANS (anthocyanidin synthase). These anthocyanidins are modified for glycosylation and methylation to produce anthocyanins through the action of UFGT (UDP glucose: flavonoid glucosyltransferase) and MT (methyl transferase) [14,15]. Earlier studies proved that the expression of the structural genes is regulated by synchronized interactions of R2R3-MYB (myeloblastosis family), WD40-repeat transcription factors, and bHLH (basic helix–loop–helix) [16]. MBW complex (R2R3-MYB, WD40, and bHLH) is very important to plants; its transcriptional regulatory system might be different between mono- and dicotyledonous plants, and all members of this MBW complex are not obligatory for anthocyanin biosynthesis in all species [14,15].
Transcriptome analysis is a dynamic way to detect leaf color genes in any plant species. Therefore, RNA-Seq (RNA sequencing) is used for profiling the transcripts with gene annotation and identification in different plant species [17,18,19]. In this study, RNA-seq was performed by using the leaf of pl mutant and its wild type (WT; Zhenong 34) at the grain filling stage. This pl mutant was developed by EMS (ethyl methane sulfonate) mutagenesis. It was phenotypically characterized in M2 population [20]. We executed transcriptomic analysis of pl vs. WT to identify the DEGs for explaining the mechanisms involved in leaf coloration in rice.

2. Results

2.1. Phenotypic Characterization

The leaf color of pl mutant was found to be purple at the grain filling stage, which was not detected in the seedling to tillering stages. The pl plants attained a purple color at the grain filling stage (Figure 1A–C), whereas WT remained green (Figure 1D,E). In pl mutant, purple coloration started firstly in the leaf blade and midrib, then gradually covered the entire leaf. Additionally, pl plants produced shorter grains with fewer numbers of grains/panicle and a lower weight per 1000 grains compared to WT.

2.2. RNA-Sequence Analysis

RNA-seq was performed using leaf samples of pl mutant and WT plants for profiling the transcriptome. Illumina technology was used to sequence the cDNA of pl mutant and WT to identify the possible DEGs and related pathways that are involved in purple leaf coloration. A total of 54,237,834 bp and 55,160,750 bp raw reads were extracted from pl and WT samples, respectively (Table 1). After eliminating the adaptor and poor-quality sequences, in total, 52,941,300 bp and 53,875,169 bp clean reads were retrieved from pl and WT samples, respectively, with an average length of 142.99 and 144.14 bp. The average GC content varied from 57.71% to 57.83% in different libraries with >97% Q30 (Table 1).
In total, 106,816,469 clean reads were attained from six cDNA libraries, and they were categorized into three different classes: uniquely mapped, multiple mapped, and total mapped. From the total clean read, 84.95–86.70 % was totally mapped, 77.99–82.21% was uniquely mapped, and 2–6% was multiple mapped (Table 2).
In the heat map of correlation coefficients of the samples of pl mutant and WT (Figure 2), pl samples had 0.98 (pl1, pl2, and pl3) average correlation coefficient while WT (WT1, WT2, and WT3) had 0.96, indicating their vicinity toward the unity.
We found 59,358 and 58,361 SNPs as well as 5,648 and 5,297 indels in pl and WT, respectively. SNP regions were separated as inter-genic and genic regions (UTR and CDS), where most of the SNPs were found in the genic regions of the pl and WT samples (Table S1). Two types of nucleotide substitutions—transversion (A/T, A/C, G/A C/G and G/T) and transition (A/G, G/A, T/C and C/T)—were found in homozygote SNPs. We found more transition than transversion substitution of SNPs in both pl and WT sequences (Table S2).

2.3. DEGs Identification

The DEGs pl vs. WT were identified and it was predicted that the genes are involved in the purple leaf color of pl mutant (Table S3). The Volcano plot shows the distribution of DEGs at the threshold level for selection. Among 609 DEGs, 513 and 96 were up- and down-regulated in pl vs. WT, respectively (Figure 3A, Table S4). In the Venn diagram, a total of 261 genes were identified as commonly shared for all samples (Figure 3B).

2.4. GO Analysis

GO (gene ontology) is a global gene function classification that categorizes genes based on their biological processes, cellular components, and molecular functions. Among the DEGs, 40.35% were classified based on their biological processes, 29.82% based on their molecular functions, and 26.32% based on their cellular functions (Figure 4 and Table S5). The majority of the genes in the biological process category were sub-categorized into metabolic process; the rest of the genes were included into the cellular process and single-organism subcategories. A significant proportion of the genes was classified based on catalytic and binding activity. The majority of this portion of genes was included in the cell and cell part followed by membrane, membrane part, and organelle sub-categories (Figure 4). The up- and down-regulated DEGs were categorized according to enriched GO terms by bar plot. A total of 23,373 genes were allotted in GO terms and only 460 (1.96%) DEGs were confirmed as annotated in the database, of which 407 genes were up- and 53 were down-regulated. The scattered plot was constructed by using up- and down-regulated DEGs, which exhibited that most of the genes are involved in metabolic process and rest of them are involved in carboxylic acid and phenylpropanoid biosynthesis processes (Figure S1 and Table S6).

2.5. Analysis of DEGs via KEGG and KOG Pathway

The KEGG (Kyoto Encyclopedia of Genes and Genomes) database was implemented to explore the pathways involved in leaf coloration of the pl mutant at the late grain filling stage. Bar-plot analysis in the KEGG database was used to identify the uni-gene positioning of metabolic pathways. Only 96 out of 3,335 assigned genes were found to be involved in four pathways: cellular processes, environmental, genetic information processing, and metabolism (Figure 5A and Table S7).
The scattered plot analysis using the top 30 KEGGs showed that the DEGs designated as ko00940 (phenylpropanoid biosynthesis), ko00360 (phenylalanine metabolism), ko00941 (flavonoid biosynthesis), and ko00400 (phenylalanine, tyrosine, and tryptophan biosynthesis) were involved in the leaf coloration-related pathway. The uni-gene positioning was explored through the bar-plot analysis, and it was found that only 96 genes were involved in metabolic, cellular, genetic, and environmental information processing pathways (Figure 5B and Table S8).
A gene interaction network diagram was generated to identify the involvement of DEGs and their contribution to specific biological phenomena. Most of the genes were found in phenylalanine metabolism, phenylpropanoid biosynthesis, flavonoid biosynthesis, and stilbenoid, diarylheptanoid with gingerol biosynthesis pathways (Figure 6 and Table S8).
Leaf coloration of pl mutant might be predicted through the analysis of enrichment. The function of the genes was grouped using the information available in the COG/KOG (Clusters of Orthologous Groups) database based on their orthologs. A total of 13,979 uni-genes were identified as non-redundant in the KOG; among them, 296 genes were differentially expressed. We classified all of the identified uni- and special DEGs into 26 sub-functional groups (Figure 7, Table S9). Maximum genes were found to involve in signal transduction system; thereafter, the rest of the DEGs were found in protein transporter, general post-translational modification, chaperon, secondary metabolite biosynthesis, transport, catabolism, and general function prediction and cell motility.

2.6. Validation of the RNA-Seq Data by qRT-PCR

RNA-seq data was validated through qRT-PCR (quantitative real-time polymerase chain reaction) considering 20 randomly selected genes related to purple color in plants (Table 3). Expression pattern of the selected genes showed consistency in sense of up- and down-regulation with RNA-seq data, indicating reliability of the transcriptomes (Table 3). The pl mutant had colored leaves with four- and eight-times higher anthocyanin content in the leaf and leaf sheath, respectively, than WT (Table S10). In pl mutant, we identified the DEGs related to anthocyanin pathway, which are enriched with MYB, bHLH, and WD40 transcription factors. We also identified 17 DEGs related to the phenylpropanoid biosynthesis rice gene’s ortholog. Among them, six DEGs are related to flavonoid biosynthesis and seven DEGs to the UDP transferase family were up-regulated significantly in pl mutant compared to WT (Table 3), which supported our RNA-seq data.

2.7. Involvement of Transcription Factors and Anthocyanin Biosynthesis in pl Mutant

The transcription factor families—MYB, WD40, and bHLH—were found to be involved in the anthocyanin biosynthetic process in pl mutant. Four DEGs of each of the MYB, bHLH, and WD40 family were identified in RNA-seq data (Table 3 and Figure 8). Candidate transcription factors were screened out through quantification of relative expression by qRT-PCR. Among them, two MYB genes (Os05g0553400 (OsMYB55) and Os08g0428200), two bHLH family genes (Os01g0196300 and Os04g0300600), and two WD40 family genes (Os11g0132700 and Os11g0610700) showed up-regulation in pl mutant. These genes might have vital roles in establishing the purple leaf color in pl mutant and could provide information on anthocyanin biosynthesis mechanisms in rice.

2.8. Expression Analysis of Anthocyanin Biosynthesis Genes in pl Mutant

We found some structural and early biosynthesis genes; OsPAL (Os02g0626100, Os02g0626400, Os04g0518400, Os05g0427400 and Os02g0627100), OsF3H (Os03g0122300), OsC4HL (Os05g0320700) and Os4CL5 (Os08g0448000) were highly upregulated in pl mutant than WT (Figure 9). Out of two MYB genes, Os05g0553400 (OsMYB55; LOC_Os05g48010) was predicted as a candidate in the long arm of chromosome 5 through bulk segregation analysis (BSA) in our previous study [20]. This predicted gene was also identified as the potential candidate by RNA-seq and showed differential expression both in RNA-seq data and qPCR analysis with a 2.5-fold higher expression in pl mutant compared with WT (Figure 8). This higher relative expression levels of genes together with regulatory transcription factor genes might contribute to higher anthocyanin content in pl plants.

3. Discussion

Currently, purple color in plants has achieved much attention because of its advantageous biological functions. The distribution and amount of three types of colors (chlorophylls, carotenoids, and anthocyanins/flavonoids) are responsible for various pigments that are found in different plant tissues. Chlorophyll is higher than other types of pigments in green plants, whereas in yellow-leaved plants, carotenoid is higher; in purple-, red-, and blue-leaved plants, anthocyanin is higher than the other two pigments [21]. However, this insight is obviously lacking and still not well understood. For a better understanding of the underlying reason of leaf coloring, researchers are searching for mutants with changed phenotypes. A new rice mutant, p1, with a purple leaf color at the grain filling stage was identified [20], which confirmed that a purple leaf color in pl mutant was attributed due to a single recessive gene. The p1 plants showed abnormal and disintegrated chloroplast cells with decreased chlorophyll content [20]. In this study, RNA-seq data showed variations in the transcriptome of p1 mutant from the WT of rice plants at the grain filling stage. RNA-seq helps in the detailed transcriptomic profiling, annotation, and identification of genes [17,22,23,24].
The correlation coefficient of gene expression values among the tested samples reflects the consistency of the RNA-Seq data of the samples. The correlation coefficient among the p1 samples was 0.98, whereas it was 0.96 in WT samples, indicating the consistency of the RNA-seq data among the samples because the values were in the vicinity of unity [18].
Confirmation of the expression of predicted DEGs through qRT-PCR is necessary for further functional analysis of the transcriptomes [18,25]. The qRT-PCR was conducted using 20 randomly selected genes, which were predicted as DEGs. GO analysis conferred that up-regulated DEGs are involved in phenylpropanoid biosynthesis process (Figure S1). Moreover, KEGG and pathway enrichment analyses showed the uni-genes involved in phenylpropanoid, phenylalanine metabolism, and flavonoid biosynthesis pathways as well as leaf coloration [8,21,26].
Anthocyanins are a group of flavonoids and are synthesized in the phenylpropanoid pathway, which has been studied thoroughly [26]. At the beginning of the flavonoid pathway, 4-coumaroyl CoA is catalyzed to naringenin chalcone to naringenin and to dihydroflavonol with the catalyzation by chalcone synthase (CHS), chalcone isomerase (CHI), and flavanone 3-hydroxylase (F3H), respectively. Thereafter, dihydrokaempferol (DHK) is hydroxylated by F3′H to produce dihydroquercetin, which is converted to leucocyanidin and then to cyanidin by DFR and ANS activities, respectively [9,27]. Finally, decorative anthocyanins are produced and transported to vacuole due to activities of transcription factors and transporter genes. The anthocyanin pathway has been genetically and biochemically explained in Arabidopsis. Due to knockout mutations of biosynthetic genes, less pro-anthocyanidin and anthocyanin is produced in Arabidopsis seeds [4,28,29]. Four genes encode the isomers of PAL, of which PAL1 and PAL2 are the preferred flavonoid pathways at nitrogen deficit and cold temperature conditions [28,30,31,32]. In our study, the DEGs of phenylpropanoid and flavonoid biosynthesis pathways were up-regulated, for example: naregenin-chalcone synthase (CHS), chalcone isomerase (CHI), flavonol synthase (FAS) and UDP-glucosyl transferase. The orthologous of PAL family genes: OsPAL1 (Os02g0626100), OsPAL2 (Os02g0626400), OsPAL4 (Os02g0627100), OsPAL6 (Os04g0518400), and OsPAL7 (Os05g0427400) were highly up-regulated in pl samples. In addition, C4H (Os05g0320700), 4CL genes (Os02g0177600 and Os08g0448000), OsF3H (Os03g0122300), anthocyanidin 5,3-O-glucosyltransferase (Os01g0734800) were also up-regulated in pl samples (Table S3).
Anthocyanins are synthesized due to the up-regulation of structural and transcription factor genes regulated by different regulatory genes. The downstream genes of anthocyanin biosynthesis are regulated by WBM complex (WD40-bHLH-MYB). The R2R3-MYB proteins contain two MYB domain, such as MYB75/PAP1, MYB90/PAP2, MYB113 and MYB114, and are regulated by anthocyanin biosynthesis in vegetative tissues [28,33,34]. In plants, different tissue attains various levels of anthocyanin due to the degree of expression of R2R3-MYB transcription factors of MBW complex [35,36,37,38], consistent with our findings. Anthocyanin accumulation in apple and purple-head Chinese cabbage are regulated by MdMYB1 and BrMYB2, respectively [39], and also the up-regulation of structural genes BrF3’H, BrDFR1, BrANS1, and BrUGTs [40]. In pear (Pyrus pyrifolia) fruit, light-induced anthocyanins are produced through the binding of PpMYB10 to the promoters of CHS and CHI genes [41]. The roles of bHLH transcription factors and WD proteins are pleiotropic [42]. Simultaneous up-regulation of MYB and bHLH with structural genes in Arabidopsis, cauliflower, and petunia is essential for anthocyanin biosynthesis [9,35,42,43].
In red-skinned pear (Pyrus comunis L.), PybHLH3 and PyMYB114 transcription interact together and promote anthocyanin accumulation [44]. Anthocyanin biosynthesis in apple is modulated by the interaction of MdMYB308L and MdbHLH33 through binding to MdDFR promoter [45].
Two MYB genes ((Os05g0553400 (OsMYB55) and Os08g0428200)), two bHLH family genes (Os01g0196300 and Os04g0300600), and two WD40 family genes (Os11g0132700 and Os11g0610700) showed up-regulation in pl mutant. In our previous study, we reported that OsMYB55 (LOC_Os05g48010) with the ID Os05g0553400 in RAP_locus might be the reason for accumulating anthocyanin in the p1 mutant [20], consistent with other studies [35,46]. Our studies corroborated similar results of RNA-seq and qPCR expression analyses or BSA techniques for predicting candidate gene and identification. We narrate that the MYB, bHLH, and WD40 transcriptional complexes regulate phenylpropanoid metabolism through binding with flavonoid biosynthesis genes for enhancing anthocyanin accumulation (Figure 10). We speculate that R2R3-MYB transcription factors might be responsible for activating anthocyanin-biosynthesis genes, whereas, bHLH activates the transcription of MYB and consequently promotes anthocyanin accumulation.
Our RNA-seq data and the relative expression of the transcription factors showed up-regulation of two genes from each of the MYB, bHLH, and WD40 families are predicted to involve in anthocyanin accumulation through coordinated transcriptional activation of anthocyanin biosynthesis genes in pl mutant. Up-regulation of these three transcription factors might be influenced by the activity of the OsMYB55 (R2R3-MYB transcription factor) coordinately in p1 mutant for anthocyanin accumulation and purple coloration.
We profiled transcriptomes of pl vs. WT plants to figure out the transcripts involved in leaf coloration of rice. Through GO enrichment, KEGG pathway and KOG/COG analyses, we found that multiple pathways are intermingled in leaf coloration in rice. Based on the literature, we presumed that the phenylpropanoid pathway might contribute more than any other pathways to leaf coloration of rice. From our transcript data, we found that up-regulation of the members of transcriptional regulation complex (MYB, bHLH, and WD40) coordinately regulates anthocyanin accumulation. Among this complex, the R2R3-MYB transcription factors activate the biosynthetic genes and take part in leaf coloration. These results could give an advanced understanding of the leaf coloration process in rice and might also provide the experimental insight for manipulating purple color traits of rice through breeding programs.

4. Materials and Experimental Methods

4.1. Plant Materials and Sample Collection

Flag leaves of cultivar Zhenong 34 and the purple leaf (pl) mutant of Oryza sativa L. ssp. indica were used for RNA sequencing. The pl mutant was induced through ethyl methane sulfonate (EMS) mutagenesis and it was identified phenotypically with purple coloration of leaf and leaf sheath at the grain filling stage in M2 population [20]. Ten flag leaves from each of the pl and WT plants at the grain filling stage were collected and pooled together. The pooled samples of pl and WT were immersed immediately in liquid nitrogen and stored in a freezer at −80 °C until RNA extraction.

4.2. Isolation of RNA, Library Construction and Sequencing

Total RNA was extracted using a total RNA extractor kit (Trizol) (B511311, Sangon, Shanghai, China) following the manufacturer’s instructions. The RNA quality and quantity were measured using agilent 2100 Bioanalyzer (Agilent Technologies, Blbd Santa Clara, CA, USA) and Nanodrop (IMPLEN, Westlake Village, CA, USA), respectively. High-grade RNA samples were sent to Sangon Biotech Co. Ltd. (Shanghai, China) for sequencing and further analyses. The mRNAs from total RNA were refined using poly-T oligo-attached magnetic bead. The first-strand cDNA was synthesized by using random hexamer primer and M-MuLV reverse transcriptase (RNase H-) and second-strand cDNA was synthesized by using DNA polymerase I and RNase H. Exonuclease/polymerase was used to change the remaining overhangs into blunted ends and the adaptors were ligated for preparing the libraries after adenylation at 3′ends of DNA fragments. Forward primer: 5′ AGATCGGAAGAGCACACGTCTGAAC 3’ and reverse primer: 5′ AGATCGGAAGAGCGTCGTGTAGGGA 3′ were used in preparing the libraries. A total of six libraries—three from each of the pl mutant and WT—were produced and Illumina HiSeqTM2000 (Illumina, San Diego, CA, USA) was used for paired-end sequencing of the constructed libraries.

4.3. RNA-Seq Data Assay

Raw reads gathered from Illumina HiSeqTM 2000 sequencing were trimmed to retrieve the clear reads. Adapter, poly-N, and poor-quality reads were eliminated during purification of the data to improve the accuracy of the analysis. The raw reads were cleaned by FASTQ [18,47]. The cleaned and purified reads were kept in FASTQ format and were mapped to the reference gene using hierarchical indexing for spliced alignment of transcripts (Hisat2; version 2.1.0) with default parameters [48]. The raw sequences are uploaded in the NCBI under accession no. prjna760250 (https://www.ncbi.nlm.nih.gov/bioproject/?term=prjna760250/, accessed on 10 September 2021). Nipponbare genome (http://rice.plantbiology.msu.edu/, accessed on 12 October 2020) was used as reference to identify the perfect match and mismatch [49]. For defining the transcript isoforms of the similar gene, RNA-seq by expectation maximization (RSEM) was employed [50]. Whole RNA-seq was conducted using the pipeline described in Figure S2.

4.4. Transcriptomic Analysis

DESeq2 (version 1.12.4) was utilized for identifying the DEGs from the samples [51]. Genes were considered as differentially expressed, where q < 0.001 and log2|fold change| > 2. Here, q is the quality factor used to control the false-positive rate whilst carrying out multiple tests, and the fold change showed the multiple of expression difference. Strict filtering was applied via independent filtering for removal of the expression values which were equal to zero or NA on the mean of normalized counts and plotted on log plot. The transcript per million (TPM) was used to normalize the reads and was expressed as fragments per kilo base per million (FPKM).
DEGs were assigned to GO terms (biological functions) in the database using ClusterProfiler [52] and a hypergeometric assessment was carried out for identifying GO enriched for the genes from the reference genome. KEGG pathway and the DEGs were compared to reference genome for identifying significantly enriched metabolic pathways.

4.5. qRT PCR Analysis

For validation of RNA-sequencing data, a total of 20 genes were selected randomly from DEGs and determined the expression level using qRT PCR (quantitative real-time PCR). Primer 5 software version 5 was used for designing the gene-specific primers (Table S11). The qRT-PCR was conducted with 10 μL reaction mixture using Ex TaqII (Takara, Tokyo, Japan) in a Roche light cycler 480 (Roche, Stuttgart, Germany) and the reaction was repeated five times as replication. The reaction mixture contained 1 µL template cDNA (50 ng), 2 µL ddH2O, 5 µL 2× SYBR mix, and 1 µL (10 pmol) of each forward (F) and reverse (R) primers. The qRT-PCR was performed using the thermal protocol; denaturation 95 °C for 30 s; thereafter, 40 cycles of denaturation; 95 °C for 5 s, annealing; 55 °C for 20 s and extension; 72 °C for 10 s. The OsActin gene was used as an internal control for calculating the relative expression [53]. For calculating the relative gene expression, cq values were normalized through cq of OsActin following the 2−ΔΔCT method [54]. A significance test for mean differences of gene expression was performed for 5 repeats between pl and WT plants per gene following a Tukey test (p ≤ 0.05) using Minitab 18 (State College, Borough, PA, USA).

4.6. Anthocyanin Content Measurement

Total anthocyanin in leaf and leaf sheath of pl and WT was estimated by using spectro-phtotometer UV-1800 (Shanghai, China) following the protocol [8]. Frozen leaf samples were grinded in a mortar and pestle placed inside liquid nitrogen and anthocyanin was extracted using acidic alcohol (1% HCl and 95% ethanol). Total anthocyanin was measured as: anthocyanins = (A530 − 0.25 × A657) × FW − 1; here, A530 = absorbance at 530 nm, A657 = absorbance at 657 nm, and FW = weight of fresh leaf (g). Each sample was repeated 3 times in quantification.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/ijms22189787/s1.

Author Contributions

D.A., R.P. and C.S. designed and developed the experiments. D.A. and R.X. performed the main experiments and drafted the manuscript. C.S. identified and provided the populations of the pl mutant for genetic study and gene mapping. Y.Z., Y.F. and R.P. helped to carry out the experiments, and collect and analyze some data. U.K.N. analyzed the data and edited the manuscript. Y.G. and R.P. helped to revise the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This study was supported by grants from the Leading Innovative and Entrepreneur Team Introduction Program of Zhejiang (2019R01002), Zhejiang Provincial Natural Science Foundation of China (LQ19C130004), Zhejiang University Student Research Practice Program (P2020052) and a grant from the Science and Technology Office of Zhejiang Province (2012C12901-2, 2016C32G2010016 and 2016C02050-6).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data will be available in the NCBI accession PRJNA760250 along with this manuscript upon publication.

Acknowledgments

The authors acknowledge Ran Qin, in the Department of Agronomy, Zhejiang University, Hangzhou, China for her technical support.

Conflicts of Interest

The authors declared no conflict of interest.

References

  1. Liu, X.; Sun, X.; Wang, W.; Ding, H.; Liu, W.; Li, G.; Jiang, M.; Zhu, C.; Yao, F. Fine mapping of Pa-6 gene for purple apiculus in rice. J. Plant Biol. 2012, 55, 218–225. [Google Scholar] [CrossRef]
  2. Xudong, C.L.Q.Q.Z.; Zhenmin, Z.D.M.S.X. Breeding of a Photo-thermo Sensitive Genie Male Sterile Indica Rice Zhongzi S with a Purple-leaf Marker and the Hoterosis of Its Hybrid Rice Produced with It. Acta Agron. Sin. 1999, 25, 44–49. [Google Scholar]
  3. Zhentao, T.; Qingzhong, X. Observation on Flowering Habit of Photo—Thermoperiod Sensitive Genic Sterile Rice Lines with a Pale—Green Leaf Marker. Fen Zi Zhi Wu Yu Zhong 2004, 2, 241–245. [Google Scholar]
  4. Winkel-Shirley, B. Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology. Plant Physiol. 2001, 126, 485–493. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Bradshaw, H.D.; Schemske, D.W. Allele substitution at a flower colour locus produces a pollinator shift in monkeyflowers. Nature 2003, 426, 176–178. [Google Scholar] [CrossRef]
  6. Castellarin, S.D.; Pfeiffer, A.; Sivilotti, P.; Degan, M.; Peterlunger, E.; Di Gaspero, G. Transcriptional regulation of anthocyanin biosynthesis in ripening fruits of grapevine under seasonal water deficit. Plant Cell Environ. 2007, 30, 1381–1399. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Lev-Yadun, S.; Gould, K.S. Role of Anthocyanins in Plant Defence. In Anthocyanins; Springer: Berlin/Heidelberg, Germany, 2008; pp. 22–28. [Google Scholar] [CrossRef]
  8. Goswami, G.; Nath, U.K.; Park, J.I.; Hossain, M.R.; Biswas, M.K.; Kim, H.T.; Kim, H.R.; Nou, I.S. Transcriptional regulation of anthocyanin biosynthesis in a high-anthocyanin resynthesized Brassica napus cultivar. J. Biol. Res. 2018, 25, 1–15. [Google Scholar] [CrossRef] [PubMed]
  9. Zhang, Y.; Hu, Z.; Zhu, M.; Zhu, Z.; Wang, Z.; Tian, S.; Chen, G. Anthocyanin accumulation and molecular analysis of correlated genes in purple kohlrabi (Brassica oleracea var. gongylodes L.). J. Agric. Food Chem. 2015, 63, 4160–4169. [Google Scholar] [CrossRef]
  10. Reddy, V.S.; Dash, S.; Reddy, A.R. Anthocyanin pathway in rice (Oryza sativa L.): Identification of a mutant showing dominant inhibition of anthocyanins in leaf and accumulation of proanthocyanidins in pericarp. Theor. Appl. Genet. 1995, 91, 301–312. [Google Scholar] [CrossRef]
  11. Petrussa, E.; Braidot, E.; Zancani, M.; Peresson, C.; Bertolini, A.; Patui, S.; Vianello, A. Plant flavonoids-biosynthesis, transport and involvement in stress responses. Int. J. Mol. Sci. 2013, 14, 14950–14973. [Google Scholar] [CrossRef] [PubMed]
  12. Dixon, D.P.; Skipsey, M.; Edwards, R. Roles for glutathione transferases in plant secondary metabolism. Phytochemistry 2010, 71, 338–350. [Google Scholar] [CrossRef]
  13. Tanaka, Y.; Sasaki, N.; Ohmiya, A. Biosynthesis of plant pigments: Anthocyanins, betalains and carotenoids. Plant J. 2008, 54, 733–749. [Google Scholar] [CrossRef]
  14. Petroni, K.; Tonelli, C. Recent advances on the regulation of anthocyanin synthesis in reproductive organs. Plant Sci. 2011, 181, 219–229. [Google Scholar] [CrossRef] [PubMed]
  15. Hichri, I.; Barrieu, F.; Bogs, J.; Kappel, C.; Delrot, S.; Lauvergeat, V. Recent advances in the transcriptional regulation of the flavonoid biosynthetic pathway. J. Exp. Bot. 2011, 62, 2465–2483. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Koes, R.; Verweij, W.; Quattrocchio, F. Flavonoids: A colorful model for the regulation and evolution of biochemical pathways. Trends Plant Sci. 2005, 10, 236–242. [Google Scholar] [CrossRef]
  17. Mortazavi, A.; Williams, B.A.; McCue, K.; Schaeffer, L.; Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 2008, 5, 621–628. [Google Scholar] [CrossRef]
  18. Wang, J.; Islam, F.; Li, L.; Long, M.; Yang, C.; Jin, X.; Ali, B.; Mao, B.; Zhou, W. Complementary RNA-sequencing based transcriptomics and iTRAQ proteomics reveal the mechanism of the alleviation of quinclorac stress by salicylic acid in oryza sativa ssp. japonica. Int. J. Mol. Sci. 2017, 18, 1975. [Google Scholar] [CrossRef] [PubMed]
  19. Haas, B.J.; Zody, M.C. Advancing RNA-Seq analysis. Nat. Biotechnol. 2010, 28, 421–423. [Google Scholar] [CrossRef]
  20. Akhter, D.; Qin, R.; Nath, U.K.; Eshag, J.; Jin, X.; Shi, C. A rice gene, OsPL, encoding a MYB family transcription factor confers anthocyanin synthesis, heat stress response and hormonal signaling. Gene 2019, 699, 62–72. [Google Scholar] [CrossRef] [PubMed]
  21. Zhuang, W.; Wang, H.; Liu, T.; Wang, T.; Zhang, F.; Shu, X.; Zhai, H.; Wang, Z. Integrated physiological and genomic analysis reveals structural variations and expression patterns of candidate genes for colored- and green-leaf poplar. Sci. Rep. 2019, 9, 11150. [Google Scholar] [CrossRef]
  22. Yang, S.S.; Tu, Z.J.; Cheung, F.; Xu, W.W.; Lamb, J.A.F.S.; Jung, H.J.G.; Vance, C.P.; Gronwald, J.W. Using RNA-Seq for gene identification, polymorphism detection and transcript profiling in two alfalfa genotypes with divergent cell wall composition in stems. BMC Genom. 2011, 12, 1–19. [Google Scholar] [CrossRef] [Green Version]
  23. Wang, Z.; Gerstein, M.; Snyder, M. Nrg2484-1. Nat. Rev. Genet. 2009, 10, 57–63. [Google Scholar] [CrossRef]
  24. Akhter, D.; Qin, R.; Nath, U.K.; Eshag, J.; Jin, X.; Shi, C. Transcriptional profile corroborates that bml mutant plays likely role in premature leaf senescence of rice (Oryza sativa L.). Int. J. Mol. Sci. 2019, 20, 1708. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Ma, Z.; Tan, Y.; Cui, G.; Feng, Y.; Cui, Q.; Song, X. Transcriptome and gene expression analysis of DHA producer Aurantiochytrium under low temperature conditions. Sci. Rep. 2015, 5, 1–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Feng, F.; Li, M.; Ma, F.; Cheng, L. Phenylpropanoid metabolites and expression of key genes involved inanthocyanin biosynthesis in the shaded peel of apple fruit in response to sun exposure. Plant Physiol. Biochem. 2013, 69, 54–61. [Google Scholar] [CrossRef] [PubMed]
  27. Bogs, J.; Ebadi, A.; McDavid, D.; Robinson, S.P. Identification of the flavonoid hydroxylases from grapevine and their regulation dining fruit development. Plant Physiol. 2006, 140, 279–291. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Shi, M.-Z.; Xie, D.-Y. Biosynthesis and Metabolic Engineering of Anthocyanins in Arabidopsis thaliana. Recent Pat. Biotechnol. 2014, 8, 47–60. [Google Scholar] [CrossRef] [Green Version]
  29. Shirley, B.W.; Kubasek, W.L.; Storz, G.; Bruggemann, E.; Koornneef, M.; Ausubel, F.M.; Goodman, H.M. Analysis of Arabidopsis mutants deficient in flavonoid biosynthesis. Plant J. 1995, 8, 659–671. [Google Scholar] [CrossRef] [PubMed]
  30. Huang, J.; Gu, M.; Lai, Z.; Fan, B.; Shi, K.; Zhou, Y.H.; Yu, J.Q.; Chen, Z. Functional analysis of the Arabidopsis PAL gene family in plant growth, development, and response to environmental stress. Plant Physiol. 2010, 153, 1526–1538. [Google Scholar] [CrossRef] [Green Version]
  31. Rohde, A.; Morreel, K.; Ralph, J.; Goeminne, G.; Hostyn, V.; De Rycke, R.; Kushnir, S.; Van Doorsselaere, J.; Joseleau, J.P.; Vuylsteke, M.; et al. Molecular phenotyping of the pal1 and pal2 mutants of Arabidopsis thaliana reveals far-reaching consequences on phenylpropanoid, amino acid, and carbohydrate metabolism. Plant Cell 2004, 16, 2749–2771. [Google Scholar] [CrossRef] [Green Version]
  32. Olsen, K.M.; Lea, U.S.; Slimestad, R.; Verheul, M.; Lillo, C. Differential expression of four Arabidopsis PAL genes; PAL1 and PAL2 have functional specialization in abiotic environmental-triggered flavonoid synthesis. J. Plant Physiol. 2008, 165, 1491–1499. [Google Scholar] [CrossRef]
  33. Stracke, R.; Werber, M.; Weisshaar, B. The R2R3-MYB gene family in Arabidopsis thaliana. Curr. Opin. Plant Biol. 2001, 4, 447–456. [Google Scholar] [CrossRef]
  34. Dubos, C.; Stracke, R.; Grotewold, E.; Weisshaar, B.; Martin, C.; Lepiniec, L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010, 15, 573–581. [Google Scholar] [CrossRef] [PubMed]
  35. Chiu, L.W.; Zhou, X.; Burke, S.; Wu, X.; Prior, R.L.; Li, L. The purple cauliflower arises from activation of a MYB transcription factor. Plant Physiol. 2010, 154, 1470–1480. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Feng, S.; Wang, Y.; Yang, S.; Xu, Y.; Chen, X. Anthocyanin biosynthesis in pears is regulated by a R2R3-MYB transcription factor PyMYB10. Planta 2010, 232, 245–255. [Google Scholar] [CrossRef]
  37. Mathews, H.; Clendennen, S.K.; Caldwell, C.G.; Liu, X.L.; Connors, K.; Matheis, N.; Schuster, D.K.; Menasco, D.J.; Wagoner, W.; Lightner, J.; et al. Activation tagging in tomato identifies a transcriptional regulator of anthocyanin biosynthesis, modification, and transport. Plant Cell 2003, 15, 1689–1703. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Borevitz, J.O.; Xia, Y.; Blount, J.; Dixon, R.A.; Lamb, C. Activation tagging identifies a conserved MYB regulator of phenylpropanoid biosynthesis. Plant Cell 2000, 12, 2383–2393. [Google Scholar] [CrossRef] [Green Version]
  39. An, J.P.; Zhang, X.W.; You, C.X.; Bi, S.Q.; Wang, X.F.; Hao, Y.J. MdWRKY40 promotes wounding-induced anthocyanin biosynthesis in association with MdMYB1 and undergoes MdBT2-mediated degradation. New Phytol. 2019, 224, 380–395. [Google Scholar] [CrossRef]
  40. He, Q.; Wu, J.; Xue, Y.; Zhao, W.; Li, R.; Zhang, L. The novel gene BrMYB2, located on chromosome A07, with a short intron 1 controls the purple-head trait of Chinese cabbage (Brassica rapa L.). Hortic. Res. 2020, 7. [Google Scholar] [CrossRef] [PubMed]
  41. Bai, S.; Tao, R.; Tang, Y.; Yin, L.; Ma, Y.; Ni, J.; Yan, X.; Yang, Q.; Wu, Z.; Zeng, Y.; et al. BBX16, a B-box protein, positively regulates light-induced anthocyanin accumulation by activating MYB10 in red pear. Plant Biotechnol. J. 2019, 17, 1985–1997. [Google Scholar] [CrossRef] [Green Version]
  42. Gonzalez, A.; Zhao, M.; Leavitt, J.M.; Lloyd, A.M. Regulation of the anthocyanin biosynthetic pathway by the TTG1/bHLH/Myb transcriptional complex in Arabidopsis seedlings. Plant J. 2008, 53, 814–827. [Google Scholar] [CrossRef] [PubMed]
  43. Spelt, C.; Quattrocchio, F.; Mol, J.N.M.; Koes, R. Anthocyanin1 of Petunia encodes a basic helix-loop-helix protein that directly activates transcription of structural anthocyanin genes. Plant Cell 2000, 12, 1619–1631. [Google Scholar] [CrossRef] [Green Version]
  44. Li, C.; Wu, J.; Hu, K.D.; Wei, S.W.; Sun, H.Y.; Hu, L.Y.; Han, Z.; Yao, G.F.; Zhang, H. PyWRKY26 and PybHLH3 cotargeted the PyMYB114 promoter to regulate anthocyanin biosynthesis and transport in red-skinned pears. Hortic. Res. 2020, 7. [Google Scholar] [CrossRef] [Green Version]
  45. An, J.P.; Wang, X.F.; Zhang, X.W.; Xu, H.F.; Bi, S.Q.; You, C.X.; Hao, Y.J. An apple MYB transcription factor regulates cold tolerance and anthocyanin accumulation and undergoes MIEL1-mediated degradation. Plant Biotechnol. J. 2020, 18, 337–353. [Google Scholar] [CrossRef] [Green Version]
  46. Mano, H.; Ogasawara, F.; Sato, K.; Higo, H.; Minobe, Y. Isolation of a regulatory gene of anthocyanin biosynthesis in tuberous roots of purple-fleshed sweet potato. Plant Physiol. 2007, 143, 1252–1268. [Google Scholar] [CrossRef] [Green Version]
  47. Cock, P.J.A.; Fields, C.J.; Goto, N.; Heuer, M.L.; Rice, P.M. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2009, 38, 1767–1771. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Li, B.; Colin, N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 1–16. [Google Scholar] [CrossRef] [Green Version]
  51. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 1–21. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. ClusterProfiler: An R package for comparing biological themes among gene clusters. Omics J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef] [PubMed]
  53. Zeng, D.D.; Qin, R.; Li, M.; Alamin, M.; Jin, X.L.; Liu, Y.; Shi, C.H. The ferredoxin-dependent glutamate synthase (OsFd-GOGAT) participates in leaf senescence and the nitrogen remobilization in rice. Mol. Genet. Genom. 2017, 292, 385–395. [Google Scholar] [CrossRef]
  54. 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]
Ijms 22 09787 g001 550
Figure 1. Phenotypic features of WT and pl mutant rice plants. Plants at (A) seedling stage, bar: 2 cm; (B) tillering stage, bar: 10 cm; (C) grain filling stage, bar: 10 cm; (D) pl plants showing purple color at field condition; (E) magnification of flag leaf from the top of the plant, bar = 5 cm.
Figure 1. Phenotypic features of WT and pl mutant rice plants. Plants at (A) seedling stage, bar: 2 cm; (B) tillering stage, bar: 10 cm; (C) grain filling stage, bar: 10 cm; (D) pl plants showing purple color at field condition; (E) magnification of flag leaf from the top of the plant, bar = 5 cm.
Ijms 22 09787 g001
Ijms 22 09787 g002 550
Figure 2. Correlation coefficient using inter-sample RNA-seq values (normalized reads) presented as heat map. Color blocks indicate the correlation coefficient values. The gray color indicates a lower correlation and the orange color represents a higher correlation coefficient between samples.
Figure 2. Correlation coefficient using inter-sample RNA-seq values (normalized reads) presented as heat map. Color blocks indicate the correlation coefficient values. The gray color indicates a lower correlation and the orange color represents a higher correlation coefficient between samples.
Ijms 22 09787 g002
Ijms 22 09787 g003 550
Figure 3. Volcano plot and Venn diagram of DEGs. (A) DEGs are presented by volcano map for pl vs. WT. Horizontal axis shows fold-change value (log2(B/A)) of the differential expression of the genes in different sample groups, whereas the vertical axis shows the level of significance with p value depicting the changes of gene expression. Every dot in the figure denotes a gene: green as down-regulation, red as up-regulation, and black as non-differentiated genes. (B) DEGs for pl vs. WT at grain filling stage are plotted in Venn diagram. Comparisons between different groups are characterized by different colors; numbers in the figure signify the number of shared/specific DEGs. The overlapping regions show the number of shared DEGs by various sets and the non-overlapped regions show the number of unique DEGs.
Figure 3. Volcano plot and Venn diagram of DEGs. (A) DEGs are presented by volcano map for pl vs. WT. Horizontal axis shows fold-change value (log2(B/A)) of the differential expression of the genes in different sample groups, whereas the vertical axis shows the level of significance with p value depicting the changes of gene expression. Every dot in the figure denotes a gene: green as down-regulation, red as up-regulation, and black as non-differentiated genes. (B) DEGs for pl vs. WT at grain filling stage are plotted in Venn diagram. Comparisons between different groups are characterized by different colors; numbers in the figure signify the number of shared/specific DEGs. The overlapping regions show the number of shared DEGs by various sets and the non-overlapped regions show the number of unique DEGs.
Ijms 22 09787 g003
Ijms 22 09787 g004 550
Figure 4. GO analysis of the functional enrichments of down- and up-regulated DEGs between pl vs. WT at late grain filling stage. Here, horizontal axis represents functional category and vertical axis represents the gene number in the category (right) and the percentage of the total genes on the annotation (left). Various colors signify various categories. Light color on the axis and histogram represents the differentially expressed genes, and dark color represents total number of genes in the respective categories.
Figure 4. GO analysis of the functional enrichments of down- and up-regulated DEGs between pl vs. WT at late grain filling stage. Here, horizontal axis represents functional category and vertical axis represents the gene number in the category (right) and the percentage of the total genes on the annotation (left). Various colors signify various categories. Light color on the axis and histogram represents the differentially expressed genes, and dark color represents total number of genes in the respective categories.
Ijms 22 09787 g004
Ijms 22 09787 g005 550
Figure 5. KEGG analysis of DEGs for pl vs. WT by bar plot and scattered diagram. (A) DEGs analysis by bar plot. Horizontal and vertical axes represent functional category and number of genes in the category (right) and percentage of total genes on the annotation (left), respectively. Colors represent different categories of KEGG enrichment. Light color in the histogram represents the differentially expressed genes and dark color represents total number of genes. (B) Scattered plot represents the important enrichment functions. Horizontal axis represents the rich factor function, and the vertical axis represents functional annotation. Top 30 GOs were used for analyzing the Q value, where size is shown by dot color and the smaller Q value is close to red.
Figure 5. KEGG analysis of DEGs for pl vs. WT by bar plot and scattered diagram. (A) DEGs analysis by bar plot. Horizontal and vertical axes represent functional category and number of genes in the category (right) and percentage of total genes on the annotation (left), respectively. Colors represent different categories of KEGG enrichment. Light color in the histogram represents the differentially expressed genes and dark color represents total number of genes. (B) Scattered plot represents the important enrichment functions. Horizontal axis represents the rich factor function, and the vertical axis represents functional annotation. Top 30 GOs were used for analyzing the Q value, where size is shown by dot color and the smaller Q value is close to red.
Ijms 22 09787 g005
Ijms 22 09787 g006 550
Figure 6. Significantly enriched functions of gene interaction network in KEGG. Square, circular, and edges nodes represent the function information, genes, and correlation between genes and their functions, respectively. The node size is proportional to the degree of connectivity of the nodes; therefore, more edges are connected to the larger node. The color of the circular node represents degree of difference in gene expression. In the log2FC value, green and red color represent down- and upward adjustment, respectively, and the intensity of color indicates degree of difference between up- and downward adjustments. The color of the square node represents the p value of enrichment program of the function. High enrichment corresponds to low p value with dark color. Square node with comparatively larger area indicates a greater number of DEGs involved in the biological phenomena. DEGs related to this function are drawn by taking only top 10 functions with the highest degree of enrichment.
Figure 6. Significantly enriched functions of gene interaction network in KEGG. Square, circular, and edges nodes represent the function information, genes, and correlation between genes and their functions, respectively. The node size is proportional to the degree of connectivity of the nodes; therefore, more edges are connected to the larger node. The color of the circular node represents degree of difference in gene expression. In the log2FC value, green and red color represent down- and upward adjustment, respectively, and the intensity of color indicates degree of difference between up- and downward adjustments. The color of the square node represents the p value of enrichment program of the function. High enrichment corresponds to low p value with dark color. Square node with comparatively larger area indicates a greater number of DEGs involved in the biological phenomena. DEGs related to this function are drawn by taking only top 10 functions with the highest degree of enrichment.
Ijms 22 09787 g006
Ijms 22 09787 g007 550
Figure 7. COG/KOG analysis of the DEGs of pl vs. WT. Horizontal axis represents the functional category and number of genes (right) and vertical axis represents percentage of the total genes and the total number of DEGs (left). Different categories of DEGs are represented with a wide range of color. Light color in the histogram represents the DEGs and dark color represents total gene number of that category.
Figure 7. COG/KOG analysis of the DEGs of pl vs. WT. Horizontal axis represents the functional category and number of genes (right) and vertical axis represents percentage of the total genes and the total number of DEGs (left). Different categories of DEGs are represented with a wide range of color. Light color in the histogram represents the DEGs and dark color represents total gene number of that category.
Ijms 22 09787 g007
Ijms 22 09787 g008 550
Figure 8. Expression of the transcription factor genes in anthocyanin biosynthesis in pl leaf as compared to the WT leaves. * indicates significant difference at p ≤ 0.05.
Figure 8. Expression of the transcription factor genes in anthocyanin biosynthesis in pl leaf as compared to the WT leaves. * indicates significant difference at p ≤ 0.05.
Ijms 22 09787 g008
Ijms 22 09787 g009 550
Figure 9. Relative expression patterns of the early and late-structural genes of anthocyanin biosynthesis in leaf samples of pl and WT. Colored scales indicate the levels of gene expression.
Figure 9. Relative expression patterns of the early and late-structural genes of anthocyanin biosynthesis in leaf samples of pl and WT. Colored scales indicate the levels of gene expression.
Ijms 22 09787 g009
Ijms 22 09787 g010 550
Figure 10. Proposed model for summarizing the transcriptional complex; MYB, bHLH, and WD40 proteins coordinately regulate anthocyanin accumulation through binding to the anthocyanin biosynthetic structural genes.
Figure 10. Proposed model for summarizing the transcriptional complex; MYB, bHLH, and WD40 proteins coordinately regulate anthocyanin accumulation through binding to the anthocyanin biosynthetic structural genes.
Ijms 22 09787 g010
Table
Table 1. Summary of reads of the samples retrieved from Illumina sequencing.
Table 1. Summary of reads of the samples retrieved from Illumina sequencing.
Sample NameRaw ReadsAverage Raw Read Length (bp)Clean ReadsAverage Clean Read Length (bp)Error (%)Q20
(%)		Q30
(%)		GC Content
(%)
WT55,160,75015053,875,169144.140.0099.2097.0657.83
pl54,237,83415052,941,300142.990.0099.1897.0357.71
Table
Table 2. Statistical summary of alignment results of RNA-seq used in mapping.
Table 2. Statistical summary of alignment results of RNA-seq used in mapping.
Sample NameTotal Clean Reads in MappingClean Read Mapped (%)Multiple Mapped (%)Uniquely Mapped (%)Non-Splice Reads (%)Splice Reads (%)Reads Mapped in Proper Pairs (%)
WT45,766,95684.952.7482.2155.1327.0881.31
pl45,900,10786.706.0677.9951.7226.2883.06
Table
Table 3. Relative gene expression analysis of 20 randomly selected genes for making the comparison of pl and WT group conducted by RNA- seq and qRT PCR.
Table 3. Relative gene expression analysis of 20 randomly selected genes for making the comparison of pl and WT group conducted by RNA- seq and qRT PCR.
Gene IDGene Name *pl vs. WTGene DescriptionPathway/Function Involved
RNA-seq log2 fold changeqRT PCR fold change
Os02g0626100OsPAL13.192868323.01 ± 0.01Similar to Phenylalanine ammonia-lyaseAnthocyanin Biosynthesis
Os02g0626400OsPAL22.109962311.93 ± 0.03Phenylalanine ammonia-lyasePhenylpropanoid metabolism, Anthocyanin Biosynthesis
Os04g0518400OsPAL65.027380523.10 ± 0.19Similar to Phenylalanine ammonia-lyasePhenylpropanoid metabolism, Anthocyanin Biosynthesis
Os05g0427400OsPAL73.038833431.25 ± 0.10Similar to Phenylalanine ammonia-lyasePhenylpropanoid metabolism
Os02g0177600Os4CL31.828286481.19 ± 0.104-coumarate: coenzyme A ligasePhenylpropanoid metabolism
Os08g0448000Os4CL54.871148542.08 ± 0.144-coumarate: coenzyme A ligase 5Phenylpropanoid metabolism
Os03g0122300OsF3H4.764684872.37 ± 0.24flavonol synthaseFlavonoid biosynthesis, Anthocyanin Biosynthesis
Os08g0441500 4.48865033.93 ± 0.015cinnamoyl-CoA reductasePhenylpropanoid metabolism, Anthocyanin Biosynthesis
Os05g0320700 1.82828641.65 ± 0.31cinnamate 4-hydroxylase 2Phenylpropanoid metabolism, Anthocyanin Biosynthesis
Os10g0512400OsCAld5H1; F5H13.16127752.57 ± 0.15coniferaldehyde 5-hydroxylase 1, ferulate 5-hydroxylase 1Lignin biosynthesis
Os01g0196300DPF; OsbHLH0255.60103493.65 ± 0.020basic helix-loop-helix protein 025
Os01g0838350Phosphate-limitation inducible gene 1−1.0446037−0.59 ± 0.09Conserved hypothetical protein
Os01g0734800 1.1200711.65 ± 0.10anthocyanidin 5,3-O-glucosyltransferaseAnthocyanin Biosynthesis
Os11g0116300 2.702188282.25 ± 0.02chalcone isomeraseAnthocyanin Biosynthesis
Os02g0611800 2.381978851.26 ± 0.08Hydroxycinnamoyltransferase 2Phenylpropanoid and Flavonoid biosynthesis
Os02g0194700OsLOX1−2.9658063−0.68 ± 0.06Lipoxygenase 1
Os11g0610700OsWD40-1901.29405961.71 ± 0.10WD40 repeat-like domain containing protein
Os07g0543100 −2.9658063−0.30 ± 0.01Beta-amylase 1
Os02g0738100 3.247568521.19 ± 0.13basic helix-loop-helix protein
Os04g0662600 4.4920931.66 ± 0.07Flavanone 3-hydroxylase 1Flavonoid biosynthesis
* Gene IDs are arranged following the gene order of the anthocyanin biosynthetic pathway.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Xu, R.; Pan, R.; Zhang, Y.; Feng, Y.; Nath, U.K.; Gan, Y.; Shi, C.; Akhter, D. RNA-Seq-Based Profiling of pl Mutant Reveals Transcriptional Regulation of Anthocyanin Biosynthesis in Rice (Oryza sativa L.). Int. J. Mol. Sci. 2021, 22, 9787. https://doi.org/10.3390/ijms22189787

AMA Style

Xu R, Pan R, Zhang Y, Feng Y, Nath UK, Gan Y, Shi C, Akhter D. RNA-Seq-Based Profiling of pl Mutant Reveals Transcriptional Regulation of Anthocyanin Biosynthesis in Rice (Oryza sativa L.). International Journal of Molecular Sciences. 2021; 22(18):9787. https://doi.org/10.3390/ijms22189787

Chicago/Turabian Style

Xu, Ruonan, Ronghui Pan, Yuchan Zhang, Yanlei Feng, Ujjal Kumar Nath, Yinbo Gan, Chunhai Shi, and Delara Akhter. 2021. "RNA-Seq-Based Profiling of pl Mutant Reveals Transcriptional Regulation of Anthocyanin Biosynthesis in Rice (Oryza sativa L.)" International Journal of Molecular Sciences 22, no. 18: 9787. https://doi.org/10.3390/ijms22189787

APA Style

Xu, R., Pan, R., Zhang, Y., Feng, Y., Nath, U. K., Gan, Y., Shi, C., & Akhter, D. (2021). RNA-Seq-Based Profiling of pl Mutant Reveals Transcriptional Regulation of Anthocyanin Biosynthesis in Rice (Oryza sativa L.). International Journal of Molecular Sciences, 22(18), 9787. https://doi.org/10.3390/ijms22189787

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