Genome-wide identification of protein phosphatase 2C family members in Glycyrrhiza uralensis Fisch. and their response to abscisic acid and polyethylene glycol stress

The protein phosphatase 2C (PP2C) gene family, which is known to participate in cellular processes, is one of the most conserved plant-specific gene families that regulate signal transduction in eukaryotic organisms. However, the PP2C gene family of Glycyrrhiza uralensis Fisch. has not yet been systematically characterized. In this study, 76 PP2C genes of G. uralensis (GuPP2Cs) were identified. Phylogenetic analysis divided these GuPP2Cs into 12 subfamilies. The gene structure and protein motifs of the GuPP2Cs were then analysed. The analysis of the transcriptomic data of the GuPP2Cs revealed that nine genes were highly sensitive to abscisic acid (ABA) stress and might respond to ABA treatment. These nine genes were validated via quantitative real-time polymerase chain reaction. This study provides a foundation for future studies on the functional PP2C genes of G. uralensis and the improvement of the cultivation and quality of G. uralensis.


Introduction
Abscisic acid (ABA), one of the six major hormones in plants, acts as an endogenous messenger of plant responses to biotic and abiotic stresses [1], such as drought stress [2], cold stress [3], salt stress [4] and pathogen infection [5]. Therefore, revealing the genes involved in the ABA signalling pathway is essential for improving plant quality and yield. Numerous proteins, especially ABA receptors, that are associated with the ABA signalling pathway and ABA stress responses have been identified. Type 2C protein phosphatase (PP2C), which is known as an interaction protein of the ABA receptor, plays an important role in response to various stresses, especially drought [6]. When the ABA signal is missing, the PP2C gene inactivates SnRK2 kinase through phosphorylation; in response to ABA signalling, the receptor PYR/PYL/RCAR binds to ABA and inhibits the activity of PP2C. SnRK2 kinase is then activated so that it can phosphorylate downstream factors to achieve signal transmission and cell response to ABA stress. In this manner, the ABA signalling pathway can be transmitted downwards, and the drought resistance of plants can be improved [7,8].
A total of 80 (13 subfamilies) and 78 (11 subfamilies) PP2C family members have been annotated in the genomes of Arabidopsis thaliana and Oryza sativa, respectively [9]. Several PP2Cs in plants have been found to be negative regulators within the ABAmediated signalling network [10][11][12][13][14]. The PP2CA subfamily of A. thaliana has been identified as a key negative regulator in ABA signalling [15]. BdPP2CA6, a PP2CA gene of Brachypodium distachyon, reportedly plays a positive regulatory role in the ABA and stress signalling pathways [16]. Other subfamilies of PP2C, such as subfamilies PP2CB, PP2CC, PP2CD and PP2CE, have also been reported [17][18][19]. As key regulators of ABA responses, PP2Cs are central to understanding the integrative network of ABA signalling.
Glycyrrhiza uralensis is the main source of liquorice, which has great medical and economic value. Its roots and rhizomes are widely used as herbal medicine and as the raw materials of traditional Chinese medicine. It is also extensively utilized as food and chemical © 2022 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
additives [20]. Pharmacological studies have shown that G. uralensis has anti-inflammatory, antiviral, antimicrobial, anticancer, immunomodulatory, gastroprotective, hepatoprotective, neuroprotective and cardioprotective properties [21]. Glycyrrhizic acid and liquiritin are the main secondary metabolites in G. uralensis. Previous studies have demonstrated that moderate drought or ABA stress can promote the accumulation of these active components, thereby improving the quality of G. uralensis [22,23]. Therefore, the PP2C genes of G. uralensis must be identified to understand the mechanism through which the accumulation of secondary metabolites is regulated under stresses, as well as to improve the quality of G. uralensis.
The availability of the whole G. uralensis genome provides convenience to the research on the PP2C genes of G. uralensis [24]. In this study, we identified PP2C genes by using the genomic data of G. uralensis. We then investigated the gene families of G. uralensis via the phylogenetic analysis of the PP2C genes of A. thaliana. The transcriptome of G. uralensis under treatment with different concentrations of ABA and a specific concentration of polyethylene glycol (PEG) and at different time points was sequenced to study the expression patterns of PP2C genes. Nine GuPP2C genes responding to ABA stress were screened by using transcriptomic data and verified via quantitative real-time polymerase chain reaction (qRT-PCR). The results of this work lay a foundation for further studies on the molecular and regulatory mechanisms of ABA signalling to secondary metabolites and pave the way for improving the cultivation and quality of G. uralensis.

