The Pentatricopeptide Repeat Gene Family in Salvia miltiorrhiza: Genome-Wide Characterization and Expression Analysis

The pentatricopeptide repeat (PPR) gene family is one of the largest gene families in plants and plays important roles in posttranscriptional regulation. In this study, we combined whole genome sequencing and transcriptomes to systematically investigate PPRs in Salvia miltiorrhiza, which is a well-known material of traditional Chinese medicine and an emerging model system for medicinal plant studies. Among 562 identified SmPPRs, 299 belong to the P subfamily while the others belong to the PLS subfamily. The majority of SmPPRs have only one exon and are localized in the mitochondrion or chloroplast. As many as 546 SmPPRs were expressed in at least one tissue and exhibited differential expression patterns, which indicates they likely play a variety of functions in S. miltiorrhiza. Up to 349 SmPPRs were salicylic acid-responsive and 183 SmPPRs were yeast extract and Ag+-responsive, which indicates these genes might be involved in S. miltiorrhiza defense stresses and secondary metabolism. Furthermore, 23 salicylic acid-responsive SmPPRs were co-expressed with phenolic acid biosynthetic enzyme genes only while 16 yeast extract and Ag+-responsive SmPPRs were co-expressed with tanshinone biosynthetic enzyme genes only. Two SmPPRs were co-expressed with both phenolic acid and tanshinone biosynthetic enzyme genes. The results provide a useful platform for further investigating the roles of PPRs in S. miltiorrhiza.

In this study, 562 SmPPRs were genome-wide identified and characterized in S. miltiorrhiza. Different expression patterns were investigated. The majority of elicitor-responsive SmPPRs and subsets of SmPPRs co-expressed with phenolic acid and tanshinone biosynthetic genes were shown. Our results provide some references to further functional studies of this family gene in S. miltiorrhiza.

Genome-Wide Identification and Bioinformatic Analysis of 562 SmPPRs
Through BLAST analysis of Arabidopsis thaliana PPRs against the current assembly of the S. miltiorrhiza genome and subsequent gene prediction of the retrieved genomic DNA sequences, a total of 562 SmPPR genes were identified. They were named SmPPR1-SmPPR562, respectively (Supplementary Table S1). These sequence data have been submitted to the GenBank databases under an accession number MH004461-MH005022. The number of PPR genes in S. miltiorrhiza is comparable with that in other plants such as Arabidopsis, rice, maize, and foxtail millet [16][17][18][19]. It indicates that the identified SmPPRs represent an almost complete set of PPRs in S. miltiorrhiza even though it may be not a fully complete set. The number of introns in SmPPR coding regions varied from 0 to 9 with about 88% (495/562) containing no intron, about 10% (57/562) containing one intron, and only 2% (10/562) more than one (See Supplementary Table S1). It was accorded with the previous studies showing that a noticeable feature of PPR genes is that most of them do not contain introns within the coding sequence [16][17][18][19][20][21]. The deduced SmPPR proteins have amino acid numbers widely ranging from 243 to 1486, theoretical pI values ranging between 4.47 and 9.69, and predicted molecular weights ranging from 29.98 kDa to 168.15 kDa (see Supplementary Table S1), which suggests the divergence of SmPPRs.
Generally, PPR proteins have 2-26 tandem arrays of a degenerate 35 amino acid repeat motif and are split into two major subfamilies including P and PLS based on the feature of motifs. P subfamily proteins have the typical P motif while PLS subfamily proteins contain the P motif and P motif-derived variants (the short (S) and the long (L) motifs). Furthermore, based on the C-terminal motifs, the PLS subfamily could be further divided into the P-L-S, E, E + , and DYW subgroups [16]. Motif analysis using the HMMER3.0 package showed that 299 of the 562 identified SmPPRs belonged to the P subfamily while the other 263 were members of the PLS subfamily. Among the 263 PLS SmPPRs, 11 were members of the P-L-S subgroup, 107 belonged to the E subgroup, 53 were included in the E + subgroup, and the other 92 were members of the DYW subgroup (see Figure 2). We failed to construct a phylogenetic tree via the neighbor-joining method, which was consistent with previous studies that found an attempt to relate PPR genes with a phylogenetic tree was not meaningful since the PLS has multiple motifs that are homologous within and among them and the remaining sequences are too divergent to be aligned correctly [20].
The subcellular localization prediction showed that SmPPRs were widely located in chloroplast, mitochondrion, secretory pathway, and other locations (see Supplementary Table S1) with the majority (410/562) to be localized in the mitochondrion or chloroplast.

