Comparative Analysis of Transcriptome Proﬁles Reveals the Mechanisms in the Difference of Low Potassium Tolerance among Cultivated and Tibetan Wild Barleys

: Potassium (K) deﬁciency is a bottleneck for crop production. Thus, developing low K (LK)-tolerant crop cultivars to relieve the issue is extremely urgent. Our previous studies had found that Tibetan annual wild barley accessions showed a higher LK tolerance than the cultivated barley. In this study, RNA-sequencing was performed on three barley genotypes, wild (XZ153, LK tolerance; XZ141, LK sensitivity) and cultivated (ZD9, LK sensitivity) barley genotypes, to compare the transcriptome proﬁles of their shoots at two time points after LK stress. In total, 4832 genes displayed differential expression at 48 h and 15 d among three genotypes after K stress treatment, with XZ153 having much more differentially expressed genes (DEGs) at 48 h than 15 d, but it was the opposite in ZD9. Meanwhile, GO annotation analysis and KEGG pathway enrichment were implemented on 555 and 814 LK tolerance-associated DEGs at 48 h and 15 d after LK stress, respectively. Three barley genotypes differed signiﬁcantly in transcriptional level after LK treatment. The high tolerance in wild genotype XZ153 could be attributed to many factors, mainly including K channels, Ca 2+ signaling pathway, ethylene biosynthesis process, TCA cycle, glycolysis, pentose phosphate pathway, and photosynthesis. Furthermore, some candidate genes identiﬁed in this study may be used to improve the LK tolerance of barley. and further supplemented the previous studies on LK response in barley. Based on the results, we constructed and improved the possible regulatory networks of LK response in barley. Furthermore, some candidate genes identified in this study may be used for the genetic improvement of LK tolerance in barley, but further study is needed. Therefore, based on the obtained results, we will make functional assumptions about the candidate gene. For example, when the target gene is a transporter, we will focus on its transport capacity with the help of yeast. In brief, we will next carry out the overexpression and gene editing of target genes in barley for gene function analysis, and it will eventually lay the foundation for LK toler- ance breeding of barley and other crops.


