Genome-Wide Identification, Expression and Interaction Analyses of PP2C Family Genes in Chenopodium quinoa

1. Introduction

Protein phosphorylation and dephosphorylation are the main forms of reversible post-translational modifications, which control the important regulatory mechanisms of many biological processes by regulating the localization, conformation, stability, and activity of substrate proteins in eukaryotes [1]. The phosphorylation state of proteins is dynamically controlled by protein kinase (PK) and protein phosphatase (PP), where protein kinase transfers the phosphate group of donor ATP to the side chain of receptor protein, while protein phosphatase dephosphorylates phosphoprotein. According to their mechanism of catalysis, substrate specificity and specific response to inhibitors, eukaryotic PPs can be divided into protein tyrosine (Tyr) phosphatases (PTP), phosphoprotein phosphatase (PPP), metallo-dependent protein phosphatase (PPM), and aspartate (Asp)-dependent phosphatase [2]. The PTP family includes Tyr specific phosphatases (PTPs) and dual-specificity phosphatase (DsPTP) that dephosphorylates serine (Ser), Thr (threonine), and Tyr phosphoresidue. The PPP and PPM families are Ser/Thr-specific phosphatases (STPs) [3]. The PPP family consists of seven members: PP1, PP2A, PP2B, and PP4/5/6/7 [2]. Although PP1, PP2A, and PP2B share a common catalytic domain, these enzymes exhibit the greatest differences in their N- and C-terminus without catalytic function, as well as the formation of various intact enzymes through their different regulatory subunits. PP2A holoenzyme is a heterotrimer consisting of A (scaffold), B (regulatory) and C (catalytic) subunits [4]. The PPM family includes protein phosphatase 2C (PP2C), pyruvate dehydrogenase phosphatase, and other magnesium (Mg²⁺)/manganese (Mn²⁺)-dependent STPs [5].

PP2Cs widely exist in prokaryote and eukaryote, are evolutionarily conservative, and significantly regulate stress signal pathways [6]. The relatively conserved catalytic domain in eukaryotic PP2C protein is located at the N- or C-terminus, while the region of non-catalytic domain is not highly conserved and has diverse amino acid sequences with different functions. The non catalytic domain region is important for defining the function of PP2C members, as it contains sequence motifs and/or transmembrane regions related to cellular signaling, including those that interact with protein kinases [7,8]. PP2Cs, as negative regulators of protein kinase cascades activated in different processes, participate in regulating signaling pathways. In fission yeast, genetic evidence has shown that PP2Cs are involved in negative regulation of osmotic sensing signals transmitted through the Wis1-MAPK cascade [9]. In budding yeast, two PP2Cs, PTC1 and PTC3, are negative regulators of the PBS2-HOG1 MAPK pathway. In the yeast HOG pathway, four types of PP2C phosphatases Ptc1-Ptc4 dephosphorylated differently two activated phosphorylation sites of Pbs2 MAP2K [10]. In humans, PP2Calpha negatively regulates the stress-responsive MAPK cascades through dephosphorylation and inactivation of MKK6, SEK1, and MAPK (p38) [11], while PP2Cbeta dephosphorylates and inactivates the MAPKK kinase TAK1 to negatively regulate the TAK1 stress-signaling pathway [12] and PP2Cepsilon associates stably with TAK1 and dephosphorylates TAK1 to inhibit the TAK1 signaling pathway [13]. PpABI1A and PpABI2B of group A PP2C are directly involved in ABA response, acting downstream of ABA activated kinase and regulating ABA induced genes in the moss P. patens [14].

