*Article* **Combination of Genomics, Transcriptomics Identifies Candidate Loci Related to Cold Tolerance in Dongxiang Wild Rice**

**Dianwen Wang †, Yulong Xiao †, Hongping Chen, Cheng Huang, Ping Chen, Dazhou Chen, Wei Deng and Jilin Wang \***

> Rice National Engineering Research Center (Nanchang), Rice Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang 330200, China

**\*** Correspondence: wangjilin1982@163.com; Tel.: +86-133-8753-2293

† These authors contributed equally to this study.

**Abstract:** Rice, a cold-sensitive crop, is a staple food for more than 50% of the world's population. Low temperature severely compromises the growth of rice and challenges China's food safety. Dongxiang wild rice (DXWR) is the most northerly common wild rice in China and has strong cold tolerance, but the genetic basis of its cold tolerance is still unclear. Here, we report quantitative trait loci (QTLs) analysis for seedling cold tolerance (SCT) using a high-density single nucleotide polymorphism linkage map in the backcross recombinant inbred lines that were derived from a cross of DXWR, and an indica cultivar, GZX49. A total of 10 putative QTLs were identified for SCT under 4 ◦C cold treatment, each explaining 2.0–6.8% of the phenotypic variation in this population. Furthermore, transcriptome sequencing of DXWR seedlings before and after cold treatment was performed, and 898 and 3413 differentially expressed genes (DEGs) relative to 0 h in cold-tolerant for 4 h and 12 h were identified, respectively. Gene ontology and Kyoto encyclopedia of genes and genomes (KEGG) analysis were performed on these DEGs. Using transcriptome data and genetic linkage analysis, combined with qRT-PCR, sequence comparison, and bioinformatics, *LOC\_Os08g04840* was putatively identified as a candidate gene for the major effect locus *qSCT8*. These findings provided insights into the genetic basis of SCT for the improvement of cold stress potential in rice breeding programs.

**Keywords:** Dongxiang wild rice; seedling cold tolerance; quantitative trait locus; transcriptomics; differentially expressed genes

#### **1. Introduction**

Rice (*Oryza sativa* L.) is the staple food for half of the world's population. Increased rice yield has become an important issue for the global economy and food security due to rapid population growth and a reduction in arable land. The growth and development of rice are sensitive to temperature fluctuations, and the optimum temperature for rice cultivation is 25–30 ◦C [1]. Double-cropping early rice often encounters cold weather in late spring, and suffers from chilling stress during its seedling stage, causing chlorosis, reduction of growth rate and tillering, and low seedling vigor [1–4]. Thus, improving cold tolerance (CT) in cultivars to promote high and stable rice yield has been a major goal of rice breeding [5].

CT in rice is a quantitative trait that is controlled by quantitative trait loci (QTL) and is largely influenced by the environment. Numerous loci were studied for cold tolerance by QTL mapping and association analysis in rice, including more than 80 QTLs related to seedling cold tolerance (SCT) [4,6–16]. Among these loci, only a few QTLs have been thoroughly researched and cloned using map-based cloning strategies, and most of the functional mechanisms remain largely unknown [16]. *COLD1* encodes a Ca2+ signaling regulator that interacts with *qCTS9*, encoding a novel Brassinosteroid Insensitive-1 protein,

**Citation:** Wang, D.; Xiao, Y.; Chen, H.; Huang, C.; Chen, P.; Chen, D.; Deng, W.; Wang, J. Combination of Genomics, Transcriptomics Identifies Candidate Loci Related to Cold Tolerance in Dongxiang Wild Rice. *Plants* **2022**, *11*, 2329. https:// doi.org/10.3390/plants11182329

Academic Editors: Lixi Jiang, Mingxun Chen and Yuan Guo

Received: 27 July 2022 Accepted: 2 September 2022 Published: 6 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and co-regulates cold tolerance in rice [17,18]. Functional interaction between *bZIP73* and *bZIP71* can significantly enhance cold tolerance in rice seedlings [19]. The substitution of a single nucleotide at position 343 from A to G in *qPSR10* can effectively enhance cold tolerance in the rice seedling stage [20]. *HAN1* encodes a biologically active jasmonyl-Lisoleucine (JA-Ile) to the inactive form 12-hydroxy-JA-Ile (12OH-JA-Ile) oxidase, which is mediated by fine-tuning JA cold tolerance of rice [21]. *OsWRKY115* (*qCT7*) transcription factor positively regulates the cold tolerance of rice seedlings [22]. In this context, a large effort is still required to explore multiple beneficial alleles for the improvement of cold tolerance in rice.

The development of DNA sequencing technology has supported the possibility of identifying functional genes at the genome-wide level [23–25]. RNA-sequencing (RNA-seq), characterizing the genome-wide gene expression profile, is an efficient means to detect stimuli-responsive genes at the genome-wide level [26]. Transcriptome analysis of seedlings before and after cold treatment showed that there are many genes involved in response to cold stress at the seedling stage, such as transport, metabolism, signal transduction, and transcriptional regulation [27–29]. In particular, comparative transcriptome analysis of germplasms with distinct cold tolerance phenotypes provides important clues for revealing key genes and mechanisms controlling cold tolerance in rice [30–32]. Moreover, the strategy combining transcriptome analysis with QTL mapping was explored to accelerate the identification of candidate genes, including *LOC\_Os07g22494*, for a CT-related QTL [33–35]. Jiangxi Dongxiang wild rice (DXWR) is the northernmost wild rice in China (28◦14 N), and its seedlings are extremely cold-tolerant [36]. Several QTLs and candidate genes associated with cold tolerance have been identified in DXWR [14,37,38]. However, the genetic basis for the regulation of cold tolerance in DXWR remains unclear.

In this present study, a high-density single nucleotide polymorphism (SNP) linkage map was constructed in the backcross recombinant inbred line (BRIL) population, which was derived from a backcross of DXWR and GZX49 to detect QTLs associated with SCT. At the same time, we cold-treated DXWR seedlings and performed transcriptome sequencing to analyze differentially expressed genes. For a novel major effect, *qSCT8* candidate gene analysis was performed, and *LOC\_Os08g04840* encoding an MYB transcription factor was used as the *qSCT8* candidate gene. These findings provide new insights into the genetic basis of SCT in rice.

#### **2. Results**

#### *2.1. Development and Construction of High-Density Bin Map of BRILs*

To construct the backcross recombinant inbred lines (BRILs), the F1 plant derived from a cross between Ganzaoxian49 (GZX49) and Dongxiang wild rice (DXWR) was backcrossed as a receptor with GZX49. The BC1F1 seeds were harvested and single-seed descent to obtain a BRIL population comprising 132 lines (Figure 1A). The population was genotyped using the genotyping-by-sequencing (GBS) method, and generated a total of 70.76 Gb of raw data, approximately 536.09 Mb each line, with an overall effective mean depth coverage of 1.97-fold. Based on the Rice Genome Annotation Project (version 7), a total of 55,036 SNPs were detected in the population and were evenly distributed throughout the genome (Figure 1B, Table S1). According to the high-density genotype of the BRIL population, the recombination breakpoint of each line is determined. In addition, a bin map was constructed based on the recombination breakpoints of each line. The bin map contains 1658 recombination bins distributed on 12 chromosomes, and the average size of the bin is 220.4 kb, with a size range of 12.8 kb to 4.97 Mb. The genetic linkage map was further constructed using a bin map, the total length of the linkage map was 2341.9 cM, and the average interval between adjacent bins was 1.44 cM (Figure 1C, Table S1).

**Figure 1.** Development and construction of high-density bin map of BRILs. (**A**) Flowchart of developing backcross recombinant inbred line (BRIL) population. ⊗ represents self-cross. (**B**) Polymorphic SNPs between DXWR and GZX49 distributed on chromosomes. (**C**) Genetic linkage map constructed by population BRILs.

#### *2.2. QTL analysis of Seedling Cold Tolerance*

Seedling cold tolerance (SCT) was measured for the two parental lines and the BRILs at the seedling stage. The seedlings grew to trifoliate (Figure 2A), were treated at a low temperature of 4 ◦C for 84 h in the incubator, and then resumed growth at 28 ◦C for 7 days; the results showed that the survival rate of DXWR was still 100% after cold treatment, while the survival rate of GZX49 was almost zero (Figure 2B). This indicated that DXWR had strong SCT compared with GZX49. The survival rate of seedlings after cold treatment in the BRIL population showed a continuous distribution from 0 to 100% (Figure 2C,D). These results suggested SCT was in quantitative inheritance controlled by polygenes or QTLs, and there was a large genetic variation for SCT among the BRILs.

