Transcriptomic Proﬁling of Pomegranate Provides Insights into Salt Tolerance

: Pomegranate ( Punica granatum L.) is widely grown in arid and semi-arid soils, with constant soil salinization. To elucidate its molecular responses to salt stress on mRNA levels, we constructed 18 cDNA libraries of pomegranate roots and leaves from 0 (controls), 3, and 6 days after 200 mM NaCl treatment. In total, we obtained 34,047 genes by mapping to genome, and then identiﬁed 2255 DEGs (di ﬀ erentially expressed genes), including 1080 up-regulated and 1175 down-regulated genes. We found that the expression pattern of most DEGs were tissue-speciﬁc and time-speciﬁc. Among root DEGs, genes associated with cell wall organization and transmembrane transport were suppressed, and most of metabolism-related genes were over-represented. In leaves, 41.29% of DEGs were ﬁrst suppressed and then recovered, including ions / metal ions binding-related genes. Also, ion transport and oxidation-reduction process were restricted. We found many DEGs involved in ABA, Ca 2 + -related and MAPK signal transduction pathways, such as ABA-receptors, Ca 2 + -sensors, MAPK cascades, TFs, and downstream functional genes coding for HSPs, LEAs, AQPs and PODs. Fifteen genes were selected to conﬁrm the RNA-seq data using qRT-PCR. Our study not only illuminated pomegranate molecular responses to salinity, but also provided references for selecting salt-tolerant genes in pomegranate breeding processes. methodology, C.L. and Y.Z.; formal analysis, C.L., Y.Z., J.W., and X.Z.; investigation, C.L., Y.Z., X.Z., and J.W.; writing—original draft preparation, C.L.; writing—review and editing, C.L., X.Z., M.G., and Z.Y.; supervision, M.G. and Z.Y.; Funding acquisition, Z.Y.


Sequence Assembly
Low quality reads and reads containing adapter and ploy-N were filtered by NGS QC ToolKit [26] from raw reads. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. The clean data from all 18 libraries were separately mapped to the pomegranate genome assembly (ASM220158v1) using HISAT2 software [27]. StringTie [28] was used to construct and identify both known and novel transcripts from HISAT2 alignment results. StringTie was also used to count the reads numbers mapped to each gene. After that, FPKM (Fragments Per Kilobase of exon per million fragments Mapped) were calculated based on the length of the gene and reads count mapped to this gene [27].

