Genome-Wide Comprehensive Analysis the Molecular Phylogenetic Evaluation and Tissue-Specific Expression of SABATH Gene Family in Salvia miltiorrhiza

The plant SABATH gene family is a group of O-methyltransferases (O-MTs), which belongs to the S-adenosyl-l-methionine-dependent methyltransferases (SAM-MTs). The resulting reaction products of SABATH genes play an important role in various processes of plant development. In this study, a total of 30 SABATH genes were detected in Salvia miltiorrhiza, which is an important medicinal plant, widely used to treat cardiovascular disease. Multiple sequence alignment and phylogenetic analyses showed that SmSABATH genes could be classified into three groups. The ratios of non-synonymous (Ka) and synonymous (Ks) substitution rates of 11 pairs paralogous of SmSABATH genes revealed that the SmSABATH genes had gone through purifying selection. Positive selection analyses using site models and branch-site models indicated that SmSABATH genes had undergone selective pressure for adaptive evolution. Functional divergence analyses suggested that the SmSABATH subgroup genes were divergent in terms of functions and positive selection sites that contributed to a functional divergence among the subgroups that were detected. Tissue-specific expression showed that the SABATH gene family in S. miltiorrhiza was primarily expressed in stems and leaves.


Introduction
Methylation is a ubiquitous reaction that takes place in bacteria, fungi, plants, and mammals. The process of methylation is catalyzed by S-adenosyl-L-methionine-dependent methyltransferases (SAM-MTs), and involves the transfer of the methyl group of S-adenosyl-L-methionine (SAM) to carbon, nitrogen, oxygen, or sulfur atoms, and modifies DNA, RNA, proteins, or small molecules with the formation of corresponding methylated products and S-adenosyl-L-homocysteine (SAH) [1]. Enzymatic methylation of hydroxyl and carboxyl moieties are catalyzed by O-methyltransferases (O-MTs) [2], of which there are three defined types, via protein X-ray crystallography [3][4][5]. Type 1 O-MTs exclusively methylate oxygen atoms of the hydroxyl moieties of phenylpropanoid-based compounds [6] and type 2 O-MTs are specific to phenylpropanoid esters of the coenzyme A, and are found in all lignin-producing plants [3]. The type 3 O-MTs specifically methylate carboxyl groups of

