Transcriptome Analysis of Rice Seedling Roots in Response to Potassium Deficiency

Rice is one of the most important food crops in the world, and its growth, development, yield, and grain quality are susceptible to a deficiency of the macronutrient potassium (K+). The molecular mechanism for K+ deficiency tolerance remains poorly understood. In this study, K+ deficient conditions were employed to investigate the resulting changes in the transcriptome of rice seedling roots. Using ribonucleic acid sequencing (RNA-Seq) and analysis, a total of 805 differentially expressed genes were obtained, of which 536 genes were upregulated and 269 were downregulated. Gene functional classification showed that the expression of genes involved in nutrient transport, protein kinases, transcription processes, and plant hormones were particularly altered in the roots. Although these changes were significant, the expression of most genes remained constant even in K+-deficient conditions. Interestingly, when our RNA-Seq results were compared to public microarray data, we found that most of the genes that were differentially expressed in low K+ conditions also exhibited changes in expression in other environmental stress conditions.

a theoretical foundation for the selection and cultivation of rice cultivars with low potassium tolerance. With the development of molecular techniques, more research on the molecular genetic mechanisms of potassium nutritional activities in plants is currently being conducted. Many genes for potassium transporters, potassium ion channels, and related regulatory proteins have been cloned to study K + flux and protein localization. In addition, forward and reverse genetic analyses have become important molecular tools for observing the regulation of potassium nutritional activities (absorption and transport) in plants 1,12 . Although genes related to efficient nutritional activity have already been cloned for functional and regulatory analyses and have applications in the improvement of potassium nutrition in crops, investigations focused on the functional genomics and molecular genetic regulatory networks for high nutritional efficiency in plants are still ongoing.
Transcriptomics has already been broadly adopted in functional genomics research. With the rapid development of sequencing technology, the application of RNA-Seq has become an important new method for gene expression and transcriptome analysis 13 . When compared to gene microarray technology, RNA-Seq has an even greater sensitivity and resolution for gene expression research and transcriptional profiling 14,15 . In addition, it has advantages for identifying new transcripts, single nucleotide polymorphisms, gene structures, and splice variants 16,17 . Currently, this technology is already broadly utilized in the study of differential gene expression during plant responses to various biological and nonbiological stressors [18][19][20][21][22][23] .
We previously screened a group of rice germplasm accessions that are tolerant of low potassium during the seedling phase 24 . Additional experiments have confirmed that multiple genes regulate low potassium tolerance in rice and a partial quantitative trait locus for low potassium stress was identified 25 . However, there are few data on the transcriptomics of the response to low potassium stress in rice. The present study used Nipponbare rice as source material and Illumina RNA-Seq technology combined with real-time fluorescent quantitative PCR (QPCR) to analyze and compare the rice root transcriptome under low and normal potassium conditions, and to identify signal transduction pathways and regulatory networks related to low potassium tolerance.