Identification and analysis of the PP2C genes of G. uralensis.
The amino acid sequences of each subfamily (A to L) in A. thaliana were downloaded from TAIR (https://www. arabidopsis.org/) [25] in March, 2021. The amino acids of G. uralensis with E-values of 1×10 −5 were aligned by using the protein basic local alignment search tool (BLASTP) [26]. PP2C gene sequences were extracted from the genome of G. uralensis [24]. The coding sequence (CDS), DNA and amino acid information of PP2Cs was controlled by using APOLLO [27]. The PFMA website (http://pfam.xfam.org/search/batch) was applied to predict the domains of PP2Cs [28]. The ExPaSy website (https://web.expasy.org/protparam/) was utilized to determine the properties of PP2Cs, including the number of amino acids, molecular weight (MW), theoretical isoelectric point (pI), instability index (II), aliphatic index (AI) and grand average of hydropathicity (GRAV). Gene Structure Display Server 2.0 (GSDS, http://gsds.cbi.pku.edu.cn/) was used to show the structural features of the PP2C genes of G. uralensis [29]. MEME (https://meme-suite.org/) was used to analyse the motifs of PP2Cs [30].

Phylogenetic analysis of PP2C genes of G. uralensis
The full-length PP2C protein sequences of G. uralensis and A. thaliana [9] were aligned by using MAFFT [31], and a maximum likelihood (ML) tree was constructed in accordance with the LG + F+R5 model and 1000 bootstraps by using IQ-TREE [32] to confirm the evolutionary relationships among the members of the PP2C gene family.

Transcriptomic sequencing and expression analysis of GuPP2C genes in response to ABA stress
The genome of G. uralensis was downloaded from http://ngs-data-archive.psc.riken.jp/Gur-genome/dow nload.pl [24]. A total of 48 2-year-old G. uralensis samples were obtained from Beijing Medical Botanical Garden in Beijing City and identified by Professor Yulin Lin from the Institute of Medicinal Plant Development. G. uralensis plants with good and consistent growth were chosen. The leaves of G. uralensis were subjected to different stress conditions, including 10, 25, 50, and 100 mg/L ABA and 20% PEG4000. The leaves were collected at 0, 3, 6 and 12 h [33] after ABA and 20% PEG4000 treatments. The samples were immediately frozen in liquid nitrogen and stored in a freezer at −80°C for total RNA extraction and transcriptome sequencing. Three samples were taken from each treatment as biological repeats. Total RNA was isolated by using RNAprep Pure Plant Plus Kit (TIANGEN, DP441, China) by following the manufacturer's instructions. Transcriptomic data were obtained with high-throughput sequencing technology. The heat maps of transcriptomic data were drawn by using TBtools software [34]. Complementary DNA was obtained with the M-MLV reverse transcriptase synthesis system (Promega, Madison, WI, USA) by following the manufacturer's protocols. Nine GuPP2C genes with relatively large differences in expression were screened for reverse qRT-PCR verification. The actin gene from the genome of G. uralensis was selected as the internal reference gene. The specific primers of GuPP2C genes were designed by the Primer5 programme (PREMIER Biosoft International, CA, USA). The primers of the nine genes are listed in Table S1. The amplification conditions were 95°C for 3 min, 95°C for 3 s, 57°C for 30 s and 72°C for 15 s for 40 cycles; 65°C for 5 s and 95°C for 0.5 s.

Statistical analysis
The results were calculated on the basis of three biological replicates (three technical replicates per biological replicate). Statistical analyses and significant correlations were performed by using SPSS (IBM, USA).

Genome-wide identification and phylogenetic analysis of the PP2C genes of G. uralensis
A total of 76 PP2C genes were detected from the genomic sequence database of G. uralensis. All of these 76 GuPP2Cs were included to construct a ML tree, and 26 AtPP2Cs from A. thaliana were used as references ( Figure 1). In accordance with the AtPP2C subfamily, the GuPP2Cs were divided into 12 subfamilies, namely, PP2CA, B, C, D, E, F1, F2, G, H, I, J and L, which contained 11, 2, 3, 15, 11, 7, 4, 8, 5, 2, 1 and 3 GuPP2C genes, respectively. The PP2CA, PP2CD and PP2CE subfamilies had the largest gene numbers, whereas the PP2CJ subfamily had only one gene. In contrast to A. thaliana, G. uralensis did not contain PP2CK genes. Four genes (GuPP2CX01-04) clustered into an independent branch that did not belong to the 12 subfamilies. As shown in Figure 1, the PP2CC and PP2CD subfamilies and the PP2CE and PP2CL subfamilies constituted a monophyletic sister branch with 100% bootstrap support, suggesting that the PP2CC and PP2CD subfamilies and the PP2CE and PP2CL subfamilies had a close evolutionary relationship.

Analysis of the gene structures and protein motifs of the PP2C genes of G. uralensis
The lengths of the GuPP2Cs varied from 797 bp to 13 598 bp with the CDSs ranging from 633 bp to 3 264 bp in length (Table S2)  34.2%) had II values of less than 40 and were considered as stable, and the others were considered as unstable. The AI was between 66.28 and 94.91. All of the GuPP2C proteins were hydrophilic given their negative GRAVY values.
The gene structure and protein motifs of the PP2C genes of G. uralensis were analysed. The results of gene structure analysis showed remarkable differences, which were mainly reflected in the numbers and lengths of introns (Figure 2a), amongst these subfamilies. The majority of these subfamilies, except for the PP2CL and PP2CA subfamilies, had similar numbers of exons. The PP2CE, PP2CI, PP2CB, PP2CC and PP2CG subfamilies had high genetic similarity. In GuPP2Cs, the numbers of introns ranged from 1 to 18, and the numbers of exons ranged from 2 to 19. The numbers of introns and exons greatly varied amongst the 12 subfamilies. The GuPP2CJ01 gene had the largest number of introns (18) and exons (19). Consistent with their evolutionary relationship depicted by the phylogenetic tree, most members in the same subfamily shared a similar intron/exon structure and gene length.
Twenty-six protein motifs were identified in the GuPP2C proteins (Figure 2b). Their sequences are presented in Table 1. Certain differences were observed in motif composition amongst these subfamilies. Some protein motifs were widespread amongst the GuPP2C proteins, whereas some were specific to only one or two subfamilies. Motifs 1, 2, 3 and 4 widely existed in over 90% of the GuPP2C protein sequences and were present in all of the subfamilies. However, 12 motifs were specific to a certain subfamily. Four motifs (18, 22, 23 and 26) were distributed in the PP2CH subfamily only; three motifs (14, 15 and 17) existed in the PP2CE subfamily only; two motifs (6 and 17) were present in the PP2CD subfamily only; two motifs (19 and 20) were found in the PP2CF1 subfamily only and motif 21 existed in the PP2CA subfamily only. Motifs 5 and 8 were specific to the PP2CC and PP2CD subfamilies; motif 10 was specific to the PP2CD and PP2CH subfamilies and motif 25 was specific to the PP2CA and PP2CD subfamilies. Amongst the different gene members of the same subfamilies, the main motifs were similar to each other. The composition of some subfamilies, such as PP2CC and PP2CD, was very conservative as reflected by their identical protein motifs, except for individual members. The motifs of the PP2CA subfamily greatly varied. The different motif distributions in the protein sequences may serve as a reference for studying the divergence of gene functions in different subfamilies. Moreover, the highly similar motif distribution of the proteins in the same subfamily supported the close evolutionary relationship of these proteins [9].

