Transcriptomic Analysis of Rice Plants Overexpressing PsGAPDH in Response to Salinity Stress

In plants, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) is a main enzyme in the glycolytic pathway. It plays an essential role in glycerolipid metabolism and response to various stresses. To examine the function of PsGAPDH (Pleurotus sajor-caju GAPDH) in response to abiotic stress, we generated transgenic rice plants with single-copy/intergenic/homozygous overexpression PsGAPDH (PsGAPDH-OX) and investigated their responses to salinity stress. Seedling growth and germination rates of PsGAPDH-OX were significantly increased under salt stress conditions compared to those of the wild type. To elucidate the role of PsGAPDH-OX in salt stress tolerance of rice, an Illumina HiSeq 2000 platform was used to analyze transcriptome profiles of leaves under salt stress. Analysis results of sequencing data showed that 1124 transcripts were differentially expressed. Using the list of differentially expressed genes (DEGs), functional enrichment analyses of DEGs such as Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were performed. KEGG pathway enrichment analysis revealed that unigenes exhibiting differential expression were involved in starch and sucrose metabolism. Interestingly, trehalose-6-phosphate synthase (TPS) genes, of which expression was enhanced by abiotic stress, showed a significant difference in PsGAPDH-OX. Findings of this study suggest that PsGAPDH plays a role in the adaptation of rice plants to salt stress.


Introduction
Plants continuously face environmental stressors, including drought, high salinity, and extreme temperatures. These stressors affect both biomasses and grain yields of crops. Salt stress is a particularly important abiotic stress that can seriously affect the development, growth, and productivity of plants. Severe stress may threaten their survival. Rice (Oryza sativa) is a major staple food crop in the world. It is also a model plant for the genomics research of monocotyledons [1]. Salinity is one of the most disastrous abiotic stresses for rice. Salt-affected soils currently account for approximately 20% of global agricultural production [2]. Recently, many studies have reported molecular and cellular mechanisms such as ERF1, OsMYB91, and OsGLYII-2 involved in the tolerance of rice to To construct an overexpression vector, full-length cDNA of PsGAPDH was cloned into a plant expression vector pPZP, which was constructed by introducing CaMV35S promoter and PinII terminator containing Bar as a plant selectable marker ( Figure 1A) [19,30]. The 35S::PsGAPDH construct was introduced into an Agrobacterium tumefaciens strain LB4404 to generate transgenic rice [19]. Rice plants were transformed by Agrobacterium infection as described previously [31]. Transgenic rice plants were regenerated from transformed calli on selection medium containing 4 mg/L phosphinothricin and 500 mg/L cefotaxime and subsequently grown in a greenhouse.
All tissue samples collected were immediately frozen in liquid nitrogen and stored at −70 °C.
To construct an overexpression vector, full-length cDNA of PsGAPDH was cloned into a plant expression vector pPZP, which was constructed by introducing CaMV35S promoter and PinII terminator containing Bar as a plant selectable marker ( Figure 1A) [19,30]. The 35S::PsGAPDH construct was introduced into an Agrobacterium tumefaciens strain LB4404 to generate transgenic rice [19]. Rice plants were transformed by Agrobacterium infection as described previously [31]. Transgenic rice plants were regenerated from transformed calli on selection medium containing 4 mg/L phosphinothricin and 500 mg/L cefotaxime and subsequently grown in a greenhouse.

Genomic DNA Extraction and TaqMan Copy Number Analysis
Genomic DNA of transgenic plants was extracted using an IncloneTM Genomic DNA Prep Kit (Inclone, Korea). Copy number assay was performed on an Applied Biosystems StepOnePlusTM (Applied Biosystems, Foster City, CA, USA) with a TaqMan ® Gene Expression Master Mix (Applied Biosystems) kit. Primers and probes of the predesigned TaqMan ® copy number assay for rice tubulin α-1 chain gene (AK102560) as an endogenous control were used. For the transgene, primers and probes were specifically