Results
Transcriptome sequencing in rice roots identified differentially expressed genes sensitive to K + deficiency in rice. An Illumina HiSeq ™ 2500 was used to conduct high-throughput transcriptome analysis of control and mixed low potassium stress-treated rice root samples at three time points; 11,432,307 and 15,120,257 clean reads were obtained, respectively, and accounted for over 98% of the total sequences (Table S1). FASTQC was used to assess sequence quality in both control and treatment samples ( Figure S2). Transcriptome analysis was conducted using the complete genomic sequence of Nipponbare rice as a reference. The length of the longest gene obtained was over 4,000 bp; the lengths of most genes were distributed between 1,000 and 2,500 bp (Table S3). Most gene sequences were well matched with the reference genome (Table S4). The majority of the matched genes had a coverage of 60% or more ( Figure S5 and Fig. 1a).
The reads per kilobase of transcript per million mapped reads calculation method was employed to compare transcriptome differences between the control and the treated group. This study yielded a total of 805 differentially expressed genes, of which 536 genes were upregulated and 269 genes were downregulated (Table S6 and Fig. 1b). A threshold of two-fold or above was used as a standard for evaluating differential gene expression. Among these, the expression level of the most upregulated gene was 52-fold that of the control, and the expression level of the most downregulated gene was 79-fold below that of the control (Table S6). There were 98 downregulated genes with a three-fold or more decrease in expression level, and there were 219 upregulated genes with a three-fold or more expression level increase (Table S6). Among genes with a five-fold or higher difference in expression, 58 genes were upregulated and 47 genes were downregulated. Among genes with a 10-fold or higher difference in expression, 19 genes were downregulated and 14 genes were upregulated (Table S6). Functional classification of differentially expressed genes identified in low K + conditions. To better understand the function of the genes affected by low K + , we used GoPipe to conduct GO classification on the differentially expressed genes. A total of 12,414 predicted proteins were matched with 81,604 GO terms. Among the molecular functions, genes that encoded proteins with binding and catalytic activity accounted for a large portion. Additionally, genes encoding proteins localized in organelles also accounted for a large number of those that were differentially expressed under K + deficiency. Many of these genes participate in cell growth and metabolism, and their expression was stimulated by adverse environmental conditions (Fig. 2a). GO classification analysis of the differentially expressed genes revealed that these genes involved in a very large portion of the whole life process (Fig. 2b,c).
To further understand the functions of the differentially expressed genes and the related biological processes they participate in, GO and KEGG enrichment analysis were conducted (Fig. 3). GO enrichment showed that genes with functions related to catalytic activity, binding activity, and nutrient storage were up or downregulated under conditions of K + stress. The biological processes that these genes primarily participate in involve metabolism, growth and development, multicellular processes, and signal transduction processes (Fig. 3a). KEGG analysis also showed that the differentially expressed genes were primarily involved in metabolic pathways and signal transduction (Fig. 3b). Enrichment analysis also showed that the expression of genes involved in replication and transcriptional pathways are strongly affected, implying that K + stress also affects the normal biological functions of root cells.
Real-time fluorescent quantitative analysis confirms differential expression of genes in response to low K + . Transcriptome analysis results show that the expression of many genes changed. To further validate the accuracy of these expression differences, 24 genes among the differentially expressed genes were selected for fluorescent quantitative validation according to their expression pattern (Table S7). Of these, 14 genes were upregulated and 11 genes were downregulated (Fig. 4). Analysis results showed that there was good correlation between the RNA-Seq and QPCR experiments (R 2 = 0.9647), and that the changes in expression of most differentially expressed genes agreed with the results of the quantitative analysis.
Gene families are differentially expressed under low K + . Rice has many gene families, and members of the same family all have very important effects upon the same biological processes. To further understand the gene families whose expression was affected by low K + stress, we identified families represented in our list of differentially expressed genes (Table S8). Results showed that members from gene families of core histones H2A/ H2B/H3/H4, peroxidases, protein kinases, and transmembrane amino acid transporter proteins closely related to ion transport accounted for many of the genes differentially expressed under low K + conditions. This implies that genes in these families were the most affected by low K + stress, and also shows that these gene families may be important for plant adaptability to low K + stress conditions (Fig. 5).
Identification of transcription factors that respond to low K + . High-throughput sequencing is a more sensitive assay for gene expression when compared to microarrays; it can capture the expression of less abundant genes, such as transcription factors. Further analysis of differentially expressed genes in rice roots showed that there were 47 transcription factors belonging to 22 transcription factor families ( Table 1). Of these,  heat stress transcription factors and APETALA 2/ethylene-responsive element binding factor family members were the most numerous. The heat stress transcription factors family is related to various stress responses in plants, and the ethylene response factor family primarily affects DNA binding. Other transcription factors are involved in biological processes, including signal transduction and defense against adverse environmental conditions.

Construction of protein interaction networks responsive to K + stress.
To further understand the relationships between the differentially expressed genes in more detail and to search for possible critical pathways, protein-protein interaction analysis was conducted on the genes identified in our sequencing experiment ( Figure S9). The gene with the most abundant connection points was the Putative activator of 90-kDa heat shock protein adenosinetriphosphatase homolog 1, which was downregulated in this study ( Figure S9). Because life processes were highly related with energy, the identification of many differentially expressed genes involved in energy-related process seems quite reasonable. Results showed that upregulated genes were involved in transcription, implying that K + stress affects the transcription process and the role of K + involved in transcription needs further investigation.
Microarray data analysis confirms genes related to K + deficiency are also involved in other environmental stress responses. To further understand the behavior of the differentially expressed genes identified in this experiment in other adverse environments and stresses, microarray data from public databases were used to analyze their expression profiles under conditions of pathogen assault 32 , drought, salt, and low temperature stress 33 . Results showed that most genes that were significantly differentially expressed in low K + were also responsive to other adverse environmental conditions, suggesting that genes related to K + stress are not limited to ion transport, but also play important roles in plant responses to several adverse environmental conditions (Figures S10, S11 and S12).