Expression Patterns of SmPPRs in Different Tissues
It has been shown that many PPR proteins are involved in plant developmental processes [23][24][25][26][27][28][29][30][31][32][33]. In order to preliminarily know the role of SmPPRs in S. miltiorrhiza growth and development, the expression of SmPPRs in roots, stems, leaves, and flowers of S. miltiorrhiza was analyzed transcriptome-wide of which 546 were expressed (see Supplementary Table S2). Among the 546 expressed SmPPRs, 217 exhibited tissue-specific expressions (see Figure 3a). These include 142 expressed mainly in leaves, 29 in flowers, 27 in stems, and 19 in roots. Among these tissue-specific expression SmPPRs, six belong to PLS subgroup, 16 belong to E + subgroup, 44 belong to DYW subgroup, 46 belong to E subgroup, and the remaining 105 belong to P subgroup (see Supplementary  Table S1). These genes probably play tissue-specific roles. The other 329 were highly expressed in at least two tissues (see Figure 3b) of which five belong to the PLS subgroup, 37 belong to the E + subgroup, 47 belong to the DYW subgroup, 60 belong to the E subgroup, and the remaining 180 belong to the P subgroup (see Supplementary Table S1). It indicates that they have more ubiquitous roles in S. miltiorrhiza. A total of 16 SmPPRs had RPKM (reads per kilobase of exon model per million mapped reads) values less than one in all tissues analyzed. They are pseudogenes or expressed only at specific developmental stages or under special conditions [19].

Expression Patterns of SmPPRs in Different Tissues
It has been shown that many PPR proteins are involved in plant developmental processes [23][24][25][26][27][28][29][30][31][32][33]. In order to preliminarily know the role of SmPPRs in S. miltiorrhiza growth and development, the expression of SmPPRs in roots, stems, leaves, and flowers of S. miltiorrhiza was analyzed transcriptome-wide of which 546 were expressed (see Supplementary Table S2). Among the 546 expressed SmPPRs, 217 exhibited tissue-specific expressions (see Figure 3a). These include 142 expressed mainly in leaves, 29 in flowers, 27 in stems, and 19 in roots. Among these tissue-specific expression SmPPRs, six belong to PLS subgroup, 16 belong to E + subgroup, 44 belong to DYW subgroup, 46 belong to E subgroup, and the remaining 105 belong to P subgroup (see Supplementary  Table S1). These genes probably play tissue-specific roles. The other 329 were highly expressed in at least two tissues (see Figure 3b) of which five belong to the PLS subgroup, 37 belong to the E + subgroup, 47 belong to the DYW subgroup, 60 belong to the E subgroup, and the remaining 180 belong to the P subgroup (see Supplementary Table S1). It indicates that they have more ubiquitous roles in S. miltiorrhiza. A total of 16 SmPPRs had RPKM (reads per kilobase of exon model per million mapped reads) values less than one in all tissues analyzed. They are pseudogenes or expressed only at specific developmental stages or under special conditions [19].

