Genome-wide identification, evolution and expression analysis of the PP2C gene family in ice plant

As sessile organisms, plants have developed complex mechanisms in response to stress. Protein reversible phosphorylation is one of the most important regulatory strategies. Protein kinases phosphorylate serine (Ser), threonine (Thr), and tyrosine (Tyr), while protein phosphatases eliminate phosphate groups^[1,2]. According to their substrate specificity, protein phosphatases can be divided into three categories: Ser/Thr phosphatases, phosphotyrosine phosphatase, and dual specificity phosphatases^[2,3]. PP2C (Protein phosphatases 2C), the metal ion-dependent protein phosphatase, belongs to Ser/Thr phosphatases. PP2Cs, which merged early in evolution as distinct protein phosphatases, are the largest family of protein phosphatases, and are identified in archaea, bacteria, fungi, plants, and animals^[4,5]. To date, PP2Cs have been reported in many higher plants including Arabidopsis (80), rice (90), wheat (95), peanut (178), tomato (92), maize (102), cucumber (56), cotton (181 in Gossypium hirsultum, 87 in Gossypium arboreum, 147 in Gossypium barbadense, 99 in Gossypium raimondii), sweet cherry (76), Rosa chinensis Old Blush (69), soybean (126), quinoa (117), strawberry (56 in Fragaria vesca, 228 in Fragaria ananassa), etc^[1,2,5−15]. PP2Cs are highly conserved in evolution. They are divided into 13 clades (A−M) and some single branch members based on phylogenetic analysis^[4,7]. Most of them have a conserved catalytic region at the C-terminus and a poorly conserved region at the N-terminus. The diversity of the N-terminus indicates that they have different functions in the signal transduction mechanism^[1,2,4].

PP2Cs play vital roles in the plant response to abiotic environmental pressures^[1,4]. Clade A PP2Cs participate in abscisic acid (ABA) signaling by physically interacting with SnRK2s, which inhibits activation of the downstream transcription factors^[4,16−18]. When ABA accumulates, it binds to PYR1, which in turn binds to and inhibits PP2Cs. Then, the downstream SnRK2 kinases are released to active downstream target genes. Thus, clade A PP2Cs function as a switch of the ABA signaling network. This regulation often happens under osmotic stress conditions caused by drought and salinity^[4,15−19]. However, transcription analysis demonstrates that many PP2Cs are up or down expressed under salt treatment, indicating other clades of PP2Cs are also involved in stress response^[6,20]. Several PP2C B, and F repress the MAPK signal pathways activated by drought^[4,21−23]. AtPP2C49 (G clade) represses AtHKT1.1 to regulate Na⁺ distribution in Arabidopsis after salt treatment^[24]. AtPP2C.D6 and D7 interact with SOS1 and inhibit its Na⁺/H⁺ antiporter activity in unstressed conditions. Under salt stress, SCaBP8 interacts with the PP2Cs suppresses their phosphatase activity, and simultaneously releases PP2C.D6 from the plasma membrane^[25]. AtPP2CG1 is fully expressed under salt stress and positively regulates salt tolerance in an ABA-dependent manner^[26]. Over-expression of ZmPP2C2 (clade B) in tobacco enhances cold stress tolerance by displaying higher germination speed and rate, higher antioxidant enzyme activities, and lower cold-induced electrolyte leakage and malondialdehyde contents^[27]. A PP2C D clade gene, AtFGT2, affects the heat stress memory. Fgt2 mutants are defective in the memory of heat stress^[28]. Taken together, PP2Cs are important for plant abiotic stress responses.