Identification and Functional Annotation of DEGs
All genes were compared against various protein databases by BLASTX, including the Nr (NCBI non-redundant protein sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins), and Swiss-Prot (A manually annotated and reviewed protein sequence database) with an E-value cut-off of 10 −5 . Then genes with the best BLAST hit (the highest score) were chosen along with their protein functional annotations. DESeq and Q-value were employed to evaluate the differential expression genes (DEGs) between controls and treatments. The false discovery rate (FDR) and log 2 FC (log of fold change) were calculated for all genes, and only transcripts with FDR < 0.05 and|log 2 (fold change)|≥1 were considered as DEGs. To annotate the DEGs with gene ontology (GO) terms, the Nr BLAST results were imported into the Blast2GO program [29]. To examine the expression Agronomy 2020, 10, 44 4 of 17 patterns of DEGs, the expression data from roots and leaves (T0, T1 and T2) were normalized to 0, log 2 (T1/T0), log 2 (T2/T0), and then clustered by Short Time-series Expression Miner software (STEM) separately, using a FDR correction method and p-value ≤ 0.05 as the cutoff [30]. The results of GO annotations in each pattern were enriched and refined using TBtools v0.6652 (Toolbox for Biologists, https://github.com/CJ-Chen/TBtools). We used KOBAS software [31] to test the enrichment of DEGs in KEGG pathways.

Validation of RNA-Seq by qRT-PCR
We selected 15 DEGs as experimental validation by quantitative real-time PCR (qRT-PCR), which were performed with three biological and three technical replicates for each cDNA sample. The primers for these genes were designed with NCBI primer-BLAST (Table S1). Reverse transcription was conducted with the GoScriptTM Reverse Transcription System (Promega, Madison, WI, USA). We conducted qRT-PCR in a 7500 fast Real-Time PCR system (Applied Biosystems, CA, MA, USA), and analyzed results with the ∆∆Ct method, and Pgactin (F: ATCCTCCGTCTTGACCTTG, R: TGTCCGTCAGGCAACTCAT) gene was used as a reference gene. Each reaction was carried out in a final volume of 20 µL, containing 7.5 µL of ddH 2 O, 10 µL of SYBR Green PCR master mix, 0.5 µL of each gene-specific primer and 2 µL of diluted cDNA. The PCR thermal cycling conditions were as follows: 95 • C for 10 min; 40 cycles of 95 • C for 5 s, 60 • C for 30 s and 72 • C for 30 s. Data were collected during the extension step: 95 • C for 15 s, 60 • C for 1 min, 95 • C for 30 s and 60 • C for 15 s.

Overview of Sequencing and Mapping
To obtain a global overview of the salt-induced changes at whole-transcriptomic scale in pomegranate, we constructed 18 cDNA libraries from the roots and leaves of 0 (controls), 3 and 6 days after 200 mM NaCl treatment. A high correlation (R 2 > 0.80) between biological replicates was observed for all treatments ( Figure S2), which indicated that the biological replicates were reliable in this study. In totally, sequencing libraries yielded 519.47 million reads with 150 bp for both paired ends. After adapter removal and refining, we obtained 155.84 Gb of clean data, and Q30 percentage were all above 90.08% (Table 1). The ratios of reads mapping to the pomegranate genome were high, with values ranging from 93.61% to 95.33%. Then, the unique-mapped reads were 91.89-93.79% and the mapped sequences in genome exon regions were 62.88-69.32% (Table 1). Finally, we assembled 29,226 transcripts from the 18 cDNA libraries of pomegranate roots and leaves under salt stress. To splice the genes more completely and accurately, the StringTie software was employed to reconstruct the transcripts, and then retrieved 34,047 assemblies, including 5396 novel transcripts. Finally, 26,444 genes (including 2421 new genes) were assembled from annotated genes using 7 different databases (see methods), and can be used as a reference of pomegranate genome annotation (Table S2).

Identification and Annotation of DEGs
A total of 2255 DEGs were identified from salt treatments with FDR < 0.05 and|log 2 (fold change)|≥1, including 1080 up-regulated and 1175 down-regulated genes, and more DEGs were found in roots (1623 genes) than in leaves (632 genes) (Table S2). To provide insights into the underlying functions of pomegranate transcripts under salt stress, DEGs were annotated after they were compared to well-studied sequences in the GO, Nr, Swiss-Prot, KEGG, and KOG databases (Table S2). In this study, most of DEGs in T1R were distributed in T subgroup (Signal transduction mechanisms), E subgroup (Amino acid transport and metabolism) and G subgroup (Carbohydrate transport and metabolism) (Figure 1a). DEGs of T2R and T1L samples were mostly concentrated in R subgroup (General function prediction only), Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms), and K subgroup (Transcription) (Figure 1b,c). The DEGs in T2L were distributed in O subgroup (Posttranslational modification, protein turnover and chaperones), P (Inorganic ion transport and metabolisms) and Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms) ( Figure 1d). These subgroups were close to the metabolism, such as DNA transcription, cell division, protein modification and energy conversion.

Identification and Annotation of DEGs
A total of 2255 DEGs were identified from salt treatments with FDR < 0.05 and|log2(fold change)|≥1, including 1080 up-regulated and 1175 down-regulated genes, and more DEGs were found in roots (1623 genes) than in leaves (632 genes) (Table S2). To provide insights into the underlying functions of pomegranate transcripts under salt stress, DEGs were annotated after they were compared to well-studied sequences in the GO, Nr, Swiss-Prot, KEGG, and KOG databases (Table S2). In this study, most of DEGs in T1R were distributed in T subgroup (Signal transduction mechanisms), E subgroup (Amino acid transport and metabolism) and G subgroup (Carbohydrate transport and metabolism) (Figure 1a). DEGs of T2R and T1L samples were mostly concentrated in R subgroup (General function prediction only), Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms), and K subgroup (Transcription) (Figure 1b,c). The DEGs in T2L were distributed in O subgroup (Posttranslational modification, protein turnover and chaperones), P (Inorganic ion transport and metabolisms) and Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms) (Figure 1d). These subgroups were close to the metabolism, such as DNA transcription, cell division, protein modification and energy conversion. The DEGs of roots and leaves were also analyzed by KEGG pathways enrichment. These pathways were mainly classified into 5 categories ( Figure S3) and showed the top 20 pathways (Figure 2). The majority of pathways in pomegranate roots and leaves were associated with The DEGs of roots and leaves were also analyzed by KEGG pathways enrichment. These pathways were mainly classified into 5 categories ( Figure S3) and showed the top 20 pathways (Figure 2). The majority of pathways in pomegranate roots and leaves were associated with metabolisms, such as global and overview maps, amino acid metabolism, carbohydrate metabolism, lipid metabolism, and energy metabolism, etc. ( Figure S3). The metabolism-related genes were significantly enriched in metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), cysteine and methionine metabolism (ko00270), phenylpropanoid biosynthesis (ko00940), starch and sucrose metabolism (ko00500), amino sugar and nucleotide sugar metabolism (ko005200) and protein processing in endoplasmic reticulum (ko04141) pathways. Also, there were some DEGs involved in environmental adaptation (plant-pathogen interaction, ko04626) and signal transduction (plant hormone signal transduction, ko04075) pathways ( Figure 2; Figure S3). Genes clustered in genetic information processing were enriched in folding, sorting and degradation, replication and repair, transcription and translation, including protein processing in endoplasmic reticulum (ko04141), DNA replication (ko03030), RNA polymerase (ko03020), RNA transport (ko03013) and RNA degradation (ko03018). In addition, genes involved in the photosynthetic pathways, such as carotenoid biosynthesis (ko00906), porphyrin and chlorophyll metabolism (ko00860) and photosynthesis-antenna proteins (ko00196) were suppressed by salinity.