Discussion
While K + is important for plant growth and development, the amount of K + that plants can absorb directly from the soil is very limited. This is not only because the amount of free and readily useable K + in the soil is low, but it is also due to an increased grain yield and cropping index, which depletes K + from the soil and leads to insufficiently supplemented soils over time. According to statistical analyses, one-quarter of the soil in the global tropical and subtropical zones is in a state of K + deficiency 34 . This forces plants in these regions to grow in low K + environments for long periods. To adapt to this environment, plants must take various actions to promote K + absorption. To effectively absorb K + in environments with low potassium availability, plants initiate two important mechanisms to improve the efficiency of K + uptake 35 . In addition, to transport potassium more effectively, K + transport also needs to be more efficient. Nutrient content and availability in the soil complicate K + uptake, as does how plants access and utilize the nutrient; this suggests that regulation at the molecular and cellular levels may also be highly complex. All of these activities require the participation of a series of genes and may also require certain organs to change shape (Figure S13) 1 .
Low K + stress affects the growth and development of rice. This study shows that low K + stress during the seedling phase results in increased Malonic Dialdehyde (MDA) content, MDA is an important indicator for membrane lipid peroxidation, which indicates membrane system damage, suggesting an oxidative stress aggravation in  roots. Correspondingly, the activity of some antioxidases also increased ( Figure S14). As the primary plant organ that directly contacts the soil, the roots sense and absorb K + . Analyses of mechanisms related to K + uptake and roots are integral to the understanding of the molecular regulation induced by low K + conditions and provide a model for addressing plant responses to K + soil deficiency. Plant roots are responsible for detecting external K + concentrations and their fluctuations during plant growth. Although no potassium receptor has been identified, research has shown that shaker transporter proteins are involved in the detection and absorption of K + in rice 48 . Significant changes in the expression of the shaker transporter proteins were not detected in the present study, which may be due to differences between their expression in the roots and in other parts of the plant. Protein kinases, particularly pyruvate kinases 36 , play an important role in intracellular detection of K + . The present study found two similar pyruvate kinases, one upregulated and the other downregulated. In addition, protein kinases also have K + channel regulating activity 1 . These genes should be closely related to the regulation of K + channels.
We identified genes clearly annotated as encoding protein kinases; 21 genes were upregulated and six were downregulated. In addition, the upregulated protein kinases included mitogen-activated protein kinase pathway genes, implying that these kinases may be responsible for K + signal transduction. Previous research had always considered calcium and reactive oxygen species (ROS) as important secondary messengers of K + signal transduction 1 ; however, kinase signal transduction functions have not been studied in detail. Plant hormones are required in each of the growth phases of plants. In low K + stress conditions, hormones similarly play an important role. In this study, genes related to jasmonic acid, auxin, and ethylene were the most abundant genes involved hormone-related processes, which is consistent with previous reports 37 . Two genes in the ASR gene family in rice, OsASR4 and OsASR6, were detected in this study and were both downregulated. Previous research has shown that this gene family has different effects on different hormones and has a certain enhanced effect under adverse environmental conditions 38 . Additionally, the genes in this family are related to aluminum stress responses 39 , implying that there may be crosstalk between heavy metal stress and potassium stress response pathways in plants. However, previous research has concentrated on interaction between K + and other ions such as Na + and NH 4 + and has not focused on the relationship between heavy metals and potassium ions. Pyrrolysine gene family members, which are possible receptors for abscisic acid 40 , were regulated somewhat differently; OsPYL6 was downregulated and OsPYL9 was upregulated. These two genes can significantly increase drought and cold tolerance in rice 40 , suggesting that there may be a close relationship between K + stress and other nonbiological stressors. The OsHAK1 gene is the rice homolog of the Arabidopsis HAK gene and belongs to the KT/KUP/HAK gene family, which is closely related to potassium transport. New research by Chen et al. shows that just as the Arabidopsis HAK gene plays important roles in potassium transport, the OsHAK1 gene does so regardless of whether potassium conditions are low in rice 41 . In their study, the OsHAK1 gene was significantly upregulated under K + stress conditions 41 , but in this study, OsHAK1 was among the downregulated genes, possibly because of differences in sample mixing and the time of sample collection ( Figure S14). Gene expression studies under long-term K + stress conditions would benefit from the study of how the expression of genes that are essential for K + transport changes over time and from gaining a deeper understanding of how gene expression changes in rice under K + stress.
ROS have been an important component in research on K + stress 1 . The present study identified many peroxidases that were upregulated, indicating that ROS are important in plant responses to low potassium. However, there are currently very few reports on the interactions between ROS and K + . The level of participation of ROS in signal transduction and K + transport is not clear and requires further research. In addition, some of these peroxidases are also related to other biological stressors, and we observed that the gene expression for some proteins related to disease was upregulated. Thus, the relationship between K + stress and biological stressors warrants further investigation. Genes related to ethylene in this study were upregulated, showing that ethylene has a positive function in the response to K + stress. In addition, a large number of genes related to auxin was also upregulated, indicating the possibility of hormonal participation in plant morphology during K + stress. Genes related to rice root morphology had differential expression levels, implying that rice root morphology is also susceptible to K + stress.
The growth of rice is regulated by many environmental factors. This study showed that the effects of K + stress on rice are not only limited to kinase and transporter proteins but also have a very large effect through changes in the expression of genes encoding heat shock proteins and other transcription factor families, implying that K + stress is closely related to temperature. Ion transport in the plant body inevitably involves many transmembrane proteins. We found that the upregulation of genes for transmembrane proteins were significant for three members of the transforming growth factor beta receptor family. However, the role of this gene family in K + stress in rice has not been documented and awaits further experimental confirmation.
With the continued advancement of functional genomics research in rice, more K + transport proteins and their functions, such as OsHAK, OsAKT, and OsHKT, have been confirmed [42][43][44] . Recent research has found that there are many other ions that promote or inhibit K + absorption. Lan et al. found that the high-affinity K + transporter in rice can also mediate Ca 2+ transport 45 . Ma et al. utilized gene microarrays to analyze gene expression in rice under long-term and short-term K + stress and found that many groups of genes, especially protein kinases and metal ion transporter family, were significantly upregulated 46 . Although genes related to efficient nutritional activity have already been cloned and characterized, further investigation of the functional genomics and molecular genetic regulatory networks, and their possible application in the improvement of potassium nutritional activities in crops plants, should be a focus of future research 1,12,47 .
This study employed high-throughput sequencing using mixed materials with different processing times for the analysis of the rice root transcriptome under potassium stress conditions. The results of the study revealed that a large number of genes exhibit clear changes in expression in response to low K + growing conditions. Among the large number of genes differentially expressed in this study, many have unknown functions (Fig. 6), implying that research on the transcriptional changes in rice caused by K + stress is a rich area for future investigations.