**Figure 2.** Phenotypic identification of seedling for cold tolerance. (**A**,**B**) represent the phenotypes of GZX49 and DXWR before and after cold treatment, respectively. (**C**) Phenotypes of some lines of the BRIL population after cold treatment. The scale bar of (**A**–**C**) means 5 cm. (**D**) Frequency distribution of seedling survival rate after cold treatment. Arrows indicate the means of parental lines GZX49 and DXWR.

The ridge regression analysis for the QTL detection was performed in the BRIL with the bin genotypes. A total of 10 QTLs for SCT were identified and found to explain 48.9% of the phenotypic variance, distributed on all chromosomes except chromosomes 5, 6, and 10 (Figure 3A, Table S2). The majority (9/10) of the QTLs for SCT had a positive effect, indicating that alleles from DXWR increased SCT (Table S2). The phenotypic variation explained (PVE) by each QTL ranged from 2.0% to 6.8% (Table S2), confirming that the SCT was a complex quantitative trait controlled by multiple genes. Among these QTLs, 8 QTLs correspond to bin intervals less than 500 kb, and 5 QTLs correspond to bin intervals that were smaller (less than 200 kb) (Table S2). The detection of QTLs in relatively small intervals was of great significance for further fine mapping and candidate gene mapping of these QTLs. Based on the mapping interval of two QTLs, *qSCT1* and *qSCT2*, it was found that these two QTLs contained the two reported SCT-related genes *OsMYB3R-2* [39] and *OsTPP1* [40], respectively. Using qRT-PCR technology to detect the expression of *OsMYB3R-2* under 4 ◦C of treatment, we found that the expression of *OsMYB3R-2* was significantly up-regulated in DXWR but not significantly changed in GZX49 (Figure 3B). Sequence comparison of *OsMYB3R-2* between DXWR and GZX49 revealed in its promoter region (2-kb upstream of the predicted start codon), 17 SNPs and 4 Indels (Figure 3C). Thus, *OsMYB3R-2* was a possible candidate gene for *qSCT1*. Similarly, qRT-PCR technology and sequence comparison were conducted to screen out the degree of response of *OsTPP1* expression to cold stress and the causal polymorphisms of *OsTPP1* (Figure 3D,E). Thus, *OsTPP1* is a possible candidate gene for *qSCT2*. Additionally, these results indicate that the QTLs are entirely effective.

**Figure 3.** QTL analysis of the BRIL population. (**A**) Manhattan plots of the loci for seedling cold tolerance (SCT). The *x*-axis represents single nucleotide polymorphism (SNP) along each numbered chromosome; the *y*-axis represents the negative logarithm of the *p*-value (-log10 p) for the SNP association. Horizontal dashed lines in the plots indicate the declaration thresholds. (**B**,**D**) Expression levels of the *OsMYB3R-2* and *OsTPP1* in GZX49 and DXWR after cold stress measured by qRT-PCR, respectively. The results were statistically analyzed using Student's *t*-test (\*\* *p* < 0.01, \*\*\* *p* < 0.005). Transcription levels relative to 0 h, which was set to 1, are presented as the mean and SE of triplicates. *LOC\_Os03g13170* (Ubiquitin) is the control gene. (**C**,**E**) Sequence comparison of *OsMYB3R-2* and *OsTPP1* among GZX49 and DXWR. The vertical bars and triangle represent SNPs and nucleotide deletion, respectively.

#### *2.3. Transcriptome Analysis of Differentially Expressed Genes in DXWR under Cold Stress*

To identify the cold stress-responsive genes contained in DXWR, DXWR seedlings were treated at 4 ◦C, and transcriptome sequencing of aboveground tissues treated for 0 h, 4 h, and 12 h, respectively, was performed. There were 898 and 3413 differentially expressed genes (DEGs) relative to 0 h at 4 ◦C for 4 h and 12 h, respectively (Figure 4A). Among them, 291 and 607 genes were up- and down-regulated at 4 h, respectively, while

1294 and 2119 genes were up- and down-regulated at 12 h, respectively, by 4 ◦C treatment (Figure 4A). Additionally, there were 1047 DEGs between 4 h and 12 h of treatment at 4 ◦C, including 610 up-regulated genes and 437 down-regulated genes (Figure 4A). Among the DEGs at 4 h, 93.4% (273/291) of the up-regulated genes and 86.7% (526/607) of the downregulated genes were repeatedly detected in the DEGs of 12 h, respectively (Figure 4B,C). In addition, among the differentially expressed genes at 12 h, 79% (1021/1293) of up-regulated genes and 75.2% (1593/2119) of down-regulated genes had no detectable differences at 4 h (Figure 4B,C). This indicated that the response mechanism of DXWR to cold stress was mainly at 12 h of cold treatment.

**Figure 4.** Transcriptome analysis of the genetic mechanism of DXWR in response to cold stress. (**A**) Differentially expressed gene statistics. (**B**,**C**) are Venn diagrams of up-regulated and downregulated genes, respectively. (**D**,**E**) are gene ontology (GO) functional enrichment histogram and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment histogram, respectively. 0\_4, 0\_12, and 4\_12 represent the comparison of DXWR seedling treatment at 4 ◦C for 0 h and 4 h, 0 h and 12 h, and 4 h and 12 h, respectively. PGI, FAE, CSWB, SSM, ASNSM, GM, GPM, LAM, ALAM, NM, PB, FB, PST, PPER, and PPI in subfigure (**E**) represent pentose and glucuronate interconversions, fatty acid elongation, cutin suberine and wax biosynthesis, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycerolipid metabolism, gycerophospholipid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism, nitrogen metabolism, phenylpropanoid biosynthesis, flavonoid biosynthesis, plant hormone signal transduction, protein processing in endoplasmic reticulum, and plant–pathogen interaction, respectively.

Furthermore, the DEGs were classified according to gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) terms. DEGs at 4 h and 12 h were enriched in 14 GO terms, mainly containing mitochondrion, cytosol, chloroplast, plastid, and DNA binding (Figure 4D). KEGG analysis showed that DEGs at 4 h and 12 h displayed enrichment in diverse metabolic pathways, including photosynthesis metabolism, phenylalanine metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and glycerolipid metabolism (Figure 4E). In addition, the DEGs were classified as plant–pathogen interaction, plant hormone signal transduction, protein processing in the endoplasmic reticulum, and MAPK signaling pathway (Figure 4E).

#### *2.4. Candidate Gene Analysis for qSCT8*

To identify the possible novel genes for STC within the QTLs, *qSCT8* was thoroughly analyzed because its *p*-value is the most significant. *qSCT8* was first identified within a 638 kb interval. Using Rice Genome Annotation Project (uga.edu), it was found that the interval of *qSCT8* contains about 90 genes. According to the transcriptome sequencing data of DXWR, it was found that among the 90 genes, 18 genes were detected to be expressed at the seedling stage, of which only *LOC\_Os08g04840*, encodes an MYB family transcription factor, was up-regulated after cold treatment (Figure 5A). The expression of *LOC\_Os08g04840* during 4 ◦C cold stress in DXWR and GZX49 was further analyzed by qRT-PCR technology, and it was found that the expression of *LOC\_Os08g04840* did not change significantly in DXWR and GZX49 after 4 h of treatment compared with 0 h, but it was significantly up-regulated by 2-fold after 12 h of treatment. It was up-regulated by 2-fold, and finally, DXWR was up-regulated by 3.5-fold, and GZX49 was only 2.3-fold after 24 h of treatment (Figure 5B). This indicated that the response intensity of *LOC\_Os08g04840* in DXWR was higher than that of GZX49 in the late 4 ◦C cold treatment. Furthermore, sequence comparison revealed that *LOC\_Os08g04840* had substantial variation among rice varieties. In its promoter region (2-kb upstream of the predicted start codon), 4 SNPs were found between DXWR and GZX49. In the coding region, there were no nucleotide variations. Thus, *LOC\_Os08g04840* was putatively a possible candidate gene for *qSCT8*.

**Figure 5.** *qSCT8* candidate gene analysis. (**A**) Clustering heat maps of the relative expression levels of genes within the qSCT8 localization interval determined using RNA-seq data. Standard scores (Z-scores) were used as the numerical signs to evaluate the standard deviations from the mean of the corresponding samples. (**B**) Expression levels of the *LOC\_Os08g04840* in GZX49 and DXWR after cold stress measured by qRT-PCR. The results were statistically analyzed using Student's *t*-test (\*\*\* *p* < 0.005). (**C**) Sequence comparison of *LOC\_Os08g04840* among GZX49 and DXWR. Vertical bars represent SNPs.