Expression Patterns of SmPPRs in Different Tissues
It has been shown that many PPR proteins are involved in plant developmental processes [23][24][25][26][27][28][29][30][31][32][33]. In order to preliminarily know the role of SmPPRs in S. miltiorrhiza growth and development, the expression of SmPPRs in roots, stems, leaves, and flowers of S. miltiorrhiza was analyzed transcriptome-wide of which 546 were expressed (see Supplementary Table S2). Among the 546 expressed SmPPRs, 217 exhibited tissue-specific expressions (see Figure 3a). These include 142 expressed mainly in leaves, 29 in flowers, 27 in stems, and 19 in roots. Among these tissue-specific expression SmPPRs, six belong to PLS subgroup, 16 belong to E + subgroup, 44 belong to DYW subgroup, 46 belong to E subgroup, and the remaining 105 belong to P subgroup (see Supplementary  Table S1). These genes probably play tissue-specific roles. The other 329 were highly expressed in at least two tissues (see Figure 3b) of which five belong to the PLS subgroup, 37 belong to the E + subgroup, 47 belong to the DYW subgroup, 60 belong to the E subgroup, and the remaining 180 belong to the P subgroup (see Supplementary Table S1). It indicates that they have more ubiquitous roles in S. miltiorrhiza. A total of 16 SmPPRs had RPKM (reads per kilobase of exon model per million mapped reads) values less than one in all tissues analyzed. They are pseudogenes or expressed only at specific developmental stages or under special conditions [19].  To confirm the results from RNA-seq data, eight genes were selected for qRT-PCR analysis. It includes SmPPR195 and SmPPR554 mainly expressed in roots, SmPPR158 and SmPPR401 mainly expressed in stems, SmPPR312 and SmPPR420 mainly expressed in leaves, and SmPPR271 and SmPPR278 mainly expressed in flowers (see Figure 4). The results showed that these genes had similar trends between the results from qRT-PCR analysis and RNA-seq data, which validates the RNA-sequence results.
Molecules 2018, 23, x FOR PEER REVIEW 5 of 12 To confirm the results from RNA-seq data, eight genes were selected for qRT-PCR analysis. It includes SmPPR195 and SmPPR554 mainly expressed in roots, SmPPR158 and SmPPR401 mainly expressed in stems, SmPPR312 and SmPPR420 mainly expressed in leaves, and SmPPR271 and SmPPR278 mainly expressed in flowers (see Figure 4). The results showed that these genes had similar trends between the results from qRT-PCR analysis and RNA-seq data, which validates the RNA-sequence results.

Expression of SmPPRs in Response to Salicylic Acid, Yeast Extract, and Ag + Treatments
PPR genes play significant roles in plant response to stress [34][35][36][37]. However, no information is available for SmPPRs. In this study, to investigate the roles of SmPPRs, we mapped RNA-sequence data of S. miltiorrhiza suspension cells treated with or without salicylic acid for 0 hours, 2 hours, and 8 hours to SmPPRs using the SOAP 2.0 software. A total of 558 SmPPRs were found to be expressed (see Supplementary Table S3). Compared with the levels in non-treated control, 349 were differentially expressed at least at one time point (see Supplementary Table S3) of which 186 belong to the P subfamily and the other 163 belong to the PLS subfamily (see Supplementary Table S1). Among them, 62 were up-regulated and 223 were down-regulated with the regulation significant in at least a time-point (see Figure 5a,b). The other 64 were up-regulated at a time-point and downregulated at the other (see Figure 5c). The results suggest that over 62% (349/562) of SmPPRs are salicylic acid-responsive of which the majority were down-regulated after treatment.