Expression pattern analysis of GuPP2C genes under different stress conditions
The gene expression patterns of GuPP2Cs under different ABA and PEG stress conditions were investigated. The transcriptome from the leaves of G. uralensis were sequenced. Their fragment per kilobase of exon model per million mapped read (FPKM) values were calculated (Table S3). The results of differential expression analysis showed that 27.6% (n = 21) of the genes had a low expression (FPKM < 10), and included six genes with silent expression (FPKM < 1). By contrast, 72.4% (n = 55) of the genes had high expression (FPKM ≥ 10) (Figure 3). The expression levels of 18 GuPP2C genes increased under ABA stress. The average expression levels of GuPP2CA08, GuPP2CA12, GuPP2CA13, GuPP2CB02, GuPP2CC01, GuPP2CE04, GuPP2CE10, GuPP2CF14 and GuPP2CG02 in the ABA stress group were at least two times higher than those in the control group, indicating that these genes were more sensitive to ABA stress than other genes. These genes may respond to ABA treatment and can be used as candidate genes for further research. In addition, the expression levels of 22 GuPP2Cs were decreased, and those of 36 GuPP2Cs were not influenced by environmental stresses. Nine candidate genes (GuPP2CA08, GuPP2CA12, GuPP2CA13, GuPP2CB02, GuPP2CC01, GuPP2CE04, GuPP2CE10, GuPP2CF14 and GuPP2CG02) with increased sensitivity to ABA stress and that might respond to ABA treatment were screened by calculating FPKM values. The expression patterns of the nine candidate genes after treatment with 10, 25, 50 and 100 mg/L ABA for 3, 6 and 12 h were validated by qRT-PCR. The qRT-PCR results showed that the expression levels of six genes, i.e. GuPP2CA08, GuPP2CA12, GuPP2CE04, GuPP2CE10, GuPP2CF14 and GuPP2CG02, were significantly changed after ABA treatment (Figure 4). The expression patterns of GuPP2CA08, GuPP2CE10 and GuPP2CG02 were similar. Overall, the expression of these three genes significantly increased after 3 and 6 h of ABA treatment and then sharply decreased after 12 h of treatment. The differences amongst the three genes were as follows: GuPP2CA08 expression showed no significant change after 12 h of treatment; GuPP2CE10 expression significantly increased after 12 h of treatment with 50 and 100 mg/L ABA and GuPP2CG02 expression significantly decreased after 12 h of treatment. The gene expression levels of GuPP2CE04 and GuPP2CF14 were significantly increased under ABA treatment at all concentrations and treatment times. However, the expression of GuPP2CF14 significantly decreased under treatment with 50 mg/L ABA for 6 h. Compared with that under the control treatment, the expression of GuPP2CA12 was significantly reduced under ABA treatment at all concentrations and treatment times but not under treatment with 25 mg/L ABA for 6 h.
Significant correlations between the expression values obtained through qRT-PCR and transcriptome analysis were identified ( Table 2). The expression values of GuPP2CA08, GuPP2CE10 and GuPP2CG02 in different datasets were significantly correlated, whereas the others were not. Poor correlations may be due to sample differences.