Conclusions
Transcriptome sequencing was applied to perform transcriptomic analysis of mixed rice root samples at different times under K + stress conditions. Physiological experiments showed that low potassium treatment during the seedling phase resulted in increased Superoxide Dismutase (SOD) and Peroxidase (POD) activity and increased MDA content. Transcriptome analysis yielded a total of 805 differentially expressed genes, including genes related to plant hormones, peroxidases, and a series of genes related to the detection and transport of K + . The differentially expressed genes were related to root growth and development and included transcription factors from many different gene families, implicating a complex mechanism for the response to low potassium stress conditions in rice. Changes in rice growth and biological pathways under low potassium conditions are highly important for understanding K + detection and transport in rice varieties that are tolerant to low potassium.

Material and Methods
Plant materials and growth conditions. A hydroponic experiment was conducted in a plant growth chamber set to a photoperiod of 14 h light/10 h dark. Seeds of Nipponbare rice were sterilized with 2% H 2 O 2 for 30 min, rinsed with distilled water three times, and then incubated for two days in the dark at 26 °C. The seedlings were cultured and transplanted as previously described 25 . The nutrient solution was prepared according to an established protocol 26  O, and 3.6 × 10 −2 mM FeCl 3 · 6H 2 O. The pH value of the culture solution was adjusted to 5.1 using 1 M HCl or NaOH solution as required every two days. Ten days after transplanting the seedlings to the full-strength nutrient solution, rice seedlings were subjected to low-K treatment. The potassium concentration was adjusted to 0.01 mM (low-K treatment) or 1 mM for the control. The experiment was conducted as a split-plot design with three independent replicates.