Genomic DNA Extraction and TaqMan Copy Number Analysis
Genomic DNA of transgenic plants was extracted using an IncloneTM Genomic DNA Prep Kit (Inclone, Korea). Copy number assay was performed on an Applied Biosystems StepOnePlusTM (Applied Biosystems, Foster City, CA, USA) with a TaqMan ® Gene Expression Master Mix (Applied Biosystems) kit. Primers and probes of the predesigned TaqMan ® copy number assay for rice tubulin α-1 chain gene (AK102560) as an endogenous control were used. For the transgene, primers and probes were specifically designed for the terminator of the nos gene. Gene-specific primers used were NOS-F 5 -GCATGACGTTATTTATGAGATGGGTTT-3 and NOS-R 5 -TGCGCGCTA TATTTTGTTTT-CTATCG-3 . NOS-probe 5 -TAGAGTCCCGCA ATTAT-3 was also used for the nos gene. PCR conditions were as follows: 10 min at 95 • C; followed by 40 cycles of 95 • C for 15 s, 60 • C for 1 min, and 72 • C for 1 min. In the reaction plate, each sample was measured in triplicate. To calculate the copy number of a target gene, relative quantitation analysis of genomic DNA target was analyzed based on real-time PCR data using Applied Biosystems CopyCaller ® Software v2.0 (Applied Biosystems) according to the manufacturer's instructions.

Transient Expression Analysis in Protoplast
The cDNA sequence of PsGAPDH was cloned with a GFP fusion protein (PsGAPDH: GFP) into a pJJ2485-GFP expression vector using a gateway recombination system (Invitrogen, Life technologies, Carlsbad, CA, USA). Protoplast transformation was performed as previously described [32]. Transformed protoplasts were observed with a confocal microscope (LSM 510 META, Carl Zeiss, Jena, Germany). Autofluorescence of chlorophyll was used as a chloroplast marker.

Salt Stress Treatment
Salt stress treatment for transgenic rice plants was performed using a published method [33]. For the phenotypic analysis of PsGAPDH-OX plants in response to salt stress, sterilized rice seeds were transferred into 1/2 Murashige and 1/2 MS medium containing 100 mM NaCl or 1/2 MS liquid medium containing 200 mM NaCl. Seeds were incubated at 28 • C for 5 days or 14 days. Germination and seedling growth of transgenic plants and wild-type seeds sown on MS agar plates at each nutrient level was monitored. Three replicates with 15 seeds in each replicate were used. All data were analyzed by t-tests. Significant differences are indicated by asterisks (*, p < 0.05; **, p < 0.01). All the data from control and treatments were subjected to statistical analysis using SPSS 20.0.
For RNA sequence analysis of PsGAPDH-OX in response to salt stress, two-week-old seedlings were transferred to soil culture in a greenhouse. After 7 weeks, all plants were subjected to treatment with 200 mM NaCl (salt stress) for 3 days. RNA sequence analyses of three biological replicates for each wild-type and transgenic plants were then performed.

RNA Extraction, cDNA Library Construction, and Sequencing
Total RNAs were extracted from rice leaves using an RNeasy ® Plant mini kit (QIAGEN, Hilden, Germany) according to the manufacturer's protocol. RNA quality was determined with a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Only samples with an RNA integrity number >8 were used for library preparation. Each paired-end cDNA library was constructed using a TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA).