Expression of SmPPRs in Response to Salicylic Acid, Yeast Extract, and Ag + Treatments
PPR genes play significant roles in plant response to stress [34][35][36][37]. However, no information is available for SmPPRs. In this study, to investigate the roles of SmPPRs, we mapped RNA-sequence data of S. miltiorrhiza suspension cells treated with or without salicylic acid for 0 h, 2 h, and 8 h to SmPPRs using the SOAP 2.0 software. A total of 558 SmPPRs were found to be expressed (see Supplementary Table S3). Compared with the levels in non-treated control, 349 were differentially expressed at least at one time point (see Supplementary Table S3) of which 186 belong to the P subfamily and the other 163 belong to the PLS subfamily (see Supplementary Table S1). Among them, 62 were up-regulated and 223 were down-regulated with the regulation significant in at least a time-point (see Figure 5a,b). The other 64 were up-regulated at a time-point and down-regulated at the other (see Figure 5c). The results suggest that over 62% (349/562) of SmPPRs are salicylic acid-responsive of which the majority were down-regulated after treatment. To confirm the results from RNA-seq data, eight genes were selected for qRT-PCR analysis. It includes SmPPR195 and SmPPR554 mainly expressed in roots, SmPPR158 and SmPPR401 mainly expressed in stems, SmPPR312 and SmPPR420 mainly expressed in leaves, and SmPPR271 and SmPPR278 mainly expressed in flowers (see Figure 4). The results showed that these genes had similar trends between the results from qRT-PCR analysis and RNA-seq data, which validates the RNA-sequence results.

Expression of SmPPRs in Response to Salicylic Acid, Yeast Extract, and Ag + Treatments
PPR genes play significant roles in plant response to stress [34][35][36][37]. However, no information is available for SmPPRs. In this study, to investigate the roles of SmPPRs, we mapped RNA-sequence data of S. miltiorrhiza suspension cells treated with or without salicylic acid for 0 hours, 2 hours, and 8 hours to SmPPRs using the SOAP 2.0 software. A total of 558 SmPPRs were found to be expressed (see Supplementary Table S3). Compared with the levels in non-treated control, 349 were differentially expressed at least at one time point (see Supplementary Table S3) of which 186 belong to the P subfamily and the other 163 belong to the PLS subfamily (see Supplementary Table S1). Among them, 62 were up-regulated and 223 were down-regulated with the regulation significant in at least a time-point (see Figure 5a,b). The other 64 were up-regulated at a time-point and downregulated at the other (see Figure 5c). The results suggest that over 62% (349/562) of SmPPRs are salicylic acid-responsive of which the majority were down-regulated after treatment.  Similarly, we mapped RNA-sequence data of S. miltiorrhiza hairy roots treated with or without yeast extract and Ag + to SmPPRs for 0 h, 12 h, 24 h, and 36 h using the SOAP 2.0 software. A total of 241 SmPPRs had RPKM value greater than one in at least a time-point (see Supplementary Table S4). Compared to the levels in non-treated control, 183 SmPPRs were differentially expressed at least at one time point (see Supplementary Table S4) of which 131 belong to the P subfamily and the other 52 belong to the PLS subfamily (see Supplementary Table S1). Among them, 27 were up-regulated with the regulation to be statistically significant at least at one-time point, 83 were down-regulated, and 73 were fluctuated (see Figure 6a-c). The results suggest that over 32% (183/562) of SmPPRs are developed with yeast extract and Ag + -responsive. Similarly, we mapped RNA-sequence data of S. miltiorrhiza hairy roots treated with or without yeast extract and Ag + to SmPPRs for 0 hours, 12 hours, 24 hours, and 36 hours using the SOAP 2.0 software. A total of 241 SmPPRs had RPKM value greater than one in at least a time-point (see Supplementary Table S4). Compared to the levels in non-treated control, 183 SmPPRs were differentially expressed at least at one time point (see Supplementary Table S4) of which 131 belong to the P subfamily and the other 52 belong to the PLS subfamily (see Supplementary Table S1). Among them, 27 were up-regulated with the regulation to be statistically significant at least at onetime point, 83 were down-regulated, and 73 were fluctuated (see Figure 6a-c). The results suggest that over 32% (183/562) of SmPPRs are developed with yeast extract and Ag + -responsive.  Tables S3 and S4). In poplar, on the basis of genome-wide transcriptomic analysis, 154 of the PtrPPR genes were induced by biotic and abiotic treatments of which 11 were chosen for verification by qRT-PCR [21]. This suggested that transcriptomic analysis was feasible. In addition, 14 of the 16 SmPPRs with RPKM less than one in roots, stems, leaves, and flowers of normal plants were expressed in suspension cells treated with salicylic acid and/or hairy roots treated with yeast extract and Ag + . It suggests that these genes are only expressed at specific developmental stages or under special conditions.