Ka and Ks Calculation
The paralogs for the SmSABATH genes were inferred from the phylogenetic tree that was explained in Section 2.3. Non-synonymous (Ka) and synonymous (Ks) substitution rates, and the Ka/Ks ratio of each paralogous gene pair, were determined by PAL2NAL program (http://www.bork.embl.de/pal2nal/) [41]. Meanwhile, the Ka/Ks ratios for all of the paralogous genes were calculated with a sliding window of 20 aa [42].

Tests of Positive Selection
To determine whether the SmSABATH gene family exhibited evidence of positive selection under the site models and branch-site models [43], the codeml program in PAML v4.9a was applied to test the hypothesis of positive selection. We reconstructed the phylogenetic tree with the amino acid sequences of SmSABATH under the model of JTT + I + G in MrBayes [35,36]. The selected model also used the program ProtTest [37]. In the site model, M0 (one ratio), M3 (discrete), M1a (neutral), M2a (selection), M7 (beta), and M8 (beta & ω) were applied to the alignments, and the variation in the ω parameter among sites were detected using LRT (likelihood ratio test) for M0 vs. M3, M1a vs. M2a, and M7 vs. M8. The branch-site model [44] was used to compare the Ka/Ks ratio between the branches. The positive selection amino acid sites of SmSABATH were tested by the improved branch-site model [44]. The branches that were tested for positive selection were used as the foreground, while all of the other branches on the tree were used as the background. Furthermore, the ratio of non-synonymous and synonymous substitution rates for each branch was calculated under the Null Model and Alternative Model. In the Null Model, the omega was set as 1; in the Alternative Model the omega was set as >1. The positive selection sites were detected by comparing the significance level between the Null Model and Alternative Model, using LRT. If LRT suggested the presence of codons under positive selection on the foreground branch, the codon was probably from the site class of positive selection [45]. Posterior probabilities (Qks) were estimated with the Bayes Empirical Bayes (BEB) method [46].

Estimation of Functional Divergence
An analysis of the functional divergence between the SmSABATH genes of subgroups was performed using DIVERGE version 3.0 [47]. The method could estimate significant changes in the site-specific shifts, based on maximum likelihood procedures [48]. The estimation was based on the neighbor-joining tree, which was reconstructed with SmSABATH amino acid sequences using MEGA 6.0 [49], and the coefficients of Type-I and Type-II functional divergences (θ I and θ II ) between two clusters were calculated. The coefficients of Type-I and Type-II functional divergences (θ I and θ II ) that were greater than 0 indicated that site-specific altered selective constraints, and a radical shift in amino acid physiochemical properties occurred after gene duplication or speciation [50]. The posterior probabilities (Qks) of amino acid sites that were responsible for functional divergence could also be estimated by this program. A large posterior probability (Qk) represented a high possibility that the evolutionary rate or the radical change in the amino acid property of a site was different between two clusters [50]. Additionally, Qk > 0.8 was empirically used as the cutoff in the identification of Type-I functional divergence-related residues between gene groups. Meanwhile, Type-II functional divergence divergence-related residues were identified with Qk > 1.0 as the cutoff [48].

Plant Materials, RNA Extraction and Real-Time qPCR
Roots, stems, leaves, and flowers were collected from tow-year-old, field-grown S. miltiorrhiza Bunge plants from Shangluo County, Shaanxi Province, China, and were stored in liquid nitrogen until use. Total RNA was extracted from the S. miltiorrhiza tissues of roots, stems, leaves, and flowers, using the Quick RNA Isolation Kit (Huayueyang, Beijing, China). RNA quantity was determined using a NanoDrop 2000C Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The first-strand cDNA was synthesized using Prime-Script RT Master Mix (TaKaRa, Beijing, China), according to the manufacture's protocol. Real-time qPCR was performed using a Light Cycler 96 Instrument (Roche, Basel, Switzerland). The reaction mixture contained 10 µL of 2× SYBR Premix Ex Taq II (Takara, Beijing, China), 20 ng of first-strand cDNA, and 500 nM each of sense and antisense primers. Initial thermal-cycling at 95 • C for 30 s was followed by 45 cycles of 95 • C for 10 s and 60 • C for 30 s. Relative expression was calculated by the 2 −∆∆Ct method [51]. All real-time qPCRs were repeated in three biological and three technical replicates. The relative expressions were analyzed as means ± standard deviation (SD). The lengths of the amplicons were between 100 bp and 250 bp. According to Yang et al.'s method [52] for screening the reference genes for real-time qPCR analysis in various tissues of S. miltiorrhiza, Smβ-actin (DQ243702) was selected as a reference gene. Gene-specific primers were designed using Premier 5.0 and are listed in Table S2.

Sequence Feature of SABATH Genes in S. miltiorrhiza
In order to identify SABATH genes in the S. miltiorrhiza, tBLASTn analyses against the S. miltiorrhiza genome was performed using AtSABATH amino acid sequences as queries. With the BLAST search, a total of 30 SmSABATH genes were detected. The gene lengths of SmSABATH varied from 1010 bp (SMil_00022020) to 5649 bp (SMil_00022342) (Table S3), and the lengths of the SmSABATH cDNAs and proteins varied from 561 bp and 186 aa (SMil_00022020) to 1293 bp and 430 aa (SMil_00023670) ( Table S3). The molecular weights of the predicted proteins ranged from 20.01 kDa (SMil_00022020) to 47.70 kDa (SMil_00023670) (Table S3), and the theoretical isoelectric points were predicted to range from 5.03 (SMil_00021702) to 9.58 (SMil_00008156) (Table S3).

Phylogenetic Analysis of SABATH Gene Family in S. miltiorrhiza
To detect the evolutionary relationship of the SABATH gene family in S. miltiorrhiza with the members of SABATHs from other species, an un-rooted phylogenetic tree was constructed using the SmSABATH, AtSABATH and known SABATHs from other species amino acid sequences ( Figure 1). Based on the phylogenetic tree, all of the SABATHs were divided into three major groups.
The bootstrap values for all of the subgroups were high, suggesting that genes in the same subgroup might share a similar origin (Figure 1). Nearly all of the kindred plant clustered together, and the same function SABATHs were also classed into the same clade ( Figure 1). In other words, the SABATHs with higher homology maybe have the same attributes and same function. Thus, it suggested that we could infer the function of unknown SABATHs according to the clustering situation. SMil_00008666 and SMil_00025720, SMil_00001154 and SMil_00022343 in Group A; SMil_00007297 and SMil_00010605, SMil_00007747 and SMil_00026995, SMil_00007772 and SMil_00021640, SMil_00023670 and SMil_00028890 in Group B; SMil_00008156 and SMil_00021702, SMil_00028867, and SMil_00030124 in Group C (Figure 1). In addition, one pairs of orthologs genes of SMil_00017556 and At1g19640 regarding to the SABATH gene family in S. miltiorrhiza were identified (Figure 1), which may have the same function.

Gene Structure Analysis of SABATH Genes in S. miltiorrhiza
The structural features of each of the SABATH genes in groups are shown in Figure 2. Structural analyses of all of the SABATH genes in S. miltiorrhiza revealed that the number of exons varied from two to five. There were no intronless genes. The average exon numbers in the groups were two or three. We also found that genes in the same group had similar gene structures (Figure 2), which indicated that the genes in same group may have similar functions. Otherwise, symmetric exons were the exons that have the same splicing phase at both ends. The excess of symmetric exons and phase 0 introns are likely to facilitate exon shuffling, recombinational fusion, and protein domain exchange [55,56]. According to the 92 exons analyzed herein, 22 exons were symmetric with phase 0 introns, and no exons were symmetric with phase 1 and 2 introns. Among the 62 introns of the SABATH genes, 53 were phase 0, six were phase 1, and three were phase 2. Therefore, our analyses of the gene structures indicated diversity amongst the SABATH genes in S. miltiorrhiza. For the SABATHs in S. miltiorrhiza, all of the SmSABATH genes were also divided into three major groups. The fourteen and five SmSABATHs clustered together (Group A and Group C), while the other eleven SmSABATHs were divided among other species (Group B) ( Figure 1). This suggested that the members in Group B may be divergent with other groups in function. Another main objective of this phylogenetic study was to identify putative orthologous and paralogous. Paralogs usually display different functions while orthologs may retain the same function [53,54]. According to the phylogenetic tree (Figure 1), among the 30 SmSABATH genes, 11 pairs of paralogous were identified from the SABATH gene family in S. miltiorrhiza including SMil_00003309 and SMil_00022342, SMil_00020191 and SMil_00020192, SMil_00003310 and SMil_00022020, SMil_00008666 and SMil_00025720, SMil_00001154 and SMil_00022343 in Group A; SMil_00007297 and SMil_00010605, SMil_00007747 and SMil_00026995, SMil_00007772 and SMil_00021640, SMil_00023670 and SMil_00028890 in Group B; SMil_00008156 and SMil_00021702, SMil_00028867, and SMil_00030124 in Group C ( Figure 1). In addition, one pairs of orthologs genes of SMil_00017556 and At1g19640 regarding to the SABATH gene family in S. miltiorrhiza were identified (Figure 1), which may have the same function.

Gene Structure Analysis of SABATH Genes in S. miltiorrhiza
The structural features of each of the SABATH genes in groups are shown in Figure 2. Structural analyses of all of the SABATH genes in S. miltiorrhiza revealed that the number of exons varied from two to five. There were no intronless genes. The average exon numbers in the groups were two or three. We also found that genes in the same group had similar gene structures (Figure 2), which indicated that the genes in same group may have similar functions. Otherwise, symmetric exons were the exons that have the same splicing phase at both ends. The excess of symmetric exons and phase 0 introns are likely to facilitate exon shuffling, recombinational fusion, and protein domain exchange [55,56]. According to the 92 exons analyzed herein, 22 exons were symmetric with phase 0 introns, and no exons were symmetric with phase 1 and 2 introns. Among the 62 introns of the SABATH genes, 53 were phase 0, six were phase 1, and three were phase 2. Therefore, our analyses of the gene structures indicated diversity amongst the SABATH genes in S. miltiorrhiza.

Analysis of Conserved Domains and Motifs
The SABATH domains that are contained in all of the deduced SABATH proteins in S. miltiorrhiza were identified using the Pfam program, and the results showed that all members of the SmSABATH gene family contained a conserved domain, which was conserved among O-methyltransferases [57,58]. Then, the conserved domain sequences were aligned using software DNAMAN; most members of SmSABATH proteins contained binding sites (the motifs I and III) of SAM (S-adenosyl-L-methionine) ( Figure 3), a well-known methyl donor in plant cells [11]. Using the MEME program, we identified 13 conserved motifs in SmSABATH amino acid sequences (Table S4). The lengths of the motifs varied from 14 to 50 amino acids, and the number of motifs in each SmSABATH ranged from 4 to 11. The frequency of all 13 conserved motifs in SmSABATH proteins varied from 7 to 28 (Table S4). Some motifs, such as motifs 1, 2, 3, 4, and 8, were diffusely distributed among SmSABATH proteins (Table S5). Among the 13 conserved motifs, motifs 1, 2, 3 and 4 were located in the SmSABATH conserved domain, while the other nine motifs were located outside the conserved domain. Otherwise, motif 4 matched motif I while motif 2 matched motif III (Figure 3). Based on the groups indicated by the phylogenetic tree, most SmSABATH proteins in the same group had similar motifs compositions and the orders of the motifs were very similar in each group (Table S5), which suggested that the proteins in the same group might have similar functions in plant development. Some specific motifs were located in the proteins of specific groups; for instance,

Analysis of Conserved Domains and Motifs
The SABATH domains that are contained in all of the deduced SABATH proteins in S. miltiorrhiza were identified using the Pfam program, and the results showed that all members of the SmSABATH gene family contained a conserved domain, which was conserved among O-methyltransferases [57,58]. Then, the conserved domain sequences were aligned using software DNAMAN; most members of SmSABATH proteins contained binding sites (the motifs I and III) of SAM (S-adenosyl-L-methionine) (Figure 3), a well-known methyl donor in plant cells [11]. Using the MEME program, we identified 13 conserved motifs in SmSABATH amino acid sequences (Table S4). The lengths of the motifs varied from 14 to 50 amino acids, and the number of motifs in each SmSABATH ranged from 4 to 11. The frequency of all 13 conserved motifs in SmSABATH proteins varied from 7 to 28 (Table S4). Some motifs, such as motifs 1, 2, 3, 4, and 8, were diffusely distributed among SmSABATH proteins (Table S5). Among the 13 conserved motifs, motifs 1, 2, 3 and 4 were located in the SmSABATH conserved domain, while the other nine motifs were located outside the conserved domain. Otherwise, motif 4 matched motif I while motif 2 matched motif III (Figure 3). Based on the groups indicated by the phylogenetic tree, most SmSABATH proteins in the same group had similar motifs compositions and the orders of the motifs were very similar in each group (Table S5), which suggested that the proteins in the same group might have similar functions in plant development. Some specific motifs were located in the proteins of specific groups; for instance, motif 7 was specific to Groups A and B (Table  S5). This indicated that the SmSABATH proteins may exhibit functional divergence in different groups.

Driving Forces for Genetic Divergence
Gene duplication is an important event for gene family expansion and functional diversity during evolution [59]. To detect whether Darwinian positive selection was involved in the driving of gene divergence after duplication, the Ka/Ks ratios of non-synonymous and synonymous substitution rates were calculated using the CDS of paralogous SmSABATH. Generally, a Ka/Ks ratio < 1 indicates a negative or purifying selection, a ratio = 1 indicates a neutral evolution and a ratio > 1 indicates a positive selection [60]. According to the paralogous genes determined from the phylogenetic tree, there were 11 pairs of paralogs in SmSABATH genes (Figure 1). In other words, more than 73% of SmSABATH genes appeared to be duplicated. This suggested that most SmSABATH genes have happened functional diversity and gene family expansion during evolution. The Ka/Ks ratios for all of the 11 SmSABATH paralogous pairs were <1 (Table S6), suggesting that the SmSABATH genes have experienced purifying selection pressure. Meanwhile, we also calculated the Ka/Ks ratios for all of the paralogous genes with a sliding window of 20 aa. The ratios of Ka/Ks that were higher than 1 in the regions of all paralogous genes indicated that these regions had gone through positive selections. However, the proportion of such regions were few. Most of the regions

Driving Forces for Genetic Divergence
Gene duplication is an important event for gene family expansion and functional diversity during evolution [59]. To detect whether Darwinian positive selection was involved in the driving of gene divergence after duplication, the Ka/Ks ratios of non-synonymous and synonymous substitution rates were calculated using the CDS of paralogous SmSABATH. Generally, a Ka/Ks ratio < 1 indicates a negative or purifying selection, a ratio = 1 indicates a neutral evolution and a ratio > 1 indicates a positive selection [60]. According to the paralogous genes determined from the phylogenetic tree, there were 11 pairs of paralogs in SmSABATH genes (Figure 1). In other words, more than 73% of SmSABATH genes appeared to be duplicated. This suggested that most SmSABATH genes have happened functional diversity and gene family expansion during evolution. The Ka/Ks ratios for all of the 11 SmSABATH paralogous pairs were <1 (Table S6), suggesting that the SmSABATH genes have experienced purifying selection pressure. Meanwhile, we also calculated the Ka/Ks ratios for all of the paralogous genes with a sliding window of 20 aa. The ratios of Ka/Ks that were higher than 1 in the regions of all paralogous genes indicated that these regions had gone through positive selections. However, the proportion of such regions were few. Most of the regions showed that the ratios of Ka/Ks were less than 1 in paralogous genes, and this also indicated that the SABATH genes in S. miltiorrhiza had undergone purifying selection (Figure 4).

Positive Selection on SmSABATH Genes
In order to preliminarily examine the evolutionary mechanisms of the SmSABATH gene family, we tested the hypothesis of positive selection of the genes using site models and branch-site models in PAML program [43,46] base on the phylogenetic tree ( Figure S1).
In site models, there were six codon substitution models [43], including M0, M1a, M2a, M3, M7, and M8, which were applied to the alignments and these models assume variation in ω among sites. The results of these models' parameter estimates, log likelihood, and the LRT tests are shown in Table 1. To examine how dN/dS (non-synonymous/synonymous) ratios differed among codon positions, models M0 and M3 were compared. M0 (one ratio) assumes that the different sites have the same evolution rate, while M3 (discrete) assumes a general discrete distribution with three site classes (p0, p1, p2) [48]. The log likelihood of M0 for SmSABATH sequences was ι = −24,100.940606, with an estimate of ω = 0.32227. The log likelihood of M3 was ι = −23,631.868500, with an estimate of

Positive Selection on SmSABATH Genes
In order to preliminarily examine the evolutionary mechanisms of the SmSABATH gene family, we tested the hypothesis of positive selection of the genes using site models and branch-site models in PAML program [43,46] base on the phylogenetic tree ( Figure S1).
In site models, there were six codon substitution models [43], including M0, M1a, M2a, M3, M7, and M8, which were applied to the alignments and these models assume variation in ω among sites. The results of these models' parameter estimates, log likelihood, and the LRT tests are shown in Table 1. To examine how dN/dS (non-synonymous/synonymous) ratios differed among codon positions, models M0 and M3 were compared. M0 (one ratio) assumes that the different sites have the same evolution rate, while M3 (discrete) assumes a general discrete distribution with three site classes (p 0 , p 1 , p 2 ) [48]. The log likelihood of M0 for SmSABATH sequences was ι = −24,100.940606, with an estimate of ω = 0.32227. The log likelihood of M3 was ι = −23,631.868500, with an estimate of ω 0 = 0.06316, ω 1 = 0.25310, and ω 2 = 0.65803 (Table 1). These results indicated that all of the codons were under purifying selection. Additionally, the value of twice the log likelihood difference (2∆lnL) between M3 and M0 was 938.1442, which was strongly statistically significant (p < 0.01) and suggested that M3 was better than M0. Therefore, the results indicated that different sites bare different selection pressures and also indicated fluctuations in the overall level of selective constraints.
The codon substitution models of M2a and M8 allow for positive selection, while the models of M1a and M7 hypothesize a nearly neutral selection [48]. We could test whether positive selection promoted divergence between genes through comparing the M2a vs. M1a and M8 vs. M7, respectively. The log likelihood of M1a and M2a for SmSABATH genes was ι = −23,766.855086. The value of 2∆lnL between M1a and M2a was 0; it was not statistically significant, and no sites were positively selected at a level of 95% (Table 1). The log likelihood of M7 and M8 for SmSABATH genes was ι = −23,622.773935 and ι = −23,622.775486, respectively. The value of 2∆lnL between M7 and M8 was close to 0 (Table 1); it was also not statistically significant. Although six sites were detected, they were not positively selected at a level of 95% (Table 1). In both cases, no significant evidence of positive selection was found. Branch-site models were designed to detect positive selection that was affecting a few sites along particular lineages, and allowed ω ratios to simultaneously vary among sites and branches [46]. The parameter estimates for branches under positive selection are list in Table 2. When Group A was set as the foreground branch, the value of 2∆lnL between the Null and Alternative models was close to 0, and no positive sites were found ( Table 2). In addition, three sites were found when Group B was set as the foreground branch, but they were not positively selected at a level of 95%, and the value of 2∆lnL between Null and Alternative models was 5.653506 (p = 0.017) ( Table 2). When Group C was set as the foreground branch, the value of 2∆lnL between Null and Alternative models was 12.07159 (p < 0.01) ( Table 2). A total of four sites were found but only one site in Group C was positively selected at a level of 95% (Table 2). In other words, only one positive site and one lineage group were found to be under positive selection. Moreover, the results suggested that different SmSABATH lineages may have different evolutionary rates. Group A evolution seemed to be more conservative; Group B might be confronted with positive Darwinian selection, but no positive sites were detected; and, Group C could be confronted with strong positive Darwinian selection, since significant positive sites were detected at the 0.05 significance level (Table 2).

Functional Divergence Analysis (FDA) of SmSABATH Proteins
With the program DIVERGE 3.0, the shifted evolutionary rates and altered amino acid properties after gene duplication could be evaluated [50,61]. The Type-I functional divergence (θ I ) was based on evolutionary rate [50] and the Type-II functional divergence (θ II ) was based on differences in the biochemical properties of amino acids [61]. According to the neighbor-joining tree, the SmSABATH amino acid sequences were also divided into three major clusters (Cluster A, Cluster B, and Cluster C) ( Figure S2). Posterior probability (Qk) was carried out to estimate Type-I and Type-II between the SmSABATH clusters.
Using DIVERGE, we found that all of the coefficients for Type-I functional divergence (θ I ) were greater than zero in the three group pairs, Group A vs. Group B, Group A vs. Group C, and Group B vs. Group C, and ranged from 0.286372 to 0.683167 (Table S7). This suggested that certain amino acid sites may have experienced significant site-specific changes between these group pairs, leading to a subgroup-specific functional evolution after their diversification. We also found that Type-I functional divergence (θ I ) of Group A vs. Group B and Group B vs. Group C were statistically significant (p = 0.016727 and p = 0.015877, respectively) (Table S7), while Group A vs. Group C was not statistically significant (p = 0.198275). This indicated that the members in Group B may be different from those of other groups in terms of function. This also verifies the phylogenetic analyses that the members in Group B have divergent with those of other groups in function. Some residues in this group probably play important roles in the functional divergence of SmSABATH genes. The coefficients for Type-II functional divergence (θ II ) in two group pairs, Group A vs. Group B, Group A vs. Group C, were greater than zero (0.041610 and 0.117169, respectively), which indicated a radical shift in amino acid properties, while the coefficients in Group B vs. Group C was less than zero (Table S8). However, the Type-II coefficients were not statistically significant among the three group pairs (p = 0.468158, p = 0.370624, and p = 0.312775, respectively) (Table S8), which suggested that most residues in the SmSABATH family should not experience obvious physical and chemical property changes.
To extensively reduce positive false, we used Qk > 0.8 and 1.0 as the cutoff to identify Type-I and Type-II functional divergence-related positive selection sites between gene groups, respectively. Using the detailed posterior probabilities analysis, we found that the distribution and the number of positive selection sites for functional divergence in group pairs were different. In Type-I functional divergence, when Qk > 0.8, all three of the group pairs contained positive selection sites (Table S7). Meanwhile, in Type-II functional divergence, when Qk > 1.0, we found that only one group pair (Group A vs. Group C) contained three positive selection sites. No positive selection sites were found in other tow group pairs (Table S8). This suggested that these sites probably play an important part in the functional divergence of SmSABATH during the evolutionary process. A detailed distribution of site-specific predictions for Type-I and Type-II functional divergence of SmSABATH between groups are showed in Figures 5 and 6.

Tissue-Specific Expression of SABATH Gene Family in S. miltiorrhiza
In order to preliminarily detect the tissue-specific expression of the SABATH gene family in S. miltiorrhiza, we analyzed the expression of SmSABATH in the roots, stems, leaves, and flowers of S. miltiorrhiza plants. All of the 30 SmSABATHs that were identified exhibited differential expression patterns (Figure 7). Among the 30 SmSABATHs, three (10.0%) showed predominant expression in roots, seven (23.3%) mainly expressed in leaves, 13 (43.3%) were primarily expressed in stems, and two (6.7%) were in flowers. In addition, five (16.7%) SmSABATHs were mainly expressed in at least two tissues analyzed, indicating that these genes may play a more ubiquitous role in S. miltiorrhiza. On the whole, the SABATH gene family in S. miltiorrhiza was primarily expressed in stems and leaves. This suggested that the SABATH gene family is related to the plant defense against insects that gnaw at the leaves and stems of plants. Once a plant is damaged from the outside environment, the products of SABATH genes, such as methyl jasmonate or methyl salicylate, can enhance the level of plant resistance [11].

Tissue-Specific Expression of SABATH Gene Family in S. miltiorrhiza
In order to preliminarily detect the tissue-specific expression of the SABATH gene family in S. miltiorrhiza, we analyzed the expression of SmSABATH in the roots, stems, leaves, and flowers of S. miltiorrhiza plants. All of the 30 SmSABATHs that were identified exhibited differential expression patterns (Figure 7). Among the 30 SmSABATHs, three (10.0%) showed predominant expression in roots, seven (23.3%) mainly expressed in leaves, 13 (43.3%) were primarily expressed in stems, and two (6.7%) were in flowers. In addition, five (16.7%) SmSABATHs were mainly expressed in at least two tissues analyzed, indicating that these genes may play a more ubiquitous role in S. miltiorrhiza. On the whole, the SABATH gene family in S. miltiorrhiza was primarily expressed in stems and leaves. This suggested that the SABATH gene family is related to the plant defense against insects that gnaw at the leaves and stems of plants. Once a plant is damaged from the outside environment, the products of SABATH genes, such as methyl jasmonate or methyl salicylate, can enhance the level of plant resistance [11]. Site-specific profile for predicting critical amino acid residues responsible for type-II functional divergence between groups of SmSABATH. The X-axis represents the locations of sites. The Y-axis represents the probabilities of each group. The red line indicates a cutoff = 1.0.

Tissue-Specific Expression of SABATH Gene Family in S. miltiorrhiza
In order to preliminarily detect the tissue-specific expression of the SABATH gene family in S. miltiorrhiza, we analyzed the expression of SmSABATH in the roots, stems, leaves, and flowers of S. miltiorrhiza plants. All of the 30 SmSABATHs that were identified exhibited differential expression patterns (Figure 7). Among the 30 SmSABATHs, three (10.0%) showed predominant expression in roots, seven (23.3%) mainly expressed in leaves, 13 (43.3%) were primarily expressed in stems, and two (6.7%) were in flowers. In addition, five (16.7%) SmSABATHs were mainly expressed in at least two tissues analyzed, indicating that these genes may play a more ubiquitous role in S. miltiorrhiza. On the whole, the SABATH gene family in S. miltiorrhiza was primarily expressed in stems and leaves. This suggested that the SABATH gene family is related to the plant defense against insects that gnaw at the leaves and stems of plants. Once a plant is damaged from the outside environment, the products of SABATH genes, such as methyl jasmonate or methyl salicylate, can enhance the level of plant resistance [11].

Figure 7.
Tissue-specific expression of the SABATH gene family in S. miltiorrhiza. The expression level of SmSABATHs was analyzed by the 2 −ΔΔCt method. Y-axis indicates the relative expression levels. X-axis indicates different tissues. Transcript levels in roots were arbitrarily set to 1, and the levels in other tissues were determined relative to this. Error bars represent standard deviations of mean value from three biological and three technical replicates. Analysis of variance was calculated using SPSS. p < 0.05 was considered to be statistically significant.

Conclusions
SABATH genes are ubiquitous in higher plants and play important roles in various processes of plant development. Further studies on this family could not only illustrate the SABATH genes' vital functions in the developmental processes of higher plants, but also elucidate the evolutionary relationships between different species. In this study, we identified 30 members of SABATH genes in the S. miltiorrhiza genome database. The identified SmSABATH genes were characterized using a comprehensive approach, including gene structure analyses, SABATH domain characterization, phylogenetic analyses, conserved motif identification, positive selection analyses, functional divergence analyses, and tissue-specific expression. We showed that 30 SmSABATH genes could be divided into three groups, according to the phylogenetic tree. SABATH genes in S. miltiorrhiza have experienced strong purifying selection pressures, since the Ka/Ks ratios for all 11 paralogous were <1. A total of 13 conserved motifs were identified, of which some group-specific motifs might be attributable to the functional divergence of SABATH genes in S. miltiorrhiza. Functional divergence analyses also showed that the SmSABATH genes have diverged in terms of function. Positive selection analyses with site model and branch-site model showed that SABATH genes in S. miltiorrhiza experienced positive selection. Tissue-specific expression showed that the SABATH Figure 7. Tissue-specific expression of the SABATH gene family in S. miltiorrhiza. The expression level of SmSABATHs was analyzed by the 2 −∆∆Ct method. Y-axis indicates the relative expression levels. X-axis indicates different tissues. Transcript levels in roots were arbitrarily set to 1, and the levels in other tissues were determined relative to this. Error bars represent standard deviations of mean value from three biological and three technical replicates. Analysis of variance was calculated using SPSS. p < 0.05 was considered to be statistically significant.

Conclusions
SABATH genes are ubiquitous in higher plants and play important roles in various processes of plant development. Further studies on this family could not only illustrate the SABATH genes' vital functions in the developmental processes of higher plants, but also elucidate the evolutionary relationships between different species. In this study, we identified 30 members of SABATH genes in the S. miltiorrhiza genome database. The identified SmSABATH genes were characterized using a comprehensive approach, including gene structure analyses, SABATH domain characterization, phylogenetic analyses, conserved motif identification, positive selection analyses, functional divergence analyses, and tissue-specific expression. We showed that 30 SmSABATH genes could be divided into three groups, according to the phylogenetic tree. SABATH genes in S. miltiorrhiza have experienced strong purifying selection pressures, since the Ka/Ks ratios for all 11 paralogous were <1. A total of 13 conserved motifs were identified, of which some group-specific motifs might be attributable to the functional divergence of SABATH genes in S. miltiorrhiza. Functional divergence analyses also showed that the SmSABATH genes have diverged in terms of function. Positive selection analyses with site model and branch-site model showed that SABATH genes in S. miltiorrhiza experienced positive selection. Tissue-specific expression showed that the SABATH gene family in S. miltiorrhiza is primarily expressed in stems and leaves. These results provide abundant information regarding SmSABATH and are useful in further studying SABATH gene functions in S. miltiorrhiza.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/8/12/365/s1. Table S1: List of the SABATH genes from other species. Table S2: Primers for qRT-PCR. Table S3: Gene features of SmSABATH. Table S4: Normal expression sequences of 13 motifs identified in 30 SmSABATH proteins. Table S5: The SABATH protein motif diagram of S. miltiorrhiza. Table S6: Ka/Ks and divergence analysis of SABATH paralogous in S. miltiorrhiza. Table S7: The coefficient of Type-I functional divergence (θ I ) from pairwise comparisons between SmSABATH groups. Table S8: The coefficient of Type-II functional divergence (θ II ) from pairwise comparisons between SmSABATH groups. Figure S1. Phylogenetic tree was reconstructed using the Bayesian inference method under the JTT + I + G model with the SmSABATH amino acid sequences. Figure S2. The neighbor-joining phylogenetic tree was reconstructed with SmSABATH amino acid sequences using MEGA 6.0.