The Expressional Patterns of DEGs
Numbers of DEGs increased in roots with the duration of salt-stress, which were opposite in leaves. The genes were both with little overlap in roots and leaves at two points of stress time ( Figure  3a). Only a small portion (2.0%) of DEGs (20 up-regulated and 24 down-regulated) shared the common expression tendency between roots and leaves, and even less (1.2%) DEGs showed an utterly opposite tendency between two tissues (3) genes were up-regulated in leaves and down-regulated in roots, 23 genes were down-regulated in leaves but up-regulated in roots) ( Figure  3b). The remaining majority of DEGs were exclusively up-regulated or down-regulated in either tissue. Almost 85.9% of up-regulated and 53.4% of down-regulated DEGs were activated in roots

The Expressional Patterns of DEGs
Numbers of DEGs increased in roots with the duration of salt-stress, which were opposite in leaves. The genes were both with little overlap in roots and leaves at two points of stress time (Figure 3a). Only a small portion (2.0%) of DEGs (20 up-regulated and 24 down-regulated) shared the common expression tendency between roots and leaves, and even less (1.2%) DEGs showed an utterly opposite tendency between two tissues (3) genes were up-regulated in leaves and down-regulated in roots, 23 genes were down-regulated in leaves but up-regulated in roots) (Figure 3b). The remaining majority of DEGs were exclusively up-regulated or down-regulated in either tissue. Almost 85.9%  Next, DEGs within these profiles were subjected to GO-term enrichment analysis. They were classified into three main categories, including cellular component, biological process, and molecular function, and the top 30 sub-categories (if more than 30, p-value < 0.05) of each category were shown (Figure 4). In roots, genes in up-regulated pattern of Profile 4 were just inducted in T2R. They were enriched in chitinase activity, chitin binding, oxidoreductase activity, metal ion, and cation binding under molecular function; chitin, glucosamine-containing compound, aminoglycan, amino sugar and cell wall macromolecule metabolic process under biological process (Figure 4c; Table S3). Gens were continuously suppressed with the extension of stress in Profile 7, and they were enriched in endopeptidase activity, peptidase activity, hydrolase activity, and catalytic activity under molecular function and proteolysis and metabolic process under biological process ( Figure  4d; Table S3). The top subcategories of down-regulated pattern (Profile 3) were extracellular region (GO:0005576), cell wall (GO:0005618) and external encapsulating structure (GO:0030312) under The STEM was employed to analyze the expressional patterns of DEGs. All the 1,623 DEGs in roots were clustered into 8 profiles, whereby 1,340 DEGs were significantly clustered into 4 profiles with p-value ≤ 2 × 10 −5 (Figure 4a). The two up-regulated patterns were Profile 4 (31.55%, 512 DEGs) and Profile 7 (15.16%, 246 DEGs), and two down-regulated patterns were Profile 3 (28.59%, 464 DEGs) and Profiles 0 (7.27%, 118 DEGs). All the 632 DEGs in leaf samples were also clustered into 8 profiles and then 2 enriched profiles, including one first down-regulated and then the up-regulated pattern of Profile 2 (41.29%, 261 DEGs, p-value = 3 × 10 −32 ), and one down-regulated pattern of Profile 1 (16.61%, 105 DEGs, p-value = 6 × 10 −17 ) (Figure 4b).
oxidoreductase activity, metal ion, and cation binding under molecular function. The top subcategories in down-regulated pattern of Profile 1 were coenzyme binding, cofactor binding and oxidoreductase activity under molecular function, and oxidation-reduction process, cation and ion transport under biological process. We identified genes coding ion transporters, such as cation/calcium exchanger, zinc transporter, ammonium transporter 1 member, copper-transporting ATPase, potassium transporter (Figure 4h; Table S5).  Next, DEGs within these profiles were subjected to GO-term enrichment analysis. They were classified into three main categories, including cellular component, biological process, and molecular function, and the top 30 sub-categories (if more than 30, p-value < 0.05) of each category were shown ( Figure 4). In roots, genes in up-regulated pattern of Profile 4 were just inducted in T2R. They were enriched in chitinase activity, chitin binding, oxidoreductase activity, metal ion, and cation binding under molecular function; chitin, glucosamine-containing compound, aminoglycan, amino sugar and cell wall macromolecule metabolic process under biological process (Figure 4c; Table S3). Gens were continuously suppressed with the extension of stress in Profile 7, and they were enriched in endopeptidase activity, peptidase activity, hydrolase activity, and catalytic activity under molecular function and proteolysis and metabolic process under biological process (Figure 4d; Table S3). The top subcategories of down-regulated pattern (Profile 3) were extracellular region (GO:0005576), cell wall (GO:0005618) and external encapsulating structure (GO:0030312) under cellular component, and lipid biosynthetic process (GO:0008610), cell wall organization (GO:0071555), and external encapsulating structure organization (GO:0045229) under biological process. (Figure 4e; Table S4). Many DEGs coding cell wall components were mostly suppressed in T2R samples, including four extensins, two pectin acetylesterase, three polygalacturonase, and three xyloglucan endotransglucosylase/hydrolase Agronomy 2020, 10, 44 9 of 17 proteins. At the same time, transmembrane transport (GO:0055085) and oxidation-reduction process (GO:0055114) in Profile 0 were inhibited under salt stress. Many genes were identified in these biological process, such as adenine/guanine permease, aquaporin, cationic amino acid transporter, sugar carrier, and zinc transporter (Figure 4f; Table S4). Profile 2 showed a first down-and then up-regulated pattern in leaves, there were many subcategories, such as proteolysis and oxidation-reduction process under biological process, oxidoreductase activity, metal ion, and cation binding under molecular function. The top subcategories in down-regulated pattern of Profile 1 were coenzyme binding, cofactor binding and oxidoreductase activity under molecular function, and oxidation-reduction process, cation and ion transport under biological process. We identified genes coding ion transporters, such as cation/calcium exchanger, zinc transporter, ammonium transporter 1 member, copper-transporting ATPase, potassium transporter ( Figure 4h; Table S5).