Co-Expression Pattern Analysis of SmPPRs with Phenolic Acid or Tanshinone Biosynthetic Genes
Phenolic acids are a class of bioactive compounds in S. miltiorrhiza. Previous studies have shown that salicylic acid affects the expression of genes involved in phenolic acid biosynthesis and leads to the accumulation of these compounds in S. miltiorrhiza [7,9]. In this study, we found that more than 62% of the identified 562 SmPPRs responded to salicylic acid treatment. In order to gain insight into the relationship between SmPPRs and phenolic acid biosynthetic genes, we investigated the coexpression patterns of 349 salicylic acid-responsive SmPPRs and phenolic acid biosynthetic genes including SmPAL1, SmC4H1, Sm4CL1, SmTAT1, SmHPPR1, SmRAS1, and SmCYP98A14 [13]. As a result, 25 co-expressed SmPPRs were identified (see Figure 7a).
Tanshinones are the other class of bioactive substances in S. miltiorrhiza. Many tanshinone biosynthetic enzyme genes responded to yeast extract and Ag + treatment [6,44]. In this study, we identified 183 yeast extract and Ag + -responsive SmPPRs. The relationship between SmPPRs and tanshinone biosynthetic genes were investigated. It results in the identification of 18 co-expressed SmPPRs [13] (see Figure 7b).  Tables S3 and S4). In poplar, on the basis of genome-wide transcriptomic analysis, 154 of the PtrPPR genes were induced by biotic and abiotic treatments of which 11 were chosen for verification by qRT-PCR [21]. This suggested that transcriptomic analysis was feasible. In addition, 14 of the 16 SmPPRs with RPKM less than one in roots, stems, leaves, and flowers of normal plants were expressed in suspension cells treated with salicylic acid and/or hairy roots treated with yeast extract and Ag + . It suggests that these genes are only expressed at specific developmental stages or under special conditions.

Co-Expression Pattern Analysis of SmPPRs with Phenolic Acid or Tanshinone Biosynthetic Genes
Phenolic acids are a class of bioactive compounds in S. miltiorrhiza. Previous studies have shown that salicylic acid affects the expression of genes involved in phenolic acid biosynthesis and leads to the accumulation of these compounds in S. miltiorrhiza [7,9]. In this study, we found that more than 62% of the identified 562 SmPPRs responded to salicylic acid treatment. In order to gain insight into the relationship between SmPPRs and phenolic acid biosynthetic genes, we investigated the co-expression patterns of 349 salicylic acid-responsive SmPPRs and phenolic acid biosynthetic genes including SmPAL1, SmC4H1, Sm4CL1, SmTAT1, SmHPPR1, SmRAS1, and SmCYP98A14 [13]. As a result, 25 co-expressed SmPPRs were identified (see Figure 7a).
Tanshinones are the other class of bioactive substances in S. miltiorrhiza. Many tanshinone biosynthetic enzyme genes responded to yeast extract and Ag + treatment [6,44]. In this study, we identified 183 yeast extract and Ag + -responsive SmPPRs. The relationship between SmPPRs and tanshinone biosynthetic genes were investigated. It results in the identification of 18 co-expressed SmPPRs [13] (see Figure 7b). Co-expression pattern analysis showed that 23 SmPPRs were co-expressed with phenolic acid biosynthetic enzyme genes only, 16 SmPPRs were co-expressed with tanshinone biosynthetic enzyme genes only, and two SmPPRs were co-expressed with both phenolic acid and tanshinone biosynthetic enzyme genes (see Figure 7).