#### **3. Discussion**

In recent years, direct seeding of rice has gained popularity due to its time and labor savings and low input. However, the transplanting of seedlings was carried out in a greenhouse, so the genes related to cold tolerance in the seedling stage of rice were lost due to domestication, which leads to major problems such as poor chlorosis, reduced growth rate and tillering, and low seedling vigor in direct-seeded rice when it encounters a cold spring [1,2,41]. Therefore, cultivating rice varieties with strong cold tolerance at the seedling stage is a necessary condition for the application of direct seeding cultivation of rice. DXWR survives heavy snow and can naturally overwinter in its habitat (Figure 6A). The cold tolerance of DXWR is stronger than the main variety LongJing31 in northeast China (Figure 6B). However, most of the identified cold tolerance genes are from japonica rice, which limits the development of cold tolerance breeding in cultivated rice [4,16,18,42–44]. Therefore, the identification of SCT-related genes and analysis of the genetic basis of cold tolerance in DXWR has important scientific value for improving the low-temperature seedling vigor ability of rice.

**Figure 6.** Cold tolerance phenotype of DXWR. (**A**) The phenotype of DXWR encountering heavy snow weather in its habitat. (**B**) Phenotypes of DXWR, Longjing 31, and GZX49 after 4 days of treatment at 0 ◦C at the seedling stage. LJ31 stands for LongJing31. The scale bar of (**B**) means 5 cm.

In this study, a high-generation genetic population was constructed using DXWR/GZX49, and a high-density linkage map was constructed using this population, which contained 1658 bins with an average physical interval of 220.4 kb (Table S1). Using this population, 10 QTLs affecting seedling survival after 78 h of cold treatment at 4 ◦C were identified, and these QTLs explained 2.06–6.8% of the phenotypic variation range. By comparing the positions of these QTLs, 6 QTLs overlapped with the reported QTLs related to cold tolerance (Table S2). The location of *qSCT1* was very close to that of the *QTL3* [44], and the identified cold tolerance-related gene *OsMYB3R-2* was included in this mapping interval [39]. The *qSCT2* mapping region contained the identified cold tolerance gene *OsTPP1* [40]. *qSCT4* was mapped near *qCTB4-1*, and two cold tolerance-related genes, *CTB4a* and *CTB2*, have been cloned in this interval [45,46]. *qSCT7* and *qSCT9* contain intervals that overlap with the cloned cold tolerance genes *qCT7* [22] and *qCTS-9* [17], respectively. Two QTL *qSCT11.1* and *qSCT11.2* were identified on chromosome 11, which were detected as *qCTS11-2* and *qCTS11-4*, respectively, in previous reports [43]. The results of the above comparison reflected the accuracy of this study and also illustrated the complexity of genetic regulation of cold tolerance in the rice seedling stage. In addition, this study showed that DXWR was extremely cold tolerant due to its inclusion of multiple cold-tolerant genes. These results not only reinforced those of previous studies but also reflected the complexity of developing strong cold-tolerant rice varieties.

Many rice cold tolerance genes have been identified using reverse genetics, which is important for understanding the regulatory mechanisms of rice cold tolerance [16]. However, only cold tolerance genes with natural allelic variation can be directly used for breeding improvement. Cold stress could induce the up-regulated expression of *OsMYB3R-2*; meanwhile, overexpression of *OsMYB3R-2* in *Oryza sativa japonica* 'Zhonghua 10' could significantly enhance the cold tolerance of its seedlings [39]. Low-temperature stress could activate the activity of *OsMAPK3*, and the active *OsMAPK3* could phosphorylate *OsbHLH002* and make it accumulate without ubiquitinase degradation. *OsbHLH002* can promote the expression of *OsTPP1*, thereby increasing trehalose content and resistance to low-temperature damage [40]. Our study found that the *OsMYB3R-2DXWR* and *OSTPP1DXWR* allele was significantly more up-regulated than *OsMYB3R-2GZX49* and *OSTPP1GZX49* under cold stress conditions (Figure 3B,D), indicating that *OsMYB3R-2DXWR* and *OSTPP1DXWR* may be a natural allelic variation with strong cold tolerance, which can be directly used for breeding improvement. This further indicates that DXWR contains many excellent cold tolerance alleles, which can provide important genetic resources for the improvement of rice cold tolerance breeding.

In order to further analyze the cold tolerance characteristics of Dongxiang wild rice at the level of gene expression, transcriptomes were performed on DXWR during cold treatment in this study, and more than 3000 cold stress response genes were identified (Figure 4A–C). GO terms and KEGG pathways analysis were performed on these genes, and it was found that these genes have been involved in multiple biological processes (Figure 4D,E). This further illustrated the complexity of DXWR resistance to cold stress.

Combined QTL-mapping and transcriptomes analysis, compared to the approaches for identifying candidate genes using only traditional QTL-mapping or high-throughput expression profiling, takes less time, reduces labor costs, and increases the selection veracity of the target regions or candidate genes [35,47]. In this study, a novel cold-tolerant QTL, *qSCT8*, was first identified within a 638 kb interval (Table S2 and Figure 3A). The genes *LOC\_Os08g04840*, coding an MYB family transcription, that responded to cold stress in the *qSCT8* mapping region were screened out by transcriptome, and qRT-PCR was used to further verify *LOC\_Os08g04840* which showed differences in expression before and after cold treatment (Figure 5A,B). In particular, the two G/A variant at the -350 and -354 site upstream of the start codon may cause MADS combined element (AAAAAAAAAGAAAG) defects in *LOC\_Os08g04840*GZX49 (Figure 5C) [48]. Transcription factors of MYB and MADS family have been reported to regulate cold tolerance [39,49–52]. The interaction mode of MADS—MYB might play an important role in the signal transmission of cold tolerance in DXWR. This presumption might provide a new idea for the follow-up study on the genetic mechanism of cold tolerance in the DXWR seedling stage.

#### **4. Materials and Methods**

#### *4.1. Plant Materials*

A backcross recombinant inbred line (BRIL) population was developed by single-seed descent from a backcross (BC1F1) of Dongxiang wild rice (DXWR) as donor and *indica* cultivar GZX49 as the recurrent parent. This population consisted of 132 lines, which were backcrossed for 1 generation with GZX49 as the parent, and then selfed for 8 generations (Figure 1A). The BRIL population was grown at the experimental field of Jiangxi Academy of Agricultural Sciences in 2019 at Nanchang (28.57 N, 115.9 E), China.

#### *4.2. Evaluation of Cold Tolerance at the Seedling Stage*

The harvested BRIL seeds were incubated at 40 ◦C for approximately 36 h to break dormancy and soaked in deionized water at 30 ◦C for approximately 60 h for germination. A total of 30 germinated seeds from each line were selected and sown in soil in pots and cultivated under a 12 h light/12 h dark cycle at 28 ◦C/26 ◦C with 80% humidity. When most of the seedlings grew to trifoliate, weak seedlings were removed, and the rest were treated at a low temperature of 4 ◦C for 84 h in the incubator and then resumed growth at 28 ◦C for 7 days. Survival rates were determined after 14-day treatments by counting the surviving plants (leaf contains about 40% green part) and the dead plants. All lines were subjected to three independent replicates.

#### *4.3. DNA Extraction and SNP Genotyping*

Genomic DNA extraction was carried out using the plant genomic DNA extraction Kit (TIANGEN, Beijing, China), following the manufacturer's instructions. RNase A was then added to digest RNA. The quality and concentrations of DNA were detected using a NanoDrop 2000 (Thermo Fisher Scientifc, Waltham, MA, USA), while DNA integrity was examined by electrophoresis on 1% agarose gels. To prepare the reduced representation libraries for sequencing, the GBS protocol was carried out according to the method reported by Elshire et al. [53]. In brief, the genomic DNA was first digested by restriction enzymes. In this case, *EcoRI* and *MseI* were selected to effectively reduce genome complexity. Barcode adapters were designed and modified according to the standard Illumina adapter design for paired-end read libraries. The ligation reaction was incubated for 1 h at 22 ◦C with T4 DNA ligase (Thermo Scientific, Madison, WI, USA) and inactivated at 65 ◦C for 20 min. The ligation products from each sample were pooled in a single tube, and the products were amplified with 10 cycles of PCR. The amplified library was purified using a QIA quick PCR purification kit (Qiagen, Hilden, Germany), quantified on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), and finally sequenced on an Illumina HiSeq3000 instrument (Illumina, San Diego, CA, USA), which generated 150 bp paired-end reads.

