# Genome-Wide Identification, Evolution, and Transcriptional Profiling ofPP2C Gene Family inBrassica rapa

Authors: Nadeem Khan; Han Ke; Chun-mei Hu; Emal Naseri; Muhammad Salman Haider; Aliya Ayaz; Waleed Amjad Khan; Jianjun Wang; Xilin Hou
Journal: BioMed Research International
 (2019)

---

## Abstract

The type 2C protein which belongs to the major group of protein phosphatases (PP2C) plays a vital role in abscisic acid (ABA) signaling and signal transductions processes. In the present study, 131PP2C genes were identified in total inBrassica rapa and categorized into thirteen subgroups based on their phylogenetic relationships. TheseB. rapa PP2C are structurally conserved based on amino acid sequence alignment, phylogenetic analysis, and conserved domains. Moreover, we utilized previously reported RNA-sequence data on various tissues (root, stem, leaf, flower, and silique), which suggests overlapping expression pattern in 29 paralogous gene pairs. The qRT-PCR validation of 15 paralogous gene pairs depicts distinct expression patterns in response to various abiotic stresses, such as heat, cold, ABA, and drought. Interestingly, stress-responsiveBraPP2C candidate genes were also identified, suggesting their significance in stress-tolerance mechanism inB. rapa. The evolutionary analysis for 15 paralogous gene pairs suggested that only three pairs have the positive selection and remaining were purifying in nature. The presented results of this study hasten our understanding of the molecular evolution of thePP2C gene family inB. rapa. Thus, it will be ultimately helping in future research for facilitating the functional characterization ofBraPP2C genes in developing the abiotic stress tolerant plants.


---

## Body

## 1. Introduction

In cellular signaling, the reversible protein phosphorylation catalyzed by protein kinases (PKs) and phosphatases (PPs) is known to participate in critical processes [1]. In plants, the reversible protein phosphorylation is one of the critical modification processes, which helps in the regulation of some important physiological and biochemical reactions. The PPs provide modulations of protein phosphoregulation by reversing the action of PKs. The PPs can be categorized into two major classes, such as protein tyrosine phosphatase (PTPs) and protein serine/threonine phosphatases (PSPs) based on the substrate specificity [2]. Moreover, most of the PSPs are further subcategorized into two groups; group one includes protein phosphatase 1 (PP1), PP2A, PP2B, PP4, PP5, and PP6, while group two contains protein phosphatase M (PPM), including the PP2C and also pyruvate dehydrogenase phosphatases [2].From prokaryotes to higher eukaryotes, the PP2C are evolutionarily conserved and mostly found in archaea, bacteria, fungi, plants, and animals [3, 4]. Furthermore, it acts as a negative modulator of PKs cascades activated and then implicated in regulating stress-signaling pathways. Like PKs, PPs also share a crucial role under various stress conditions and might be pivotal at different developmental stages of plants and signal transduction processes. Particularly inArabidopsisandBrassica raparecently, several groups of the researcher are engaged by deciphering the involvement of different kinases against biotic and abiotic stress networks such as MAPKs [5, 6], CDPKS [7, 8], and CIPKs [9, 10]. For instance, temperature (high and low) stresses are the major environmental factor that limits crop productivity in vegetable production [11, 12]. Generally, during cold treatment in which plants adjust their metabolism are considered as cold acclimation. To date, several studies highlighted the importance of PP2C that may play pivotal roles in various processes, including both biotic-abiotic stress factors and plant development [13, 14]. In plants, PP2C members have been described as a regulator of desiccation tolerance and act as a negative regulator in the ABA-signaling pathway [15–17]. For example, in moss, the only twoPpABI1A andPpABI1Bbelonging to subfamily A of PP2C are known to induce desiccation tolerance in vegetable and are directly involved in ABA signaling [18]. In addition, drought stress considerably affects the crop productivity with an ultimate decline of photosynthetic assimilates, due to the osmotic stress-imposed constraints. Plants can utilize adaptive strategies against drought stress, such as escape, avoidance, and tolerance [19]. The validation of maizeZmPP2C-A10 by transgenic studies has confirmed its negative regulation in drought stress tolerance [20]. In particular, the twoPP2C genes encoding ABA-insensitive mutants (abi1-1 and abi2-1) are known to partake in the various physiological processes after exposure to abiotic stimuli, including salt, drought, and freezing [21–24]. Plant development and freezing tolerance were accelerated by downregulation ofPP2CA genes and the antisense mediation [25]. Hence, for a better understanding of protein phosphatases and their functional significance, the identification of PP2C provides a stepping stone in stress-signaling pathways.Brassica rapa, a diploid species, share a complex history with model plantArabidopsis and has experienced two main duplication events such as (WGD α and β) and additional one whole-genome triplication event (WGT γ) 13–17 million years ago (MYA) [26, 27]. Thus, theB. rapadue to WGT event has experienced considerable fractionation in the genome (i.e., duplicate gene loss); as a result, it presents to us an opportunity forPP2C genes to study its evolutionary implications. In this study, all thePP2C genes through various approaches were isolated fromArabidopsis andB. rapaand their phylogenetic relationships were conducted to categorize them into multiple subgroups. In addition, we comparatively analyzedPP2C genes from ten representative plant species, including one of the earliest angiosperm plantAmborella trichopoda, Populus trichocarpa, Vitis vinifera, Capsella rubella, Citrus sinensis, Carica papaya, Solanum lycopersicum, Fragaria vesca, Arabidopsis lyrata,and a bryophyte (Physcomitrella patens), the earliest sequenced land plant [28].We further compared thePP2C genes betweenB. rapaandArabidopsis,to explore and identify both shared and specific subgroups. Following the gene structure organization analysis, conserved protein motifs, cis-elements, and interaction network, we traced paralogous gene pairs and their evolutionary divergence that likely resulted in the expansion of thePP2C gene family. To shed a light of some criticalBraPP2C genes on their associated indigenous functional role were further exposed to heat, cold, ABA-signaling, and drought stress conditions. Additionally, the transcriptional profiling of theBraPP2Cgenes for various tissues and qRT-PCR analysis of 15 paralogous gene pairs were analyzed and compared. In both eukaryotes and prokaryotes, the PPs are vital in both diverse signaling pathways and kinase-counteracting components. Indeed, for phosphatases, most of the functionally characterized literatures have been studied in Arabidopsis, and minimal studies have been done in crop plants. Therefore, to our knowledge, this is the first systemic reported genome-wide and transcriptional profiling of thePP2C genes inB. rapa, and it is affirmative to carry out a comprehensive study to understand the regulation of phosphatases inB. rapaduring stress and development.


## 2. Materials and Methods

### 2.1. Data Resources for Sequence Retrieval