Discussion
In this study, 562 SmPPRs were identified genome-widely in S. miltiorrhiza of which 299 members belong to the P subfamily and the others were members of the PLS subfamily. Among the 263 members of the PLS subfamily, 11 were members of the P-L-S subgroup, 107 belonged to the E subgroup, 53 were included in the E + subgroup, and the other 92 were members of the DYW subgroup (see Supplementary Table S1). However, in Arabidopsis, 251 members are from the P subfamily and 199 members are from the PLS subfamily including six members of the P-L-S subgroup, 47 members of the E subgroup, 59 members of the E + subgroup, and 87 members of the DYW subgroup [16]. Comparative analysis of PPR numbers in each subfamily showed that S. miltiorrhiza had more P, P-L-S, and E PPRs than Arabidopsis and the number of P-L-S SmPPRs was about twice that of P-L-S AtPPRs and the number of E SmPPRs was more than twice that of P-L-S AtPPRs, which suggests the member of these subfamilies expanded in S. miltiorrhiza.
It has been reported that some PPRs play important roles in plant developmental processes including cytoplasmic male sterility and fertility restoration, embryogenesis, and seed development [45]. In this study, the tissue-specific expression patterns of SmPPRs in roots, stems, leaves, and flowers in S. miltiorrhiza were investigated based on transcriptome data, which shows that the majority of SmPPRs were expressed ubiquitously with the exception of a few genes expressed in specific tissues (see Figure 3, Supplementary Table S2). It is in line with previous studies in maize [19], which suggests that SmPPRs are multifunctional and are involved in a wide range of biological processes.
More research has shown that PPRs are involved in plant responses to stresses and secondary metabolism [34][35][36][37]42,43]. However, the function of many PPRs is still unknown and there is a lack of knowledge about SmPPRs. Previous studies suggested that salicylic acid could inhibit the activity of plasma membrane H + -ATPase, could make cell generate oxidative stress, could promote the synthesis of phenolic acids, and up-regulate the expression of PAL, TAT, and RAS in S. miltiorrhiza cell culture [7,9]. Ag + could trigger a burst of reactive oxygen species (ROS) and yeast extract could improve the antioxidant defense of S. miltiorrhiza hairy roots [6,11]. The combined use of yeast extract and Ag + was more effective than single yeast extract or Ag + in the hairy root cultures of S. miltiorrhiza Co-expression pattern analysis showed that 23 SmPPRs were co-expressed with phenolic acid biosynthetic enzyme genes only, 16 SmPPRs were co-expressed with tanshinone biosynthetic enzyme genes only, and two SmPPRs were co-expressed with both phenolic acid and tanshinone biosynthetic enzyme genes (see Figure 7).