The sequencing reads were aligned to the Nipponbare genome sequence (IRGSPv7) using SOAP2 software (version 2.20) [54]. SNP calling was performed with realSFS based on the Bayesian estimation of site frequency at every site. All SNPs were filtered using a Practical Extraction and Report Language (PERL) script based on the following criteria: loci with >50% missing data and minor allele frequency less than 5% in the population. In the end, we obtained 55,036 SNPs. Based on SNP genotyping, the bin was defined by a unique overlapping recombination segment across the BRILs according to a previously reported approach [55]. A bin without breakpoints was generated using the R/qtl package function *fill.geno* with the "argmax" method. The high-density bin map was constructed as previously described [56]. Briefly, the sliding window approach was adopted to evaluate a group of consecutive SNPs for genotyping and determination of recombination breakpoints along the chromosomes of each individual. Blocks with a length less than 250 kb in which the number of sequenced SNPs was fewer than 5 were masked as missing data to avoid false double recombination. Genotypes of bins for regions at the transitions between two different genotype blocks were imputed using the R/qtl package. The genetic linkage map was constructed using the R/qtl package function *est.map* with the Haldane map method [57]. Finally, we used 48,339 SNPs to construct bin maps with 1658 bins. The genotype data of the BRIL population are available on request.

#### *4.4. QTL Analysis*

The QTL analysis of the phenotypic data with bin maps with 1658 bins in the BILs was performed using the linear ridge regression method to reduce the multicollinearity among markers, as described previously [58]. A significance level of *p* < 0.005 was set as the threshold in the BILs to declare the presence of a putative QTL in a given bin. If several adjacent bins showed *p*-values lower than the threshold, the QTL was tentatively located in the bin (peak bin) with the lowest *p-*value [58]. The phenotypic variance explained by each QTL was decomposed using the "relaimpo" package of R ("lmg" function). QTL nomenclature followed the principles suggested by a previous report [59].