Differential Expressed Gene (DEG) Analysis
Differentially expressed genes (DEGs) were identified using an edgeR Bioconductor package based on a Generalized Linear Model (GLM) used for RNA-seq data analysis, considering gene expression as a negative binomial [37]. The EdgeR-robust method [38] was used to reduce the effect of outlier genes. Based on quantified gene expression, relationships among samples and stages were investigated using a multidimensional scaling (MDS) method. A heatmap was generated to visualize DEGs which included the 500 most expressed genes (TMM upper value) across samples using R. A volcano plot was obtained to visualize differential expression between samples. Significant differences (Log-fold change > 2.0 or Log-fold change < −2.0, p-value < 0.05) between two different conditions were illustrated in the volcano plot. All identified proteins were annotated with the Os-Nipponbare-Reference-IRGSP-1.0 database (https://rapdb.dna.affrc.go.jp accessed on 26 March 2021). An adjusted p-value (FDR) < 0.05 was used as the significance cutoff for differentially expressed genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis
DEGs were further processed with DAVID 6.8 Functional Annotation Tool (http: //david.abcc.ncifcrf.gov/ accessed on 26 March 2021) for term enrichment analysis [39]. Results were filtered using a Fisher's exact statistic methodology as previously described [39].

Quantitative Real-Time PCR (qRT-PCR)
Transcriptome sequencing data were validated by qRT-PCR. Briefly, RNAs (1 µg) were reverse-transcribed to cDNAs using an AmfiRevert cDNA synthesis kit (Gendepot, USA). Synthesized cDNAs were used as templates for qRT-PCR. Primer pairs specific for the amplification of target genes are shown in Supplementary Table S4. Relative quantification was performed to calculate expression levels of target genes in different treatments using the 2 −∆∆Ct method. The expression level of OsActin1 was used for the normalization of real-time PCR results. All data are expressed as the mean ± SD from three independent experiments.

Generation of Rice Plants Overexpressing PsGAPDH
The PsGAPDH gene from Pleurotus sajor-caju (Oyster mushroom) is highly induced by salt, drought, cold, and heat stress conditions [16,17]. Interestingly, it has been reported that transgenic potato plants constitutively expressing the PsGAPDH gene can increase their salt tolerance under salt stress conditions [17]. Thus, rice plants overexpressing PsGAPDH (PsGAPDH-OX) were generated by introducing a construct containing the fulllength cDNA of PsGAPDH under the control of a CaMV35S promoter ( Figure 1A). A total of 149 T0 generation lines of PsGAPDH-OX rice plants were generated. TaqMan copy number assay was used to select single-copy lines from PsGAPDH-OX plants ( Figure  1B). Flanking sequence tag (FST) analysis of T0 plants was carried out only for those transgenic plants in which a single copy T-DNA insertion was confirmed. Transgenic plants with the T-DNA inserted in an intergenic region are less likely to be affected by other genes. Hence, single copy/intergenic transgenic plants were selected using FST analysis (Supplementary Table S1) [19]. We obtained five single copy/intergenic transgenic plants from 21 lines of transformants. PsGAPDH-OX (#5, #6, #10, and #17) plants with a singlecopy and an intergenic insertion were selected and separated to produce homozygous T3 generation. Transcriptional levels of the PsGAPDH gene in these four lines of singlecopy/intergenic/homozygous PsGAPDH-OX plants were determined using qRT-PCR. Results of qRT-PCR revealed that PsGAPDH was overexpressed in transgenic rice plants. PsGAPDH-OX #10 line showed a relatively low level of expression compared to other transgenic rice lines ( Figure 1C). Single-copy/intergenic/homozygous PsGAPDH-OX (#5, #6, and #17) plants were then selected and used in subsequent experiments.

Subcellular Localization of PsGAPDH
GAPDH catalyzes the conversion of glyceraldehyde-3-phosphate to 1,3-bisphospho glycerate. Both cytosolic (GAPCs) and plastidial (GAPCps) GAPDH activities have been described in plants. To determine the subcellular localization of PsGAPDH expression, the distribution of PsGAPDH in maize and Arabidopsis protoplasts was examined. PsGAPDH-GFP proteins were transiently expressed in protoplasts isolated from leaves using a polyethylene glycol (PEG)-mediated transfection system [32]. A fusion protein with GFP was expressed and its localization was visualized using a laser scanning confocal microscope. As shown in Figure 1D, PsGAPDH-GFP was predominantly present in the cytoplasm of Arabidopsis and maize. This result demonstrates that PsGAPDH is protein-localized in the cytoplasm.

Effect of PsGAPDH-OX on Seedling Growth under Salt Stress Condition
To assess the role of PsGAPDH in the salt stress response of rice, seeds of PsGAPDH-OX and wild-type rice were sown on MS agar plates containing different concentrations (100 mM and 200 mM) of NaCl. Their germination rates and seedling growths were monitored. There was no significant difference between the wild-type and PsGAPDH-OX cultured in MS medium without NaCl ( Figure 2). However, when they were cultured in MS medium containing 100 mM or 200 mM NaCl for 8 days, PsGAPDH-OX #5, #6, and #17 rice lines showed different germination rates from that of the wild-type. In MS medium containing 100 mM NaCl, the germination rate of wild-type rice seeds was 68%, while those of PsGAPDH-OX #5, #6, and #17 were 95%, 76%, and 85%, respectively ( Figure 2). In MS medium containing 200 mM NaCl, the germination rate of the wild-type rice seed was 49%, while those of PsGAPDH-OX #5, #6, and #17 were 60%, 67%, and 49%, respectively. After 14 days of culturing, PsGAPDH-OX #6 seedlings showed significantly better growth performance than the wild-type (p < 0.05) on MS media supplemented with 100 mM or 200 mM NaCl ( Figure 3). These results suggested that PsGAPDH-OX lines were resistant to salt stress.
were then selected and used in subsequent experiments.

Subcellular Localization of PsGAPDH
GAPDH catalyzes the conversion of glyceraldehyde-3-phosphate to 1,3-bisphospho glycerate. Both cytosolic (GAPCs) and plastidial (GAPCps) GAPDH activities have been described in plants. To determine the subcellular localization of PsGAPDH expression, the distribution of PsGAPDH in maize and Arabidopsis protoplasts was examined PsGAPDH-GFP proteins were transiently expressed in protoplasts isolated from leave using a polyethylene glycol (PEG)-mediated transfection system [32]. A fusion protein with GFP was expressed and its localization was visualized using a laser scanning confo cal microscope. As shown in Figure 1D, PsGAPDH-GFP was predominantly present in the cytoplasm of Arabidopsis and maize. This result demonstrates that PsGAPDH is pro tein-localized in the cytoplasm.

Effect of PsGAPDH-OX on Seedling Growth under Salt Stress Condition
To assess the role of PsGAPDH in the salt stress response of rice, seeds of PsGAPDH OX and wild-type rice were sown on MS agar plates containing different concentration (100 mM and 200 mM) of NaCl. Their germination rates and seedling growths were mon itored. There was no significant difference between the wild-type and PsGAPDH-OX cul tured in MS medium without NaCl ( Figure 2). However, when they were cultured in MS medium containing 100 mM or 200 mM NaCl for 8 days, PsGAPDH-OX #5, #6, and #17 rice lines showed different germination rates from that of the wild-type. In MS medium containing 100 mM NaCl, the germination rate of wild-type rice seeds was 68%, while those of PsGAPDH-OX #5, #6, and #17 were 95%, 76%, and 85%, respectively ( Figure 2) In MS medium containing 200 mM NaCl, the germination rate of the wild-type rice seed was 49%, while those of PsGAPDH-OX #5, #6, and #17 were 60%, 67%, and 49%, respec tively. After 14 days of culturing, PsGAPDH-OX #6 seedlings showed significantly bette growth performance than the wild-type (p < 0.05) on MS media supplemented with 100 mM or 200 mM NaCl (Figure 3). These results suggested that PsGAPDH-OX lines were resistant to salt stress.

Statistical Test for DEG
Based on quantified gene expression, relationships among samples and stages were investigated using the MDS method. Multi-dimensional plot indicated that samples could be distinguished by transgenic rice (Figure 4). After checking distances of samples, statistical tests were performed to identify DEGs between control and transgenic samples. Total 1124 genes expressed significant differences that up-regulated and down-regulated genes were 603 and 521, respectively ( Figure 5).
The top 10 down-regulated DEGs and top 10 up-regulated DEGs matched with Oryzabase gene symbols and Oryzabase gene name synonym(s) with an adjusted p-value < 0.05 (Table 1). Typically, many F-box protein genes were up-or down-regulated. A total of 396 F-box genes were included in 1124 DEGs (Supplementary Table S3). The heatmap of the 500 most expressed genes (TMM upper value) across samples showed that control samples were clustered with each other (Figure 6). A total of 1124 DEGs were selected based on a significance cutoff of FDR < 0.05.

Statistical Test for DEG
Based on quantified gene expression, relationships among samples and stages were investigated using the MDS method. Multi-dimensional plot indicated that samples could be distinguished by transgenic rice (Figure 4). After checking distances of samples, statistical tests were performed to identify DEGs between control and transgenic samples. Total 1124 genes expressed significant differences that up-regulated and down-regulated genes were 603 and 521, respectively ( Figure 5).
The top 10 down-regulated DEGs and top 10 up-regulated DEGs matched with Oryzabase gene symbols and Oryzabase gene name synonym(s) with an adjusted p-value < 0.05 (Table 1). Typically, many F-box protein genes were up-or down-regulated. A total of 396 F-box genes were included in 1124 DEGs (Supplementary Table S3). The heatmap of the 500 most expressed genes (TMM upper value) across samples showed that control samples were clustered with each other (Figure 6). A total of 1124 DEGs were selected based on a significance cutoff of FDR < 0.05. Genes 2021, 12, x FOR PEER REVIEW 8 of 15

GO and KEGG PATHWAY Analysis
Using the list of 1124 DEGs, functional enrichment of DEGs such as GO terms and the KEGG pathway were investigated ( Table 2). In KEGG pathways of DEGs, metabolism of starch and sucrose was significantly enriched (Figure 7). Trehalose-6-phosphate (T6P), an intermediate of trehalose biosynthesis, is an essential signal metabolite in plants. It is linked to growth and development according to carbon status. In this study, the trehalose-6phosphate synthase (TPS) gene for transferring Gluc-6P to T6P was down-regulated, whereas the trehalose-6-phosphate phosphatase (TPP) gene for transferring T6P to trehalose was upregulated (Table 3). After TPS converts glucose-6-phosphate and UDP-glucose into T6P, T6P is then dephosphorylated into trehalose by TPP.    Table 2. Functional classification of differentially expressed proteins in rice samples. Proteins are classified into three main categories by GO analysis: biological process, cellular component, and molecular function. Column headers indicate the number of genes in a category, the percentage of a specific category of genes in that category, and the name of genes in that category. KEGG analysis showed one metabolism category. Results were filtered with a p-value < 0.05 as a cutoff.  We identified the metabolism pathway using DAVID.
Expression levels of OsTPS1 and OsTPS5 genes were decreased in PsGAPDH-OX plants compared to those in the wild-type. Furthermore, the transcript levels of the TPP genes were significantly induced ( Figure 8). Results of qRT-PCR correlated well with transcription levels estimated from RNA sequencing data, thus supporting the involvement of these genes in the defense against salinity stress in transgenic rice plants. Additionally, these results indicated that low levels of T6P and high levels of trehalose in PsGAPDH-OX rice plants might play a role in the resistance of transgenic plants against abiotic stress. PsGAPDH-OX plants compared to those in the wild-type. Furthermore, the transcript levels of the TPP genes were significantly induced ( Figure 8). Results of qRT-PCR correlated well with transcription levels estimated from RNA sequencing data, thus supporting the involvement of these genes in the defense against salinity stress in transgenic rice plants.
Additionally, these results indicated that low levels of T6P and high levels of trehalose in PsGAPDH-OX rice plants might play a role in the resistance of transgenic plants against abiotic stress. Figure 8. Expression levels of starch and sucrose-associated genes in wild-type and PsGAPDH-OX plants under NaCl stress conditions. Quantitative real-time PCR data were analyzed using the 2 −∆∆CT method with OsActin1 gene as an internal control. Values are presented as means ± SD of three independent measurements. Different lowercase letters indicate significant differences (ANOVA with Tukey's HSD, p < 0.05).

Discussion
GAPDH genes are commonly used as internal controls for the relative quantitation of gene expression [40]. However, a recent study has revealed that GAPDH displays unacceptably high expression variability, thus limiting its use as an internal control [41]. Indeed, several GAPDH genes have been characterized in relation to abiotic stresses. In particular, transgenic rice plants overexpressing OsGAPC3 show significantly increased tolerance to salt stress, both during and after germination [40]. In the present study, we found that PsGAPDH-OX rice plants were more adaptive to salt stress than wild-type ones at germination and seedling growth stages (Figures 2 and 3). Especially, PsGAPDH-OX #6 was the most tolerant against salt stress. This result suggests that GAPDH activity is required for the acclimation of rice plants to salt stress in the early growth stage. GAPDH is a highly conserved glycolytic enzyme that plays an important role in carbon metabolism Figure 8. Expression levels of starch and sucrose-associated genes in wild-type and PsGAPDH-OX plants under NaCl stress conditions. Quantitative real-time PCR data were analyzed using the 2 −∆∆CT method with OsActin1 gene as an internal control. Values are presented as means ± SD of three independent measurements. Different lowercase letters indicate significant differences (ANOVA with Tukey's HSD, p < 0.05).

Discussion
GAPDH genes are commonly used as internal controls for the relative quantitation of gene expression [40]. However, a recent study has revealed that GAPDH displays unacceptably high expression variability, thus limiting its use as an internal control [41]. Indeed, several GAPDH genes have been characterized in relation to abiotic stresses. In particular, transgenic rice plants overexpressing OsGAPC3 show significantly increased tolerance to salt stress, both during and after germination [40]. In the present study, we found that PsGAPDH-OX rice plants were more adaptive to salt stress than wild-type ones at germination and seedling growth stages (Figures 2 and 3). Especially, PsGAPDH-OX #6 was the most tolerant against salt stress. This result suggests that GAPDH activity is required for the acclimation of rice plants to salt stress in the early growth stage. GAPDH is a highly conserved glycolytic enzyme that plays an important role in carbon metabolism pathways. Glycolysis is a principal metabolic pathway that provides precursors for biosynthetic reactions and takes part in the stress adaptation of plants [12][13][14]. Recent metabolomics analyses have demonstrated that long-term salt stress can alter cellular metabolic processes, including glycolysis, gluconeogenesis, and sucrose metabolism [13,40]. Plants can accumulate osmotic substances such as proline, trehalose, and sucrose to adapt to salt stress, because they can induce metabolic processes associated with osmolyte production [13,40].
Our RNA sequencing analysis data indicated that overexpression of PsGAPDH caused transcriptomic changes for a considerable number of genes involved in various physiological processes of rice. Those genes found to be responsive to salt stress in rice plants overexpressing PsGAPDH in the present study have previously been reported to be involved in trehalose biosynthesis processes such as trehalose-6-phosphate synthase (TPS) and trehalose-6-phosphate phosphatase (TPP). Among various osmolytes, trehalose (α-D-glucopyranosyl-α-D-glucopyranoside), a non-reducing disaccharide, is a principal compatible solute and a potential signal metabolite in plants. In plants, there is only one pathway for trehalose biosynthesis. It consists of a two-step process involving TPS and TPP that synthesize and subsequently dephosphorylate trehalose-6-phosphate (T6P) as an intermediate. Previous studies have shown that genetic modification of the TPS/TPP gene can increase stress resistance in different species [42,43]. Heterologous expression of bacterial or yeast trehalose biosynthesis gene in tobacco, Arabidopsis, rice, and potato can increase their stress tolerance [44]. Overexpression of OsTPP1 can improve the resistance to salt and cold. It can also trigger the expression of a series of abiotic response genes, thus increasing the stress tolerance of rice [45]. Previous studies have shown that some TPS/TPP genes are directly involved in stress tolerance by improving trehalose contents in several plants [46].

Conclusions
Our results suggest that the PsGAPDH gene might participate in salt tolerance through coordination with these responsive genes. Hence, we speculate that PsGAPDH is possibly involved in various cellular processes in addition to the metabolism of starch and sucrose. Taken together, our findings suggest that PsGAPDH plays an important role in the metabolism of starch and sucrose during the adaptation of plants to salt stress.

Conflicts of Interest:
The authors declare no conflict of interest.