Discussion
In this study, 562 SmPPRs were identified genome-widely in S. miltiorrhiza of which 299 members belong to the P subfamily and the others were members of the PLS subfamily. Among the 263 members of the PLS subfamily, 11 were members of the P-L-S subgroup, 107 belonged to the E subgroup, 53 were included in the E + subgroup, and the other 92 were members of the DYW subgroup (see Figure 2, Supplementary Table S1). However, in Arabidopsis, 251 members are from the P subfamily and 199 members are from the PLS subfamily including six members of the P-L-S subgroup, 47 members of the E subgroup, 59 members of the E + subgroup, and 87 members of the DYW subgroup [16]. Comparative analysis of PPR numbers in each subfamily showed that S. miltiorrhiza had more P, P-L-S, and E PPRs than Arabidopsis and the number of P-L-S SmPPRs was about twice that of P-L-S AtPPRs and the number of E SmPPRs was more than twice that of P-L-S AtPPRs, which suggests the member of these subfamilies expanded in S. miltiorrhiza.
It has been reported that some PPRs play important roles in plant developmental processes including cytoplasmic male sterility and fertility restoration, embryogenesis, and seed development [45]. In this study, the tissue-specific expression patterns of SmPPRs in roots, stems, leaves, and flowers in S. miltiorrhiza were investigated based on transcriptome data, which shows that the majority of SmPPRs were expressed ubiquitously with the exception of a few genes expressed in specific tissues (see Figure 3, Supplementary Table S2). It is in line with previous studies in maize [19], which suggests that SmPPRs are multifunctional and are involved in a wide range of biological processes.
More research has shown that PPRs are involved in plant responses to stresses and secondary metabolism [34][35][36][37]42,43]. However, the function of many PPRs is still unknown and there is a lack of knowledge about SmPPRs. Previous studies suggested that salicylic acid could inhibit the activity of plasma membrane H + -ATPase, could make cell generate oxidative stress, could promote the synthesis of phenolic acids, and up-regulate the expression of PAL, TAT, and RAS in S. miltiorrhiza cell culture [7,9]. Ag + could trigger a burst of reactive oxygen species (ROS) and yeast extract could improve the antioxidant defense of S. miltiorrhiza hairy roots [6,11]. The combined use of yeast extract and Ag + was more effective than single yeast extract or Ag + in the hairy root cultures of S. miltiorrhiza [44,45]. They could promote the accumulation of tanshinones in S. miltiorrhiza and stimulate the activities of HMGR and DXS enzymes [6]. These results suggest that salicylic acid, yeast extract, and Ag + are effective elicitors for bioactive compound production and related gene expression in S. miltiorrhiza. In this study, 349 SmPPRs were responsive to salicylic acid treatment and 183 SmPPRs were responsive to yeast extract and Ag + treatment, which indicates the importance of SmPPRs in S. miltiorrhiza stress responses and secondary metabolites. Co-expression analysis has previously been described to be a suitable tool for gene function identification [46]. Various genes such as Arabidopsis starch metabolism-related genes [47], rice cell-wall formation-related genes [48], and potato starch biosynthesis-associated transcription factor genes [49], were identified based on co-expression analysis. Gene co-expression analysis showed that 25 salicylic acid-responsive SmPPRs were co-expressed with phenolic acid biosynthetic genes and 18 yeast extract and Ag + -responsive SmPPRs were co-expressed with tanshinone biosynthetic genes (see Figure 7). It includes two co-expressed with both phenolic acid and tanshinone biosynthetic enzyme genes, which brings the total number of co-expressed SmPPRs to 41. These SmPPRs could be involved in phenolic acid and/or tanshinone biosynthesis in S. miltiorrhiza. Subcellular localization prediction showed these SmPPR proteins mainly located in mitochondrion and chloroplast (see Supplementary Table S1). They could be related with the respiratory cytochrome pathway or the photosynthesis of chloroplasts to regulate the biosynthesis of phenolic acid and tanshinone [42,43]. They also could be through the mitochondrial electron transport system or hormonal responses to regulate plant environmental stresses -responsive [35][36][37]. Further physiological and biochemical analysis will be performed to verify the functions of these SmPPRs in the future.

Plant Materials
Salvia miltiorrhiza Bunge (line 993) with whole genome sequence available was grown in a field nursery at the Institute of Medicinal Plant Development (Beijing, China). Roots, stems, leaves, and flowers were collected from two-year-old plants and stored in liquid nitrogen until use.

