Genome-Wide Identification, Classification, and Expression Divergence of Glutathione-Transferase Family in Brassica rapa under Multiple Hormone Treatments

The GSTs is one of the most important multifunctional protein families which has been playing a crucial role in the different aspects of plant growth. This extensive study about GSTs may establish a solid foundation for the brief functional analysis of BraGSTs in future. In this study, a total of 75 genes were identified in B. rapa. Phylogenetic analysis characterized them into eight different subclasses, while Tau and Phi subclasses were the most numerous. The exon-intron structure and the motif composition of BraGSTs were exhibited accordingly to their subclasses. Notably, we also investigated 15 tandem paralogous pairs of genes, which highlighted that all the pairs were purifying in nature as their synonymous values were lower than 1.00. Duplication analysis indicated that about 45.33% of genes mainly occurred through tandem duplication in B. rapa. Predominately, the tandem cluster of genes in subclass Tau was greater than the other subclasses. Furthermore, among eight multiple hormonal treatments (ABA, GA, BR, ETH, IAA, IBA, NPA, and JA), most number of BraGSTs was activated by NPA, BR, and ABA treatments. This analysis has provided comprehensive information about GSTs family which may assist in elucidating their exact functions in B. rapa.


Introduction
The plant glutathione transferases, formally glutathione-stransferase (GSTs; EC 2.5. 1.18), are known for their diversity, being multifunctional proteins, and being widely distributed in most of the organisms. To date, the role of GST proteins was considered to be crucial in multiple plant functions such as herbicide detoxification, plant developmental processes, signal transduction, oxidative damage, heavy metals toxicity, and other several key biotic and abiotic factors [1,2]. It has been reported that GSTs participated in the endogenous developmental process, for example, in maize [3], petunia [4], Arabidopsis [5], grapevine, and cyclamen [6]; they were involved in flavonoid binding. The versatility of the GSTs was further extended through previous studies, which have proved their involvement in the regulation of genes. It increases the transcript level through exposure to biotic and abiotic factors including hormones, specifically auxins, ethylene, salicylic acid (SA), abscisic acid (ABA), and methyl jasmonate (MeJA) [7,8].
The GSTs superfamily have shown variation in their structure and sequences, although they share structural homology based on considered codomain thioredoxin/ glutaredoxin-like N-terminal and a larger C-terminal. Thus, variability in nature leads to different hydrophobic substrate specificities in plant GSTs [9]. On the other hand, as compared to mammals plant, GSTs possess a larger cleft for co-substrate, as a result it accepted much more and larger diverse substrates [10,11]. The GSTs can be categorized into eight different subclasses based on sequence similarities of amino acid and genes structure such as Phi, Tau, Lambda, Zeta, Theta, Lambda, dehydroascorbate reductase (DHAR), -subunit of translation elongation factor (EF1G), and tetrachlorohydroquinone dehalogenase (TCHQD). Subclasses of GSTs like Phi, Tau, Lambda, and DHAR are considered as plant-specific [12] while Phi and Tau are the most numerous 2 BioMed Research International and highly functioning subclasses for their involvement in the detoxification of xenobiotics [13]. Predominantly, most of the GSTs were characterized and investigated in several plant species belonging to Phi and Tau subclass. The overexpression of different genes member of Phi and Tau subclass has been reported to increase the herbicide tolerance [14]; salinity and oxidative stress [15,16]; and chilling stress [17]. Meanwhile, in Sorghum, 62.5% of genes from subclass Tau have shown significant responses against multiple abiotic factors such as cold, drought, and salinity stresses [7]. However, others subclasses of GSTs such as Lambda and DHAR are mainly involved in redox and thiol transfer reactions rather than xenobiotics [18,19]. DHAR subclass plays important role in stress resistance and is involved in the ascorbate-GSH recycling reactions [20][21][22]. Although, two subclasses theta and zeta mainly helps in the function as GSH-dependent peroxidase, isomerase, and act as counterparts in the mammalian system [23,24]. The EFIG class is basically based on two domains, a typical GST and EF1G domain, which mainly function as GSH peroxidase [25]. GSTs stimulate the scavenging pathways and help them to reduce the effect of toxic material in plant tissues [26]. In addition, glutathione is a specific tripeptide, produced by the combination of glutamic acid, cysteine, and glycine, and referred to as gamma-glutamyl-cysteinylglycine [27]. It is present in body fluid (blood) and cells [28]. On the other hand, it is also located in cytosol and other cell organelles such as endoplasmic reticulum, nucleus, and mitochondria [29,30]. In particular, glutathione plays a key role in various metabolic pathways that focused on improving the quality of nutrients. Its functions are also linked with the controlling the gene expression patterns, signal transduction, synthesis of DNA, and protein [28]. It has been reported recently that allyl isothiocyanate elevated the expression of GST genes by inducing the oxidative stress in A. thaliana [31]. The above studies have provided a brief background of GSTs and their importance in plant cell structure and function. However, their associated endogenous role has still need to be identified in B. rapa, against such multiple hormone treatments, and thus is largely obscure.
The B. rapa genome (Chiifu-401-42) has recently been sequenced and assembled [32], which elucidated its close relationships with A. thaliana and after its divergence, it has experienced a whole genome triplication event (WGT) [33,34]. In this study, a total of 75 GSTs family members were identified in B. rapa genome. Here, we conducted a systematic study of GST genes family in B. rapa, to identify the characterization and phylogenetic relationships with A. thaliana and rice along with collinear correlation between B. rapa and A. thaliana. Interaction networks and expression divergence of BraGSTs under multiple hormonal treatments were also investigated. As a result, this study provides a foundation in understanding the mechanism of GSTs family in B. rapa and valuable information for further investigation of hormonal stress responsive genes.