So far, PP2C family genes have been studied in many plants. In the study of Arabidopsis, six out of nine members (ABI1, ABI2, HAB1, HAB2, AHG1 and AHG3) of the A subfamily have been verified as ABA co-receptors, negatively regulating the ABA signaling pathway [15]. The MAPK phosphatase AP2C1 of the PP2C subfamily B interacts with MPK3, MPK4, and MPK6 to control their activity [16]. The expression of AP2C1 and the accumulation of AP2C1 protein are strongly and locally enhanced at the induction site of the syncytium, indicating that AP2C1 acts as a negative regulatory factor for MAPK (MPK3, MPK4, and MPK6) to ensure inhibition of MAPK activation in the developing syncytium [17,18]. POLTERGEIST (POL) and POL-LIKE1 (PLL1) encode the related protein phosphatases 2C of the PP2C subfamily C, which are crucial for the establishment of shoot and root meristem tissues during embryogenesis and the maintenance of stem cell pools during post-embryonic development in Arabidopsis [19,20]. AtPP2C of subfamily D participates in the response to saline and alkali stresses [21]. In subfamily E, AtPP2C-6-6 interacts with histone acetyl transferase AtGCN5 to participate in transpiration by regulating stomatal signaling [22]. In subfamily F, WIN2 is involved in modulating plant defense by interacting with the bacterial effector HopW1-1 [23]. Likewise, the protein phosphatase homologue 1 (PPH1) of unclustered PP2Cs involves maintaining efficient photosynthesis through dephosphorylation of Lhcb1 and Lhcb2 in plants [24]. In rice, OsPP2C09 (Os01g62760) of subfamily A PP2C interacts with RING-H2 type E3 ligase OsRF1 to participate in salt tolerance of rice [25]. SAL1 (sensitive to aluminum 1) encodes PP2C. D phosphatase located on the plasma membrane, interacts with PM H⁺- ATPase and inhibits its activity, participating in rice aluminum resistance [26]. The subfamily F PP2C phosphatase ZmPP84 participates in regulating drought stress responses by dephosphorylating ZmMEK1 to inhibit its kinase activity in maize [27]. Likewise, the subfamily B PP2C phosphatase ZmPP2C26 can dephosphorylate ZmMAPK3 and ZmMAPK7, participating in negative regulation of drought tolerance in maize [28]. All these researches indicate that PP2Cs have multiple functions and are worthy of further research.

Quinoa (Chenopodium quinoa Wild.) is a cultivated crop with high nutritional value in the world. At the same time, it exhibits strong resistance to various soil and climatic conditions such as high salinity, drought, and cold, which allows quinoa to be planted on marginal land. PP2C is a multifunctional gene that regulates plant growth, development, and stress response, and has been studied in many plants. However, genome-wide study on the PP2C gene family has not yet been performed in quinoa. This study conducted a comprehensive genome-wide analysis of the CqPP2C gene family in quinoa, including gene identification, chromosomal localization, phylogenetic relationships, gene structures, conserved motifs and domains, gene duplication analysis, cis-acting elements analysis, and relative expression of CqPP2C genes. In addition, we analyzed the possible interactions between subfamily A CqPP2Cs and CqSnRK2s, both of which respond to abscisic acid and abiotic stress. These results will provide important information for understanding the mechanisms of PP2C in abiotic stress signal transduction.

2. Results

2.1. Identification and Basic Information of PP2C Genes in Quinoa

A total of 121 PP2C-coding candidate genes were identified via BLASTP and HMM searching in quinoa (Chenopodium quinoa). By using the CDD program with default settings, it was found that 4 of the 121 candidate PP2Cs did not contain PP2C catalytic domains. Therefore, 117 genes in quinoa were identified as members of the PP2C family and labeled as CqPP2C1 to CqPP2C117 based on their order on chromosomes. The information of gene name, gene ID, chromosome location, number of amino acid (aa), molecular weight (Mw), isoelectric point (pI), instability index, hydrophilic coefficient, and subcellular localization prediction of 117 PP2C proteins were analyzed (Table S1). The lengths of proteins ranged from 110 aa residues (CqPP2C94) to 1501 aa residues (CqPP2C101). The Mw ranged from 11495.04kDa (CqPP2C94) to 167616.03kDa (CqPP2C101) and pI varied from 4.08 (CqPP2C8) to 9.44 (CqPP2C70). The result of the instability index indicated that 62.3% CqPP2Cs were unstable proteins. Except for CqPP2C7, CqPP2C94 and CqPP2C108, all other CqPP2C showed GRAVY below zero, indicating that these proteins are hydrophilic. The results of subcellular localization prediction showed that most of the quinoa PP2Cs might be located in the cytoplasm, chloroplast or nucleus. In addition, only CqPP2C6 might be located in the plasma membrane, CqPP2C7 and CqPP2C83 might be located in the mitochondria, CqPP2C51, CqPP2C72, CqPP2C73 and CqPP2C99 might be located in the endoplasmic reticulum, and CqPP2C57 might be located in the vacuole membrane. These results indicated that CqPP2C proteins were randomly distributed in cells and played a role in various environments.

2.2. Phylogenetic Analysis of CqPP2C Genes

To study the phylogenetic relationships between PP2C genes in quinoa and Arabidopsis, we used the maximum likelihood method to construct a phylogenetic tree based on the alignments of 80 PP2C protein sequences in Arabidopsis and 117 in quinoa (Figure 1).