Genome-Wide Identification of SmPPRs
Arabidopsis PPR (AtPPR) protein sequences were downloaded from TAIR (http://www. arabidopsis.org/). BLAST analysis of AtPPR proteins against the genome database of S. miltiorrhiza (line 993) [50] was carried out using tBLASTn. An e-value cut-off of e −10 was applied. Gene models of SmPPRs were predicted based on the alignments between the retrieved DNA sequences and PPR proteins from other plant species via BLASTx. The predicted gene models were examined and comparatively analyzed with the genome database of the other S. miltiorrhiza line (http://www.herbalgenome.cn). The domain information of SmPPR proteins was analyzed using the software package HMMER 3.0 (http://hmmer.janelia.org/software/archive). The HMMbuild program in the HMMER 3.0 package was used for producing the matrices specific to each type of SmPPRs. Motifs were searched using the HMMsearch program in the HMMER 3.0 package. The gene structures were analyzed using the Gene Structure Display Server (GSDS 2.0, http://gsds.cbi.pku.edu.cn/index.php). Subcellular localization of SmPPR proteins was predicted with TargetP version 1.1 (http://www.cbs.dtu.dk/ services/TargetP/). The molecular weight (MW) and theoretical isoelectric point (pI) were predicted via the compute pI/MW tool on the ExPASy server (http://web.expasy.org/compute_pi/).

Expression Analysis of SmPPRs
Expression patterns of SmPPR genes in different organs were analyzed using the transcriptome datasets of roots, stems, leaves, and flowers of S. miltiorrhiza downloaded from GenBank (SRP051564, SRP028388). The changes of SmPPR expression levels in response to stress treatments were analyzed using the transcriptome datasets of suspension cells treated with salicylic acid (SA, 0.16 mM) (SRX1423774) for 0 h, 2 h, and 8 h and hairy roots treated with yeast extract (100 µg/mL) and Ag + (30 µM) (SRR924662) for 0 h, 12 h, 24 h, and 36 h. RNA-sequence reads were mapped to SmPPRs using SOAP 2.0 [51] and analyzed as described previously [12]. Co-expression pattern analysis of SmPPR genes with phenolic acid and tanshinone biosynthetic genes in S. miltiorrhiza was analyzed using the software package of Co-expression Pattern Clustering Analysis in the BMKCloud cloud server (http://www.biocloud.net/). The abundance of transcripts was estimated by RPKM. Genes with RPKM value greater than one in at least a tissue was considered to be expressed. To draw heat maps, the value of RPKM for each gene was normalized between −1.0 and 1.0 using the R package version 3.4.1 and analyzed for differential expression by the Fisher's exact test. p < 0.05 was considered as differentially expressed.

Quantitative Real-Time Reverse Transcription-PCR (qRT-PCR)
Gene expression was analyzed using qRT-PCR. Total RNA was extracted from prepared samples using the Plant RNA EASYspin Plus (Aidlab biotech, Beijing, China). RNA integrity was analyzed on a 1.2% agarose gel and its quantity was determined using a NanoDrop 2000C Spectrophotometer (Thermo Scientific, Waltham, MA, USA). Reverse transcription was carried out using the PrimeScript™ RT reagent kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China). PCR was performed using the SYBR ® Premix Ex Taq™ II (Tli RNaseH Plus, Shiga, Japan) (TaKaRa) as described previously [12]. Gene-specific primers for SmPPRs were listed in Supplementary Table S5. The length of amplicons was between 80 bp and 250 bp. SmUBQ10 was used as an internal control. qRT-PCRs were performed in three biological replicates with three technical replicates each. The expression level in stems was set to one and the levels in other tissues were given relative to this. The relative expression levels of genes were calculated by the 2 −∆∆Ct method. ANOVA (analysis of variance) was calculated using SPSS (Version 19.0, IBM, Chicago, IL, USA). p < 0.05 was considered statistically significant.

Conclusions
In this study, bioinformatic analysis for PPR gene family in S. miltiorrhiza was performed. A total of 562 PPR protein genes were identified in S. miltiorrhiza genome. Gene structure and classification showed important features for this family. Tissue-specific expression analysis indicated that SmPPRs have more ubiquitous roles in S. miltiorrhiza. The expression pattern in response to different elicitors revealed SmPPRs might be involved in stress defense and secondary metabolites. Furthermore, some SmPPRs co-expressed with phenolic acid and tanshinone biosynthetic genes were shown. The results provide valuable information for future studies on characterizing the biological functions of PPR protein genes in S. miltiorrhiza.

Supplementary Materials:
The following are available online. Table S1: Sequence features of the PPR gene family members in S. miltiorrhiza. Table S2: The expression levels of SmPPR genes in different tissues (RPKM). Table S3: RPKM value for SmPPR expression in S. miltiorrhiza suspension cell treated with or without salicylic acid. Table S4: RPKM value for SmPPR expression in hairy roots treated with or without yeast extract and Ag + . Table S5: Primers used for qRT-PCR.