ABA Signaling Pathway
The ABA signaling pathway mainly consists of three protein classes: ABA receptors (PYR/PYL/RCAR), type 2C protein phosphatase (PP2Cs) and sucrose non-fermenting1-related protein kinase 2 (SnRK2s). In our study, three ABA receptors, PYLs, were significantly down-regulated ( Figure 5), while 13 PP2Cs have different expressional patterns in pomegranate roots and leaves under salt stress. Among these genes, seven PP2Cs were up-regulated in roots or leaves, five PP2Cs down-regulated, and another gene up-regulated in roots and down-regulated in leaves ( Figure 5). In contrast, just one SnRK2 was identified among the DEGs, and it was up-regulated in leaves, which allowed the accumulation of phosphorylated downstream ABFs and activation of ABA-response genes.

ABA Signaling Pathway
The ABA signaling pathway mainly consists of three protein classes: ABA receptors (PYR/PYL/RCAR), type 2C protein phosphatase (PP2Cs) and sucrose non-fermenting1-related protein kinase 2 (SnRK2s). In our study, three ABA receptors, PYLs, were significantly down-regulated ( Figure 5), while 13 PP2Cs have different expressional patterns in pomegranate roots and leaves under salt stress. Among these genes, seven PP2Cs were up-regulated in roots or leaves, five PP2Cs down-regulated, and another gene up-regulated in roots and down-regulated in leaves ( Figure 5). In contrast, just one SnRK2 was identified among the DEGs, and it was up-regulated in leaves, which allowed the accumulation of phosphorylated downstream ABFs and activation of ABA-response genes.