Mesembryanthemum crystallinum (ice plant) is a halophyte that can grow under high salt conditions^[29−32]. Mature ice plants are capable of growing in soil that contains 450 mM NaCl, which is higher than the concentration found in seawater^[29]. On MS medium, the root growth of ice plant continues in 140 mM NaCl, while the root growth of Arabidopsis thaliana is completely inhibited^[29]. Moreover, ice plant switches from C3 photosynthesis to crassulacean acid metabolism (CAM) under high salinity stress, which is related to the key enzyme-phosphoenolpyruvate carboxylase (McrPEPC/08G01316) induction in leaf^[30,31]. The C3-CAM switch allows ice plants to improve their salt tolerance by high efficiency of water use and CO₂ fixation^[32]. In 2022, a high-quality ice plant genome sequence was reported, providing the data resource for comparative and functional genomic studies of ice plants on salt tolerance and C3-CAM transitions^[31].

Since PP2Cs play crucial roles in plant stress response and ice plant can be adapt to extreme environments, we conducted a genome-wide investigation of the PP2C gene family in ice plant. In total, 59 McrPP2Cs were identified and phylogenetic evolution, gene replication, selection pressure, and cis-regulatory elements were analyzed using bioinformatic methods. In addition, transcriptomic data was used to study McrPP2Cs expression patterns under different NaCl treatments. The results will lay the foundation for understanding the potential roles of McrPP2Cs under salt stress and why ice plant has a strong salt tolerance.