The results showed that 103 CqPP2C proteins were divided into 13 subfamilies (A-L), including A (14), B (8), C (8), D (14), E (17), F1 (8), F2 (6), G (9), H (6), I (4), J (2), K (3), L (4) (Table S2). In addition, quinoa has seven separate branches. Among them, CqPP2C1, CqPP2C53, CqPP2C70 and AT1G18030 formed a separate branch; CqPP2C60, CqPP2C79 and AT2G40860 formed a separate branch; CqPP2C55, CqPP2C83 and AT4G27800 formed a separate branch; CqPP2C107, CqPP2C54 and AT3G23360 grouped into a separate branch; CqPP2C93 and AT4G11040 formed a separate branch; CqPP2C21, CqPP2C68 and CqPP2C109 grouped into a separate branch; CqPP2C99 formed a separate branch.

2.3. Gene Structural and Conserved Domain Analyses of CqPP2Cs

We analyzed the exon/intron structure patterns of CqPP2C genes and protein conserved motifs. In quinoa, the number of exons of CqPP2Cs varied from 1 to 18. Among them, CqPP2C93 in unclustered PP2Cs and CqPP2C94 and CqPP2C108 in subfamily A only had one exon without intron, while CqPP2C72 in subfamily J had 18 exons and 17 introns, and had the longest intron. Generally, most genes in the same subfamily shared a similar exon/intron structure (Figure 2A).

Using MEME motif search tool, twenty motifs were identified for CqPP2C proteins. As shown in (Figure 2B), the number of motifs ranged from 3 to 11 with 8–50 residues in all CqPP2C proteins. Among them, 107 CqPP2C proteins all contained motif 2. Each group had specific motifs, except for common motifs. For example, motifs 8 and 9 existed in group D but not in other groups, motifs 13 and 17 only existed in groups C and D, motifs 18 only existed in group H, and motifs 19 only exist in group F1. Proteins in the same group had similar motif distribution patterns, with the G family having exactly the same motif distribution, indicating that CqPP2C members in the same cluster might have similar functions.

We used the NCBI CDD/SPARCLE database to predict conservative structural domains (Figure 2C). Most of CqPP2Cs have a PP2Cc domain, with the exception of the nine members in G subfamily, the three members in K subfamily and CqPP2C21. The seven members of the G subfamily (CqPP2C2, 9, 24, 25, 67, 71, 78) have PP2C_C super family domain and CqPP2C10 and CqPP2C38 in G subfamily have a PLN03145 domain belonging to the PP2C super family. CqPP2C8, 47, 95 in K subfamily and CqPP2C21 have a PP2Cc superfamily domain (Serine/threonine phosphatases, family 2C, catalytic domain). Additionally, domains FHA-PP2C70-like, PKC-like superfamily, GUB_WAK_bind, MDR superfamily, ZnF-BED, CAP-ED and 2A194 superfamily also appear in quinoa PP2C protein sequences.

As is well known, the A subfamily proteins (PP2CAs) of PP2Cs are involved in controlling abscisic acid (ABA) signaling and responding to various abiotic stresses, and have a negative regulatory effect on plant growth and development. To further investigate their biological functions, further comparison was made between the subfamily A protein of PP2Cs in quinoa and the reported PP2C proteins in Arabidopsis. The catalytic domain of PP2C proteins contain 11 conserved motifs, in which 5 conserved residues participate in Mg²⁺/Mn²⁺ coordination. The multiple alignments results of the 14 CqPP2C (CqPP2CA) and 9 AtPP2C (AtPP2CA) in subfamily A indicated that not all CqPP2CA members contained all the 11 conserved motifs (Figure S1). It was found that CqPP2C94, CqPP2C108, and CqPP2C44 were partial deletions of the N-terminus or/and C-terminus of the PP2C catalytic domain, which might lead to the elimination of some important motifs and loss of functions. Five sites responsible for Mg²⁺/Mn²⁺ coordination were found within all CqPP2C catalytic domains: [xxxD], [DGxG], [CGD], [DG] and [xxDN] (C-cysteine; D-aspartic acid; G-glycine; N-asparagine) (Figure S1). Subfamily A proteins of PP2C (PP2CA) had several residues involved in their phosphatase activity in the catalytic domain. Among them, the critical active-site residues [Arg138, Glu142, Asp143, Asp177, Gly178, His179, Asp347 and Asp413 in ABI1] had been conserved. Similarly, the ABA-sensing tryptophan [Trp385 (W385) in HAB1] had been conserved. Whereas, the Arg residue [Arg505 (R505) in HAB1] that mediates interaction between the ABA box and HAB1 [29] showed less conservativeness. The well-described residues responsible for ABI1-PYL1 interaction [30] had also been conserved (Figure S1). These results indicated that the structures of PP2CA proteins were similar, especially within the highly conserved catalytic domain, despite greater changes in the N-terminal region (Figure S1).

2.4. Chromosomal Location and Duplication of CqPP2C Genes