RNA-Seq sampling and RNA isolation.
For RNA-Seq sampling, seedlings at the three-leaf-stage were exposed to low-K + stress (0.01 mM) for 24 h, 48 h, and 84 h, respectively. The roots of 12 seedlings were collected and pooled together, immediately frozen in liquid nitrogen, and stored at −80 °C. RNA was isolated according to the instructions of the miRNeasy mini kit (QIAGEN, Germany). RNA sample quantity and purity were tested to confirm they met the protocol requirements. Subsequently, equal amounts of total RNA from 24 h, 48 h, and 84 h were mixed and stored. RNA extraction. TRIZOL reagent (Invitrogen, USA) was used for RNA extraction per the manufacturer's instructions. Extracted RNA was stored at −80 °C until it was used for transcriptome sequencing and real-time fluorescent quantitative PCR validation. Transcriptome sequencing. Total RNA was treated with 10 U DNaseI (Ambion, USA) at 37 °C, and then mRNA was purified using the Micropoly(A) Purist mRNA purification kit (Ambion). First-strand complementary deoxyribonucleic acid (cDNA) synthesis was performed using GsuI-oligo(dT) (Invitrogen) as the primer for reverse transcription, 10 µg mRNA as template, and Superscript II reverse transcriptase (Invitrogen) for elongation at 42 °C. Next, NaIO 4 (Sigma, USA) was applied to oxidize 5′ mRNA cap structures and conjugate them to biotin. Dynal M280 magnetic beads (Invitrogen) were used to screen for biotin-conjugated mRNA/cDNA, and cDNA was released by alkaline hydrolysis. Next, DNA ligase (TaKaRa, Japan) was employed to add an adapter to the 5′ end of the first-strand cDNA, and ExTaq polymerase (TaKaRa) was used to synthesize the second strand of cDNA. Finally, the GsuI restriction enzyme was used to cleave the poly(A) tail and the 5′ end adapter. Synthesized cDNA was broken into 300-500-bp fragments using a sonicator (Fisher) and then purified with Ampure beads (Agencourt, USA). The purified cDNA was used to produce a library with the TruSeq TM DNA Sample Prep Kit -Set A (Illumina, USA), followed by amplification with the TruSeq PE Cluster Kit (Illumina). Finally, the sequencing reaction was performed on an Illumina Hiseq 2500.
Sequencing data processing. Clean reads were obtained using the FASTX-Toolkit software (http://hannonlab.cshl.edu/fastx_toolkit/) and the following data processing procedure and statistical analysis: 1) FASTX Clipper tool was used to remove linker sequences from the reads; 2) FASTX Quality Filter tool was used to remove unconfirmed bases from the reads in the 3′ to 5′ direction until the first confirmed base was reached; 3) After the FASTX Quality Filter removed successive low-quality bases, if the length of clean reads was less than 40 bp, the sequence and its complement were deleted; 4) A local script was used to complement pair end reads (using pair-end sequences).
Clean read mapping was performed with a provided reference genome (http://plants.ensembl.org/Oryza_ sativa/Info/Index). Statistical analysis of mapping data was performed as follows. Gene Ontology (GO) analysis was performed using GoPipe 27 . Predicted proteins were first compared with the SwissProt and TrEMBL 28 databases using the Protein Basic Local Alignment Search Tool and an e-value < 1e-5. Next, the GoPipe program was Scientific RepoRts | 7: 5523 | DOI:10.1038/s41598-017-05887-9 used to compare results based on gene2go (National Center for Biotechnology Information), and GO data of predicted proteins were obtained.
Predicted proteins were compared using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database with bidirectional BLAST (http://www.genome.jp/tools/kaas/help.html), and the KEGG ORTHOLOGY number of predicted proteins was obtained. Based on the KEGG ORTHOLOGY number, data for the metabolic pathways predicted proteins participated in were obtained.
Statistical analysis of genes with read numbers from two samples was converted into reads per kilobase of transcript per million mapped reads 29 . Finally, the microarray-plot-based method with random sampling model) in the DEGseq program package 30 was adopted to calculate the difference in expression abundance between the two samples for each gene. False discovery rate values less than 0.001 were considered to be significantly different.
Construction of protein-protein interaction networks. Medusa software 31 was used for data display based on the comparison of the protein-protein interaction data from rice (http://arabidopsisreactome.org/download/all_interactions.html).

QPCR validation.
To further evaluate the accuracy of the sequencing data, 30 genes with significant changes in their expression were selected for investigation. Fluorescent quantitative analysis of gene expression from rice roots grown in K + stress conditions and correlation analysis of transcriptome data were conducted. Primers and probes for quantitative analysis used in this study were designed using the online tool at https://www.genscript. com/ssl-bin/app/primer.