Ca 2+ -Related Signaling Pathways
Thirty DEGs involved in Ca 2+ -related signaling pathway were identified, these transcripts coded function proteins included three Ca 2+ -ATPases (ACAs), four cation/H + antiporters (CAXs), two cation/calcium exchangers (CCXs), four glutamate receptor (GLRs), ten calcium-binding proteins (CaM/CMLs), two CBL-interacting protein kinases (CIPKs) and five calcium-dependent protein kinases (CDPKs) ( Figure 6). Expectedly, two ACAs were up-regulated in roots (T2R), and one was down-regulated in leaves ( Figure 6). Ten CAXs expressed differently, including two CAXs located on vacuoles were both up-regulated in T2R. Two CIPKs were significantly up-regulated in leaves, but no obvious change was observed in roots. The majority of CaM/CMLs were up-regulated in roots and down-regulated in leaves. The CDPKs were up-or down-regulated in roots and/or leaves. Also, two DEGs, that were involved in the MAPK signaling pathway were identified, and both were down-regulated in leaves ( Figure 6).

Ca 2+ -Related Signaling Pathways
Thirty DEGs involved in Ca 2+ -related signaling pathway were identified, these transcripts coded function proteins included three Ca 2+ -ATPases (ACAs), four cation/H + antiporters (CAXs), two cation/calcium exchangers (CCXs), four glutamate receptor (GLRs), ten calcium-binding proteins (CaM/CMLs), two CBL-interacting protein kinases (CIPKs) and five calcium-dependent protein kinases (CDPKs) ( Figure 6). Expectedly, two ACAs were up-regulated in roots (T2R), and one was down-regulated in leaves ( Figure 6). Ten CAXs expressed differently, including two CAXs located on vacuoles were both up-regulated in T2R. Two CIPKs were significantly up-regulated in leaves, but no obvious change was observed in roots. The majority of CaM/CMLs were up-regulated in roots and down-regulated in leaves. The CDPKs were up-or down-regulated in roots and/or leaves. Also, two DEGs, that were involved in the MAPK signaling pathway were identified, and both were down-regulated in leaves ( Figure 6).