In total, 113 CqPP2Cs were localized to the 18 chromosomes of quinoa and only four CqPP2Cs could not be mapped to any chromosome and were therefore assigned to chromosome zero. According to their order on chromosomes, they were named CqPP2C1–CqPP2C117 (Table S1 and Figure 3.

The largest number of CqPP2C genes were localized to chromosome 07 (17 CqPP2Cs), while chromosome 09, chromosome 13 and chromosome 18 had the smallest number of CqPP2Cs (only three CqPP2Cs). Gene duplication caused by polyploidization or replication related segments and tandem duplication is the main mechanism for producing new genes, which contributes to the gene family expansion in the plant kingdom [31]. In this study, we found that there were no tandem duplication gene pairs, but 61 pairs of paralogous CqPP2C genes were involved in segment duplication events, indicating that segment duplication was the driving force for the expansion of the quinoa PP2C gene family. Among them, 2 CqPP2C genes pairs, including CqPP2C30/32/64/116 and CqPP2C29/74/101/113, were found to have four copies. 3 CqPP2C genes pairs, including CqPP2C12/13/36, CqPP2C31/63/117 and CqPP2C44/56/94, were found to have three copies, other genes pairs contained two copies. There is only one copy of the remaining 18 CqPP2C genes in quinoa (Figure 3). These results suggested that gene loss might also occur in the PP2C gene family of quinoa, resulting in the loss of some homologous copies. Similar findings were also observed in the CqWRKY and CqNAC gene family in quinoa [32]. The ratio of Ka/Ks was less than 1 (Table S4), indicating that purification selection was the main evolutionary mode for all duplicated CqPP2C genes.

2.5. Cis-acting Elements Analysis

To further investigate the function of CqPP2C genes, 2000 bp upstream promoter sequences of CqPP2C genes were analyzed with PlantCARE. In addition to common cis-elements such as TATA boxes and CAAT boxes, other cis-elements were related to abiotic stress responses, light, hormones, plant growth and development, and other regulatory stresses (Figure 4).

Among them, light response elements accounted for the largest proportion and were found in almost all CqPP2C genes. Hormone responsive elements included abscisic acid, auxin, gibberellin, salicylic acid, and methyl jasmonate responsive elements. Abiotic stress responsive elements included drought, anaerobic induction, low temperature, and defense and stress responsive elements. Obviously, there were many cis-elements related to plant abiotic stress in the promoter region of the CqPP2C gene. The promoter of 89 out of 117 CqPP2C genes had ABA-responsive elements (ABRE), and CqPP2C110 contained twelve ABA responsive elements, with the largest number. These indicated that CqPP2C genes played an important role in abiotic stress responses through the ABA signaling pathway

2.6. Expression of CqPP2C Genes in Different Quinoa Tissues

To investigate the expression of CqPP2C family genes in different quinoa tissues including apical meristems, leaves petioles, flowers and immature seeds, dry seeds, stems, seedling, internode stems, inflorescence, leaves, fruit of white sweet quinoa, flowers of white sweet quinoa, flowers of yellow bitter quinoa, and fruit of yellow bitter quinoa, a heatmap was constructed using previously published RNA-seq data. As shown in Figure 5, CqPP2C family genes had different expression levels in different tissues.

In general, most CqPP2C genes displayed an extremely broad expression range, and several CqPP2C genes showed high expression in at least one or more tissues. We identified three genes with leaves petioles-specific expression (CqPP2C55, CqPP2C59, and CqPP2C83). CqPP2C12 was highly expressed in apical meristems and CqPP2C51, CqPP2C61, and CqPP2C111 accumulated in the flower and immature seeds. In addition, nine CqPP2C genes were abundant in internode stems and nineteen genes were found to be highly expressed in seedling. We also found subfamily A genes of PP2C (PP2CA) were abundant in dry seeds and fruit of white sweet quinoa except CqPP2C5, CqPP2C22, and CqPP2C105. This analysis identified candidate CqPP2C genes that might play important roles in the development of different organs.

2.7. Expression Patterns of CqPP2C Genes under Stress Conditions

The prediction of cis-acting elements had indicated that CqPP2C genes might be involved in the response to drought, cold, and NaCl stress. In addition, many studies had shown that PP2C gene expression in different plant species was regulated by abiotic stress and hormone treatment [33]. In this study, we analyzed the expression levels of CqPP2C genes in root and shoot under different abiotic stresses using transcriptome data (Figure 6).

It was observed that the expression of PP2C genes in shoot was lower than that in root in control. Some genes (CqPP2C30, CqPP2C55, CqPP2C82, and CqPP2C83) in shoot and genes (CqPP2C5, CqPP2C19, CqPP2C22, CqPP2C32, CqPP2C40, CqPP2C43, CqPP2C44, CqPP2C56, CqPP2C64, CqPP2C90, CqPP2C104, CqPP2C109, and CqPP2C116) in root were significantly increased under drought stress compared with control. Some of these genes (CqPP2C2, CqPP2C6, CqPP2C51, CqPP2C52, CqPP2C60 and CqPP2C105) in shoot and genes (CqPP2C7, CqPP2C16, CqPP2C33, CqPP2C34, CqPP2C36, CqPP2C42, CqPP2C62, CqPP2C63, CqPP2C80, CqPP2C81, CqPP2C100, CqPP2C110, CqPP2C112, CqPP2C114, and CqPP2C117) in root were significantly increased under heat stress compared with control. Some of these genes (CqPP2C9, CqPP2C24, and CqPP2C101) in shoot and genes (CqPP2C13, CqPP2C14, CqPP2C17, CqPP2C28, CqPP2C31, CqPP2C48, CqPP2C59, CqPP2C61, CqPP2C68, CqPP2C99, and CqPP2C107) in root were significantly increased under low phosphorus stress compared with control. Only CqPP2C71, CqPP2C73, CqPP2C84, and CqPP2C89 in root were significantly increased under salt stress compared with control. The expression of three genes (CqPP2C42, CqPP2C43, and CqPP2C46) in shoot and three genes (CqPP2C52, CqPP2C82, and CqPP2C83) in root were significantly lower than that in control under drought stress. The expression of three genes (CqPP2C4, CqPP2C39, and CqPP2C47) in shoot were significantly lower than that in control under low phosphorus. In conclusion, we found that most of the CqPP2C genes were highly or moderately expressed in root and shoot in different treatments, which indicated that most of the CqPP2C genes responded to abiotic stress.

It is well known that the subfamily A PP2Cs in rice and A. thaliana were transcriptionally regulated in abiotic stress responses dependent on ABA signaling pathways. The expression of several members of CqPP2C genes in subfamily A was examined by qRT-PCR under drought and salt stress in root and shoot. Six genes from subfamily A (CqPP2C5, CqPP2C30, CqPP2C44, CqPP2C104, CqPP2C105, and CqPP2C116) were randomly selected for this analysis (Figure 7).

The results showed that all six CqPP2C genes had varying degrees of response to these two stresses. Under drought treatment, six genes were up-regulated in both root and shoot, while under salt treatment, these six genes were only up-regulated in the root. In the shoot, only CqPP2C5 was up-regulated during salt treatment, while the other five genes were down-regulated. Some results were consistent with the analysis of the transcriptome data (Figure 6). In summary, the expression patterns of subfamily A CqPP2C genes indicated that these genes were responsive to abiotic stress.

2.8. Protein Interaction between Subfamily A CqPP2Cs and CqSnRK2s

Subfamily A PP2Cs (PP2CAs) had been shown to be negative regulators of the ABA signaling pathway. In this study, we studied the interactions between subfamily A CqPP2Cs and CqSnRK2s using yeast two-hybrid assay. Five CqSnRK2 members had been isolated previously in our lab, including three members of subclass II CqSnRK2s [CqSnRK2.1 (AUR62027801), CqSnRK2.4 (AUR62003254), and CqSnRK2.9 (AUR62007175)] and two members of subclass III CqSnRK2s [CqSnRK2.6 (AUR62003840) and CqSnRK2.11 (AUR62011423)]. We cloned the six subfamily A PP2C genes (CqPP2C5, CqPP2C30, CqPP2C44, CqPP2C104, CqPP2C105, and CqPP2C116) that had been detected for gene expression using qRT-PCR (Figure 7). As shown in Figure 8, CqPP2C104 and CqPP2C105, the homolog of AtABI1/2, could interact with CqSnRK2.1 and CqSnRK2.4, while CqPP2C104 also strongly interacted with CqSnRK2.9 but weakly interacted with CqSnRK2.11. CqPP2C5, CqPP2C30, and CqPP2C44 strongly interacted with CqSnRK2.11 and weakly interacted with CqSnRK2.1 and CqSnRK2.4.

In addition, the homologue CqPP2C30 of AtAHG3 strongly interacted with CqSnRK2.6, while the homologue CqPP2C116 of AtHAI1/2/3 only weakly interacted with CqSnRK2.1 and CqSnRK2.4. The results showed that the six subfamily A CqPP2Cs exhibited complex interactions with the five CqSnRK2s.

3. Discussion

Plant PP2Cs play a vital role in regulating many key biological processes involved in development and response to stresses [34,35]. In current study, we comprehensively analyzed the CqPP2C genes in Chenopodium quinoa. A total of 117 CqPP2C genes were identified. Compared with A. thaliana (80), rice (78) [36], B. distachyon (86) [6], M. truncatula (94) [37], cucumber (56) [38], tomato (56) [39], and maize (97) [40], the amount of PP2C in quinoa was much more. Although the genome sizes of higher plants such as rice and Arabidopsis are comparable to those of lower plants such as green algae (Chlamydomonas reinhardtii), lycophyte (Selaginella moellendorffii) and moss (Physcomitrella patens), there are only 10 PP2C genes in green algae, and 50 PP2C genes in lycophyte and moss, while 78 in rice and 80 in Arabidopsis [41]. This indicated that there were differences in the expansion of the PP2C genes among different species, which may be related to the evolution of plants from unicellular organisms to multicellular organisms.

According to the phylogenetic tree, the PP2C genes of quinoa were divided into 13 subfamilies (Figure 1), consistent with the PP2C groups in A. thalinan, O. sativa [36], B. distachyon [6], and cucumber [38]. In phylogenetic analysis, different PP2C groups of quinoa and Arabidopsis were arranged together to form a common branch, indicating that PP2C had sequence conservation and similar evolutionary lineages. Phylogenetic analysis can identify homologous genes from different species to predict gene function. In the A subfamily, the CqPP2C18, CqPP2C19, CqPP2C104 and CqPP2C105 protein was homologous with AtABI1 (AT4G26080) and AtABI2 (AT5G57050), indicating that these four CqPP2C proteins might be involved abiotic stress in plants and were believed to have a negative regulatory effect on ABA signaling [42,43]. Similarly, the AtPP2CF1 (AT3G05640) protein in the E subfamily can activate cell proliferation and expansion, as well as accelerate inflorescence growth, and its homologous CqPP2C41 and CqPP2C88 may have the same function [44].

The exon/intron structure of genes and protein conserved motifs are important markers of the evolutionary relationship of family genes. Accordingly, we analyzed gene structure and protein conserved motifs of CqPP2Cs (Figure 2). The results showed that CqPP2Cs within the same subfamily had a similar exon/intron structure, with some exceptions, which might be due to different reasons. Previous studies on Brachypodium distachyon, Fragaria vesca, and Fragaria ananassa had shown that there were many PP2C genes with intron deletion [6,45], and similar results had been found in quinoa research. Twenty conserved motifs were identified. As shown in the Figure 2B, CqPP2Cs in the same subfamily exhibited similar motif distribution. When analyzing conservative domains, in addition to the main PP2C phosphatase domain, we also found 14 other domains. KAPP (kinase-associated protein phosphatase) is an Arabidopsis PP2C that contains FHA (forkhead associated domain) at the N-terminus of its kinase interaction region, which is crucial for connecting to phosphorylated target proteins and thus facilitates signal transduction [46]. Therefore, studying the important functions of CqPP2Cs carrying these special structural domains would be of interest.

The spatial expression analysis of CqPP2C genes revealed tissue-specific expression patterns. The abundant presence of the PP2C subfamily A gene in the dry seeds and fruit of white sweet quinoa indicated that these genes are involved in ABA mediated seed development, dormancy, and germination [47]. Cis-acting elements are important regulators of resistance to various stresses and hormone response in plant development. The analysis of promoter region of CqPP2C genes showed that ABREs and DREs (Drought response elements) elements were abundant in most CqPP2C gene promoter regions, indicating that CqPP2C genes may be involved in the ABA signaling pathway acting on drought or salt stress resistance. Further expression analysis showed that most CqPP2C genes were responsive to drought, salt, heat stress, and Pi starvation. Many studies have shown that subfamily A PP2Cs are involved in the ABA signaling pathway and play an important role in plant responses to abiotic stress [48,49]. We identified fourteen PP2C genes belonging to the A subfamily in quinoa through sequence alignment and evolutionary analysis. We found that most members of the A subfamily in quinoa were significantly up-regulated or down-regulated under drought and/or salt stress, which was consistent with reports from other plants. AtABI1 and AtABI2 have been identified as the main components of the ABA signaling pathway in response to stresses [41,48,49]. CqPP2C64 and CqPP2C116, which are homologs of AtHAI PP2Cs, were significantly up-regulated under drought stress. AtHAI PP2Cs had unique drought resistance functions in Arabidopsis. The HAI PP2C mutant reduced the expression of several defense related genes under low water potential, but increased the expression of abiotic stress-related genes encoding late embryogenesis abundant proteins and dehydratin, as well as increased the accumulation of proline and osmoregulatory solutes [49]. Likewise, BdPP2Cs in B. distachyon and MtPP2Cs in M. truncatula from subfamily A were also up-regulated by cold, heat, drought, salt or H₂O₂ treatment [6,37]. In P. euphratica, ABA has a moderate inducing effect on PeHAB1, while drought stress has a significant inducing effect on PeHAB1 [50]. Otherwise, the expression of TaPP2C59 in wheat and most FvPP2Cs in Fragaria vesca were significantly down-regulated under drought and high salt stress, suggesting that these genes play a negative regulatory role [45]. Most members of subfamily A AtPP2Cs have been identified as negative regulators of ABA signaling. The expression of these genes was highly inducible under ABA and abiotic stresses [36,51]. The induction of PP2C gene expression may be the ABA desensitization mechanism that regulates ABA signaling and maintains plant homeostasis. Therefore, ABA and abiotic stresses upregulate genes encoding both positive and negative effectors of its signaling network.

Similar to the research results of other plants, our study in quinoa also indicated that some CqPP2Cs from subfamilies other than subfamily A was induced by abiotic stresses. In the B subfamily, the expression of CqPP2C63 and CqPP2C117 was highly induced under heat stress, and the expression of CqPP2C28 and CqPPS2C31 was significantly altered after Pi starvation treatment in root. The expression of CqPP2C52 in shoot and the expression of CqPP2C7, CqPP2C36, and CqPPS2C31 in root were highly induced under heat stress. A study had shown that almost all members of subfamily D in Arabidopsis and soybean contain heat stress response elements (HSE) in their promoters, and subfamily D genes in wheat respond to heat treatment [21,52]. These results indicate that CqPP2Cs can play important functions under different stresses, and their detailed roles still need further exploration.

In Arabidopsis, members of subfamily A PP2Cs can interact with both ABA receptor PYLs and subclass III SnRK2s, mediating the ABA signaling pathway to regulate seed germination and response to abiotic stress [53,54]. Subclass III AtSnRK2 (AtSnRK2.2, AtSnRK2.3, AtSnRK2.6) proteins always interact with subfamily A PP2Cs and are inactivated by direct dephosphorylation of subfamily A PP2Cs [53,55]. In our study, results of yeast two-hybrid assay showed that all six CqPP2Cs interacted with one or two members of CqSnRK2s. Subfamily A CqPP2Cs interacted not only with subclass III CqSnRK2s, but also with subclass II CqSnRK2s. This result was consistent with studies in rice, B. distachyon, and wheat. In rice, OsSAPK2 was classified as Class 2b (subclass II) and could interact with OsPP2C30 [56]. In B. distachyon, group A BdPP2C could interact with subclass II BdSnRK2.1 [57], and group A TaPP2C interacted with subclass II TaSnRK2s in wheat [52]. This suggested that subfamily A CqPP2Cs were essential for ABA signal transduction in quinoa, where subfamily A CqPP2Cs bound to CqSnRK2s in different combinations and intensities to respond to various biological processes and stresses.

4. Materials and Methods

4.1. Identification of the PP2C Gene Family in Quinoa

This experiment used BLASTP and Hidden Markov Model (HMM) [58] to identify the PP2C genes in quinoa. The sequences of 80 AtPP2C proteins were downloaded from the Arabidopsis information resource website (https://www.arabidopsis.org/). The quinoa genome data (http://www.cbrc.kaust.edu.sa/chenopodiumdb/) was used for BLASTP and the retrieval threshold was set as E-value < E⁻¹⁰. To further screen candidate genes, we compared the protein sequences of Arabidopsis PP2C using TBtools. Then, the NCBI-CDD databases (https://www.ncbi.nlm.nih.gov/cdd/) and Pfam (http://pfam.xfam.org/search#tabview=tab1) [59] were used for domain identification of candidate gene. We manually deleted candidate genes without specific domains of the PP2C (registration numbers PF00481, PF07830, and PF13672) [52]. The physicochemical properties were analyzed based on ExPASy (https://web.expasy.org/protparam/) [60]. Online software WoLF PSORT (http://www.genscript.com/wolf-psort.html) predicted the subcellular localization of CqPP2Cs [61].

4.2. Evolutionary Relationship of the PP2C Gene Family

Validated quinoa PP2C protein sequences and Arabidopsis PP2C protein sequences (AtPP2Cs) were used to establish an evolutionary relationship. A total of 197 amino acid sequences were used for this analysis. Multiple protein sequences were aligned using the MEGA 7, and the final comparison result was constructed by neighbor-joining (NJ), and the bootstrap value was set to 1000 [36,38]. The image is beautified by iTOL (https://itol.embl.de/).

4.3. Gene Structure and Protein Conserved Motif Analysis

Conserved motifs of CqPP2Cs were determined using MEME (http://meme-suite.org/tools/meme) for all CqPP2C sequences. The number of motifs was set to 20, and other parameters were defaults [62]. The gene structure was extracted from the genome annotation gff3 files. The results were displayed using TBtools. DNAMAN v.8.0 software was used for CqPP2Cs amino acid multiplex sequence alignment.

4.4. Chromosomal Location and Gene Duplication Analysis

The synteny analysis within the quinoa genome based on all-vs-all BLASTP alignments was analyzed in TBtools using the plugin MCScanX method [63]. The images of duplication events and synteny analysis were visualized using TBtools software [64]. The synonymous (Ks) and nonsynonymous (Ka) substitution rates of homologue genes in quinoa were calculated using TBtools software [64].

4.5. Analysis of Cis-acting Elements in the Promoter Regions

The promoter sequence (2000 bp upstream of the putative genes ATG) was extracted using TBtools. The promoter sequence was then submitted to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for cis-element prediction [65].

4.6. Analysis of CqPP2C Gene Expression Patterns

The quinoa RNA-seq data from different tissues (No: PRJNA394651) and different treatments (No: PRJNA306026) were downloaded from Bioproject database (http://www.ncbi.nlm.nih.gov/sra). RNA-seq data in TPM (transcripts per million reads) is normalized and performed log₂ conversion. The heatmap of PP2C gene expression was visualized using TBtools software.

4.7. Quinoa Treatment and RNA Extraction

The quinoa material used in this experiment was “YT077”. The quinoa seeds were placed in the sterilized nutrient soil and cultured in a greenhouse (average temperature 22 °C, 70-75% relative humidity, 16 h/8 h light/dark regime). To determine the expression of quinoa PP2C under salt and drought treatment. One-month-old quinoa seedlings with consistent growth were selected for salt stress (300 mM NaCl) and drought stress (20% PEG6000). The leaves and roots of seedlings were taken 0 h, 3 h, 6 h, 12 h, 24 h and 48 h after treatment. The collected samples were frozen quickly with liquid nitrogen and placed at −80 °C until further usage. Three biological replicates for each treatment. TransZol Up Plus RNA Kit (Transgen, China) was used to extract total RNA from collected samples.

4.8. qRT-PCR Analysis

The primers were designed using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) (Table S5), synthesized by Qingdao BGI. The above preserved RNA was reverse-transcribed into single-stranded cDNA using the TransScript^® One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen, China). The cDNA product obtained by reverse transcription was used as a template, and the SYBR Green I dye method was used for real-time PCR, and the quinoa Tubulin gene was used as the internal reference gene. The reaction system is 20 μL and consists of 2 μL of cDNA, 2 μL of forward and reverse primers primers, 10 μL of SYBR, and 6 μL of ddH₂O. The reaction procedure is as follows: pre-denaturation at 95 °C for 30 s; Denaturation at 95 °C for 10 s; Annealing at 57 °C for 15 s; Extension at 72 °C for 30 s and set the number of cycles to 40 times. The relative expression of genes was calculated by using the 2^-ΔΔct method [66].

4.9. Yeast Two-hybrid Assays

The primers (Table S6) for the cloning of CqPP2Cs and CqSnRK2s were designed using the Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/). Six CqPP2Cs in group A and five CqSnRK2s were amplified from the quinoa cDNA. The CqPP2C and CqSnRK2 genes were cloned into pGADT7 and pGBKT7 vectors, respectively. Yeast two-hybrid analysis was performed using yeast strain AH109 according to the manufacturer’s protocol (Clontech, USA) for at least three times. Positive transformants screened from SD medium lacking leucine and tryptophan (SD/-Leu/-Trp) are then transferred to SD medium lacking leucine, tryptophan, and histidine (SD/-Leu/-Trp/-His) for further screening. The 3-amino-1, 2, 4-triazole (3-AT) was used to eliminate the background yeast growth on SD medium lacking leucine, tryptophan, and histidine (SD/-Leu/-Trp/-His).

5. Conclusions

In present work, a total of putative 117 CqPP2C genes were identified and divided into 13 subfamilies. The chromosome localization, physical and chemical feature predictions, phylogenetic analysis, gene structure, conserved motif and domain analysis were thoroughly investigated. The results of collinearity and selection pressure analysis indicated that CqPP2C genes underwent amplification of segmental duplication and purification selection during evolution. In addition, the expression of CqPP2Cs in various tissues under different abiotic treatments was analyzed using RNA-seq data. CqPP2C genes were involved in regulating development and stress responses of quinoa. qRT-PCR results showed that six CqPP2C genes in subfamily A were up-regulated or down-regulated under salt and drought treatments. Additionally, Yeast two-hybrid assays showed that subfamily A CqPP2Cs interacted with CqSnRK2s in different combinations and intensities to respond to various biological processes and stresses. Overall, our results provide new insights and a basis for further understanding the roles of CqPP2C family in the regulation of abiotic stress response.