Introduction
Potassium (K) is crucial for plants and has important physiological functions [1]; however, available K is quite limited in most soil for plants [2]. Thus, K deficiency poses a severe limitation for crop production in the world. The best and most effective way to relieve K deficiency is to develop low K (LK)-tolerant crop varieties. It is well documented that LK tolerance differs greatly among plant species and genotypes within a species [2], so it is possible for us to improve the traits through genetic manipulation. Moreover, it is imperative to reveal the underlying plant mechanisms responding to LK stress.
At 48 h and 15 d after LK stress treatment, shoots of seedlings in all treatments and replicates [total 36, 3 genotypes × 2 treatments × 2 sampling times × 3 biological replications] were sampled for RNA-Seq. Additionally, the selection of the two sampling time points (48 h and 15 d) was in reference and summary to the team's previous research on the transcriptome of barley in response to LK stress [11,12]. Additionally, total RNA was extracted from the harvested shoots using Trizol reagent (Invitrogen, Carlsbad, CA, USA), then the kit NEBNext Poly(A) mRNA Magnetic Isolation Module, NEB#E7490 (New England Biolabs, Inc., Ipswich, MA, USA), was used to isolate mRNA from 1 µg of total RNA. Additionally, the cDNA library was constructed from the purified mRNA according to the protocol provided with the NEBNext Ultra RNA Library Prep Kit for Illumina, NEB #E7530 (New England Biolabs, Inc., Ipswich, MA, USA), and were verified by Qubit 2.0 and Agilent 2100. Additionally, the effective concentration of the cDNA library was measured by Quantitative PCR. The library was sequenced by using an Illumina HiSeqXten platform with a paired-end sequencing length of 150 bp (PE150) at Biomarker Technologies (Beijing, China). All the above steps were conducted according to Ye et al. [11]. Lastly, pure PCR products were obtained, and library quality was further assessed.
Raw data of fasta format were firstly processed after a series of steps, including cluster generation, library sequencing, and generation of paired-end reads. During this process, clean data were procured by abandoning invalid reads, which was the basis of the downstream analyses. Additionally, the accession number of raw data was PRJNA832317 (http:// www.ncbi.nlm.nih.gov/bioproject/832317, accessed on 29 April 2022). Next, the clean reads were mapped to a reference genome and gene set (IPK, 160517_Hv_IBSC_PGSB_r1_CDS_ HighConf_REPR_annotation.fasta). We only analyzed and annotated the perfect match reads or mismatch reads. Then, we aligned RNA-Seq reads to the barley reference genomes on TopHat (http://tophat.cbcb.umd.edu/, accessed on 8 April 2022) [21], and identified splice junctions between exons [22]. reads to the barley reference genomes on TopHat (http://tophat.cbcb.umd.edu/, accessed on) [21], and identified splice junctions between exons [22].

DEG Identification, GO Annotation, and KEGG Pathway Analysis
In this study, we used the FPKM to figure the library's normalized expression data [23]. A standard of p-value ≤ 0.01 and absolute log2 (fold change) ≥ 1 through DEseq2 was applied to identify DEGs [24]. GO and KEGG enrichment analysis were carried out on the BMK Cloud platform (www.biocloud.net, accessed on October 5, 2021) using GOseq R packages [25] and KOBAS software [26], respectively.

Statistical Analysis
Significance of the differences among treatments and genotypes are measured by Duncan's Multiple Range Test on SPSS statistical software, and the difference at p < 0.05 means it is significant.

Evaluation of RNA-Seq Reads and Mapping Results
A total of 299.92Gb clean reads were obtained; each sample's clean data reached up to 6.25 Gb, and the sequence alignment efficiency ranged from 79.38 % to 85.37 %. The obtained reads could be classified into two parts, multiple and unique mapped reads (Table S1). Meanwhile, we got a total of 34,719 expressed genes from all samples, which could provide favorable conditions for subsequent expression profiling analysis.

DEGs Identification, GO Function, and KEGG Analysis
A total of 4832 genes had differential expression at 48 h and 15 d after LK stress in three genotypes (Tables S3 and S4, Figure 2). In detail, 2628 and 3246 DEGs were identified at 48 h and 15 d, respectively. XZ153 had more DEGs at 48 h than 15 d, but it was the opposite in ZD9 (Figure 2). Nearly the same amount of DEGs was identified at two time points in XZ141 ( Figure 2). Thus, the gene expression pattern of XZ153 differed from the other two genotypes (XZ141 and ZD9), and it is necessary to conduct a further analysis among these three genotypes in response to LK stress.

DEGs Identification, GO Function, and KEGG Analysis
A total of 4832 genes had differential expression at 48 h and 15 d after LK stress in three genotypes (Tables S3 and S4, Figure 2). In detail, 2628 and 3246 DEGs were identified at 48 h and 15 d, respectively. XZ153 had more DEGs at 48 h than 15 d, but it was the opposite in ZD9 (Figure 2). Nearly the same amount of DEGs was identified at two time points in XZ141 ( Figure 2). Thus, the gene expression pattern of XZ153 differed from the other two genotypes (XZ141 and ZD9), and it is necessary to conduct a further analysis among these three genotypes in response to LK stress. Then, we focused on these DEGs, which were up-regulated (unchanged) in XZ153, but down-regulated/unchanged (down-regulated) in XZ141 and ZD9. Consequently, 555 and 814 DEGs met the screening criteria for further study at 48 h and 15 d after LK stress, respectively (Tables S5 and S6).
Through GO annotation analysis, these DEGs were divided into 39 and 43 groups at two time points, which were attached to three main classifications: cellular component ( (Figure 3, Tables S7 and S8). For KEGG pathway analysis, 101 and 221 enzymes were matched to 53 and 57 KEGG pathways at 48 h and 15 d after LK stress, respectively, including DNA replication, ribosome, RNA transport, ribosome biogenesis in eukaryotes, purine metabolism, mismatch repair, ubiquitin-mediated proteolysis, plant hormone signal transduction, nucleotide excision repair, homologous recombination, and so on (Figure 4, Tables S9 and S10). Then, we focused on these DEGs, which were up-regulated (unchanged) in XZ153, but down-regulated/unchanged (down-regulated) in XZ141 and ZD9. Consequently, 555 and 814 DEGs met the screening criteria for further study at 48 h and 15 d after LK stress, respectively (Tables S5 and S6).
Through GO annotation analysis, these DEGs were divided into 39 and 43 groups at two time points, which were attached to three main classifications: cellular component ( (Figure 3, Tables S7 and S8). For KEGG pathway analysis, 101 and 221 enzymes were matched to 53 and 57 KEGG pathways at 48 h and 15 d after LK stress, respectively, including DNA replication, ribosome, RNA transport, ribosome biogenesis in eukaryotes, purine metabolism, mismatch repair, ubiquitin-mediated proteolysis, plant hormone signal transduction, nucleotide excision repair, homologous recombination, and so on (Figure 4, Tables S9 and S10).

DEGs Involved in K Transposition and Ca 2+ Signaling Pathway
Here, we paid attention to these DEGs participating in K transposition in response to LK stress and identified two important K channels (Table 1). These two channels belong to shaker-type channels, one being the outward-rectifying K channel GORK, controlling K release from the guard cell, and the other being the inward-rectifying K channel KAT. Meanwhile, three Ca 2+ sensors involved in the Ca 2+ signaling pathway were found (Table 1), including two calmodulin-like proteins (CML1 and CML29) and one calcium-dependent protein kinase (CDPK1). More importantly, their expression differed greatly among the three genotypes, with XZ153 having higher expression level (Table 1).

DEGs Involved in the Ethylene Biosynthesis Process
In addition to the K transport system and Ca 2+ signaling pathway, the DEGs participating in the SAM cycle of the ethylene biosynthesis process were also identified after LK stress. Four DEGs encoded three enzymes: S-adenosylmethionine decarboxylase, 1-aminocyclopropane-1-carboxylate oxidase, and DNA (cytosine-5)-methyltransferase (Table 2, Figure 5). Moreover, four AP2-like ethylene-responsive transcription factors were identified after LK stress, including two DEGs expressed at both time points (Table 3).

The DEGs Involved in the Four Interrelated Metabolic Pathways
In this study, we identified DEGs encoding some enzymes participating in four interrelated metabolic pathways, named the TCA cycle, pentose phosphate pathway (PPP), glycolysis, and gluconeogenesis (Table 4). In detail, ATP-citrate synthase (CS), dihydrolipoyl dehydrogenase 1(DLD), and pyruvate decarboxylase 2 (PDC) were involved in the TCA cycle; phosphogluconate dehydrogenase (PGDH) participated in the PPP; pyruvate kinase (PK), phosphoglycerate kinase (PGK), and phosphoenolpyruvate carboxylase kinase (PPCK) participated in the glycolysis pathway; and 6-phosphofructokinase was associated with the gluconeogenesis pathway (Table 4, Figure 6).

DEGs Associated with Photosynthesis
Except for the above DEGs, the DEGs involved in photosynthesis were also identified at 48 h and 15d after LK stress (Table 5, Figure 7), including ATP synthase, photosystem II reaction center protein, light-harvesting complex II chlorophyll a/b binding protein, NAD(P)H dehydrogenase (quinone) FQR1-like, and NADH-plastoquinone oxidoreductase (Table 5). Moreover, three photosystem II reaction center proteins, i.e., Psb K, Psb I, and Psb Z, were detected at 15 d, and highly expressed in XZ153. LHC II chlorophyll a/b binding protein also showed the highest expression in XZ153 at 48 h after LK stress (

DEGs Associated with Photosynthesis
Except for the above DEGs, the DEGs involved in photosynthesis were also identified at 48 h and 15d after LK stress (Table 5, Figure 7), including ATP synthase, photosystem II reaction center protein, light-harvesting complex II chlorophyll a/b binding protein, NAD(P)H dehydrogenase (quinone) FQR1-like, and NADH-plastoquinone oxidoreductase (Table 5). Moreover, three photosystem II reaction center proteins, i.e., Psb K, Psb I, and Psb Z, were detected at 15 d, and highly expressed in XZ153. LHC II chlorophyll a/b binding protein also showed the highest expression in XZ153 at 48 h after LK stress (Table 5).

K Channels and Ca 2+ Signaling Pathway Contribute to LK Stress Response in XZ153
Previous research had found many K transporters and channels in plants participated in K uptake and mobilization [27]. Most of them could be induced or up-regulated by LK stress [28]. In the current study, two kinds of shaker-type channels, GORK and KAT, were up-regulated after LK stress only in XZ153 (Table 1). Moreover, functional analysis showed that KAT1, KAT2, and GORK were basically expressed in the guard cells, and synergistically mediated the inward or outward potassium ion flow [29][30][31][32]. This result was consistent with our recent finding that GORK was up-regulated in the XZ153′ leaves under LK stress [11]. Moreover, Ca 2+ is an important second messenger responding to various biotic and abiotic stresses in plants, including nutrient deficiency [33][34][35]. Calcium signals are sensed, decoded, and conducted by calcium sensors. Additionally, Ca 2+ sensors mainly include CBLs, CDPKs, CaMs, and CML [36,37]. In this study, two CML (CML1 and CML29) and one CDPK (CDPK1) were identified, being only up-regulated in XZ153 after LK treatment (Table1). In short, the tolerant genotype XZ153 could better respond to LK stress than XZ141 and ZD9 by K channels and calcium signaling pathways.

Ethylene Biosynthesis Process May Benefit to LK Stress Response
Ethylene is the first plant hormone in response to K deficiency [38,39]. In this study, 1-aminocyclopropane-1-carboxylate synthase oxidase (ACO), S-adenosylmethionine decarboxylase, and DNA (cytosine-5)-methyltransferase encoded by four DEGs involved in

K Channels and Ca 2+ Signaling Pathway Contribute to LK Stress Response in XZ153
Previous research had found many K transporters and channels in plants participated in K uptake and mobilization [27]. Most of them could be induced or up-regulated by LK stress [28]. In the current study, two kinds of shaker-type channels, GORK and KAT, were up-regulated after LK stress only in XZ153 (Table 1). Moreover, functional analysis showed that KAT1, KAT2, and GORK were basically expressed in the guard cells, and synergistically mediated the inward or outward potassium ion flow [29][30][31][32]. This result was consistent with our recent finding that GORK was up-regulated in the XZ153 leaves under LK stress [11]. Moreover, Ca 2+ is an important second messenger responding to various biotic and abiotic stresses in plants, including nutrient deficiency [33][34][35]. Calcium signals are sensed, decoded, and conducted by calcium sensors. Additionally, Ca 2+ sensors mainly include CBLs, CDPKs, CaMs, and CML [36,37]. In this study, two CML (CML1 and CML29) and one CDPK (CDPK1) were identified, being only up-regulated in XZ153 after LK treatment (Table 1). In short, the tolerant genotype XZ153 could better respond to LK stress than XZ141 and ZD9 by K channels and calcium signaling pathways.

Ethylene Biosynthesis Process May Benefit to LK Stress Response
Ethylene is the first plant hormone in response to K deficiency [38,39]. In this study, 1-aminocyclopropane-1-carboxylate synthase oxidase (ACO), S-adenosylmethionine decarboxylase, and DNA (cytosine-5)-methyltransferase encoded by four DEGs involved in ethylene biosynthesis process were identified (Table 2, Figure 5). Additionally, 90% of SAM is used for transmethylation [40]. Interestingly, the expression of two MTs [DNA (cytosine-5)-methyltransferase] was down-regulated in XZ141 and ZD9, but unchanged in XZ153 (Table 2). Additionally, the AP2 transcription factor was a factor responding to LK stress [41]. Coincidently, many AP2-like ethylene-responsive transcription factors were up-regulated/unchanged in XZ153, while down-regulated in the other two barley genotypes in this study ( Table 3). The expression was consistent with previous studies in the roots and leaves of barley after LK stress [11,12]. Therefore, the DEGs identified in this study contribute to the high LK tolerance in XZ153, and importantly complemented the ethylene biosynthesis pathway response to LK stress in barley constructed by predecessors [11,12,42,43].

The TCA Cycle, Glycolysis, and PPP Are Closely Related to LK Tolerant Stress
The TCA cycle is a vital pathway of respiration in plants [44]. Here, we identified three kinds of DEGs involved in the TCA cycle: dihydrolipoyl dehydrogenase (DLD), citrate synthase (CS), and pyruvate decarboxylase ( Table 4). In plants, dihydrolipoyl dehydrogenase (DLD) is an important component of multi-enzyme complexes of pyruvate dehydrogenase complex [45]. DLD catalyzes the irreversible reaction of pyruvate to acetyl CoA. In this study, dihydrolipoyl dehydrogenase 1 was highly up-regulated only in XZ153 (Table 4). Interestingly, a previous study found that over-expression of DLD improved photosynthesis, and as a result, biomass was increased [46]. Thus, the activity of DLD in TCA cycle may have a great influence on plant photosynthesis. Additionally, citrate synthase (CS) is a pivotal and rate-limiting enzyme catalyzing oxaloacetate and acetyl-CoA to produce citrate [47][48][49]. Here, CS was up-regulated in XZ153 under LK stress (Table 4). Furthermore, a study found some miRNAs and their targets could mutually corroborate other enzymes in TCA cycle after LK stress in barley [42], and the DEGs identified in this study just encoded some of these enzymes. Therefore, it may be assumed that DLD and CS activity in the TCA cycle has a positive effect on the plants when exposed to LK stress.
Glycolysis is a process of glucose breaking down to form pyruvate [50]. Additionally, three enzymes, pyruvate kinase (PK), phosphofructokinase (PFK), and phosphoglycerate kinase (PGK) were identified after LK stress (Table 4, Figure 6). PK irreversibly catalyzes PEP into pyruvate in glycolysis and is a key rate-limiting enzyme [51]. The effect of PK also was proved in one research on miRNAs responding to LK stress [42]. Additionally, PFK and PGK catalyze a reversible reaction in glycolysis. Thus, these enzymes can regulate the glycolysis pathway [52,53]. The expression level of these three enzymes was higher in XZ153 than XZ141 and ZD9 (Table 4). Furthermore, phosphoenolpyruvate carboxylase kinase (PPCK) is an enzyme in the pyruvate carboxylation branch (gluconeogenesis), together with pyruvate carboxylase (PC), catalyzes pyruvate from oxaloacetate to phosphoenolpyruvate bypassing the "energy barrier" [54]. Coincidently, the expression patter of PPCK1 was similar with that of the above three enzymes (Table 4). All these reactions resulted in the pyruvate accumulation in XZ153. Therefore, the high expression of these enzymes in glycolysis may be a vital factor for XZ153 in response to LK stress.
The expression of 6-phosphogluconate dehydrogenase (6PGDH) was higher in XZ153 but unchanged in XZ141 and ZD9 after LK stress in this study (Table 4). Additionally, 6PGDH is a rate-limiting enzyme at the oxidation stage in PPP [55]. Interestingly, a study of miRNAs found that bdi-miR164c negatively regulated 6PGDH in PPP after LK stress [42], and this finding was confirmed by the current study. Therefor we speculated that XZ153 could maintain relatively normal condition in the pentose phosphate pathway under LK stress.
On the whole, a brief schematic diagram characterizing each pathway and relevant enzymes can be developed ( Figure 6). In short, these results complemented the previous study on miRNAs responding to LK stress [42], and three pathways might account for LK stress tolerance in XZ153.

The DEGs Involved in Photosynthesis Contribute to LK Tolerance
Photosynthesis is important for plant growth and biomass formation [56,57]. Here, the DEGs associated with ATP synthases, photosystem II reaction center proteins, NADHplastoquinone oxidoreductase, NAD(P)H dehydrogenase (quinone), and light-harvesting complex II chlorophyll a/b binding protein were identified after LK stress (Table 5, Figure 7). Additionally, the expression of these protein was mostly up-regulated in XZ153 (Table 5). Previous studies showed LK-tolerant crops could keep relatively normal chlorophyll content and photosynthetic rate when suffered from potassium deficiency [9,58]. Furthermore, four miRNAs, negatively regulating three proteins (OEEP1, PSII Q(B) protein, and PSI RC), were involved in photosynthesis and improved the photosynthesis of XZ153 after LK stress [42]. Therefore, better performance in photosynthesis may be attributed to LK tolerance in XZ153.

Conclusions
In this study, the transcriptome profiling of XZ153, XZ141, and ZD9 differing in LK tolerance were investigated under LK stress. According to the results, XZ153 had stronger LK tolerance, which could be attributed to many factors, mainly including K channels, Ca 2+ signaling pathway, ethylene biosynthesis process, TCA cycle, glycolysis, PPP, and photosynthesis. A brief diagram was developed to explain the LK-tolerant mechanism in XZ153 ( Figure 8). The research was an integrated study on LK responses at the whole transcriptome level in the leaves of wild and cultivated barleys, and further supplemented the previous studies on LK response in barley. Based on the results, we constructed and improved the possible regulatory networks of LK response in barley. Furthermore, some candidate genes identified in this study may be used for the genetic improvement of LK tolerance in barley, but further study is needed. Therefore, based on the obtained results, we will make functional assumptions about the candidate gene. For example, when the target gene is a transporter, we will focus on its transport capacity with the help of yeast. In brief, we will next carry out the overexpression and gene editing of target genes in barley for gene function analysis, and it will eventually lay the foundation for LK tolerance breeding of barley and other crops.

The DEGs Involved in Photosynthesis Contribute to LK Tolerance
Photosynthesis is important for plant growth and biomass formation [56,57]. Here, the DEGs associated with ATP synthases, photosystem II reaction center proteins, NADHplastoquinone oxidoreductase, NAD(P)H dehydrogenase (quinone), and light-harvesting complex II chlorophyll a/b binding protein were identified after LK stress (Table 5, Figure  7). Additionally, the expression of these protein was mostly up-regulated in XZ153 (Table  5). Previous studies showed LK-tolerant crops could keep relatively normal chlorophyll content and photosynthetic rate when suffered from potassium deficiency [9,58]. Furthermore, four miRNAs, negatively regulating three proteins (OEEP1, PSII Q(B) protein, and PSI RC), were involved in photosynthesis and improved the photosynthesis of XZ153 after LK stress [42]. Therefore, better performance in photosynthesis may be attributed to LK tolerance in XZ153.

Conclusions
In this study, the transcriptome profiling of XZ153, XZ141, and ZD9 differing in LK tolerance were investigated under LK stress. According to the results, XZ153 had stronger LK tolerance, which could be attributed to many factors, mainly including K channels, Ca 2+ signaling pathway, ethylene biosynthesis process, TCA cycle, glycolysis, PPP, and photosynthesis. A brief diagram was developed to explain the LK-tolerant mechanism in XZ153 ( Figure 8). The research was an integrated study on LK responses at the whole transcriptome level in the leaves of wild and cultivated barleys, and further supplemented the previous studies on LK response in barley. Based on the results, we constructed and improved the possible regulatory networks of LK response in barley. Furthermore, some candidate genes identified in this study may be used for the genetic improvement of LK tolerance in barley, but further study is needed. Therefore, based on the obtained results, we will make functional assumptions about the candidate gene. For example, when the target gene is a transporter, we will focus on its transport capacity with the help of yeast. In brief, we will next carry out the overexpression and gene editing of target genes in barley for gene function analysis, and it will eventually lay the foundation for LK tolerance breeding of barley and other crops.   Table S7: Functional classification of LK-tolerance-related DEGs after 48 h LK stress; Table S8: Functional classification of LK-tolerance-related DEGs after 15 d LK stress; Table S9: Pathway analysis of LK-tolerance-related DEGs based on the KEGG data after 48 h of LK stress; Table S10: Pathway analysis of LK-tolerance-related DEGs based on the KEGG data after 15 d of LK stress;