The data related to the Arabidopsis thaliana genome used in this study were sourced from the TAIR website (Tair10, www.arabidopsis.org). The genome and related data of Mesembryanthemum crystallinum were sourced from the database at NGDC (https://ngdc.cncb.ac.cn/). The data related to the genomes of Vitis vinifera and Oryza sativa were sourced from the JGI database (https://genome.jgi.doe.gov/portal/). The data related to the genomes of Hylocereus undatus, Amaranthus caudatus, and Beta vulgaris were sourced from the NCBI database (www.ncbi.nlm.nih.gov). The hmm file of the PP2C protein domain (PF00481) was downloaded from the Pfam database. Then, the whole genome of M. crystallinum and the five other species was identified using HMMER3.0^[33] under the condition of an E-value of 1 × 10⁻⁵. Through literature review, the sequence numbers of 80 PP2C genes of A. thaliana were obtained. The PP2C protein sequences of 80 A. thaliana were downloaded from the TAIR database (www.arabidopsis.org). The BLASTP^[34] alignment was carried out between these sequences and the genome sequences of the other six plants at an E-value of 1 × 10⁻¹⁰, and the members of the PP2C gene family (with a score > 100) were screened out. The genes obtained by taking the union of the sequences identified by HMMER3.0 and BLASTP were regarded as the preliminary members of the PP2C family. Finally, all the queryable databases were selected in the InterPro database^[35], and the genes with the annotation of the keyword 'PP2C' were screened out as the final results.

Sequence analysis and basic information of the PP2C gene family

Wolf PSORT (WoLF PSORT: Protein Subcellular Localization Prediction) was used to predict the subcellular localization of PP2C proteins in M. crystallinum and six other plants. Expasy (SIB Swiss Institute of Bioinformatics | Expasy) was employed to analyze the molecular mass, instability coefficient, hydrophilicity, aliphatic index, and other parameters of the predicted PP2C proteins in M. crystallinum and the other plants.

Chromosomal localization of PP2C family members

Based on the genome annotation files, TBtools-II Gene Location Visualize from GTF/GFF module^[36] was used to extract the chromosome length and the positional information of PP2C genes on the chromosomes from the genome annotation files of M. crystallinum and the other six plants. Then, TBtools-II was utilized to map the positions of PP2C genes on the chromosomes.

Construction of the phylogenetic tree of the PP2C gene family

The PP2C protein sequences of M. crystallinum and the other six species were subjected to multiple sequence alignment using Clustalw2 with its default parameters. Then, the Iqtree2 software^[37] was employed to construct a phylogenetic tree. The Model-Finder Plus (MFP) method was used for the automatic estimation of the optimal substitution model, which turned out to be VT + R8. The bootstrap value was set to 1000, and other parameters were set to their default values.

Identification of collinearity and gene types in the PP2C gene family

Based on the BLASTP software, a homology search was conducted on the protein sequences of M. crystallinum, H. undatus, A. caudatus, B. vulgaris, V. vinifera, A. thaliana, and O. sativa. Then, the -icl module in WGDI was used to perform collinearity analysis^[38]. The Advanced Circos in TBtools was utilized to visualize the collinearity results. Meanwhile, the Nei-Gojobori method in the -ks module of WGDI was employed to calculate the values of Ka, Ks, and Ka/Ks^[39]. The duplicate_gene_classifier in MCScan^[40] was used to count the number of different types of duplication events. Based on the neutral substitution rate 1.5 × 10⁻⁸ substitutions per site per year, the orthologous gene pairs were used for the divergence time estimation^[41].

Analysis of the gene structure, conserved motifs, and cis-acting elements of the PP2C genes

Based on the coding sequences (CDS) and genomic DNA sequences of PP2Cs, TBtools was used to draw the gene structure distribution maps. The conserved motif analysis of PP2Cs was completed through the MEME online website, with the maximum number of motifs set to 10. Then, TBtools was employed to visualize the results. 2,000 bp upstream of the PP2Cs start codon (ATG) were extracted as the promoter sequences. The PlantCARE database^[42] was utilized to predict the cis-acting elements and their potential related functions, and home-make Python was used for visualization.

Plant growth conditions

Ice plant seeds were directly planted in soil (Klasmann-Deilmann 876) after imbibed at 4 °C for 2 d. Then the plants were grown at 23 ± 1 °C under 200 μmol·m⁻²·s⁻¹ white light with a 14 h/10 h day/night cycle. After about 6 weeks, the ice plants were irrigated with 500 mM NaCl once daily. After 36 h (short time) or 14 d (long time), the second pair of mature leaves were collected for RNA sequencing and qRT-PCR examination. Plants irrigated with water were used as a control.

Expression analysis

RNA sequencing data have been previously published^[29−32] or from Prof Song research group at the North China University of Science and Technology. For qRT-PCR analysis, total RNA was isolated using RNA Easy Fast Plant Tissue Kit. First-strand cDNA was synthesized using the FastKing RT Kit (with gDNase). QRT-PCR was performed using SuperReal PreMix Plus (SYBR Green). All the above reagents were purchased from Tiangen, Tiangen Biotech (Beijing) Co., Ltd. One μL cDNA was used as a template. The qRT-PCR condition was two-step: ① 95 °C 15 min; ② 40 cycles of 95 °C 10 s, and 60 °C 30 s. FNR1 was used as an internal control. Three replications were employed and the relative expression levels were calculated using the 2^−ΔΔCᴛ method^[1]. T-test was used to determine statistical differences (** p < 0.01, * p < 0.05). A p-value < 0.05 and absolute fold change ≥ 2.0 signified differential expression.

Since the ice plant is the only sequenced plant in the Aizoaceae family, to explore the evolution of PP2C genes, we selected several closely related species and model species, including O. sativa, V. vinifera, A. caudatus, H. undatus, B. vulgaris, and A. thaliana. A total of 59 McrPP2Cs were identified from the genome of M. crystallinum, and 115, 86, 69, 87, 49, and 80 were identified in O. sativa, V. vinifera, A. caudatus, H. undatus, B. vulgaris, and A. thaliana respectively (Supplementary Tables S1−S8). The protein sequences of PP2Cs varied greatly. Their amino acid numbers ranged from 78 to 1,832. The isoelectric points were between 3.75 and 11.86. Most PP2Cs (420/545) were acidic proteins. The largest molecular mass was 206,639.9, and the smallest was 8,163.06. Among the 545 PP2Cs, stable proteins accounted for approximately 29.9% (163), and unstable proteins accounted for approximately 70.1% (382). The hydrophilicity of the proteins ranged from 0.274 to −1.248, and the aliphatic index was between 53.04 and 106.75. Subcellular localization predictions suggested that most PP2Cs may be located in the nucleus, chloroplast, or cytoplasm (Supplementary Tables S1−S7).

Chromosomal localization of the PP2C gene family

In M. crystallinum, 59 McrPP2Cs were unevenly located on nine chromosomes (Fig. 1). Chr 3 had the largest number of McrPP2Cs (11), while Chr 7 had the fewest McrPP2Cs (3). Most McrPP2Cs were concentrated near the two ends of each chromosome while fewer genes were in the centromeric regions. There was no positive correlation between the length of a chromosome and the number of McrPP2Cs (Fig. 1; Supplementary Table S1). The distribution pattern of PP2Cs in other six studied plants was similar (Supplementary Tables S2−S8).

Figure 1. 

Chromosome distribution of the PP2C gene family in M. crystallinum. McrPP2Cs are marked on nine chromosomes and the chromosome length is demonstrated in the left scale. Different chromosomal colors indicate different gene densities, with red representing the highest density and blue representing the lowest density.

Phylogenetic analysis of the PP2C gene family

To clarify the evolutionary relationships of PP2Cs, the protein sequences of PP2Cs from the seven studied species were combined to construct a phylogenetic tree. According to the topological structure, 56 McrPP2Cs were divided into 12 clades (A−L), and three McrPP2Cs not assigned to any clade (Fig. 2, Table 1). Clade F had the largest number (9), followed by clade A (7). For the PP2Cs of the other six studied species, the clades E or F were the largest in most cases. The only exception was O. sativa, in which clade D had the largest number (Fig. 2, Table 1). Each clade contained PP2Cs from all seven species and genes from the same species tended to cluster together (Fig. 2, Table 1).

Figure 2. 

Phylogenetic analysis of the PP2C gene family from M. crystallinum, H. undatus, A. caudatus, B. vulgaris, V. vinifera, A. thaliana, and O. sativa. Phylogenetic trees were constructed using the PP2C protein sequences and different colors indicate different clades. PP2C genes in M. crystallinum are marked in red.

Amplification analysis of the PP2C gene family

Gene duplication promotes the expansion of gene families. We used MCScan to classify the gene types and divided all PP2Cs into five groups: singletons, dispersed genes, proximal duplicates, tandem duplicates, and polyploidy-related genes (Table 2). WGD contributed most to gene expansion in M. crystallinum (76.3%), V. vinifera (75.6%), and O. sativa (48.2%). In contrast, dispersed duplications contributed most to gene expansion in A. thaliana (51.3%), B. vulgaris (69.4%), and A. caudatus (52.2%). In H. undatus, the proportions of dispersed genes and WGD were the same, indicating their same contributions to gene expansion. Singletons of PP2Cs were not detected among these seven species (Table 2).

Collinearity and selection pressure analysis of the PP2C gene family

To investigate the evolutionary relationships of PP2Cs among different species, collinearity analysis was performed between M. crystallinum and the other six studied species to search for orthologous genes (Fig. 3a). The results showed that there were 80 pairs of orthologous genes between M. crystallinum and H. undatus, 95 pairs between H. undatus and A. caudatus, 57 pairs between A. caudatus and B. vulgaris, 55 pairs between B. vulgaris and V. vinifera, 80 pairs between V. vinifera and A. thaliana, and 23 pairs between A. thaliana and O. sativa (Supplementary Table S9). To explore the evolutionary patterns of McrPP2Cs, collinearity analysis was also used to search for paralogous genes. Six gene pairs were identified in the 59 McrPP2Cs, named as McrPP2C2/McrPP2C4, McrPP2C11/McrPP2C12, McrPP2C15/McrPP2C23, McrPP2C19/McrPP2C29, McrPP2C47/McrPP2C51, and McrPP2C55/McrPP2C56. The collinear genes were distributed on seven of nine chromosomes, except Chr4 and Chr6 (Fig. 3b). The divergence time between the six pairs of duplicated genes was calculated based on the mutation rate of 1.5 × 10⁻⁸ mutations per synonymous site per year (Supplementary Table S10). It was estimated that the divergence events occurred around 113.9−168.8 million years ago (Mya), very close to the whole-genome triplication event (WGT, 130−150 Mya).

Figure 3. 

Collinearity analysis of the PP2C gene family. (a) The collinear genes are shown by gray lines, while the syntenic PP2C homologous gene pairs are highlighted by blue lines. The species from top to bottom are M. crystallinum, H. undatus, A. caudatus, B. vulgaris, V. vinifera, A. thaliana, and O. sativa. (b) Collinearity analysis of McrPP2Cs. The gray lines represent all collinear genes in the ice plant, and orange lines represent the McrPP2C collinearity genes.

To further study the divergence and selection in the gene duplication of PP2Cs, the Ka/Ks values of the paralogous gene pairs in M. crystallinum were evaluated (Supplementary Table S11). The Ka/Ks values of all paralogous gene pairs in M. crystallinum and in the other six studied species were less than 1, demonstrating that these genes were mainly under purifying selection during the evolutionary process.

Analysis of the gene structure and conserved motifs of the PP2C gene family

Since the diversity of exon/intron structures plays an important role in the evolution of gene families, we analyzed the exon/intron structural patterns of McrPP2Cs based on their phylogenetic relationships. The number of exons in McrPP2Cs ranged from 2 to 19, and most members of the same clade possessed similar numbers of exons/introns with varied total gene lengths (Fig. 4). For instance, all clade A McrPP2Cs have either three or four exons, but McrPP2C17 was over 10 kb, more than 2 fold than other clade A members. The variation of gene length was mainly due to the different lengths of introns, while the length of exons was similar (Fig. 4), indicating that McrPP2Cs are relatively conserved during the evolutionary process, thus resulting in similar functions.

Figure 4. 

The motif and exon-intron structure analyses of the PP2C gene family in M. crystallinum.

To identify the common motifs in McrPP2Cs, we utilized the MEME motif search tool to determine 10 conserved motifs^[2]. The number of motifs contained in McrPP2Cs ranges from 3 to 10. Some motifs exhibited a certain degree of conservation. Motif 3 could be found in all McrPP2Cs; motif 1 did not exist only in McrPP2C53; motif 6 did not existed in McrPP2C22 and McrPP2C39 (Fig. 4; Supplementary Table S12). McrPP2Cs in the same branch had similar distribution orders and numbers of motifs (Fig. 4; Supplementary Table S12).

Prediction of cis-acting elements of the PP2C gene family

To further clarify the potential functions of McrPP2Cs, the cis-acting elements in the 2-kb upstream regions were analyzed. The results showed that McrPP2Cs contained four types of cis-acting elements: light-responsive elements, hormone-responsive elements, defense and stress-responsive elements, and growth and development-related elements (Fig. 5; Supplemental Table S13). Among them, light-responsive elements were present in all McrPP2Cs and were the most abundant in number, suggesting that McrPP2Cs are closely related to light regulation in M. crystallinum. The hormone responsive elements of McrPP2Cs included abscisic acid (ABA)-responsive elements, gibberellin-responsive elements, and methyl jasmonate (MeJA)-responsive elements. Among them, ABA-responsive elements existed in most McrPP2Cs. Only 11 of 59 McrPP2Cs did not contain ABA-responsive elements. McrPP2C2 and McrPP2C11 had the largest number of ABA-responsive elements. Besides, most McrPP2Cs contained MeJA-responsive elements. McrPP2C11 had the largest number (12). The environmental pressure response elements of McrPP2Cs consisted of low temperature, drought, anaerobic induction, and stress response elements. Most McrPP2Cs contained anaerobic induction elements and more than half McrPP2Cs contained drought induction elements. In addition, the promoters of McrPP2Cs genes also contained elements related to plant growth and development (Fig. 5; Supplementary Table S13). These results suggested that McrPP2Cs may be widely involved in the stress response, as well as in plant growth and development.

Figure 5. 

Cis-elements analysis of PP2C genes in M. crystallinum. Different colors indicate different components. The vertical bar indicate the number of elements.

Expression profiling of McrPP2Cs under salt stress

Ice plant can grow under high concentrations of NaCl^[29−32]. One reason is that it can shift from C3 to CAM under salt stress^[30−32]. Since this transition is related to the McrPEPC/08G01316 induction in leaf^[29−32], we examined the expression profiling of McrPP2Cs under salt stress before and after the McrPEPC was induced. According to previous studies, we selected 14 d as a long time salt treatment, under which McrPEPC is fully induced^[32]; and we selected 36 h as a short time salt treatment, under which McrPEPC cannot be induced^[30]. Transcriptomic data confirmed that only long time salt treatment induced McrPEPC/Mc08G01316 (Supplementary Table S14). However, the expression changes of McrPP2Cs were similar under both long and short time treatments. Compared with the control, only four McrPP2Cs (McrPP2c11/Mc02G00242, McrPP2c15/Mc03G00348, McrPP2c17/Mc03G00607, McrPP2c35/Mc05G02416) were up-regulated after both long and short time salt treatments, while most McrPP2Cs did not show significant changes (Fig. 6a, b; Supplementary Tables S15, S16). Interestingly, all four up-regulated genes belonged to clade A. QRT-PCR was used to explore the expression pattern of all seven clade A McrPP2Cs and the results were consistent with the RNA-seq data (Fig. 6e).

Figure 6. 

Expression pattern analysis of PP2C genes in M. crystallinum. Expression pattern analysis of McrPP2Cs in leaf under (a) 14 d and (b) 36 h salt stress. CK indicates the corresponding control group. Grey bars mean the gene was not detected. (c) Expression pattern analysis of clade A McrPP2Cs in root under different NaCl treatments. (d) Expression pattern analysis of clade A McrPP2Cs in leaf and root under unstressed condition. (e) qRT-PCR analysis of clade A McrPP2Cs in leaf under long and short salt stress. CK L: 14 d salt stress control; NaCl L: 14 d salt stress; CK S: 36 h salt stress control; NaCl S: 36 h salt stress. Compared with corresponding control, a p-value < 0.05 and absolute fold change ≥ 2.0 signified differential expression (* p < 0.05, ** p < 0.01).

To better understand the expression pattern of clade A McrPP2Cs under salt stress, we reanalyzed the published RNA-seq data in root^[29,31]. In root, six clade A McrPP2Cs were induced under 500 mM NaCl treatment except McrPP2C34/Mc05G01962 (Fig. 6c; Supplementary Table S17). We also explored the expression of clade A McrPP2Cs in vegetable tissues using published data^[31,42]. We found that all clade A McrPP2Cs expressed higher in root than in leaf, although their expression level varied (Fig. 6d; Supplementary Table S18). The expression level of McrPP2C2/Mc01G00236 was much higher than the other six McrPP2Cs in both leaf and root (Fig. 6d; Supplementary Table S18). Primers are shown in Supplementary Table S19.

PP2Cs are widely identified in plants and play important roles in abiotic stress^[1−5]. This study provides the first comprehensive analysis of the PP2C gene family in ice plants at the whole genome level. Finally, 59 McrPP2Cs were identified and most of them were clustered at the ends of chromosomes (Fig. 1; Supplementary Table S1). This uneven distribution might be because more repetitive sequences exist in the centromere than the ends. Similar distributions were also found in in the other six studied and some reported plants, including Arabidopsis, rice, wheat, peanut, cucumber, strawberries, tomato, soybean, sweet cherry, etc (Supplementary Table S8)^[1,2,6−15], suggesting the genes uneven distribution is prevalent in higher plants. This phenomenon indicated that the recombination active chromosome arms tend to retain stress-responsive genes, while the centromere regions enriched in repetitive sequences showed selective limitation of gene retention. Phylogenetic analysis divided the McrPP2Cs into 12 clades (A-L) and genes from the same species tended to cluster together (Fig. 2). This means PP2Cs retain their core functions while undergoing differentiation at the species level, which may be driven by adaptation to unique environmental challenges. In addition, most of the genes on the same chromosome did not belong to the same subfamily in the phylogenetic tree (Figs 1, 2), suggesting that genes on the same chromosome perform different functions.

In ice plants, the largest McrPP2C clade was F, which had nine genes (Fig. 2, Table 1). This clade is also the largest in many other plants such as Arabidopsis, grape, wheat, B. vulgaris, H. undatus, etc (Fig. 2, Table 1)^[7]. It is reported that clade F PP2C negatively regulates stomata closure under drought stress in maize^[22]. Thus, we speculated that plants developed redundant clade F PP2Cs to regulate the stomata movement during stress responses. Since ice plant converts from C3 photosynthesis to CAM under high salinity stress, the stomata convert from a night-closure to a night-opening pattern^[30,31]. Therefore, the large clade F McrPP2Cs may be involved in this stomata conversion. However, the expression pattern of clade F McrPP2Cs was not changed after salt treatment in leaf (Fig. 6a, b; Supplementary Tables S15, S16), implying that there was post-transcriptional regulation on clade F McrPP2Cs.

In plants, gene duplication promotes the expansion of gene families. Whole-genome duplication (WGD) contributed the most to the expansion of PP2Cs family in M. crystallinum, V. vinifera, and O. sativa (Table 2). Our previous study shows that M. crystallinum has only experienced one whole-genome triplication (WGT) event^[31]. The differentiation time (113.9−168.8 Mya) of six pairs of paralogous McrPP2Cs was highly consistent with their WGT events (approximately 130−150 Mya), strongly supporting that WGT was a key driving force for PP2C family expansion (Table 2). Because WGD events are often accompanied by large-scale gene loss and gene recombination^[43], the number of McrPP2Cs may be underestimated. However, compared with the 23,183 annotated genes in M. crystallinum, the 59 McrPP2Cs retention is still extremely significantly higher than expected by chance (Chi-square test p-value = 8.22446E-06 < 0.05) (Table 2)^[44], suggesting the strong selection pressure towards retaining these genes, which may be attributed to the strong stress tolerance in ice plant. In contrast, the expansion of the PP2C family in A. thaliana and B. vulgaris is mainly through dispersed replication, reflecting the differences in evolutionary strategies of species with different polyploidization histories.

Collinearity analysis and a previous study demonstrated that 46.3% of the orthologous genes between M. crystallinum and H. undatus appeared before the species divergence (~77.8 Mya)^[31], indicating that the PP2C family existed in the common ancestor of M. crystallinum and H. undatus and had already formed a certain evolutionary basis before species differentiation. Two up-regulated clade A McrPP2Cs (McrPP2C11/Mc02G00242, McrPP2C15/Mc03G00348) had colinear genes, McrPP2C12/Mc02G00368 and McrPP2C23/Mc03G2217 respectively (Fig. 3b). But expression level of the colinear genes were not changed under salt stress (Fig. 6a, b). This different expression pattern may be caused by subfunctionalization or neofunctionalization after gene duplication. The genome retains redundant genes to meet complex regulatory requirements. This dynamic evolutionary strategy may be an important mechanism for salt tolerance in ice plants. The Ka/Ks values of all McrPP2C paralogous gene pairs were far less than 1 (Supplementary Table S11), indicating that they were subjected to purify selection, which might help maintain the basic functions of PP2C genes.

Gene structure and expression analysis found that differences in exon-intron structure and motif composition were not related to the expression pattern. For example, there was one more exon in McrPP2C15/Mc03G00348 compared with other clade A McrPP2Cs, but its expression pattern was similar to the other three clade A McrPP2Cs (McrPP2C11/Mc02G00242, McrPP2C17/Mc03G00607, McrPP2C35/Mc05G02416) in leaf and root under salt stress (Figs 4, 6). Promoter analysis showed that McrPP2Cs were rich in light-responsive, hormone-responsive, and drought-responsive elements. 81.4% of McrPP2Cs contained ABA-responsive elements (ABREs) (Fig. 5; Supplementary Table S13). This is consistent with their stress responsive functions. Besides, MeJA responsive elements were also widely spread in McrPP2Cs promoter (Fig. 5; Supplementary Table S13). In Arabidopsis, JA level is elevated in response to mites and necrotrophic pathogens^[45]. Therefore, McrPP2Cs may be the cross-linking of biotic and abiotic stress signals. It is worth noting that McrPP2C2/Mc01G01236 and McrPP2C11/Mc02G00242 owned the highest number of ABREs, 9 and 10 respectively (Fig. 5; Supplementary Table S13). Both of them belong to clade A McrPP2Cs. Since this clade is tightly related to ABA signal^[4], these two genes might be key regulatory factors in the ABA dependent pathway. McrPP2C11 could be induced in leaf and root by NaCl, whereas McrPP2C2 was only induced by NaCl in root (Fig. 6a−c). This may be because the expression level of McrPP2C2 was significantly higher in root compared with that in leaf (Fig. 6d). In addition, the expression level of McrPP2C2 was much higher than the other six clade A McrPP2Cs in both leaf and root (Fig. 6d) indicating it may play other roles besides ABA-dependent salt stress response in ice plants.

In our study, more than half clade A McrPP2Cs were up-regulated under salt stress (Fig. 6; Supplementary Tables S15−S18). McrPP2C11, 15, 17, 35 were induced in both leaf and root. McrPP2C2, 23 were only induced in root. The reason that more clade A McrPP2Cs were induced in root may be because the root directly contacts with the salinity. Under salt stress, ice plant shifts from C3 photosynthesis to CAM, which is related to the key enzyme McrPEPC/Mc08G01316 induction^[30,31,42]. Our results showed that four clade A McrPP2Cs were induced earlier than McrPEPC in leaf, and the induction could be maintained for a long time (Fig. 6e; Supplementary Table S14). The quick response of McrPP2Cs suggested they may be involved in McrPEPC induction. Since clade A PP2Cs participate in the ABA signaling pathway, we believe that the ABA signal plays a central role in the C3-CAM transition. But the detailed molecular mechanism is unclear and needs to be further investigated.

In this study, we conducted a genome-wide identification and analysis of the PP2Cs in ice plants. A total of 59 McrPP2Cs were identified which were unevenly distributed across nine chromosomes with clustering near chromosome ends. Phylogenetic analysis classified PP2Cs of seven plant species into 12 clades (A−L), with McrPP2Cs distributed across all clades. A whole-genome triplication event (WGT, ~130−150 Mya) was identified as the primary driver of McrPP2C expansion, with 76.3% of gene duplication events attributed to WGT. Purifying selection dominates paralogue evolution. By predicting the cis-elements in the promoter, we found that 81.4% of McrPP2Cs contained ABA-responsive elements (ABREs). Under salt treatments, four clade A McrPP2Cs (McrPP2c11/15/17/35) displayed significant up-regulation in leaf and root, while two clade A McrPP2Cs (McrPP2c2/23) were only up-regulated in root. These results provide further insights into the ice plant PP2C gene family and indicate that they are potentially involved in the response to salt stress, which gives some clues for understanding why ice plant has strong salt tolerance.