Discussion and conclusions
The PP2C genes code for the interaction proteins of the ABA receptor, which can release SnRK2 kinase activity when its own activity is inhibited [6]. ABA plays an essential role in plant response to abiotic stress. Hence, in plants, ABA content increases to improve resistance to drought, salt stress and other abiotic stresses [35]. Therefore, identifying PP2C genes is vital for improving the abiotic stress resistance of plants. The functions of PP2C genes have been characterized in other plants, such as Salvia miltiorrhiza [36], Gossypium hirsutum [37], Brassica rapa [33] and Fragaria vesca [38], aside from A. thaliana.
To date, several studies have highlighted the importance of PP2C, which may play pivotal roles in various processes, including biotic-abiotic stresses and plant development [39]. In plants, PP2C members have been described as regulators of desiccation tolerance and as negative regulators in the ABA signalling pathway [12,13] Previous studies have reported that AtPP2CG1 regulates positively against salt tolerance in A. thaliana and is induced by drought, salt or exogenous ABA treatment [40]. In A. thaliana, two PP2C gene members show different responses. Specifically, AP2C1 expression is powerfully induced by drought, wounding and cold, whereas the same stresses only slightly induce AP2C2 expression [41]. Moreover, the expression profiles of PP2C in various species after exposure to abiotic and hormone stresses have been studied [42,43]. Shazadee et al. predicted the cis-regulatory elements in upland cotton and suggested that the GhPP2C genes in G. hirsutum may be involved in the response to heat, cold, drought and NaCl stresses. They subjected cotton seedlings to four various abiotic stress treatments and selected 30 genes for qRT-PCR. Their results showed that some GhPP2C genes exhibited high transcript levels after exposure to multiple treatments and that a few were induced by one or more treatments. For example, heat stress accounted for the dominant portion of down-regulated genes. However, cold stress up-regulated 70% of the genes and decreased the transcript level of 30% of the genes. Drought stress accounted for the up-regulation and down-regulation of 53% and 47% of GhPP2C genes, respectively, and salt stress resulted in the up-regulation and downregulation of 56% and 44% of the genes, respectively. Taken together, all of the 30 genes were induced by different abiotic stresses; the diverse expression profiles of GhPP2C genes may suggest that these genes may be critical for abiotic stress responses [37]. Similarly, the diverse expression patterns of BraPP2C genes in B. rapa are suggestive of various roles after exposure to abiotic and hormone stress conditions [33]. In addition, studies have confirmed that some PP2C genes can regulate plant tolerance to drought. Yu et al. [44] reported that the wheat protein phosphatase PP2C-A10 can promote seed germination and decrease the drought tolerance of transgenic A. thaliana by interacting with the TaDOG1L1 and TaDOG1L4 genes. Miao et al. [45] showed that in rice, OsPP2C09 acts as a negative regulator of drought tolerance through ABA signalling. Under ABA treatment, the faster and higher accumulation of OsPP2C09 transcripts in roots than in shoots may increase the root-to-shoot ratio of plants under drought stress. Xiang et al. [46] found that ZmPP2C-A10 was tightly associated with drought tolerance in Zea mays. Moreover, they found that allele-338 containing a deletion of the endoplasmic reticulum stress response element in the 5'-UTR region of ZmPP2C-A10 enhanced plant drought tolerance.
Liquorice root is one of the most popular traditional Chinese medicines, and its main active components are flavonoids. G. uralensis is an important source of liquorice. Previous studies have confirmed that in plants, moderate drought or ABA stress can promote the accumulation of secondary metabolites [23], which are important components used in quality evaluation [22]. However, the characterization and functional analyses of PP2C genes in G. uralensis that respond to ABA treatment have not been reported yet. Here, we report the PP2C genes of G. uralensis. A total of 76 PP2C genes were identified in G. uralensis. The number of the PP2C genes of G. uralensis was similar to that of other higher plants, such as A. thaliana (80) [9], O. sativa (78) [9], B. distachyon (86) [42] and Setaria italica (80) [47]. The PP2C genes were divided into 12 subfamilies, namely, subfamilies GuPP2CA-L. Phylogenetic analysis revealed that most members in the same subfamily shared similar intron/exon structures. Moreover, the results suggested that the same subfamily may perform similar biological functions. The analysis of the transcriptomic data of the GuPP2C genes indicated that GuPP2CA08, GuPP2CA12, GuPP2CA13, GuPP2CB02, GuPP2CC01, GuPP2CE04, GuPP2CE10, GuPP2CF14 and GuPP2CG02 were highly sensitive to ABA stress and might respond to ABA treatment. Gene expression pattern analysis via qRT-PCR showed that the expression levels of six genes were significantly changed after ABA treatment and were correlated with ABA concentration and treatment time. Overall, the expression patterns of GuPP2CA08, GuPP2CE10 and GuPP2CG02 changed dynamically with the treatment time. Moreover, the expression levels of GuPP2CE04 and GuPP2CF14 significantly increased, and the expression of GuPP2CA12 significantly reduced.
The present study conducted the first genome-wide identification and characterization of the phylogenetic relationships, gene structures and motifs of the PP2C gene family of G. uralensis. The expression patterns of the GuPP2C genes that responded to ABA treatment were then analysed. The results of this work provide a reference for the genome-wide identification of the PP2C gene family of other species and lay a foundation for future functional research on the PP2C genes of G. uralensis.

Data availability statements
The raw data of RNA-seq is stored in NCBI (https://www. ncbi.nlm.nih.gov/). The associated BioProject and BioSample numbers are PRJNA784170 and SAMN2347 5019, respectively.

Disclosure statement
No potential conflict of interest was reported by the authors.