For identification ofPP2C genes inBrassica rapa, we utilized the Plaza 4.0 database (https://bioinformatics.psb.ugent.be/plaza/) with the help of InterPro PP2C domain “IPR001932”. TheB. rapa genome sequences were downloaded from BRAD (http://brassicadb.org/brad/) [29], andA. thaliana sequences were retrieved from TAIR (http://www.arabidopsis.org/). The sequences of others species were downloaded from Phytozome v12.1.6 (https://phytozome.jgi.doe.gov/pz/portal.html) [30]. The domains of obtained BraPP2C proteins were also further verified using NCBI-Conserved Domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) search program and SMART databases (http://smart.embl-heidelberg.de/) [31]. Those proteins which lack PP2C domain were removed from further analysis. In addition, protein sequences that were found with obvious errors in their gene length or having less than 100 lengths were eliminated.


### 2.2. Multiple Sequence Alignment and Phylogenetic Analysis

The amino acid sequences of the BraPP2C proteins were used for further investigation, and multiple sequence alignment was performed by MUSCLE [32] using MEGA 7 software with the default options [33]. The phylogenetic trees were constructed using the maximum likelihood (ML) method. To determine the reliability of resulting tree, bootstrap values of 1000 replications were performed with the Jones, Taylor, and Thornton amino acid substitution model (JTT model), while keeping the other parameters as a default.


### 2.3. Calculation of theKa/Ks for BraPP2C Paralogous Gene Pairs

TheKa/Ks ratios were calculated for BraPP2C paralogous gene pairs using MEGA 7.0 [33]. The divergence time among the paralog pairs was calculated with the following formula: T =Ks/2r, whereKs represents the synonymous substitutions per site and r is the rate of divergence. For dicotyledonous plants, specificallyB.rapa, the assumption is 1.5 synonymous substitutions per site of 108 years [34].


### 2.4. Conserved Motifs, Exon-Intron Structure Analysis, and Physicochemical Parameters of BraPP2C Proteins

Conserved motif scanning of BraPP2C proteins was carried out through local MEME Suite Version 5.0.2. For this purpose, parameter settings were calibrated as follows: maximum number of motifs 10, with a minimum width of 100 and a maximum of 120. The other parameters were set as default [35]. For the exon-intron structure, we used the Gene Structure Display Server (GSDS 2.0) (http://gsds.cbi.pku.edu.cn) [36]. The physicochemical properties of the proteins, including molecular weight (MW), isoelectronic points (pI), and GRAVY values for each gene, were calculated using the ExPASY PROTPARAM tools (http://web.expasy.org/protparam/). The subcellular localization was predicted using the WOLF PSORT (https://wolfpsort.hgc.jp/) website.


### 2.5. Cis-Elements Predictions and Interaction Network of BraPP2C

Every BraPP2C promoter sequence (selected as 2000 upstream bp) was imported in Generic File Format (GFF) file from theB.rapagenome. Then, PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [37] was utilized to identify the cis-regulatory elements for promoters of each gene. The interaction network of BraPP2C proteins was constructed with the help of STRING software (https://string-db.org/).


### 2.6. Chromosomal Location and Paralogous Gene Pairs Identification of BraPP2C

The chromosomal locationBraPP2C genes were illustrated from top to bottom concerning their position in the genome annotation using Mapchart [38]. For synteny gene analysis, the relationships were verified between the homologs ofA. thaliana and subgenomes ofB. rapa (LF, MF1, and MF2) obtained from BRAD (http://brassicadb.org/brad/searchSynteny.php). Furthermore, the paralogous genes were identified either selecting genes pair between LF1 and MF1 or MF1 and MF2 subgenomes ofBrassica rapa, respectively. Circos program was applied to demonstrate the syntenic relationships among the chromosomes ofB. rapa andA. thaliana [39].


### 2.7. Pearson Correlation Analyses

Pearson correlation (PCC) analysis was performed with the help of Excel 2013 to evaluate the PCC values of the RNA-seq and the paralogous genes that were used for qRT-PCR [40].


### 2.8. Plant Material and Treatments

In the present study, the germinated seeds of Chinese cabbage (Chiifu-401-42) were grown in plastic pots containing a mixture of soil and vermiculite (3:1). The pots were then placed in an artificial growth chamber for five weeks. The growth conditions were as follows: the temperature was set to 24/16°C, the photoperiod was 16/8 h, and the relative humidity was 65–70%. Specific treatments were provided to the seedlings as follows: for heat and cold treatments, seedlings were exposed to 38°C and 4°C, respectively. For ABA and drought stress treatment, the seedlings were cultured in a nutrient solution medium with 100μM and 6000 PEG (w/v). All treatments were carried out in continuous time intervals of 1, 6, and 12 h, respectively, with biological triplicates. After that leaf samples were quickly frozen in liquid nitrogen and stored at −80°C for further use.


### 2.9. RNA Isolation and Transcriptional Profiling of BraPP2C under Various Stresses

Total RNA was isolated from the treated frozen leaves with Trizol (Invitrogen) following the manufacturer’s instructions. RNA was reverse-transcribed into cDNA using the Primer Script RT reagent kit (TAKARA, Dalian China) according to their instructions. Specific primers were designed using Becan Designer 7.9 and are presented in TableS1. In order to check the specificity of the primers, we used the BLAST tool against theB. rapagenome for confirmation. RT-PCR was performed according to the guidelines of previous studies [41]. Relative fold expression was calculated with the comparative Ct-method. The expression patterns of allBraPP2C genes were analyzed based on a previous study [42]. Furthermore, gene expression levels were quantified by FPKM (fragments per kilobase of transcript per million fragments mapped) values, and heat maps were generated using an online omicshares tool (http://www.omicshare.com/).


## 2.1. Data Resources for Sequence Retrieval

For identification ofPP2C genes inBrassica rapa, we utilized the Plaza 4.0 database (https://bioinformatics.psb.ugent.be/plaza/) with the help of InterPro PP2C domain “IPR001932”. TheB. rapa genome sequences were downloaded from BRAD (http://brassicadb.org/brad/) [29], andA. thaliana sequences were retrieved from TAIR (http://www.arabidopsis.org/). The sequences of others species were downloaded from Phytozome v12.1.6 (https://phytozome.jgi.doe.gov/pz/portal.html) [30]. The domains of obtained BraPP2C proteins were also further verified using NCBI-Conserved Domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) search program and SMART databases (http://smart.embl-heidelberg.de/) [31]. Those proteins which lack PP2C domain were removed from further analysis. In addition, protein sequences that were found with obvious errors in their gene length or having less than 100 lengths were eliminated.


## 2.2. Multiple Sequence Alignment and Phylogenetic Analysis

The amino acid sequences of the BraPP2C proteins were used for further investigation, and multiple sequence alignment was performed by MUSCLE [32] using MEGA 7 software with the default options [33]. The phylogenetic trees were constructed using the maximum likelihood (ML) method. To determine the reliability of resulting tree, bootstrap values of 1000 replications were performed with the Jones, Taylor, and Thornton amino acid substitution model (JTT model), while keeping the other parameters as a default.


## 2.3. Calculation of theKa/Ks for BraPP2C Paralogous Gene Pairs

TheKa/Ks ratios were calculated for BraPP2C paralogous gene pairs using MEGA 7.0 [33]. The divergence time among the paralog pairs was calculated with the following formula: T =Ks/2r, whereKs represents the synonymous substitutions per site and r is the rate of divergence. For dicotyledonous plants, specificallyB.rapa, the assumption is 1.5 synonymous substitutions per site of 108 years [34].


## 2.4. Conserved Motifs, Exon-Intron Structure Analysis, and Physicochemical Parameters of BraPP2C Proteins

Conserved motif scanning of BraPP2C proteins was carried out through local MEME Suite Version 5.0.2. For this purpose, parameter settings were calibrated as follows: maximum number of motifs 10, with a minimum width of 100 and a maximum of 120. The other parameters were set as default [35]. For the exon-intron structure, we used the Gene Structure Display Server (GSDS 2.0) (http://gsds.cbi.pku.edu.cn) [36]. The physicochemical properties of the proteins, including molecular weight (MW), isoelectronic points (pI), and GRAVY values for each gene, were calculated using the ExPASY PROTPARAM tools (http://web.expasy.org/protparam/). The subcellular localization was predicted using the WOLF PSORT (https://wolfpsort.hgc.jp/) website.


## 2.5. Cis-Elements Predictions and Interaction Network of BraPP2C

Every BraPP2C promoter sequence (selected as 2000 upstream bp) was imported in Generic File Format (GFF) file from theB.rapagenome. Then, PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [37] was utilized to identify the cis-regulatory elements for promoters of each gene. The interaction network of BraPP2C proteins was constructed with the help of STRING software (https://string-db.org/).


## 2.6. Chromosomal Location and Paralogous Gene Pairs Identification of BraPP2C

The chromosomal locationBraPP2C genes were illustrated from top to bottom concerning their position in the genome annotation using Mapchart [38]. For synteny gene analysis, the relationships were verified between the homologs ofA. thaliana and subgenomes ofB. rapa (LF, MF1, and MF2) obtained from BRAD (http://brassicadb.org/brad/searchSynteny.php). Furthermore, the paralogous genes were identified either selecting genes pair between LF1 and MF1 or MF1 and MF2 subgenomes ofBrassica rapa, respectively. Circos program was applied to demonstrate the syntenic relationships among the chromosomes ofB. rapa andA. thaliana [39].


## 2.7. Pearson Correlation Analyses

Pearson correlation (PCC) analysis was performed with the help of Excel 2013 to evaluate the PCC values of the RNA-seq and the paralogous genes that were used for qRT-PCR [40].


## 2.8. Plant Material and Treatments

In the present study, the germinated seeds of Chinese cabbage (Chiifu-401-42) were grown in plastic pots containing a mixture of soil and vermiculite (3:1). The pots were then placed in an artificial growth chamber for five weeks. The growth conditions were as follows: the temperature was set to 24/16°C, the photoperiod was 16/8 h, and the relative humidity was 65–70%. Specific treatments were provided to the seedlings as follows: for heat and cold treatments, seedlings were exposed to 38°C and 4°C, respectively. For ABA and drought stress treatment, the seedlings were cultured in a nutrient solution medium with 100μM and 6000 PEG (w/v). All treatments were carried out in continuous time intervals of 1, 6, and 12 h, respectively, with biological triplicates. After that leaf samples were quickly frozen in liquid nitrogen and stored at −80°C for further use.


## 2.9. RNA Isolation and Transcriptional Profiling of BraPP2C under Various Stresses

Total RNA was isolated from the treated frozen leaves with Trizol (Invitrogen) following the manufacturer’s instructions. RNA was reverse-transcribed into cDNA using the Primer Script RT reagent kit (TAKARA, Dalian China) according to their instructions. Specific primers were designed using Becan Designer 7.9 and are presented in TableS1. In order to check the specificity of the primers, we used the BLAST tool against theB. rapagenome for confirmation. RT-PCR was performed according to the guidelines of previous studies [41]. Relative fold expression was calculated with the comparative Ct-method. The expression patterns of allBraPP2C genes were analyzed based on a previous study [42]. Furthermore, gene expression levels were quantified by FPKM (fragments per kilobase of transcript per million fragments mapped) values, and heat maps were generated using an online omicshares tool (http://www.omicshare.com/).


## 3. Results

### 3.1. Characterization ofBraPP2C Gene Family Members

In the present study, we identified 131 putativeBraPP2C genes inB. rapa genome and were denoted as BraPP2C1 to BraPP2C131, based on phylogenetic analysis and their orthologous positions withArabidopsis. The physicochemical feature of theBraPP2C proteins along with some key information of all the identified BraPP2C is provided in Table S2. The results of sequence analysis exhibited that the length of the BraPP2C proteins varied with a ranged from 102 (BraPP2C1) to 1465 (BraPP2C63) bp, with an average of 412.37 (aa). The detailed information of BraPP2C, including molecular weights (MW), theoretical isoelectric point (pI), the grand average of hydropathicity (GRAVY), and exon number ranged from 11.76 (BraPP2C27) to 160.43 (BraPP2C63) kDa, -0.777 (BraPP2C59) to 0.001 (BraPP2C41), 4.36 (BraPP2C85 to 9.37 (BraPP2C117), and 1 (BraPP2C73) to 21 (BraPP2C21), respectively. The significant variations reflect their functional diversity among BraPP2C, while the negative value of GRAVY further demonstrated that these proteins were unstable and hydrophilic in nature. Moreover, we analyzed the values of pI, GRAVY, and exon number, as shown in (Supplementary Figures 1, 2, and 3), respectively. The subcellular localization was predicated, and results proposed that most of the BraPP2C proteins were localized in mitochondria, cytoplasm, nuclei, chloroplasts, plasma membranes, endoplasmic reticulum, and vacuoles (Table S2). The gene structure analysis of BraPP2C revealed that the exon number of eachBraPP2C gene ranged dramatically from 1 to 19 (Figure S4). Majority ofBraPP2C genes consisted of only one exon; the rest of others showed significant variations, whereas BraPP2C32 possesses 19 exons, suggesting that both exon loss and gain occurred inBraPP2C gene family. Alongside, we systematically discovered the distribution of conserved motifs and their logos by online MEME server. A total of ten different consensus motifs were obtained in all of the BraPP2C proteins, and their distribution patterns are presented in Figures S5 and S6. Hence, during the evolutionary process, the BraPP2C demonstrated the extreme conservation within subgroups.


### 3.2. Chromosome Distribution and Syntenic Gene Collinearity Analysis

The physical chromosomal locations onB. rapa, across all ten chromosomes (Br01–Br10), were determined by the Mapchart software. The 126BraPP2C genes were dispersed across all ten chromosomes, ranging from 7 to 19 per chromosome, out of which 5 genes were found on the scaffold region (Figure 1). On each chromosome, the number of BraPP2C varies drastically; the largest number of BraPP2C family members was observed on chromosomes Br05 with 19 genes, followed by Br03 with 18 genes, whereas the least numbers were revealed on chromosomes Br04 and Br07; each contains 7 genes, respectively. Among these BraPP2C, five genes (BraPP2C28,BraPP2C68,BraPP2C80,BraPP2C83, andBraPP2C111) were present on the genomic scaffold. These results suggest the uneven distribution patterns of BraPP2C across ten different chromosomes. Based on a previously reported study, 24 ancestral genomic blocks (GBs) were reconstructed [27]. The color-coding for these block were slightly modified according to their position in a proposed ancestral karyotype (AK1-07) [27, 43]. The majority (18.18%) ofBraPP2C genes were clustered in the AK5 region, followed by AK1 and AK2 (16.88%), whereas minimum genes (6.49%) were clustered in AK7 region (Figure 1).Figure 1
Chromosome locations of PP2C were obtained from the GFF file and displayed using Mapchart. The ancestral karyotypes are marked in different colors.The collinear relationships of thePP2C gene pairs betweenArabidopsisandB. rapa are shown in Figure 2. Based on three subgenomes ofB. rapa, LF contains more (42.86%)BraPP2C genes than the MF1 (32.77%) and MF2 (24.37%) subgenomes (Figure S7). In addition, the divergence time of the 29 paralogous pairs was estimated with the help of MEGA7.0 software by calculating the synonymous (Ks) and nonsynonymous substitution rates (Ka). The results proposed that most of theBraPP2C gene pairs represented less than 1.00 (Ka/Ks) ratios, suggesting the purifying selection of these genes except three pairs (BraPP2C58-BraPP2C56,BraPP2C85-BraPP2C86, andBraPP2C125-BraPP2C127) showed higher than 1.00 values, indicating positive selection. To estimate the divergence among theseBraPP2C paralogous gene pairs, we calculated the values ofKs (Figure S8 and Table S3). The average divergence time among these paralogous is 12.37 MYA, suggesting that BraPP2C divergence occurred along with theArabidopsis (9.6–16.1 MYA) [6].Figure 2
Phylogenetic relationship of PP2C betweenB. rapa andA. thaliana. The phylogenetic tree was constructed by MEGA 7 using the Maximum Likelihood Method (1000 bootstrap). Genes of PP2C paralog pairs are marked with red color. The collinear correlation for all genes of PP2C is displayed betweenB. rapa andA. thaliana. The 10 Chinese cabbage chromosomes (Br01–Br10) and fiveA. thaliana chromosomes (AT01–AT05) are shown in different random colors. The illustration was drawn using Circos Software.


### 3.3. Phylogenetic Analysis ofBraPP2C Gene and Copy Number Variations/Gene Retention

To provide an overview of the evolutionary relationships among subgroup member of PP2C betweenB. rapa andA. thaliana, a comparative phylogenetic tree was constructed using MEGA7.0 software by adapting the maximum likelihood approach with 1000 bootstrap replications (Figure S9). The phylogenetic tree and domain composition further subcategorized the PP2C into thirteen subgroups: subgroup A-L and one unclassified. The results of the tree were consistent with the previously reported studies inArabidopsis and rice [44]. Moreover, most of the BraPP2C clustered together with those fromArabidopsis, while 128 out of 131BraPP2C genes are further distributed among twelve subfamilies (A-L), and the remaining three (BraPP2C35,BraPP2C60, andBraPP261)BraPP2C genes cannot be grouped into any subfamilies and were categorized as unclassified. Most of the PP2C members were found in two subgroups E and F each with 25 genes, respectively, as compared to other subgroups. The complete overview and the distribution patterns amongBrassica rapa and various ten species are described in Table 1.Table 1
The distribution ofPP2C genes family is various species.
Subgroup of PP2C BraPP2C AthPP2C AtrPP2C CrPP2C CpPP2C FvPP2C PpPP2C VvPP2C SlPP2C AlPP2C CsPP2C PtPP2C A 15 10 2 10 4 6 0 5 6 11 11 11 B 7 6 1 4 3 1 0 0 3 2 6 5 C 10 7 4 8 3 4 0 5 5 6 18 13 D 16 9 4 11 6 8 4 9 12 9 15 29 E 25 12 4 16 6 8 14 5 9 13 18 26 F 25 13 6 16 3 4 0 8 8 12 23 25 G 9 6 3 8 2 7 0 5 6 6 12 14 H 6 3 3 4 3 3 28 3 5 3 3 15 I 4 2 0 3 2 1 19 2 3 2 10 4 J 1 2 0 1 0 0 0 0 1 1 0 1 K 6 3 1 4 3 3 0 3 2 4 5 7 L 4 0 2 1 3 1 0 2 2 2 6 4 Unclassified 3 3 1 2 0 1 0 1 1 2 0 2 ∗ BraPP2C= Brassica rapa, AthPP2C= Arabidopsis thaliana, AtrPP2C= Amborella trichopoda, CrPP2C= Capsella rubella, CpPP2C= Carica papaya, FvPP2C= Fragaria vesca, PpPP2C= Physcomitrella patens, VvPP2C= Vitis vinifera, SlPP2C= Solanum lycopersicum, AlPP2C= Arabidopsis lyrata, CsPP2C= Citrus sinensis, and PtPP2C= Populus trichocarpa.Following the similar method, we also constructed a phylogenetic tree to compare the relationships of PP2C amongArabidopsis, Amborella trichopoda, Vitis vinifera, Populus trichocarpa, Capsella rubella, Citrus sinensis, Carica papaya, Solanum lycopersicum, Fragaria vesca, Arabidopsis Iyarta, and a bryophyte (Physcomitrella patens). Intriguingly, we observed diversification within subgroups of almost every species, which might indicate that different species used for tree topologies lead to slight variations (Supplementary Figures S10 and S11).The copy number of variation and gene retention inB. rapa andA. thaliana during aBrassica-specific WGT event was also investigated. Meanwhile,B. rapagenome is based on three subgenomes and shares the same diploid ancestor withA. thaliana[26]. For each gene, we identified its syntenic paralogous and orthologous pairs by utilizing theBrassica rapa database (BRAD) (Table S4). In particular, we counted the retention ofBraPP2C genes based on subgroups and by calculating the number of gene copies among three subgenomes. It was observed that the majority of genes were found either in a single, double or triple copy. Interestingly, for subgroup E, only two pairs were found with three copies and 6 each pair were identified for subgroup D and F with two copies, whereas the remaining did not show any copy variation and group F also shows the most number (25) of genes based on subgenomes comparison (Figures 3(a) and 3(b)). Our findings also revealed that gene retention exhibited almost with varied results, 15/15 subgroup A, 7/8 subgroup B, 9/9 subgroup C, 16/17 subgroup D, 24/25 subgroup E, 24/25 subgroup F, 9/9 subgroup G, 6/6 subgroup H, 4/4 subgroup I, 1/1 subgroup J, 6/7 subgroup K, 4/4 subgroup L, and 3/3 unclassified, respectively (Table S4).Figure 3
(a) and (b) The copy number variation and the ratio of PP2C among three subgenomes ofB. rapa.
(a) (b)


### 3.4. Expression Pattern of BraPP2C in Various Tissues and the Correlation Networking Analysis of Paralogous Pairs

To address the divergence and possible involvement ofBraPP2C genes inB. rapa reflected in their growth and development. The trends of transcriptional profiles across five various tissues (roots, stems, leaves, flowers, and siliques) were determined based on previously published RNA-sequence data [45]. In this study, FPKM values were selected to representBraPP2C genes expression and heat maps were constructed to demonstrate the relative transcriptional profiling of 131BraPP2C genes in various tissues. The results displayed high alterations in expression profiling among various subgroups members of PP2C inB. rapa. Among 131BraPP2C genes,BraPP2C1,BraPP2C15,BraPP2C16,BraPP2C17,BraPP2C69,BraPP2C76,BraPP2C85, andBraPP2C98 showed no expression andBraPP2C,BraPP2C18-19,BraPP2C27,BraPP2C51,BraPP2C59-60,BraPP2C91,BraPP2C97, andBraPP2C118 showed slight expression in one of any tissues (Figure 4(a); Table S5). The rest of the PP2C were expressed in at least two or more organs. However, few members of BraPP2C were specifically observed to be selectively expressed in tissue-specific clustering (Figure 4(b)), such that four of each gene in flowers and siliques were found and two of them were observed in the stem. Intriguingly, these genes indicated a possible role in organ development ofB. rapa.Figure 4
(a) and (b) Heatmap of expression profiles (in log2-based FPKM) for PP2C in the five various tissues: root, stem, leaf, flower, and silique. The expression levels are indicated by the color bar. (b). Venn diagram analysis of the tissue expression of PP2C.
(a) (b)We next investigated the expression trends among 29 paralogous pairs forPP2C genes along with calculated their Pearson’s correlation coefficient values (Figure S12A; Table S6). These paralogous pairs exhibited significant variation in various tissues. The results showed that more than 20 pairs were highly expressed in all tissues and 13 pairs with higher PCC values (>0.6). Meanwhile, four pairs (BraPP2C2_BraPP2C1, BraPP2C77_BraPP2C76, BraPP2C85_BraPP2C86, and BraPP2C104_BraPP2C102) showed no correlation values, and four pairs (BraPP2C42_BraPP2C41, BraPP2C60_BraPP2C61, BraPP2C66_BraPP2C67, and BraPP2C72_BraPP2C71) exhibited a negative correlation. Also, the clustering image did not show any variability among various tissues (Figure S12B). Intriguingly, these results indicated that a possible involvement ofBraPP2C genes during growth by showing an overlapping in expression pattern might be the basis of crosstalk for signal transductions pathways and plant improvements. Despite the divergence in the expression profiles for BraPP2C between paralogous pairs, suggesting that after duplication in the evolutionary process, some of these pairs may acquire new functions.


### 3.5. Analysis of Putative RegulatoryCis-Element in BraPP2C and Coregulatory Expression Network Analysis

Thecis-elements in promoter regions are closely associated with gene transcription, and their response to stress as gene expression is mediated through interaction betweencis-regulatory elements and its cognate transcription factors [46]. Therefore, 2.0 kb upstream sequences were downloaded fromB. rapa database and analyzed using the PlantCARE database [37]. In the promoter regions of BraPP2C, we identified a total of 7 commoncis-regulatory elements (Figure 5 and Table S7). The skn-1 motif was responsive for endosperm expression. MBS is deliberated as an MYB binding site, which is involved in drought-inducibility. DRE and LTR arecis-acting elements that are known to control functions during adverse conditions, such as dehydration, chilling, and salt stress. Few of thecis-regulatory elements were actively receptive against hormonal stresses, such as ethylene, gibberellin, auxin, and salicylic acid. Circadiancis-regulatory element regulates the circadian rhythm of the plant system. As shown in Figure 5, the participation of hormones, light, essential elements, enhancers, stress factors, circadian, and other regulatory stress factors was higher inBraPP2C genes. As a consequence, the data suggests the presence of a large number of conservedcis-regulatory elements that are crucial in mediating responses to adverse environmental stresses or stress-related hormones. Though, to prove such assumptions, a further validation step is required. Moreover, these results highlighted the stress-responsive disposition ofBraPP2C genes.Figure 5
The ratios of different cis-elements.For functional studies, gene expression profiling provides valuable clues, manyPP2Chave been implicated in plants against abiotic stress conditions, and the presence of stress-responsivecis-elements in the promoter region further suggests their possible involvement inB. rapa response to different environmental stimuli. Therefore, we investigated the transcriptional profiles of 15 PP2C paralogous gene pairs ofB. rapa subjected to cold, heat, PEG, and ABA treatments using qRT-PCR analysis (Figure 6; Table S8). Our results indicated that theBraPP2C genes showed diversity in their expression level against four various treatments. Interestingly, the expression patterns of some genes induced rapidly when exposed to ABA and drought application. In contrary, some genes were highly upregulated or downregulated by either of one or two stresses. Under heat treatment, a dominant portion of theBraPP2C genes showed highly fluctuated transcriptional profiling levels, including 30% upregulated and 70% downregulated genes. The exposure ofB. rapa to cold stress resulted in approximately 54% of upregulated genes. The transcriptional level was significantly increased under ABA and drought treatments; most of the genes were highly 57% and 64% expressed, respectively. Taken together, the diverse expression pattern ofBraPP2C genes may suggest various roles after exposure to abiotic and hormone stress conditions. To understand the correlation and coregulatory network among these paralogous pairs, we calculated Pearson’s correlation coefficient (PCC) values based on their relative expression data. For the correlation network, we make three categories concerning PCC values, such as more than 0.6 (Highly Positive), less than 0.5 but greater than 0 (Mild Positive), and negative values with (Negative) correlation (Figure 7 and Table S9). Heat and drought stress showed highly positive and closer relationship by showing (11) PCC values each, while cold stress was found with a high number of negative (6) PCC values, suggesting its different nature among the paralogous pairs ofBraPP2Cgenes.Figure 6
Expression analysis ofPP2C genes under heat, cold, ABA, and drought stress treatments inB. rapa.Figure 7
Pearson’s correlation coefficient (PCC) values under response to heat, cold, ABA, and drought treatments.Moreover, we examined the interaction network of BraPP2C proteins with the help of STRING software. The BraPP2C proteins were highly correlated with each other in the interaction network. As shown in FigureS13, most of the subgroups of BraPP2C exhibited a dense network when compared with each other. Most of the proteins were highly located in the center and only a few of them did not interact with others, which indicated the complexity of the interaction network. Therefore, this PP2C may involve in the regulation of many downstream/upstream gene factors by playing a crucial role in fundamental molecular mechanisms of plants.


## 3.1. Characterization ofBraPP2C Gene Family Members

In the present study, we identified 131 putativeBraPP2C genes inB. rapa genome and were denoted as BraPP2C1 to BraPP2C131, based on phylogenetic analysis and their orthologous positions withArabidopsis. The physicochemical feature of theBraPP2C proteins along with some key information of all the identified BraPP2C is provided in Table S2. The results of sequence analysis exhibited that the length of the BraPP2C proteins varied with a ranged from 102 (BraPP2C1) to 1465 (BraPP2C63) bp, with an average of 412.37 (aa). The detailed information of BraPP2C, including molecular weights (MW), theoretical isoelectric point (pI), the grand average of hydropathicity (GRAVY), and exon number ranged from 11.76 (BraPP2C27) to 160.43 (BraPP2C63) kDa, -0.777 (BraPP2C59) to 0.001 (BraPP2C41), 4.36 (BraPP2C85 to 9.37 (BraPP2C117), and 1 (BraPP2C73) to 21 (BraPP2C21), respectively. The significant variations reflect their functional diversity among BraPP2C, while the negative value of GRAVY further demonstrated that these proteins were unstable and hydrophilic in nature. Moreover, we analyzed the values of pI, GRAVY, and exon number, as shown in (Supplementary Figures 1, 2, and 3), respectively. The subcellular localization was predicated, and results proposed that most of the BraPP2C proteins were localized in mitochondria, cytoplasm, nuclei, chloroplasts, plasma membranes, endoplasmic reticulum, and vacuoles (Table S2). The gene structure analysis of BraPP2C revealed that the exon number of eachBraPP2C gene ranged dramatically from 1 to 19 (Figure S4). Majority ofBraPP2C genes consisted of only one exon; the rest of others showed significant variations, whereas BraPP2C32 possesses 19 exons, suggesting that both exon loss and gain occurred inBraPP2C gene family. Alongside, we systematically discovered the distribution of conserved motifs and their logos by online MEME server. A total of ten different consensus motifs were obtained in all of the BraPP2C proteins, and their distribution patterns are presented in Figures S5 and S6. Hence, during the evolutionary process, the BraPP2C demonstrated the extreme conservation within subgroups.


## 3.2. Chromosome Distribution and Syntenic Gene Collinearity Analysis

The physical chromosomal locations onB. rapa, across all ten chromosomes (Br01–Br10), were determined by the Mapchart software. The 126BraPP2C genes were dispersed across all ten chromosomes, ranging from 7 to 19 per chromosome, out of which 5 genes were found on the scaffold region (Figure 1). On each chromosome, the number of BraPP2C varies drastically; the largest number of BraPP2C family members was observed on chromosomes Br05 with 19 genes, followed by Br03 with 18 genes, whereas the least numbers were revealed on chromosomes Br04 and Br07; each contains 7 genes, respectively. Among these BraPP2C, five genes (BraPP2C28,BraPP2C68,BraPP2C80,BraPP2C83, andBraPP2C111) were present on the genomic scaffold. These results suggest the uneven distribution patterns of BraPP2C across ten different chromosomes. Based on a previously reported study, 24 ancestral genomic blocks (GBs) were reconstructed [27]. The color-coding for these block were slightly modified according to their position in a proposed ancestral karyotype (AK1-07) [27, 43]. The majority (18.18%) ofBraPP2C genes were clustered in the AK5 region, followed by AK1 and AK2 (16.88%), whereas minimum genes (6.49%) were clustered in AK7 region (Figure 1).Figure 1
Chromosome locations of PP2C were obtained from the GFF file and displayed using Mapchart. The ancestral karyotypes are marked in different colors.The collinear relationships of thePP2C gene pairs betweenArabidopsisandB. rapa are shown in Figure 2. Based on three subgenomes ofB. rapa, LF contains more (42.86%)BraPP2C genes than the MF1 (32.77%) and MF2 (24.37%) subgenomes (Figure S7). In addition, the divergence time of the 29 paralogous pairs was estimated with the help of MEGA7.0 software by calculating the synonymous (Ks) and nonsynonymous substitution rates (Ka). The results proposed that most of theBraPP2C gene pairs represented less than 1.00 (Ka/Ks) ratios, suggesting the purifying selection of these genes except three pairs (BraPP2C58-BraPP2C56,BraPP2C85-BraPP2C86, andBraPP2C125-BraPP2C127) showed higher than 1.00 values, indicating positive selection. To estimate the divergence among theseBraPP2C paralogous gene pairs, we calculated the values ofKs (Figure S8 and Table S3). The average divergence time among these paralogous is 12.37 MYA, suggesting that BraPP2C divergence occurred along with theArabidopsis (9.6–16.1 MYA) [6].Figure 2
Phylogenetic relationship of PP2C betweenB. rapa andA. thaliana. The phylogenetic tree was constructed by MEGA 7 using the Maximum Likelihood Method (1000 bootstrap). Genes of PP2C paralog pairs are marked with red color. The collinear correlation for all genes of PP2C is displayed betweenB. rapa andA. thaliana. The 10 Chinese cabbage chromosomes (Br01–Br10) and fiveA. thaliana chromosomes (AT01–AT05) are shown in different random colors. The illustration was drawn using Circos Software.


## 3.3. Phylogenetic Analysis ofBraPP2C Gene and Copy Number Variations/Gene Retention

To provide an overview of the evolutionary relationships among subgroup member of PP2C betweenB. rapa andA. thaliana, a comparative phylogenetic tree was constructed using MEGA7.0 software by adapting the maximum likelihood approach with 1000 bootstrap replications (Figure S9). The phylogenetic tree and domain composition further subcategorized the PP2C into thirteen subgroups: subgroup A-L and one unclassified. The results of the tree were consistent with the previously reported studies inArabidopsis and rice [44]. Moreover, most of the BraPP2C clustered together with those fromArabidopsis, while 128 out of 131BraPP2C genes are further distributed among twelve subfamilies (A-L), and the remaining three (BraPP2C35,BraPP2C60, andBraPP261)BraPP2C genes cannot be grouped into any subfamilies and were categorized as unclassified. Most of the PP2C members were found in two subgroups E and F each with 25 genes, respectively, as compared to other subgroups. The complete overview and the distribution patterns amongBrassica rapa and various ten species are described in Table 1.Table 1
The distribution ofPP2C genes family is various species.
Subgroup of PP2C BraPP2C AthPP2C AtrPP2C CrPP2C CpPP2C FvPP2C PpPP2C VvPP2C SlPP2C AlPP2C CsPP2C PtPP2C A 15 10 2 10 4 6 0 5 6 11 11 11 B 7 6 1 4 3 1 0 0 3 2 6 5 C 10 7 4 8 3 4 0 5 5 6 18 13 D 16 9 4 11 6 8 4 9 12 9 15 29 E 25 12 4 16 6 8 14 5 9 13 18 26 F 25 13 6 16 3 4 0 8 8 12 23 25 G 9 6 3 8 2 7 0 5 6 6 12 14 H 6 3 3 4 3 3 28 3 5 3 3 15 I 4 2 0 3 2 1 19 2 3 2 10 4 J 1 2 0 1 0 0 0 0 1 1 0 1 K 6 3 1 4 3 3 0 3 2 4 5 7 L 4 0 2 1 3 1 0 2 2 2 6 4 Unclassified 3 3 1 2 0 1 0 1 1 2 0 2 ∗ BraPP2C= Brassica rapa, AthPP2C= Arabidopsis thaliana, AtrPP2C= Amborella trichopoda, CrPP2C= Capsella rubella, CpPP2C= Carica papaya, FvPP2C= Fragaria vesca, PpPP2C= Physcomitrella patens, VvPP2C= Vitis vinifera, SlPP2C= Solanum lycopersicum, AlPP2C= Arabidopsis lyrata, CsPP2C= Citrus sinensis, and PtPP2C= Populus trichocarpa.Following the similar method, we also constructed a phylogenetic tree to compare the relationships of PP2C amongArabidopsis, Amborella trichopoda, Vitis vinifera, Populus trichocarpa, Capsella rubella, Citrus sinensis, Carica papaya, Solanum lycopersicum, Fragaria vesca, Arabidopsis Iyarta, and a bryophyte (Physcomitrella patens). Intriguingly, we observed diversification within subgroups of almost every species, which might indicate that different species used for tree topologies lead to slight variations (Supplementary Figures S10 and S11).The copy number of variation and gene retention inB. rapa andA. thaliana during aBrassica-specific WGT event was also investigated. Meanwhile,B. rapagenome is based on three subgenomes and shares the same diploid ancestor withA. thaliana[26]. For each gene, we identified its syntenic paralogous and orthologous pairs by utilizing theBrassica rapa database (BRAD) (Table S4). In particular, we counted the retention ofBraPP2C genes based on subgroups and by calculating the number of gene copies among three subgenomes. It was observed that the majority of genes were found either in a single, double or triple copy. Interestingly, for subgroup E, only two pairs were found with three copies and 6 each pair were identified for subgroup D and F with two copies, whereas the remaining did not show any copy variation and group F also shows the most number (25) of genes based on subgenomes comparison (Figures 3(a) and 3(b)). Our findings also revealed that gene retention exhibited almost with varied results, 15/15 subgroup A, 7/8 subgroup B, 9/9 subgroup C, 16/17 subgroup D, 24/25 subgroup E, 24/25 subgroup F, 9/9 subgroup G, 6/6 subgroup H, 4/4 subgroup I, 1/1 subgroup J, 6/7 subgroup K, 4/4 subgroup L, and 3/3 unclassified, respectively (Table S4).Figure 3
(a) and (b) The copy number variation and the ratio of PP2C among three subgenomes ofB. rapa.
(a) (b)


## 3.4. Expression Pattern of BraPP2C in Various Tissues and the Correlation Networking Analysis of Paralogous Pairs

To address the divergence and possible involvement ofBraPP2C genes inB. rapa reflected in their growth and development. The trends of transcriptional profiles across five various tissues (roots, stems, leaves, flowers, and siliques) were determined based on previously published RNA-sequence data [45]. In this study, FPKM values were selected to representBraPP2C genes expression and heat maps were constructed to demonstrate the relative transcriptional profiling of 131BraPP2C genes in various tissues. The results displayed high alterations in expression profiling among various subgroups members of PP2C inB. rapa. Among 131BraPP2C genes,BraPP2C1,BraPP2C15,BraPP2C16,BraPP2C17,BraPP2C69,BraPP2C76,BraPP2C85, andBraPP2C98 showed no expression andBraPP2C,BraPP2C18-19,BraPP2C27,BraPP2C51,BraPP2C59-60,BraPP2C91,BraPP2C97, andBraPP2C118 showed slight expression in one of any tissues (Figure 4(a); Table S5). The rest of the PP2C were expressed in at least two or more organs. However, few members of BraPP2C were specifically observed to be selectively expressed in tissue-specific clustering (Figure 4(b)), such that four of each gene in flowers and siliques were found and two of them were observed in the stem. Intriguingly, these genes indicated a possible role in organ development ofB. rapa.Figure 4
(a) and (b) Heatmap of expression profiles (in log2-based FPKM) for PP2C in the five various tissues: root, stem, leaf, flower, and silique. The expression levels are indicated by the color bar. (b). Venn diagram analysis of the tissue expression of PP2C.
(a) (b)We next investigated the expression trends among 29 paralogous pairs forPP2C genes along with calculated their Pearson’s correlation coefficient values (Figure S12A; Table S6). These paralogous pairs exhibited significant variation in various tissues. The results showed that more than 20 pairs were highly expressed in all tissues and 13 pairs with higher PCC values (>0.6). Meanwhile, four pairs (BraPP2C2_BraPP2C1, BraPP2C77_BraPP2C76, BraPP2C85_BraPP2C86, and BraPP2C104_BraPP2C102) showed no correlation values, and four pairs (BraPP2C42_BraPP2C41, BraPP2C60_BraPP2C61, BraPP2C66_BraPP2C67, and BraPP2C72_BraPP2C71) exhibited a negative correlation. Also, the clustering image did not show any variability among various tissues (Figure S12B). Intriguingly, these results indicated that a possible involvement ofBraPP2C genes during growth by showing an overlapping in expression pattern might be the basis of crosstalk for signal transductions pathways and plant improvements. Despite the divergence in the expression profiles for BraPP2C between paralogous pairs, suggesting that after duplication in the evolutionary process, some of these pairs may acquire new functions.


## 3.5. Analysis of Putative RegulatoryCis-Element in BraPP2C and Coregulatory Expression Network Analysis

Thecis-elements in promoter regions are closely associated with gene transcription, and their response to stress as gene expression is mediated through interaction betweencis-regulatory elements and its cognate transcription factors [46]. Therefore, 2.0 kb upstream sequences were downloaded fromB. rapa database and analyzed using the PlantCARE database [37]. In the promoter regions of BraPP2C, we identified a total of 7 commoncis-regulatory elements (Figure 5 and Table S7). The skn-1 motif was responsive for endosperm expression. MBS is deliberated as an MYB binding site, which is involved in drought-inducibility. DRE and LTR arecis-acting elements that are known to control functions during adverse conditions, such as dehydration, chilling, and salt stress. Few of thecis-regulatory elements were actively receptive against hormonal stresses, such as ethylene, gibberellin, auxin, and salicylic acid. Circadiancis-regulatory element regulates the circadian rhythm of the plant system. As shown in Figure 5, the participation of hormones, light, essential elements, enhancers, stress factors, circadian, and other regulatory stress factors was higher inBraPP2C genes. As a consequence, the data suggests the presence of a large number of conservedcis-regulatory elements that are crucial in mediating responses to adverse environmental stresses or stress-related hormones. Though, to prove such assumptions, a further validation step is required. Moreover, these results highlighted the stress-responsive disposition ofBraPP2C genes.Figure 5
The ratios of different cis-elements.For functional studies, gene expression profiling provides valuable clues, manyPP2Chave been implicated in plants against abiotic stress conditions, and the presence of stress-responsivecis-elements in the promoter region further suggests their possible involvement inB. rapa response to different environmental stimuli. Therefore, we investigated the transcriptional profiles of 15 PP2C paralogous gene pairs ofB. rapa subjected to cold, heat, PEG, and ABA treatments using qRT-PCR analysis (Figure 6; Table S8). Our results indicated that theBraPP2C genes showed diversity in their expression level against four various treatments. Interestingly, the expression patterns of some genes induced rapidly when exposed to ABA and drought application. In contrary, some genes were highly upregulated or downregulated by either of one or two stresses. Under heat treatment, a dominant portion of theBraPP2C genes showed highly fluctuated transcriptional profiling levels, including 30% upregulated and 70% downregulated genes. The exposure ofB. rapa to cold stress resulted in approximately 54% of upregulated genes. The transcriptional level was significantly increased under ABA and drought treatments; most of the genes were highly 57% and 64% expressed, respectively. Taken together, the diverse expression pattern ofBraPP2C genes may suggest various roles after exposure to abiotic and hormone stress conditions. To understand the correlation and coregulatory network among these paralogous pairs, we calculated Pearson’s correlation coefficient (PCC) values based on their relative expression data. For the correlation network, we make three categories concerning PCC values, such as more than 0.6 (Highly Positive), less than 0.5 but greater than 0 (Mild Positive), and negative values with (Negative) correlation (Figure 7 and Table S9). Heat and drought stress showed highly positive and closer relationship by showing (11) PCC values each, while cold stress was found with a high number of negative (6) PCC values, suggesting its different nature among the paralogous pairs ofBraPP2Cgenes.Figure 6
Expression analysis ofPP2C genes under heat, cold, ABA, and drought stress treatments inB. rapa.Figure 7
Pearson’s correlation coefficient (PCC) values under response to heat, cold, ABA, and drought treatments.Moreover, we examined the interaction network of BraPP2C proteins with the help of STRING software. The BraPP2C proteins were highly correlated with each other in the interaction network. As shown in FigureS13, most of the subgroups of BraPP2C exhibited a dense network when compared with each other. Most of the proteins were highly located in the center and only a few of them did not interact with others, which indicated the complexity of the interaction network. Therefore, this PP2C may involve in the regulation of many downstream/upstream gene factors by playing a crucial role in fundamental molecular mechanisms of plants.


## 4. Discussion

In the plant kingdom, thePP2C gene family is considered as one of the largest families and has been characterized in several crop species by utilizing bioinformatics strategies. Previously, various members ofPP2C genes were identified in maize [47], rice [48],Arabidopsis[49], hot pepper [50], banana [51], andBrachypodium distachyon [3]. In this study, a comprehensive genome-wide analysis was performed, including gene identification, phylogenetic relationships, evolutionary analysis, chromosomal localizations, conserved domain motifs, and gene structure organization analysis. In addition, gene expression patterns of some keyBraPP2C genes were also determined under heat, cold, ABA, and drought stresses. Herein, a total of 131BraPP2C genes were identified. ThePP2C genes were further categorized into 11, 12, or 13 subfamilies depending on the plant species according to evolutionary analysis. The proportion of PP2C members was much higher as compared to lower plants. For this reason, we further isolated thePP2C genes in ten various species and calculated their number as well, such asAmborella trichopodawhich contains 31PP2C genes,Vitis vinifera(48), Populus trichocarpa(155), Capsella rubella(88), Citrus sinensis(127), Carica papaya(38), Solanum lycopersicum(63), Fragaria vesca(47), Arabidopsis lyrata(73),andPhyscomitrella patens (65) as shown in Figure 8 and Table S10. As a consequence, a high degree of divergence was observed betweenPP2C genes from lower plants to higher plants, suggesting that it may be correlated with the complex adaptation to environmental stimuli [52]. Overall, the number ofPP2C genes has remained distinct in each investigated species, proposing that variation in evolutionary patterns is occurring, due to gradual changes in rounds of whole-genome duplication (WGD) and subsequent gene losses/gains by natural selection constraints.Figure 8
The proportion ofPP2C genes among various species.Notably, gene duplication, being the predominant influencing force for broad expansion of the gene family, upgrades its biological function and evolutionary processes [53]. Thus, to study gene family fractionation, a linkage between them and their morphotypes,B. rapa is an excellent model plant because it has experienced WGD and an additional WGT event [26, 27]. Considering the importance of evolutionary analysis, most of the angiosperms have undergone either one or multiple polyploidization events [54–56]. As a result, for duplicated genes, it provides opportunities to diverge in various evolutionary ways. Subsequently, each of these genes experienced one of the following three evolutionary fates: subfunctionalization (the ancestral function is subdivided between copies), neofunctionalization (one copy acquires a new function), or nonfunctionalization (one copy becomes unexpressed or functionless) [57]. Furthermore, we predicted the pressure of natural selection by calculating the synonymous (Ks) and nonsynonymous substitution rates (Ka) and were identified by the MEGA7.0 program. Also, during evolutionary processes and expansion of genes family, these indicators are used for the selection history. If the value ofKa/Ks is lower than 1.00, it represents purifying selection,Ka/Ks = 1 means neutral selection, andKa/Ks> 1 means positive selection [7]. In our work, we estimated the divergence time among the 29 paralogous pairs ofBraPP2C genes. Majority of the BraPP2C paralogous pairs showed less than 1.00Ka/Ks ratios, speculating the purifying selection of these paralogous pairs. On the other hand, only three pairs show more than 1.00 values, suggesting positive selection. Moreover, the divergence time of these paralogous pairs of theBraPP2C genes ranged from 0.03 to 0.84 (Ks values) with an average mean divergence of 12.37 (MYA). Notably, the divergence time ofBraPP2C paralogous gene pairs was 12.37 MYA, which further intimates that their divergence occurred during theBrassicaandArabidopsisduplication (9.6–16.1 MYA) [7].Plants productivity encounters a severe threat by maintaining the optimum crop production, due to a wide range of environmental challenges. Surrounded by various stress factors, such as heat, cold, ABA, and drought represent the major constraint for agricultural crop production. In regulating plant’ tolerance to multiple abiotic stress factor a well-known hormone, ABA plays a crucial role. However, inB.rapa, the role of the core components of ABA signaling againstPP2C genes in responding to various stress conditions is mainly obscure. Specifically, various members of PP2C group A (ABI1, ABI2, HAB1, HAB2, AHG1, and PP2CA) resulted into increase the ABA sensitivity after double and triple mutation, suggesting the diverse outcome in ABA signaling [17, 58–61]. Previous studies reported thatPP2C regulates positively against salt tolerance inArabidopsis and drought in peach [62, 63] to modulate the stress severity. InArabidopsis, two members ofPP2C genes were responded differently, such asAP2C1 expression was strongly induced by cold, drought, and wounding, butAP2C2 was slightly influenced by these treatments [64]. These findings highlighted the significance of critical members of thePP2C gene family in the model plantArabidopsis, but their specific functions may be significantly varied inB. rapa. In the present study, variousPP2C genes showed high striking transcriptional changes followed by heat, cold, ABA, and drought stresses, indicating that some of these genes might be pivotal to stress tolerance inB. rapa. Some of the recent advances in the functional dissection of PP2C candidate genes uncovered their importance in the life cycle ofArabidopsis and rice [48, 65], although their associated roles inB. rapaare mostly unnoticed. Aiming for achieving gene expression patterns in various growth phases ofBraPP2C, we utilized the previously reported RNA-sequence data and analyzed the expression profile ofPP2C genes in various tissues (root, stem, leaf, flower, and silique). Results demonstrated that most of theBraPP2C genes showed a varying response to multiple tissues and some of the genes were highly expressed in all the tissues or in some cases; it was negative. However, few genes have shown tissue-specific expression, such that four genes were commonly identified in roots and flower and two of them were expressed in the stem as well, intimating their importance in plant development. In addition, gene family expansion and the duplication types in the neofunctionalization or subfunctionalization models are mainly associated with tissue expression divergence [66–68]. InBraPP2C genes, we also explored the promoter regions for identification of common conservedcis-regulatory elements. As to further elucidate the predicated and the possible functions ofBraPP2Cgenes in transcriptional regulation, the result represented significant variation amongBraPP2C genes and was mostly responsive to both biotic-abiotic factors and plant hormones stresses.It is indispensable to perceive the fundamental molecular mechanisms of plants for the adaptability of stress tolerance and enhancement in crop yield under adverse stress conditions. Therefore, our outcome of the study provides valuable insight and will serve as a basis for exploring the pivotal role ofBraPP2C genes in the plant against abiotic stresses. It may also assist to reveal the potential functions ofBraPP2C genes in response to both abiotic and hormone stress conditions.


## 5. Conclusions

In this study, we identified a total of 131BraPP2Cgene family members by genome-wide analysis inBrassica rapagenome. A comprehensive analysis was carried out to study the phylogenetic relationships, evolutionary analysis, exon-intron organization, chromosomal localization, protein structural features, interaction network analysis, and conserved motifs. According to the phylogenetic relationships, the BraPP2C were classified into thirteen subgroups. Moreover, transcriptional profiling revealed BraPP2C candidate genes might have participated in the plant stress tolerance particularly to heat, cold, ABA, and drought stress. TheseBraPP2C genes can be utilized to facilitate by functionally characterize them, laying the foundation for elucidating their specific regulatory mechanisms and ultimately applying them in molecular breeding programs of this important vegetable crop. In summary, the integration of our findings has provided a novel insight and unique features ofBraPP2Cgenes that may play a potential role in mediating regulation of signal transduction and particularly under various stress conditions.



---

*Source: 2965035-2019-04-03.xml*