Identification of the GST Genes in Brassica rapa.
From the Brassica database (BRAD; http://brassicadb.org/brad/) [34], all the genome sequence datasets were downloaded. The A. thaliana GST sequences were retrieved from the Arabidopsis I nformation Resource database (http://www.arabidopsis.org/). The gene information of rice were obtained from http:// rice.plantbiology.msu.edu/ [35], based on the previous reported study [36]. The A. thaliana and rice were used as queries by performing BLASTP search to identify the putative GST proteins with best domain e-value cut-off 1.0. For validation, these potential sequences SMART (http://smart.embl-heidelberg.de/) and the National Center for Biotechnology Information (NCBI) database (http://www .ncbi.nlm.nih.gov/) sever were used. Finally, the physicochemical properties of the BraGSTs were analyzed through Expasy protparam (http://web.expasy.org/protparam/) and the subcellular localization was predicted by WoLF PSORT (https://www.genscript.com/wolf-psort.html).

Phylogenetic Analysis and Genomic
Organizations Prediction. The full length protein sequences of BraGSTs were aligned through MUSCLE with default parameters [37]. For each analysis, Maximum Likelihood Method was used for constructing phylogenetic tree with JTT model (Jones-Taylor and Thornton amino acid substitution), using MEGA 7 [38]. The nucleotide divergence for all the GSTs was also analyzed through MEGA 7 with Jukes-Cantor model (1000 bootstrap values).
The coding sequences and the corresponding genomic sequences of BraGSTs were predicted by online tool (http://gsds.cbi.pku.edu.cn/) [39], for exon-intron structure. For the identification of conserved motifs of the GSTs in B. rapa, MEME (Multiple Expectation-Maximization for Motif Elicitation program) online server was used (Program version 4.12.0) [39]. The following parameters were implemented, maximum number of motifs 12 and optimum motif widths 12 and 120, while other parameters were set as default.

and
Calculation. The synonymous ( ) and nonsynonymous ( ) substitution rates among the tandem pairs of BraGST genes were calculated with the help of MEGA 7 software, based on the coding sequences alignment following the Nei and Gojobori model implemented in MEGA 7 [38].
Additionally, the divergence time was calculated through the following formula: = /2 in which was taken 1.5 × 10 −8 (synonymous substitution/year) by showing the rate of divergence [40].

Promoter Analysis and Proteins Interaction Network
Prediction. We explored a GFF (generic file format) from the B. rapa genome for all the GST promoter sequences (15,00 bp upstream). For each gene, we identified the cis-regulatory element by Plant CARE database (http://bioinformatics.psb .ugent.be/webtools/plantcare/html/). Additionally, with the help of STRING software (https://string-db.org/), we constructed an interaction network among BraGST proteins.

Chromosomal Localization and the Syntenic Paralog Pairs of BraGSTs.
All the BraGST genes were mapped on the ten B. rapa chromosomes based on their respective position in BioMed Research International 3 genome annotation with the random distribution patterns. The images for all the BraGSTs were drawn through Mapchart [41]. The syntenic relationship between A. thaliana and B. rapa was explored by searching the term "syntenic genes" in the B. rapa database [42]. To present the syntenic relationship on their chromosomes, Circos program [43] was applied.
2.6. Pearson Correlation (PCC) Analysis. The PCC values for RNA-seq data and expression patterns among the tandem pairs of BraGSTs were performed on excel sheet 2013.

RNA Isolation and Expression Pattern Analysis for
BraGSTs in Five Tissues. The total RNAs were isolated from young leaves using an RNA kit (TaKaRa, Dalian, China). The RNAs were reverse-transcribed into cDNA with a Prime-Script RT reagent kit (TaKaRa, Dalian, China). For qRT-PCR, gene specific primers were designed through Becan designer 7 software and for internal control the B. rapa actin gene (Bra028615) was used (Supplementary, Table 1).
Step one plus real-time PCR system (Applied Biosystems, Carlsbad, CA) was used for the reactions. The following PCR parameters were used: 94 ∘ C for 30 s, 40 cycles at 94 ∘ C for 05 s, 60 ∘ C for 15 s, and 72 ∘ C for 10 s, following by melting curve analysis (61 cycles at 65 ∘ C for 10 s). The relative gene expressions were further calculated by the comparative Ct value method [44]. For expression analysis of GSTs in B. rapa, we utilized the five accession of Chifu-401-42 (root, stem, leaf, flower, and silique) and the RNA-seq data were generated from previous reports [32]. Hence, the gene expression values were analyzed for each tissue and the fragments per kilobase of transcript per Million fragments mapped (FPKM) were quantified first and then heat maps were generated by using online tools (http://www.omicshare.com/).

Plant Materials and Hormone Treatments. The typical
Chinese cabbage cultivar Chiffu-401-42 was used in our experiment, as this cultivar being prominently used for research studies due to the completion of its whole genomesequencing. Seeds were first treated with sodium hypochlorite (14%) and then were raised on 0.5 MS agar plates (0.7%) in dark for three-day interval at 23 ∘ C. After that, the germinated seeds were raised in a controlled environment using small pots with 3 : 1 (soil : vermiculite mixture); growth chamber was programmed with temperature of 27 ∘ C; for photoperiod, the duration of light was 16 h and dark 8 h with 60% relative humidity in the greenhouse of Nanjing Agricultural University. One-month old seedlings were transplanted at five leaf-stage into 1/2 Hoagland's solution in small plastic pots with a pH at 6.5. After ten days of acclimatization, seedling of Chinese cabbage was grown with following

Identification of GST Proteins and Their Classification
Pattern. We identified all the putative GST genes in B. rapa through systemic BLASTP search against the B. rapa database, using the query sequences of A. thaliana (55) and rice (77). This search resulted in the identification of 75 GST proteins. Subsequently, all these protein sequences were subjected to SMART and NCBI for further verification analyses, and both GST N-and C-terminal domains were confirmed in B. rapa genome. To reveal the subclasses of GSTs in B. rapa, the results were based on conserved domains in the context of proposed nomenclature for GSTs [45,46]; it can be further divided into eight different subgroups as Phi, Tau, Zeta, Theta, Lambda, DHAR, EF1G, and TCHQD. The BraGSTs according to the subclasses were designated as BraGST. The genes of subclasses were named as BraGSTF, BraGSTU, BraGSTZ, BraGSTT, BraGSTL, BraDHAR, BraEF1G, and BraTCHQD, respectively, and corresponding number for each gene was, for example, BraGSTF1. From the predicted protein sequences of all the BraGSTs, we have collected all the basic description about their length, molecular weight (MW), isoelectric points (pI), grand average of hydropathicity (GRAVY), exon number, aliphatic index, subcellular prediction, and others which are briefly described (Supplementary, Table 2). In general, the amino acids for all the BraGSTs were varied 167-566 and the MW ranged 18.85-64.24 kDa, while the highest MW and length were recorded for subclass of EFIG with an average of 52.37 kDa and 464. The pI values were ranged 4.9-9.73, which speculated that different BraGST proteins may function in different microenvironments as described in Figure 1 Table 2). Though, the GRAVY from the subclass of BraGSTZ with only one protein were hydrophobic in nature, whereas the rest of them showed hydrophilic properties (Figure 1(c) and Supplementary, Table 2). To understand plant functions, protein localization is particularly important. In our study, most of the subclasses of BraGSTs were predicted to be located in mitochondria, chloroplast, and cytoplasm, while only a small number of proteins were involved in plasma membrane, vacuole, and nucleus.

Expansion and Structural Characteristics of BraGST Genes in Brassica rapa.
To validate the GST genes relationship among B. rapa, A. thaliana, and rice, all the putative sequences were aligned with MUSCLE in MEGA 7, to generate unrooted tree with Maximum Likelihood Method (Figure 2(a)). All the classes belonging to the eight groups of BraGSTs were in clusters with their counterparts of A. thaliana and rice. It also verified the reliability of our classification of BraGSTs, based on the conserved domain with completely matched results. The classification patterns of number BraGST genes for each subclasses are described in Table 1, while, for the comparative analysis with A. thaliana and rice, the results are shown in Figure 2 were marked with different color (Figure 2(c)). As described in Figure 2(c), subclass Tau (BraGSTU) has contained the largest (37) number of genes followed by Phi (22). Our results were in accordance with previously reported studies [7,36,47]. We also investigated the sequence features of BraGSTs through MEME program, which are used for predicting the conserved motifs; at least 12 motifs were identified and named as motifs 1-12. Through MEME, we also obtained the LOGO of these motifs and presented them along with their motifs (Figure 2(c)). The BraGSTs were mainly distributed with similar patterns of motifs according to the eight subclasses. However, motifs 1-2 were found almost among all the BraGSTs, suggesting their highly conserved domain. Besides the common patterns, subclass EFIG has contained  B ra G S T U 2 2 A T G S T U 2 7  (c)   higher number of motifs than any other subclass of BraGSTs.
Noticeably, motif 8 was found to have higher number (88) of consensus sequences while motifs 9, 10, and 12 contained fewer (16) number of consensus sequences. In addition, we also compared the genomic and cDNA sequences; we found that all the subclasses showed higher than one intron, while specifically subclass Lambda was found to contain the high (14) number of intron ( Figure 2(d)). However, the structure of Lambda was varied in nature in other genes, while DHAR and Phi showed uniformity in their structure. Notably, the gene structure of all the BraGSTs was consistent according to the subclasses.

Copy Number Variation/Gene Retention and Collinearity
Analysis of BraGSTs. To determine the copy number of variation during specific whole genome triplication event (WGT), we investigated the BraGST genes in A. thaliana and B. rapa and the collinear relationships are shown in (Figure 3(a)). Surprisingly, we found up to six numbers of copy variation between two different subclasses of BraGSTs, i.e., Phi and Tau, due to high tandem array of genes (Figure 3(b)). In subclass Tau, we explored five and four copies of genes, whereas single copy of genes was also found in subclass Tau. Our findings also demonstrated different retention of genes in the subclasses of BraGSTs, such as Phi 21/21, Tau 36/38, Zeta 3/3, Theta 2/2, Lambda 3/3, DHAR 4/4, EF1G 3/3, and TCHQD 1/1, respectively, and these results revealed high number of retention genes among all the subclasses of BraGSTs (Supplementary Table 3 and Figure 3(c)). The B. rapa contains three genomes, as based on the degree of fractionation, namely, least fractionation (LF), medium fractionated (MF1), and most fractionated (MF2). During this study, we presented the ratio of all the BraGSTs in the three subgenomes, LF genome showed the highest 40.54% of genes, and noticeably MF2 showed 32.43%, which were higher than MF1 27.03% (Figure 3(d)). On the other hand, we also highlighted the ratio of these three subgenomes with nonsynteny genes of BraGSTs; interestingly the nonsynteny ortholog and LF subgenome share 28.38% each and MF1 and MF2 share an equal of 21.62% as shown in (Figure 3(e)).

Chromosomal Localization and Tandem Array Selective
Pressure Analysis of the BraGSTs in B. rapa. All the BraGSTs were positioned on the ten B. rapa chromosomes (A01-A10) with a random distribution. Chromosome A09 contained the most BraGST genes (18.92%), followed by A07 with (17.57%), whereas chromosome A01 contained the fewest genes (1.35%). Furthermore, the duplication type (tandem array) was identified through MCScanX program; about 34 genes were identified with a random distribution on A01-A10 chromosomes and were marked with red color. Specifically, we observed that most of the tandem arrays were clustered in the region of chromosome. For example, 9 tandem array genes were clustered on A07 chromosome, all of them belonged to subclass Tau, and the only gene which was on scaffold also belonged to subclass Tau (Figure 4(a)). Additionally, we also presented a relationship of the syntenic region of BraGSTs along with A. thaliana using Circos software and the tandem array of BraGST genes was marked with red color (Figure 4(b)). A total of 15 pairs of tandem array were analyzed for selective pressure using Mega 7 and the divergence time of the duplicated genes was estimated by calculating the number of synonymous substitutions ( ) and nonsynonymous substitution rates ( ). All of the tandem pairs showed less than 1.00 / ratio, which suggest the purifying selection of the genes. The values were varied 0.13-0.77 with an average divergence of 29.80 MYA ( Table 2). These results speculated that the slow evolving nature of the tandem array of BraGSTs with a variation patterns in plant functions. The purifying nature of BraGSTs suggested the maintenance function in B. rapa.

Expression Pattern of BraGSTs in Different Tissues and the Corelation Networking Analysis of Tandem Pairs.
To investigate putative roles of BraGST genes in B. rapa growth and development, we analyzed the expression patterns across five tissues (roots, stems, leaves, flowers, and siliques) using RNAseq data [32]. Most of the BraGSTs demonstrated a high alteration in expression level in tissues-specific patterns as shown in (Figure 5(a)). For example, BraGTF18 and BraGSTF19, belonging to subclass Phi, exhibited higher expression in stem and root as compared to any other tissues, indicating that they may function in stem and root development. Meanwhile, some of the genes like from subclasses Phi, Tau, and DHAR such as BraGSTF10 and BraGST22; BraGSTU8 and BraGSTU32; and BraDHAR1 showed no expression, while BraGST11, BraGSTF17, BraGSTU3, BraGSTU19, BraGSTU33, and GSTZ3 have slight expression in any tissue (Figure 5(a) and Supplementary Table 4). We also presented tissuespecific genes; interestingly 2 genes were found in root and 1 was in silique, which highlighted the tissue-specific developmental role ( Figure 5(b)).
We also investigated trends of expression patterns among 15 pairs of tandem array and their Pearson correlation. Notably, our results demonstrated that most of the tandem pairs exhibited a high expression pattern across five tissues. Intriguingly, two tandem pairs from subclass of Phi, such as (BraGSTF19 BraGSTF20), showed a higher expression peak in root and stem, which could be speculated that these pairs due to the similarity in high expression may be involved on the same pathway for root and stem development. Meanwhile, the correlation between BraGSTU4 ATGSTU4 1 4 BraGST U5 ATG STU 3 93 43 Bra GST U3 AT GS TU 2 60 19 Br aG ST U8 AT A T G S T U 9 9 9 A T G S T U 1 0 B ra G S T U 1 0 B ra G S T U 2 4 70 Figure 4: (a) Chromosome location of the BraGSTs was obtained from the GFF file and displayed by using Mapchart. The duplication type tandem array was displayed with red color. (b) Syntenic relationship between B. rapa and A. thaliana was displayed through Circos program, while tandem array was marked with red color.  these two pairs was recorded (0.742818), which also shed light on their close relationship (Figure 6(a) and Supplementary  pseudogenization. In tissue-specific clustering, we observed only one specific gene located in root and silique each, suggesting their possible role in function of root and stem development ( Figure 6(b)).
To elucidate the coregulatory network among BraGSTs, we presented a protein-protein interaction structure, in which their functional and physical properties were examined by using STRING server (Figure 7). Most of the BraGSTs showed a highly interacting network with other proteins, except BraGSTF14. The protein interaction is not necessary to have a similar relation among the same subclasses but can also exist among different subclasses of a gene family. Predominately, most of the proteins were located in the center of the network were Phi, Zeta, Theta, and EF1G, suggesting that these proteins with other BraGSTs have a more complex interaction relationship.

Expression Profiling and Coregulatory Networks of BraGST Genes in Response to Multiple Hormone Treatments.
During recent times, many researchers primarily are focusing on how to understand plant functions under various stimulation. GSTs are to be involved in various biotic, abiotic, and hormone stresses such as auxins, ABA, and ethylene as reported in previous studies [48]. Keeping in view the importance of GSTs, we explored PlantCARE (plant cis-acting regulatory element database) for the identification of motifs only existing for hormonal stresses (Supplementary Table 6). The results demonstrated that most of the BraGST genes were involved in methyl jasmonate (32.88%), gibberellin (19.32%), abscisic acid (15.93%), salicylic acid (13.22%), and auxin (11.19%) and fewest genes were found in ethylene (7.46%) (Figure 8). Consequently, these results speculated that MeJA, GA, ABA, SA, auxin, and ethylene could affect the expression level of BraGSTs. For the functional dissection of GSTs under response to hormone stresses, mainly cis-element provided indirect evidences. However, little is known about its function in B. rapa response to hormonal stresses. To access the changes and the diversity of BraGSTs, further experimental validation step is required in future.
To understand the expression profiles of BraGST genes under eight different hormonal treatments, the expression patterns for 15 paralogous pairs of tandem array were studied using qRT-PCR experiment. Heatmap was generated in response to multiple hormone treatments for transcript expression fold change as shown in (Figure 9(a) and Supplementary Table 7). Most of the BraGSTs showed a high striking expression patterns during various hormone treatments (namely, ABA, GA, BR, ETH, IAA, IBA, NPA, and JA). More genes were particularly induced by NPA than any other hormones and expressed a higher proportion of upregulation (66%) and relatively lower downregulated level (34%), followed by BR (61%) and (39%), ABA (60%), and (40%). Meanwhile, JA were found to be sensitive to treatment as most of the genes were downregulated up to 74% and only 24% were upregulated (Figure 9(b)). To understand the correlation and coregulatory network among the tandem pairs, we also calculated the PCC values based on their relative expression values. For the correlation network, we have arranged them into three categories with respect to PCC values, such as greater than 0.6 (High), less than 0.5 but greater than 0 (Medium), and negative values with (Negative) correlation (Figure 10(a) and Supplementary Table 8). BR and IBA showed a higher and closer relationship among each other with 10 PCC values each, while ETH were found with a high number of 10 negative PCC values, suggesting its contrasting nature among the tandem pairs of genes.

Identification, Phylogeny, and Gene Duplication of
BraGSTs. In the present study, a total of 75 BraGST genes were identified in B. rapa, through genome-wide identification. Based on the domain information and phylogenetic analysis of the GSTs, they were categorized into eight subclasses in which Tau subclass was the most numerous with 37 number of genes and Phi with 22. Our analysis was further validated from previous findings with 15  dominant number of genes and the variation of copy number in plants [36,[49][50][51][52]. However, other subclasses have contained smaller number of genes with 1-4 members of GSTs. Furthermore, the results of our findings also showed a higher copy number of genes involved in the two subclasses Phi and Tau, which suggested that, after following whole genome duplication (WGD) event, it had a high degree of genes retention. As a result, our findings proved the gene dosage hypothesis for keeping network stability, following polyploidy, and underrepresented in copy variants the genes encoding members were preferentially retained [53,54]. In the evolutionary process, gene duplications are necessary for new biological functions and widespread expansion of gene family [55]. To understand duplication events, here we analyzed the expansion mechanism of BraGSTs family by MCScanX program. A large number (54.67%) of segmental type duplications were identified compared to tandem (45.33%). These analyses suggested that segmental duplication contributed in the expansion of BraGSTs family. Moreover, gene duplications played a major role in diversification by altering the genetic setup to assist them in adaptation to environmental stresses [56]. Chinese cabbage showed more close evolutionary resemblance to model plant A. thaliana [57]. They also shared the identical genomic composition but their sizes differ from each other which is ∼120 Mb in Arabidopsis [58,59]. Meanwhile, genome of B. rapa exceeds to 529 Mb, which is five times larger than that of Arabidopsis [60,61]. On the other hand, the modern diploid Brassica genome carried complexity among their three subgenomes, which are indicated in some of the preliminary studies [62][63][64]. In addition, we also found the divergence of GSTs with an average of 29.80 MYA, which has suggested a low selective pressure on BraGSTs that made them duplicate late, further demonstrated the contrasting nature with different function in subclasses of GSTs as discussed earlier. Moreover, on our results' basis, significant variation in rate of divergence reflected the similarities of BraGST sequences possibly matched with A. thaliana and rice [65,66]. Additionally, all the selected tandem pairs showed lower than 1.00 synonymous values, which lie in purifying selection nature and help in the maintenance of BraGST gene functions. The results were further validated by the synonymous values which did not show any significant differences among three subgenomes of B. rapa (LF, MF1, and MF2). The comparative analysis, phylogenetic tree, and gene structure of GSTs revealed the exon-intron accordingly to their subclasses. The similarities in genes structure and the similar patterns of exon-intron were further proved by analyzing the protein highly conserved sequences with MEME. The regulation patterns for conserved motif were tight specifically linked with motifs 1 and 2, as these were determined the most prominent in the BraGSTs family.

Expression Divergence and Regulatory Network of
BraGSTs. To validate our results, we also identified cisregulatory elements in the promoter regions of BraGSTs specific to hormonal stress. We observed that the majority of GSTs were involved in the activity of methyl jasmonate. Therefore, we can speculate that the function of BraGST genes was involved in the regulatory phenomena of hormones. In addition, the expansion of this larger gene family and the duplicated genes in both models (neofunctionalization or subfunctionalization) was more associated with process of tissue-expression divergence [67][68][69]. We also examined the tissue-expression patterns of the BraGST genes; most of them were expressed across five tissues or several at least. Three genes were tissue-specific (two in roots and one in silique) while some of them showed similar patterns. That might suggested their common importance in the regulation of plant developments. The tissue-specific genes along with the tandem pairs might be contributing to plant development by acquiring new functions. Based on the results of RNAseq data for various tissues and the response of prompter sequence analysis for cis-elements, their involvement in B. rapa response to different stresses was speculated. Thus, an opportunity to study the expression profile of BraGST genes was provided as to understand gene functions; the expression profiling can provide valuable clues [70]. Furthermore, the expression profile for 15 paralogous pairs was analyzed by qRT-PCR after application of eight multiple hormonal treatments (namely, ABA, GA, BR, ETH, IAA, IBA, NPA, and JA). In general, most of the plant growth hormones act as components of abiotic stress signaling, for instance, ABA, GA, Auxin, BR, cytokinin, and JA [71]. A large variety of cellular processes are regulated by phytohormones, although these are produced in minute concertation. To communicate cellular activities in higher plants, they work as chemical messengers [72]. The ratio for majority of treatments was upregulated in our study, such as for NPA stress (66%) and BR (61%); however JA was more sensitive as a large number (74%) of genes were downregulated. These results show that BraGSTs may also function in response to multiple hormone treatments, although most of their genes functions are still unknown. However, our analysis for phylogenetics, cis-elements, and expression profiling provides a foundation for future studies on BraGST gene functions. Although, based on PCC values, BR and IBA stress were among the highest with 10 PCC values each (>0.6) which signify a close relationship. Here, we speculated that the function of gene was enhanced and expanded through   gene duplications. During the process of evolution after duplications, the divergence in the expression profile of the tandem pairs revealed that it may acquire new functions; however functional analysis will confirm and determine the pivotal role of BraGSTs. To understand regulatory network, protein-protein interactions were elucidated and most of the genes presented a close relationship among subclasses, except for few genes. All the BraGST genes displayed a very complicated correlation, suggesting that these genes are involved in several fundamental mechanisms and are further regulated by many down-/upstream genes. Taken together our results, this study may provide novel insight into the unique features and specifically the role of the BraGSTs family in eukaryotic organisms.

Conclusion
In this study, we identified 75 BraGST genes from B. rapa and focused on those involved in response to multiple hormonal treatments as GSTs are playing a crucial role in plants. The classification, phylogenetic relationship, structural composition, evolutionary characteristics, and conserved protein motif analyses were investigated. Our study has provided a deep understanding of the BraGSTs in B. rapa. The differential expression patterns of BraGSTs in various tissues and the visible tissue-specific patterns showed that these genes are playing a key role in the developmental aspects of B. rapa. Expression analysis highlighted the involvement of BraGST genes in response to multiple hormone treatments. Furthermore, a highly interacting network and the correlation among various treatments demonstrated the importance of our study. These results will lead to novel insight by facilitating into functional divergence and will provide an assessment for further studies to understand the physiological function of GSTs in response to multiple hormone treatments in B. rapa. Table 1: Sequences of the BraGST gene primers used for quantitative real-time PCR. Table 2: The basic description of BraGST genes in Brassica rapa. Table 3: Identification of BraGST syntenic genes between A. thaliana along with three subgenomes of B. rapa. Table 4: The FPKM values of BraGST genes. Table 5: Syntenic paralog pairs of BraGSTs with PC and FPKM values. Table 6: Cis-elements of BraGST genes in Brassica rapa. Table 7: Relative expression pattern of BraGST genes along with PC values with respect to hormonal stresses by qRT-PCR. Table 8: Pearson correlation coefficient of the stress-induced BraGSTs whose PC is greater than 0.5. Figure  S1: Genetic distance among different subclasses of BraGSTs. Figure S2: Genetic distance among different subclasses of BraGSTs. Figure S3: Relative shares of different family among three subgenomes of B. rapa. Figure S4: Evolutionary tree of BraGSTs and the different subclasses being displayed with random color and tree generated with MEGA 7 using 1000bootstrap replicate values. (Supplementary Materials)