Gene annotations for a given peak bin were obtained from the Rice Genome Annotation Project Database (http://rice.plantbiology.msu.edu/, accessed on 8 June 2021).

#### *4.5. RNA-Seq Analysis*

Six RNA samples, including cold treatment and controls, were used for RNA sequencing. The aerial parts of 20 seedlings of the corresponding treatment were mixed for each sample in duplicate. In total, 3 μg of RNA per sample was used as the input material for the RNA sample preparations. The sequencing libraries were prepared with the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer's

recommendations; index codes were added to attribute sequences to each sample. Clustering of the index-coded samples was performed on a cBot Cluster Generation System with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the libraries were sequenced with an Illumina sequencing platform (Hiseq 2500), and 125-bp or 150-bp paired-end reads were generated. The reads were aligned to the reference transcript sequence (http://rice.plantbiology.msu.edu/, accessed on 8 December 2021), and the number of reads covered from the start to the end of each gene was counted. We used the RSEM v1.2.31 tool to quantify gene expression levels. Due to the influence of sequencing depth and gene length, the gene expression value of RNA-seq is generally not represented by read counts, and TPM (Transcripts Per Million mapped reads) is often used as a standardized value. TPM has successively corrected the gene length and sequencing depth. The sum of the TPM of each sample is the same, and the TPM value can reflect the ratio of reads of a gene in the comparison so that the value can be directly compared between samples. Difference analysis was used to find the reasons for the differences in different samples and the degree of their influence on the differences. Statistical analysis of sample expression data was performed to screen genes with significantly different expression levels in different states. The difference analysis is divided into three steps: first, normalize the original read counts, mainly to correct the sequencing depth; then, the statistical model is used to calculate the hypothesis test probability (*p*-value); finally, the multiple hypothesis test correction is performed. Obtain the false discovery rate (FDR) value. For samples with biological replicates, DESeq2 v1.10.1 was selected for differential gene expression analysis. The above RNA-seq-related experiments were completed by Wuhan Genoseq Bioinformatics Technology Co., Ltd. The original RNA-seq data set has been deposited in NCBI (PRJNA871989).

#### *4.6. Quantitative Real-Time PCR*

Total RNAs were isolated with the TRIzol kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The RNA was treated with DNase I (Invitrogen), and approximately 3 μg of total RNA was used to synthesize first-strand cDNA using oligo (dT)18 as a primer (Promega, Shanghai, China). The trans-intron ACTIN primer (ACTIN-M) was used to detect whether the reverse transcribed cDNA still has genomic DNA, and the cDNA without genomic DNA was used for subsequent Quantitative real-time PCR. Quantitative real-time (qRT) PCR was performed using gene-specific primers and the FastStart Universal SYBR Green Master (Roche) on a real-time PCR ViiA7 system (Applied Biosystems). Genes *ubiquitin* (*LOC\_Os03g13170*) and *actin1* (*LOC\_Os03g50885*) not differentially expressed in RNAseq were used as the internal control. The relative quantification method (2−CT) was used to evaluate gene expression level [60]. As similar expression results were observed regardless of which control genes were used, the *ubiquitin* of expressions was applied for the relative expression analyses for every assayed sample. At least three biological replicates, each containing four technical, were performed for each experiment.

The primers were designed according to the Nipponbare reference genome by Primer3 (http://redb.ncpgr.cn/modules/redbtools/primer3.php, accessed on 8 April 2022). The sequences were analyzed using Sequencer 5.0 (Gene Codes Corporation, Ann Arbor, MI, USA). All primers were synthesized at Sangon Biotech (Shanghai, China) and are listed in Table S3.

#### **5. Conclusions**

In conclusion, Dongxiang wild rice (DXWR) has strong cold tolerance. A BRIL population of DXWR/GZX49 with a high-density bin map was used to identify several QTLs and potential candidate genes for seedling cold tolerance in this study. Cold stress-responsive genes in DXWR were detected by transcriptome sequencing. Moreover, using transcriptome data and genetic linkage analysis, combined with qRT-PCR, sequence comparison, and bioinformatics, the candidate gene of the major effect locus *qSCT8* was identified as

*LOC\_Os08g04840*. These findings might provide insights into the genetic basis of SCT for the improvement of cold stress potential in rice breeding programs.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11182329/s1, Table S1: Distribution of SNP and bin on chromosome; Table S2: QTLs detected for seedling cold tolerance in BRIL population; Table S3: Primers used in the study.

**Author Contributions:** D.W. and J.W. designed and conceived the research; D.W., Y.X., J.W., C.H., P.C., W.D. and H.C. conducted the experiments; J.W., D.C. and H.C. developed the population; D.W., Y.X. and J.W. analyzed data and wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by the National Natural Science Foundation of China (31960076), the Jiangxi Modern Agricultural Scientific Research Collaborative Innovation Special Project (JXXTCXB-SJJ202204), the Jiangxi Technological Innovation Guidance Program (20212BDH81023), the Project of Discovery of Favorable Genes of Wild Rice and Breeding of Green and Efficient Varieties of Jiangxi Province (20213AAF01001) and the Doctoral Startup Fund Project of Jiangxi Academy of Agricultural Sciences (20191CBS002).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article and Supplementary Material.

**Acknowledgments:** We thank Xiaohua Tu and Xianmei Wang for assistance with phenotypic evaluation. We are grateful to Chaopu Zhang and Kun Xie for assistance in the data analysis.

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

#### **References**


**Ling-Hua Chen 1, Zu-Xin Cheng 2, Ming Xu 2, Zhi-Jian Yang 2,3 and Lin-Tong Yang 4,\***

	- College of Agriculture, Fujian Agriculture and Forestry University, Fuzhou 350002, China

**Abstract:** Organic acids metabolism and nitrogen (N) metabolism in rice seedlings and the relationship between them are not fully understood. In this study, rice (*Oryza sativa* L. ssp. Indica) variety "Huanghuazhan" was used as the experimental material, and three N levels (5 mM, 1 mM, and 0 mM NH4NO3) were set by the hydroponic method for different levels of N treatment. Our results showed that the increased content of malate in rice leaves caused by reducing N level was related to the increased synthesis of malate (the activity of leaf PEPC increased)and the decreased degradation of malate (the activity of leaf NADP-ME decreased), while the increased contents of citrate and isocitrate in rice leaves caused by reducing N level might not be caused by the increased biosynthesis, but due to the decrease in degradation of citrate and isocitrate (the activities of leaf CS, ACO, and NADP-IDH decreased). The increased content of malate in rice roots caused by reducing N level might be related to the increased biosynthesis and the decreased degradation of root malate (the activities of root NAD-MDH and PEPC increased, while the activity of NADP-ME decreased). Compared to the control (5 mM NH4NO3), the increased content of citrate in rice roots caused by reducing N level might be related to the increased biosynthesis rather than the decreased degradation of citrate, due to the higher activities of CS and ACO in rice roots under 0 mM N and 1mM N treatment when compared to that of the control ones. At the same time, the increased content of isocitrate in roots was related to the increased isomerization of isocitrate (the activity of root ACO increased) and the decreased degradation of isocitrate (the activity of root NADP-IDH decreased). With the reducing N level, the activities of N metabolism-related enzymes, such as nitrate reductase (NR), glutamine synthetase (GS), and glutamate synthase (GOGAT), decreased in rice leaves and roots, resulting in the decreased contents of total free amino acids (TFAAs) and soluble proteins in rice seedlings, and finally led to the growth inhibition. Our results showed that the dynamics of organic acids metabolism caused by reducing N level were different in rice leaves and roots. In conclusion, there was a close correlation between organic acids metabolism and N metabolism in rice leaves and roots under N-limited conditions; furthermore, such a correlation was more obvious in rice leaves than that of roots.

**Keywords:** rice; nitrogen deficiency; organic acid; nitrogen metabolism; growth

iations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Citation:** Chen, L.-H.; Cheng, Z.-X.; Xu, M.; Yang, Z.-J.; Yang, L.-T. Effects of Nitrogen Deficiency on the Metabolism of Organic Acids and Amino Acids in *Oryza sativa*. *Plants* **2022**, *11*, 2576. https://doi.org/ 10.3390/plants11192576

Academic Editors: Mingxun Chen, Lixi Jiang and Yuan Guo Received: 3 September 2022 Accepted: 27 September 2022 Published: 29 September 2022 **Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil-

#### **1. Introduction**

Nitrogen (N) is a macro-nutrient for plant growth and development. It is an important component of functional molecules, such as nucleic acid, amino acid, protein, chlorophyll, and some plant hormones. It is one of the main factors limiting plant growth and crop yield [1]. N is the most abundant nutrient element in plants, accounting for about 2–4% of the dry weight of plants [2]. Although some organic forms of N, such as urea, amino

acids, and amides, can also be absorbed by plants, ammonium and nitrate are the two main forms of N obtained by plant roots from soil solutions. These N forms are scarce in natural soil. When the total N absorbed by plant roots from the rhizosphere is less than the N required for plant growth, N deficiency will occur. In the absence of artificial N fertilizer input, N deficiency is one of the important factors restricting the crop yield and quality of agricultural products. Organic acids are the intermediate products of major carbon metabolism in plant cells and participate in various biochemical pathways, such as glycolysis, the tricarboxylic acid (TCA) cycle, photorespiration, the glyoxylic acid cycle, or the photosynthetic C4 cycle [3]. The intermediate products of organic acid metabolism, such as oxaloacetic acid and 2-ketoglutarate, provide the carbon skeletons for transamination and amino acid synthesis [4–6]. It is important to study the relationship between N level and organic acid metabolism and N metabolism for understanding plant N metabolism and guiding the application of N fertilization.

Application of N fertilizer to the field can regulate the flowering period of rice and increase the number of panicles, the number of grains per panicle, and the rice yield [7], while N deficiency can lead to the degradation of photosynthetic pigments in crop leaves, the yellowing of leaves, the reduction in carbon dioxide assimilation rate, and the inhibitionof the growth and yield of crops. The phenomenon that N deficiency reduces the photosynthetic rate and affects plant growth has been reported in many crops, including tea, citrus, corn, cucumber, sunflower, tobacco, rice, and barley [2,8–11].

The carbon (C) metabolism and N assimilation are interrelated. N assimilation requires energy and a carbon skeleton, which is provided by carbohydrate metabolism. Nitrate is an important source of inorganic N in plants. Changes in exogenous nitrate concentration could affect the rapid change in plant metabolism. For example, it could induce the synthesis of nitrate-assimilating enzymes and accelerate glycolysis and the TCA cycle to convert starch into organic acids to assimilate ammonium [12]. Excess nitrate may be accumulated into vacuoles, but excessive ammonium absorbed by crop roots and ammonium produced by nitrate through the action of nitrate reductase (NR) and nitrite reductase (NiR) have toxic effects on plants. Ammonium assimilation requires an adequate carbon skeleton and promotes carbon flow into the TCA cycle [13]. Inorganic N is assimilated and incorporated into glutamine, glutamate, asparagine, and aspartic acid, which are important N carriers in plants. These primary amino acids are the main compounds in the total free amino acid pool of many plants and can be used to synthesize other amides, such as urea, other amino acids, and amines [12,13]. In this process, the glutamine synthetase/glutamate synthase (GS/GOGAT) cycle is the main pathway of nitrogen assimilation in plants. GS catalyzes the reaction of NH4 <sup>+</sup> with glutamate to produce glutamine. GOGAT subsequently catalyzes the transfer of amide groups from glutamine to 2-ketoglutarate (2-OG), generating two glutamate molecules [14]. Among them, 2-ketoglutarate, the main carbon receptor in the GS/GOGAT pathway, is synthesized through the combined action of phosphoenolpyruvate (PEP) carboxylase (PEPC, EC 4.1.1.31), pyruvate kinase (PK), pyruvate dehydrogenase (PDH), aconitase (ACO, EC 4.2.1.3), and NADP-isocitrate dehydrogenase (NADP-IDH, EC1.1.1.42) [15].

The organic-acid-metabolism-related enzymes PEPC, NAD-malate dehydrogenase (NAD-MDH, EC 1.1.1.37), NADP-malic enzyme (NADP-ME, EC 1.1.1.40), citrate synthase (CS, EC 4.1.3.7), ACO, and NADP-IDH are considered to play an important role in balancing C and N metabolism [16,17]. PEPC catalyzes irreversible β-carboxylation of PEP in the presence of HCO3 − and Mg2 <sup>+</sup> to yield oxaloacetate, which is used to synthesize malate by NAD-MDH [18,19]. The degradation of malate is mediated by cytosolic NADP-ME in plant cells [20,21]. In the cytoplasm and mitochondrion, PEP is converted to acetyl-CoA by the sequential actions of PK and PDH. Aceytl-CoA and oxaloacetate are then combined by CS to generate citrate in the mitochondrion [22]. The citrate is reversibly isomerized to isocitrate via *cis*-aconitate by ACO [23] Isocitrate is then catalyzed to 2OG by NADP-IDH [24]. It was found on the model plant *Arabidopsis thaliana* that the interaction of C-N metabolism could significantly affect the expression levels of more than 300 genes (*p* < 0.05). Among them, the expression levels of PEPC, NAD-MDH, NADP-ME, and NADP-IDH were found to be altered by N deficiency [5]. Studies have shown that plant physiological metabolism, energy, and protein synthesis are found to be significantly affected by the interaction between C/N signals [25,26]. Specific metabolites of C/N metabolism, such as glucose, sucrose, or nitrate, serve as signals to regulate genes encoding enzymes involved in many important biological processes, including photosynthesis, carbon metabolism, nitrogen metabolism, and metabolite allocation [5]. In the study of N-deficiency transcriptome in rice, Shao et al. reported that N deficiency treatment upregulated the expression level of *PEPC* in rice leaves and increased the content of malic acid in leaves [2]. Similarly, it was also found that N-deficiency treatment increased the contents of glycolic acid and succinic acid in tea leaves and citric acid, malic acid, isocitric acid, succinic acid, and fumaric acid in tea roots [9].

Rice is an important grain crop in the world. Nitrogen fertilizer plays an important role in its growth and yield [27]. Meanwhile, increasing the soil nitrogen supply can significantly improve the soil carbon pool [28]. Previous studies have shown that under N deficiency, the activities of enzymes related to N metabolism (such as GS, GOGAT, and glutamate dehydrogenase (GDH)) were affected to varying degrees in rice [14]. Therefore, in order to further understand the mode of C and N metabolism in rice under N reduction conditions, this study carried out research on organic acid metabolism and N metabolism of rice leaves and roots under low-N and N-deficiency treatments. The results of this research can be used for rice N management, high yield, and high-quality cultivation and sustainable development of farmland.

#### **2. Results**

#### *2.1. Effects of Different Nitrogen Treatments on Rice Growth and Nitrogen Content in Leaves and Roots*

Different nitrogen levels significantly affected the growth of rice seedlings (Figure 1). Compared to the control (5 mM NH4NO3), under the treatment of 0 mM NH4NO3, rice leaves turned yellow and plant growth was inhibited. After the treatment of different N levels, the fresh weight of the plants was measured. Compared to the control, both low-N (1 mM NH4NO3) and N-deficiency (0 mM NH4NO3) treatments significantly reduced the fresh weight of rice seedlings (Figure 2A). The shoot height of rice seedlings was measured. Compared to the control, N-deficiency and low-N treatment significantly reduced the shoot height of rice seedlings (Figure 2B). After oven-drying, N content was measured. The result showed that the leaf N and root N of rice seedlings were significantly reduced by either low N (1 mMNH4NO3) or N deficiency (0 mM NH4NO3) (Figure 2C,D).

**Figure 1.** Phenotypeof *O. sativa* seedlings under different N levels. Rice seedlings were transplanted to the Hogland nutrient solution containing 5 mM (control), 1 mM (low N), and 0 mM NH4NO3 (N deficiency). Seedlings were kept under a 14 h light/10 h dark regime with white photo-illumination of 150 μmol m−<sup>2</sup> s−<sup>1</sup> as well as a relative humidity of 68% in a light incubator (28 ◦C). Fifteen days after transplanting, leaf and root samples were collected.

#### *2.2. Effects of Different N Treatments on the Contents of Organic Acids, TFAAs, and Soluble Protein in Rice Leaves and Roots*

Different N levels affected the organic acid content of rice leaves and roots. The contents of malate, citrate in leaves, and malate in roots increased significantly with decreasing N supply (Figure 3A,B,F). Compared to the control, both low N and N deficiency significantly increased the contents of leaf isocitrate (Figure 3C) and citric acid in roots (Figure 3G), while the contents of isocitrate in roots were significantly higher under N deficiency than those under low nitrogen and control treatment (Figure 3H). With the reduction in N supply, the contents of leaf TFAA (Figure 3D), leaf-soluble protein (Figure 3E), and root TFAA (Figure 3I) showed a significant downward trend, while the content of root-soluble protein had no significant difference among different N treatments (Figure 3J). All the metabolites measured above were notably lower in the root than in the shoot (Figure 3). Furthermore, it can be seen from the above results that reducing the N supply increased the contents of malate, citrate, and isocitrate in rice leaves and roots, while it decreased the contents of TFAA and leaf-soluble protein in rice, which indicated that the change trend of organic acids and TFAA content in response to reduced N availability was similar between rice leaves and roots.

#### *2.3. Effects of Different NTreatments on the Activities of Enzymes Related to Organic Acid Metabolism in Rice Leaves and Roots*

Compared to the control, the activities of NAD-MDH, NADP-ME, NADP-IDH, ACO, and CS in rice leaves were significantly decreased by both low-N and N-deficiency treatments (Figure 4A–C,E,F), but their activities were not significantly different between low-N and N-deficiency treatments. With the decrease in nutrient liquid nitrogen concentration, the PEPC enzyme activity of rice leaves showed a significant upward trend (Figure 4D). In rice roots, compared to the control, low-N and N-deficiency treatments significantly increased the enzyme activities of NAD-MDH, PEPC, and CS in the root system (Figure 4G,J,L), and significantly decreased the enzyme activities of NADP-ME and NADP-IDH in rice roots (Figure 4H,I). At the same time, there was no significant difference in NAD-MDH, PEPC, CS, NADP-ME, and NADP-IDH activities between low-N and N-deficiency treatments in rice roots (Figure 4G–J,L). With the decreased N supply, root ACO showed a significant upward trend (Figure 4K). Interestingly, all the activities of enzymes related to organic acids metabolism measured above were notably lower in the root than in the shoot (Figure 4).

#### *2.4. Effects of Different N Treatments on the Activities of Enzymes Related to N Metabolism in Rice Leaves and Roots*

Compared to the control, N-deficiency and low-N treatment significantly reduced the activities of leaf NR, leaf GS, leaf GOGAT, root NR, root GS, and root GOGAT inrice seedlings (Figure 5). In addition to no significant change in enzyme activity of leaf GS between low-N and N-deficiency treatments, low-N treatment significantly decreased the activities of the above several enzymes related to nitrogen metabolism (Figure 5B). Similar to the parameters measured above, all the activities of enzymes related to N metabolism measured above were notably lower in the root than in the shoot (Figure 5).

**Figure 4.** Effects of N deficiency on the activities of enzymes related to organic acid metabolism in *O. sativa* leaves (**A**–**F**) and roots (**G**–**L**). Bars represent means ± SD (*n* = 3). Difference among the treatments was analyzed by Duncan's multiple range test. Different letters above the bar indicate a significant difference at *p* < 0.05.

**Figure 5.** Effects of N deficiency on the activities of enzymes related to N metabolism in *O. sativa* leaves(**A**–**C**) and roots (**D**–**F**). Bars represent means ± SD (*n* = 3). Difference among the treatments was analyzed by Duncan's multiple range test. Different letters above the bar indicate a significant difference at *p* < 0.05.

#### *2.5. CorrelationAnalysis and Principal Component Analysis (PCA) Loading Plots*

The contents of leaf organic acids (malate, citrate, and isocitrate) were negatively correlated with all the parameters measured in rice leaf, whereas they were positively correlated with leaf PEPC (Figure 6A). Except for leaf PEPC, the activities of leaf-organicacids-metabolism-related enzymes were all positively correlated with leaf-N-metabolismrelated parameters (leaf N, leaf-soluble protein, leaf TFAA, leaf GOGAT, leaf GS, and leaf NR) (Figure 6A). In contrast, root-N-metabolism-related parameters, such as root GOGAT, root GS, root NR, and root TFAA, were positively correlated with enzymes related to root organic acids catabolism, such as root NADP-ME and NADP-IDH, whereas they were negatively correlated with enzymes related to root organic acids biosynthesis (such as root PEPC and NADP-MDH) (Figure 6B). Interestingly, root-soluble protein had no obvious correlation with other parameters measured in this study in rice root (Figure 6B).

The principal component analysis (PCA) loading plot generated by Sigmaplot 10.0 visualized two loadings against each other to investigate the relationships between the variables. There were three or four replicates for each treatment in PCA. Thirty-two parameters including plant weight, shoot height, tissue N contents, soluble proteins, TFAAs, organic acids, and the activities of their metabolism-related-enzymes, from *O. sativa* leaves and roots, were transformed for PCA analysis (Figure 7). The first two PCs explained 89.28% of the physiological variation in response to different N-supplying levels with PC1 being accounted for by71.56% and PC2 by 17.72%. The factor loadings are listed in Table S1. The PCA result clearly showed that tissue N contents, plant weight, and shoot height were clustered tightly. Parameters related to leaf organic acids metabolism, leaf N metabolism, and root N metabolism were clustered together, whereas parameters related to root organic acids metabolism closely clustered with leaf organic acids metabolism. Furthermore, the interconnection between leaf organic acids metabolism and leaf N metabolism was closer than that between root organic acids metabolism and root N metabolism (Figure 7).

**Figure 6.** Correlation analysis of variables in rice leaves (**A**) and root (**B**) in response to different N treatments. The direction and color of the ellipse represent the positive or negative relationship between two variables. The number in the ellipse indicates the correlation coefficient between two variables.

**Figure 7.** Principal component analysis (PCA) loading plot of the parameters in *O. sativa* leaves sand roots under different N levels. Thirty-two parameters from *O. sativa* leaves and roots were transformed for PCA analysis. The first two PCs explained 89.28% of the physiological variation in response to different N-supplying levels.

#### **3. Discussion**

Nfertilizer is an important factor for increasing yield in agricultural production. Both the previous literature and results of the current study showed that N deficiency could lead to the degradation of photosynthetic pigments in rice leaves and the yellowing of leaves, which eventually led to the reduction in photosynthetic rate and affected plant growth (Figures 1 and 2) [2]. The phenomenon that N deficiency induceda reduction in photosynthetic rate and, thus, affected plant growth has been reported in many crops, such as tea, citrus, maize, cucumber, sunflower, tobacco, and barley [2,8–11]. According to the statistics of the Food and Agriculture Organization of the United Nations (FAO), the amount of agricultural nitrogen fertilizer has increased several times in recent decades (FAOSTAT) [29]. With the prominent increase in global population and the decline in arable land, the demand for fertilizer will increase [30]. However, excessive application of fertilizer will lead to eutrophication of the atmosphere and water system, and cause severe environmental problems. Therefore, it is of great significance to evaluate the physiological and biochemical reactions of crops to the nutrient in fertilizers to guide the rational application of fertilizers.

It was found on the model plant *Arabidopsis thaliana* that C/N interaction could significantly affect the differential expression of more than 300 genes (*p* < 0.05) [5]. Among them, the organic acid metabolic enzyme gene *PEPC* was considered to play an important role in balancing C and N metabolism [17]. In the study of transcriptomic analysis inrice under N-deficiency, Shao et al. reported that N deficiency upregulated the expression of the PEPC gene in rice leaves and increased the content of malate in rice leaves [2]. In this study, we found that with the decrease inN level, the contents of malate and citrate increased significantly in rice leaves and roots (Figure 3A,B,F). Compared to the control, both low N and N deficiency significantly increased the content of isocitrate in rice leaves and roots, while the content of root isocitrate was significantly higher under N-deficiency treatment than that under low-N and the control treatment (Figure 3C,G,H). Through the determination of the activities of organic-acid-metabolism-related enzymes, it was found that the increase in malate in rice leaves caused by the reduction inN supply was related to the increased biosynthesis of malate and the decreased degradation of malate, which was consistent with the increased PEPC activity and the decreased NADP-ME activity in rice leaves (Figures 4B,D and 6A). Organic acids are synthesized by photosynthetic products through glycolysis and the TCA cycle, while 2-ketoglutarategenerated by citrate degradation in the cytoplasm is mainly used for NH4 +assimilation in leaf cells [31,32]. Masumoto et al. reported that *Osppc4* was obtained by the knockdown of a leaf PEPC gene in rice and their comparative analysis of the leaf metabolome showed that *Osppc4* knockdown greatly reduced the content of organic acids, especially malate, thus inhibiting NH4 +assimilation and subsequent amino acids biosynthesis in the GS/GOGAT cycle, resulting in the slowdown of plant growth [18]. In the current study, the increased contents of citrate and isocitrate in rice leaves caused by N reduction might not be caused by the increased activity of biosynthesis, but due to the reduced degradation of citrate and isocitrate in the TCA cycle, because, compared to the control, the enzyme activities of CS, ACO, and NADP-IDH in rice leaves under nitrogen reduction treatment all showed a downward trend (Figures 4C,E,F and 6A). Similarly, it was also reported that N deficiency led to an increase in the contents of organic acids in potato shoot [33]. As the enzyme activities of NAD-MDH and PEPC increased in rice roots under N deficiency conditions, while the activities of NADP-ME enzymes decreased, the increased content of malate in rice roots caused by low-Nand N-deficiency treatments might be related to the increase in synthesis and reduction in degradation of malate (Figures 3F, 4J,H and 6B). Compared to the control, the increased content of root citrate under low N and N deficiency might be related to the increased biosynthesis of citrate rather than the decreased degradation ofcitrate, because the activities of root CS and ACO were higher in low-N and N-deficiency treatments than that of the control ones (Figures 3G and 4K,L). At the same time, the root isocitrate content and root ACO activity were higher under N-deficiency treatment

than those under control and low-N treatment ones, while theactivity of root NADP-IDH under N-deficiency and low-N treatment was higher than that of the control ones, which indicated that the increased isocitrate content in roots was related to the increase in citrate isomerization activity and the decrease in isocitrate degradation (Figures 3H, 4I,K and 6B). It can be seen from the above results that although the change mode of organic acids in rice leaves and roots under reducing N level was basically the same, the change pattern of organic-acids-metabolism-related enzymes caused by reducing N level was different between rice leaves and roots. Our correlation analysis and PCA result showed that parameters related to leaf organic acids metabolism and root organic acids metabolism were tightly clustered together, indicating that organic acids metabolism in the leaf and root were closely interconnected (Figures 6 and 7).

Under N-deficient conditions, the content of ammonium used for amino acids biosynthesis in plants decreased. Therefore, the condensation of the amino and transamination process will decrease, which will inevitably lead to a change in the content of amino acids and proteins in plants [10,25]. In the current study, except for the soluble protein in rice roots, the contents of leaf TFAA, leaf-soluble protein, and root TFAA were decreased by the decreasing N level (Figure 3D,E,I,J). This is consistent with the results obtained in *Citrus* [10,25], tomato [34], tea [9], *Japonica* rice [4], and hybrid rice [35]. NR, GS, and GOGAT are three key enzymes for N metabolism in plants [36]. After nitrate is absorbed by roots, it is first reduced to nitrite by NR in the cytoplasm, and then reduced to NH4 + by nitrite reductase (NiR) in the plastid. Through the GS/GOGAT cycle, NH4 <sup>+</sup> produced from nitrate reduction and ammonium absorbed by rice ammonium transporters *Os*AMTs are assimilated into amino acids [27]. With the decreasing N level in nutrient solution, the activities of N-metabolism-related enzymes such as NR, GS, and GOGAT in rice leaves and roots decreased (Figure 5). Xiong et al. showed N deficiency at the tillering stage and N compensation at the spike differentiation stage, and N deficiency treatment could reduce the enzyme activities of NR and GS in the double cropping hybrid rice, while the normal N application and double N compensation groups increased the activities of NR and GS to a certain extent, thereby improving the photosynthetic performance of hybrid rice leaves and rice yield [35]. N deficiency treatment reduced the content of total N, ammonium, glutamine, and TFAAs in potato roots, and reduced the expression of N-metabolism-related genes such as *AMT* (ammonium transporter gene) and *GS* [34]. Similarly, Chen et al. and Huang et al. found that with the reduced N-supplying level in nutrient solution, the activities of NR, GS, GOGAT, and other enzymes related to N metabolism in *Citrus* leaves and roots decreased, resulting in a reduction in total free amino acid and soluble protein content in seedlings, and ultimately resulting in the inhibition of plant growth [10,25]. Interestingly, our PCA result demonstrated that leaf organic acids metabolism had a closer relationship with leaf N metabolism than those between root organic acids metabolism and root N metabolism (Figure 7; Table S1).

#### **4. Materials and Methods**

#### *4.1. Plant Material and Treatment*

The plant material used in this experiment was rice (*Oryza sativa* L. ssp. Indica) variety "Huanghuazhan" selected by the Rice Research Institute of Guangdong Academy of Agricultural Sciences. Uniform seeds of rice seeds were sowed in the plastic tray containing paddy soil with the content of alkaline hydrolysis N of 143 ± 6.98 mg/kg. When the height of the aboveground part of the rice seedlings was about 11.18 cm, the rice seedlings were transplanted to the Hogland nutrient solution containing 5 mM NH4NO3 (control), 1 mM NH4NO3 (low N), and 0 mM NH4NO3 (N deficiency), and cultured in a light incubator (28 ◦C). Seedlings were kept under a 14 h light/10 h dark regime with white photo-illumination of 150 μmol m−<sup>2</sup> s−<sup>1</sup> as well as a relative humidity of 68%. The nutrient solution was changed every two days. The formula of Hoagland nutrient solution contained the following macroelements (in mM): NH4NO3, 5 mM; KH2PO4, 1 mM; MgSO4, 2 mM; microelements (in μM):H3BO3,5 mM; MnCl2, 2 μM; ZnSO4, 2; CuSO4, 0.5 μM; (NH4)6Mo7O24, 0.065 μM; FeSO4-EDTA, 20 μM. In order to better promote the growth of rice seedlings, a concentration of 0.1 μM Na2SiO3 was added to each nutrient solution. After 15 days of culture, leaf and root samples were collected to determine the fresh weight and height of shoots after the phenotype of seedlings were recorded (Figure 1). There were four replicates for fresh weight and shoot height. At the same time, the rice plants were divided into upper and lower parts with scissors and put into aluminum foil bags. After freezing in liquid nitrogen, allthe samples were stored at −80 ◦C until assayed.

#### *4.2. Determination of N Content in Rice Leaves and Roots*

The nitrogen content was determined by the Kjeldahl method [2]. The ground samples of rice leaves and roots were digested with concentrated H2SO4-H2O2, and then the N contents in rice leaves and roots were analyzed using a Foss Kjeltec 8200 nitrogen analyzer (Hilleroed, Denmark). There were four replicates for the measurement of leaf N and root N.

#### *4.3. Determination of Organic Acids, Total Free Amino Acids (TFAA),and Total Soluble Proteins in Rice Leaves and Roots*

The extraction and determination of organic acids in rice leaves and roots were carried out according to the method described by Lu et al. [37]. Organic acids were extracted by 4% perchloric acid and centrifuged to determine the contents of malate, citrate, and isocitrate. The reaction mixture of malate contained 50 mM 3-amino-1-propanol-HCl (pH = 10), 30 mM sodium glutamate-NaOH (pH = 10), 2.7 mM NAD, 1 Uglutamate-oxaloacetate transaminase (GOT, EC 2.6.1.1), 10 U NAD-malate dehydrogenase (NAD-MDH, EC 1.1.37), and 50 μL of extracted supernatant. The reaction system of citrate contained 100 mM Tris-HCl (pH = 7.6), 0.2 mM NADH, 7 U lactate dehydrogenase (LDH, EC 1.1.1.27), 14 U NAD-MDH, 0.5 U citrate lyase (EC 4.1.3.6), and an appropriate amount of extracted supernatant. The reaction system of isocitrate contained 100 mMTris-HCl (pH = 7.6), 3.3 mM MnSO4, 0.2 mM NADP, 0.1 U NADP isocitrate dehydrogenase (NADP-IDH, EC 1.1.1.42), and an appropriate amount of extracted supernatant. There were three biological replicates for each treatment.

TFAAs were determined according to the method of Li et al. [38]. The total free amino acids were determined by the ninhydrin method. About 0.1 g of root and leaf samples were ground and extracted with 1.6 mL of 10% acetic acid. After centrifugation at 12,000× *g* for 10 min at 4 ◦C, 0.9 mL of ultrapure water, 1.5 mL of hydrated ninhydrin solution, and 0.25 mL of 1% ascorbic acid were added to a 5 mL reaction tube containing 0.1 mL of supernatant, mixed well, and reacted in a boiling water bath for 15 min. After reaction, 1.25 mL of 60% ethanol was added to the mixture, and after mixing, the absorbance value at 570 nm wavelength was measured with an Ultraviolet/Visible spectrophotometer. Leucine was used to make the standard curve and the concentration of leucine for the standard curve ranged from 0 μg/mL to 2.5 μg/mL. Each treatment was performed for three biological replicates. The total soluble protein content was determined by the Coomassie brilliant blue method described by Bradford [39]. Rice samples (about 0.1 g) were extracted by grinding with 1.6 mL of 50 mM Na2HPO4-KH2PO4 (pH = 7.0) and centrifuged at 3000× *g* for 10 min at 4 ◦C. Twenty-five microliters of supernatant and 975 μLof Coomassie brilliant blue solution were mixed and allowed to react at 25 ◦C for 5 minutes, and then the absorbance at 595 nm was determined. At the same time, bovine serum albumin (BSA) was used to prepare standard curves under the same conditions, and each treatment had three biological replicates.

#### *4.4. Determination of Enzymes Related to Organic Acid Metabolism and N Metabolism in Rice Leaves and Roots*

The extraction and activity measurement of organic-acid-metabolism-related enzymes in rice leaves and roots were carried out according to the method described by Lu et al. [37]. The extraction solution contained 50 mM HEPES KOH (pH = 7.5), 10 mM MgCl2, 2 mM EDTA-Na2, 10 mM dithiothreitol (DTT), 1% (*v/v*) Triton X-100, 5% water-insoluble PVPP, 1% (*w/v*) BSA, and 30% glycerol. About 0.1 g of leaf or root samples was ground with

1.6 mL of extraction solution. The tissue homogenate (about 1.7 mL) was centrifuged at 13,000× *g* for 5 min, and the supernatant was used for enzyme activity measurement.

The assay of the enzyme related to N metabolism in rice was carried out according to the method described by Hageman et al. [40]. The extraction solution of NR contained 10 mM cysteine, 1 mM EDTA-Na2, 5% PVPP, and 25 mM phosphate buffer (pH = 8.7). Approximately 0.1 g of rice leaf or root samples was added to pre-cooled mortar containing a small amount of quartz sand and the above phosphate buffer for grinding and extraction. The homogenate was centrifuged at 4000× *g* for 15 min at 4 ◦C. The reaction mixture for NR activity contained 0.2 mL of extract supernatant, 100 mMKNO3, and 2 mMNADH. The reaction mixture was incubated in a 30 °C water bath for 30 min. At the end of the reaction, 1% (*w/v*) sulfonamide was added immediately to terminate the reaction, and then 0.02% (*w/v*) naphthylamine solution was added to the mixture for reaction at 30 ◦C for 15 min. The total volume of the reaction mixture was 2 mL. Finally, the absorbance of the mixed solution was measured at 540 nm, and a blank control was set for each sample. NaNO2 was used to make the standard curve, and the concentrations for the curve ranged from 0.01 μg/mL to 0.1 μg/mL of NO2 −.

GS and GOGAT extraction solution contained 100 mM Tris HCl (pH = 7.5), 1 mM MgCl2, 1 mM EDTA-Na2, 1mM β-mercaptoethanol, 0.3% Triton X-100, and 5% PVPP. After grinding and centrifugation at 16,000× *g* 4 ◦C for 30 min, the supernatant was taken for enzyme determination. The GS activity measurement mixture contained a mixture of 80 mM Tris HCl, 40 mM MgCl2, 40 mM sodium glutamate, 16 mM hydroxylamine hydrochloride, and 8 mM ATP. The mixed solution was measured to react for 30 minutes at 30 ◦C in a water bath. After the reaction, a chromogenic agent containing 2% TCA solution (*w/v*), 3.5% FeCl3 (*w/v*), and 1.8% HCl (*v/v*) was added to the mixture. After the reaction for 5 min, the reaction mixture was centrifuged at 10,000× *g* for 5 min, and the supernatant was taken to measure its absorbance at 540 nm. A blank control was set for each sample. The reaction system for GOGAT determination contained 100 mM2-ketoglutarate, 10 mM KCl, 3 mM NADH, and 25 mM Tris HCl (pH = 7.6), and mixed well. Finally, 20 mM glutamate was added to start the reaction, and the enzymatic kinetic reaction was measured at 340 nm for 60 s.

#### *4.5. Experimental Design and Statistic Analysis*

In this experiment, three culture pots were set for each N treatment, and each culture pot contained 18 rice seedlings. Duncan's multiple range test was used for variance analysis between different treatments. Correlation plots were generated by Origin software (Version: Pro 2022b SR1). The experimental data and PCA were processed by Microsoft Excel 2010, analyzed by SPSS 19.0, and visualized by Sigmaplot 10.0 [41].

#### **5. Conclusions**

In this study, we found that low N and N deficiency reduced rice biomass, plant height, and N content, increased the content of main organic acids, and decreased the total free amino acids and leaf-soluble proteins in *O. sativa*. Although the dynamics of the contents of organic acids in rice leaves and roots were basically the same under low N and N deficiency, the dynamics of the activities of organic-acid-metabolism-related enzymes caused by reducing N-supplying level were different in rice leaves and roots.

With the decreased N-supplying level, the activities of N-metabolism-related enzymes, such as NR, GS, and GOGAT, decreased in rice leaves and roots, resulting in the decreased contents of TFAAs and total soluble protein, and finally resulting in the inhibition of plant growth. The signal of N deficiency might not directly act on the TCA cycle, but on the organic acid biosynthesis step in the cytoplasm or the upstream transcriptional regulation level.

The dynamics of organic acids metabolism caused by N deficiency were different in rice leaves and roots. Interestingly, our PCA result demonstrated that leaf organic acids metabolism had a closer relationship with leaf N metabolism than that between root organic acids metabolism and root N metabolism. In the future study, more effort should be focused on the specific molecular and physiological mechanism underlying how N deficiency regulates organic acids metabolism in rice plants.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11192576/s1, Table S1: Factor loadings of variables of first two principal components.

**Author Contributions:** Conceptualization, L.-H.C. and L.-T.Y.; software, Z.-X.C.; methodology, M.X.; investigation, L.-H.C. and Z.-J.Y.; writing—original draft, L.-H.C.; writing—review and editing, L.-T.Y.; funding acquisition, L.-H.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Leading Project of Department of Science and Technology of Fujian Province, China (2019N0004), Scientific Research Fund for Young Teachers of Jinshan College of Fujian Agriculture and Forestry University, China (kx211006), Collaborative Education Project of Industry-University Cooperation of Ministry of Education, China (202102549001), and National New Agricultural Science Research and Reform Practice Project of China (jx211001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