The Transcription Factors (TFs)
For identification of the TFs in pomegranate transcriptome, all the mapped genes were analyzed by BLAST against the Plant Transcription Factor Database (PlantTFDB, http://planttfdb.cbi.pku.edu.cn). A total of 1346 TFs that were classified into 55 putative families (Table  S6). Among these TFs, twenty-seven TF families, including 151 genes expressed differently under salt stress compared to controls. The most abundant of differential expression TFs included NAC, ERF, MYB-related, C2H2, MYB, bHLH, GRAS, WRKY, LBD, B3 and bZIP family genes ( Table 2; Table S6).
Under NaCl stress, we found that 12 of 19 NAC genes were up-regulated, and 10 of 19 genes were down-regulated in leaves or roots of pomegranate plants. Interestingly, we also observed that 18 PgNACs significantly changed in T2R when compared to controls, and 4 genes were up-regulated or down-regulated in roots, but reversed in leaves ( Table 2; Table S6). These results suggested that many NACs were involved in salt stress, but there were different potential responding mechanisms of NAC domain genes to salinity. Thirteen MYB genes were induced by salt treatment in pomegranate ( Table 2). Among these genes, six in roots and four in leaves were up-regulated, and four in roots and eight in leaves were down-regulated. Moreover, sixteen MYB-related genes were detected in our RNA-Seq analyses, thirteen genes were down-regulated, and seven genes were up-regulated by salt stress ( Table 2; Table S6). But so far, little was reported that the MYB-related type proteins are related with the responses to salt stress.
The AP2/ERF superfamily is divided into three families: the AP2 family proteins containing two repeated AP2/ERF domains, the ERF family proteins, containing a single AP2/ERF domain, and the RAV family proteins containing a B3 domain. There were 19 PgERFs, and 3 PgAP2s significantly expressed in treatments when compared to controls ( Table 2). Fifteen ERFs were down-regulated, and all of AP2 were repressed by salinity. There were 9 up-regulated and 8 down-regulated C2H2, three up-regulated and nine down-regulated bHLH, five up-regulated and three down-regulated WRKY identified in pomegranate roots and leaves. Remarkably, nine GRASs were up-regulated in

The Transcription Factors (TFs)
For identification of the TFs in pomegranate transcriptome, all the mapped genes were analyzed by BLAST against the Plant Transcription Factor Database (PlantTFDB, http://planttfdb.cbi.pku.edu.cn). A total of 1346 TFs that were classified into 55 putative families (Table S6). Among these TFs, twenty-seven TF families, including 151 genes expressed differently under salt stress compared to controls. The most abundant of differential expression TFs included NAC, ERF, MYB-related, C2H2, MYB, bHLH, GRAS, WRKY, LBD, B3 and bZIP family genes ( Table 2; Table S6).
Under NaCl stress, we found that 12 of 19 NAC genes were up-regulated, and 10 of 19 genes were down-regulated in leaves or roots of pomegranate plants. Interestingly, we also observed that 18 PgNACs significantly changed in T2R when compared to controls, and 4 genes were up-regulated or down-regulated in roots, but reversed in leaves ( Table 2; Table S6). These results suggested that many NACs were involved in salt stress, but there were different potential responding mechanisms of NAC domain genes to salinity. Thirteen MYB genes were induced by salt treatment in pomegranate ( Table 2). Among these genes, six in roots and four in leaves were up-regulated, and four in roots and eight in leaves were down-regulated. Moreover, sixteen MYB-related genes were detected in our RNA-Seq analyses, thirteen genes were down-regulated, and seven genes were up-regulated by salt stress ( Table 2; Table S6). But so far, little was reported that the MYB-related type proteins are related with the responses to salt stress.  NAC  97  19  1  1  12  6  1  2  1  4  ERF  114  19  1  1  8  5  1  7  9  MYB_related  91  16  --4  3  1  11  3  2  C2H2  91  13  --8  3  -1  2  MYB  79  13  -1  6  4  -8  4  2  bHLH  102  11  -1  2  5  1  4  1  1  GRAS  50  9 1 Others  419  19  --8  3  -6  1  3  Total  1346  151  5  6  70  40  6  46  13  31 The AP2/ERF superfamily is divided into three families: the AP2 family proteins containing two repeated AP2/ERF domains, the ERF family proteins, containing a single AP2/ERF domain, and the RAV family proteins containing a B3 domain. There were 19 PgERFs, and 3 PgAP2s significantly expressed in treatments when compared to controls ( Table 2). Fifteen ERFs were down-regulated, and all of AP2 were repressed by salinity. There were 9 up-regulated and 8 down-regulated C2H2, three up-regulated and nine down-regulated bHLH, five up-regulated and three down-regulated WRKY identified in pomegranate roots and leaves. Remarkably, nine GRASs were up-regulated in the root, but the no-significant change in leaves ( Table 2; Table S6).

qRT-PCR Validation
To confirm the reliability of the expression levels obtained from the RNA-Seq, fifteen DEGs, including five PgERFs, three PgMYBs and seven PgNACs were selected for qRT-PCR assays. The results showed that the expression level of each transcript closely corresponded to the transcript level estimated from the sequence data with R 2 ≥ 0.85 (Figure 7), which implies reproducibility and accuracy of the RNA-Seq results.

Discussion
RNA-sequencing techniques have proven to be beneficial and economical for scanning the transcribed genes, both in model and non-model plant species. In this study, we used the RNA-Seq approach as a powerful tool to elucidate the molecular responses to salt stress in pomegranate.

Discussion
RNA-sequencing techniques have proven to be beneficial and economical for scanning the transcribed genes, both in model and non-model plant species. In this study, we used the RNA-Seq approach as a powerful tool to elucidate the molecular responses to salt stress in pomegranate. Finally, we reconstructed the transcripts and identified 5396 novel genes from the 18 cDNA libraries. There were significant differences between the percentage of total mapped reads and the mapped reads in exonic regions. Theoretically, the ratio of sequencing reads that produced from mature mRNA mapped into exonic regions could be 100%. However, a portion of 30.68-37.12% reads were also mapped into intronic and intergenic regions for the following reasons: (1) The coding and non-coding sequences in the reference genome might have been annotated incorrectly; (2) the new alternative splicing (SNPs) and novel genes resulted from the alternative mRNA splicing and transcripts reconstructing; (3) variations exist between sequencing and reference genome genotypes, which are attributed to the gaps and diversity of their sequences [32,33].
Under salt stress, plants firstly experience osmotic stress due to a disorder of water uptake. Then a process of gradual recovery due to the partially or completely re-established uptake of water occurs within days after the salt treatment [34,35]. We identified 2255 DEGs through the 18 mRNA libraries. From the third to sixth day after treatment, there were significant differences between pomegranate roots and leaves, with little overlap of DEGs responding salinity ( Figure 3). More salt-related genes were up-regulated in roots, while the majority of suppressed genes recovered with salt-treating process in leaves. In contrast, the osmotic adjustment of leaves started only after roots had reached new water equilibrium [36]. These tissue-specific and time-specific differences between roots and leaves have also been reported in Arabidopsis thaliana (L.) [37], Populus euphratica (Oliv.) [36], and Millettia pinnata (L.) [38]. Many DEGs were involved in down-regulated patterns in roots, such as cell wall organization, transmembrane transport and oxidation-reduction process, but the proteolysis and metabolic process were over-expressed. Otherwise, most of DEGs involved in oxidation-reduction process and ion transport in leaves were suppressed (Figure 4). The cell wall plays an important role in protecting plant from salt toxicity [39]. In our study, extensin, pectin acetylesterase, polygalacturonase, and xyloglucan endotransglucosylase/hydrolase proteins were identified in roots (Table S4). These proteins are important components of the cell wall, which affect plant growth by mediating cell enlargement and expansion [39]. Many DEGs involved in transmembrane transport code for various channels, carriers and pumps, and ion transporters, such as aquaporin, cationic amino acid transporter, sugar carrier, and zinc transporter (Table S4). They worked together to transport coenzymes, amino acids, carbohydrates, and ions under salt stress. The suppressed genes coding proteins, such as dehydrogenase, cytochrome P450s and peroxidase involved in oxidation-reduction process, play essential roles in responding to salinity [40,41] (Tables S4 and S5). These results indicated that salinity affected pomegranate growth and development by inhibiting the cell division and transmembrane transport, as well as slowing down the redox reactions.
On the other hand, over-represented genes in roots were enriched in metabolic process, carbohydrate metabolic process, and catabolic process, etc. These metabolism-related genes in roots were also enriched in KEGG pathways, such as metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), phenylpropanoid biosynthesis (ko00940), starch and sucrose metabolism (ko00500), amino sugar and nucleotide sugar metabolism (ko005200) (Figure 2; Table S3). The carbohydrates such as glycan, starch, sucrose, amino sugar have been demonstrated as osmolytes and energy resouces in plant responses to salt stress, especially in halophyte species [42]. The accelerated metabolisms in pomegranate might contribute to its adaption responses to salt stress. Similar expression patterns of metabolism-related genes were observed in Oryza sativa (L.) [43], Thellungiella halophila (C. A. Mey.) [42] and Helianthus tuberosus (L.) [44] under salt stress. Metal ions/cations, such as K + , Ca 2+ , Mg 2+ , Fe 2+ , Cu 2+ , and Zn 2+ , act as cofactors participating in the catalytic activity of the enzymes they bind to [44]. The up-regulated DEGs involved in metal ion/cation binding might accelerate the catalytic activity. At the same time, the ion transport process was inhibited in roots and leaves (Tables  S3 and S5). The restriction of ion transport indicated that the excess of Na + inhibited the uptake of mineral ions [4]. The results correspond to the physiological responses in other pomegranate cultivars under salt stress [45]. We suspected that pomegranate could cope with the toxic ions to some extent via decreasing uptake of ions from soil and increasing utilization of ions in cells. The toxic ions access into the leaf cells with the transpiration, but the detrimental effect is a slow process taking weeks or months, eventually cause salt toxicity in leaves [46]. Plants were only exposed to salinity stress just for six days in our research, so further study is needed to reveal the long-term effect of salinity.
ABA is a critical hormone, which regulates plant growth, development, as well as responses to environmental stresses such as salinity, drought, heat, cold, wound, and pathogen [47,48]. Once ABA is produced, ABA-bound receptors bind ABA, inhibit PP2Cs, such as ABI1 and ABI2, thereby activate SnRK2s [49]. The SnRK2s are involved in plant responses to abiotic stresses, which can phosphorylate the ABA-responsive element binding factor (ABF) and trigger the expression of ABA-responsive genes [50]. In our study, three ABA receptors, PYLs, were significantly down-regulated in roots or leaves ( Figure 5), which is consistent with results in tea (Camellia sinensis L.) [51] and grape (Vitis vinifera L.) [52]. The down-regulation of PYLs can reduce plant sensitivity to ABA in response to salt stress and help plants adapt to salinity. The expressional differences of PP2Cs between roots and leaves in pomegranate are similar to other plants, which indicate that PP2Cs may participate in response to salinity via different strategies [52,53]. PYR/PYLs may regulate SnRK2s directly or indirectly, however, whether SnRK2s responding to various abiotic stresses in an ABA-independent or ABA-dependent pathway needs further investigation.
Ca 2+ referred as a second messenger plays a very important role in many stress-related signaling transduction pathways. Under various abiotic stresses, over-accumulation of ions induce a temporary fluctuations in plant cytosolic ([Ca 2+ ] cyt ) levels [54]. The genes involved in CIPKs and CDPKs pathways, such as ACAs, CAXs, GLRs, CaM/CMLs, CBLs etc., participate in transporting and binding Ca 2+ , sensing and relaying signals in plant cell under salinity stress [55,56]. In our study, we found two ACAs and two CAXs located on vacuoles were up-regulated in T2R, and the majority of CaM/CMLs were up-regulated in roots but down-regulated in leaves ( Figure 6). Previous reports revealed a central role for CaMs in the regulation of Ca 2+ channels and pumps, like CNGCs and ACAs [57,58]. The negative effect of CaM on CNGC activity provides a direct feedback pathway to restrict the influx of Ca 2+ into plant cells [58]. The CaM stimulates the activity of these ACAs by preventing their auto-inhibition [58]. Collectively, under salt stress, pomegranate plants restrain the excess influx of Ca 2+ into root cells under salinity. Meanwhile, it coped with the excess Ca 2+ via removing Ca 2+ from the cytosol by ACAs and CAXs in plasma membrane, including efflux of excess Ca 2+ into the outer rhizosphere and/or the influx into vacuoles [59,60].
Specifically, transcription factors, such as NAC, MYB, AP2/ERF, bHLH, WRKY, GRAS family genes, regulated the expression of downstream genes to cope with various stresses [17]. Our study indicated that the TFs of pomegranate participated in the response to salinity with different patterns. For example, most of NAC and C2H2 genes were up-regulated in roots, while most of MYB and MYB-related genes were down-regulated in leaves ( Table 2; Table S6). Interestingly, nine GRASs were all up-regulated in T2R, which is mainly due to that these proteins play roles in plant development, including root development, axillary shoot development, and maintenance of the shoot apical meristem [61]. We also found 57 DEGs coding for salt-inducted effectors, including two SODs, two APXs, nineteen PODs, five LEAs, five AQPs, ten HSFs, three AKTs, and one HKT (Table S7). Two APXs were up-regulated and two SODs were down-regulated in leaves. Fourteen and nine PODs were suppressed in roots and leaves, respectively. Most of AQPs (4 of 5) were down-regulated, while all of LEAs (5 of 5) and majority of HSFs (9 of 10) were up-regulated in pomegranate tissues under salt stress. LEAs and HSPs were reported to prevent protein denaturation and maintain cell membrane fluidity under stress [62]. Then the DEGs of LEA and HSF in pomegranate contributed to mitigating salt stress. The expressions of three AKTs and one HKT were significantly down-regulated, which suggested that the influxes of K + and Na + in the cell were restricted [63].
Briefly, when pomegranate plants are exposed to high salinity, Na + enters the cell and simultaneously generates stress signals. The messengers Ca 2+ , ROSs (reactive oxygen species), and ABA transduced the signals, and then the cascades, including Ca 2+ -sensors, MAPK cascades and TFs, involved in secondary and tertiary regulatory networks were activated ( Figure S4). TFs regulated the expression of downstream functional genes, such as HSPs, LEAs, AQPs, SODs, and APXs, to cope with salt stress.

Conclusions
In summary, our research sheds light on the pomegranate molecular response mechanisms to salinity. The differentially expressed genes could also provide references for selecting salt-tolerant breeding materials in pomegranate breeding processes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/1/44/s1, Figure S1: The images of untreated and treated plants. Figure S2: Pearson correlation of all samples, Figure  S3: Pathways classification into cellular process, environmental information processing, genetic processing, metabolism and organismal system based on the KEGG analysis, Figure S4: Signaling networks in pomegranate under salt stress, Table S1: The primers used for qPCR validation of 15 DEGs, Table S2: Total numbers of DEGs assigned to databases, Table S3: DEGs involved in up-regulated patterns of pomegranate roots, Table S4: DEGs involved in down-regulated patterns of pomegranate roots, Table S5: DEGs involved in expressional patterns of pomegranate leaves, Table S6: The numbers of differently expressed TFs in pomegranate under NaCl stress,