Transcriptomics analysis of salt stress tolerance in the roots of the mangrove Avicennia officinalis

Salinity affects growth and development of plants, but mangroves exhibit exceptional salt tolerance. With direct exposure to salinity, mangrove roots possess specific adaptations to tolerate salt stress. Therefore, studying the early effects of salt on mangrove roots can help us better understand the tolerance mechanisms. Using two-month-old greenhouse-grown seedlings of the mangrove tree Avicennia officinalis subjected to NaCl treatment, we profiled gene expression changes in the roots by RNA-sequencing. Of the 6547 genes that were differentially regulated in response to salt treatment, 1404 and 5213 genes were significantly up- and down-regulated, respectively. By comparative genomics, 93 key salt tolerance-related genes were identified of which 47 were up-regulated. Upon placing all the differentially expressed genes (DEG) in known signaling pathways, it was evident that most of the DEGs involved in ethylene and auxin signaling were up-regulated while those involved in ABA signaling were down-regulated. These results imply that ABA-independent signaling pathways also play a major role in salt tolerance of A. officinalis. Further, ethylene response factors (ERFs) were abundantly expressed upon salt treatment and the Arabidopsis mutant aterf115, a homolog of AoERF114 is characterized. Overall, our results would help in understanding the possible molecular mechanism underlying salt tolerance in plants.

Identification and functional classification of DEGs. The analysis showed that 1404 unigenes were up-regulated and 5213 unigenes were down-regulated, while a large portion of the unigenes were not differentially expressed upon salt treatment in A. officinalis roots. To better understand the relevance of gene expression profile, the DEGs were grouped into six major classes based on their biological functions. About 45% of up-and 60% of down-regulated genes could not be classified based on their functions and hence were labeled as unknown. As shown in Fig. 1b and c, the major classes of genes identified were predicted to be involved in metabolic processes (up 25%, down 14%), defense and stress response (up 14%, down 13%), signal transduction (up 5%, down 7%), transport (up 5%, down 3%), transcription-related processes (4% in both) and membrane trafficking (1% in both). Among the metabolic processes class, genes involved in metabolism of glycerophospholipid, starch and sucrose, glycolysis, ether lipid, TCA cycle, oxidative phosphorylation and pyruvate were significantly regulated by salt treatment. Under defense and stress response, various genes encoding peroxidases, chaperones, cytochrome P450s, heat shock proteins, disease resistance proteins and ubiquitin-conjugated proteases were either up-or down-regulated. Catalases and glyoxylases were only up-regulated and NADH dehydrogenases,  hydroxylases, reductases, superoxide dismutase and redoxins were found to be down-regulated. Within signal transduction class, genes encoding calmodulins (CAMs), calcineurin B-like proteins (CBLs), CBL-interacting serine/threonine-protein kinases (CIPKs), LRR family proteins, mitogen activated protein kinases (MAPKs), proline-rich receptor-like protein kinases (PERKs) and serine/threonine-protein kinases were found to be differentially regulated. Genes for rac-like GTP binding proteins, ras-related proteins, serine/threonine-protein phosphatases (PP2As) and two-component response regulators were all down-regulated. About 71 genes related to various transport processes were up-regulated while 170 genes were down-regulated. Differentially expressed transporter genes are listed in Supplemental Table S2. The major classes of up-regulated transporters were ion-, sugar-and osmolyte-transporters and carriers/permeases while, ATPases and ATP-binding cassette (ABC) transporters were down-regulated. The important up-regulated genes encoding ion transporters belonged to the following families: sodium/hydrogen exchangers (NHX2 & NHX6), K + transporters (SKOR & POT13), cation/calcium exchanger (CCX3), ABC transporters, auxin efflux carrier (PIN6) and aquaporin . Similarly, some of the down-regulated transporters included K + channels and transporters (HKT1 and HAK23), vacuolar cation/proton exchangers (CAXs), plasma membrane (11) and vacuolar (17) ATPases and ABC transporters. Transcription-related processes group included transcription factors (TFs) as well as genes involved in transcription-related processes. Interestingly, ethylene response factors (ERFs), auxin response factors (ARFs), No Apical Meristem domain-containing factors (NAC2), WRKYs and basic helix-loop-helix (bHLH) TFs were up-regulated in large numbers, while TFs such as myeloblastosis (MYB), zinc finger CCCH domain-containing factors, GATA, bHLH and bZIPs were prominently down-regulated ( Table 2). Other differentially regulated TFs include general transcription factor group-E (GTEs), Trihelix TF, TGA1, heat stress TFs and MADS-box TF. A small fraction (1%) of genes related to membrane trafficking were differentially expressed. Transcripts of dynamins, snakins, vacuolar protein sorting, vesicle-associated membrane proteins and CSN4 were differentially regulated, while clathrins, syntaxin and t-SNAREs were down-regulated.
GO enrichment analysis was carried out to further clarify the biological functions of identified DEGs that were enriched in 56 GO terms. Significantly enriched terms under biological processes are; translation, response to cadmium ion, oxidation-reduction processes, response to salt stress, response to stimulus and metabolic process (Fig. 1d). In total, 2628 DEGs were enriched in 122 KEGG pathways, which include 42 metabolic pathways (q-value ≤ 0.05) (Supplemental Table S1). Abundantly enriched biosynthetic pathways include biosynthesis of secondary metabolites (467 genes), phenylpropanoids (53 genes), unsaturated fatty acids (25 genes), valine, leucine, isoleucine (22 genes) and flavonoids (21 genes).

Experimental validation of DEGs.
To assess the reliability of our RNA-sequencing based approach to identify salt-responsive genes in A. officinalis roots, we monitored expression of DEGs by quantitative real time PCR (qRT-PCR) analysis. From the 75 DEGs tested, about 68 DEGs (~90%) showed general agreement with their differential expression determined by RNA-seq ( Fig. S3a and b), suggesting the reliability of the transcriptome profiling data. However, qRT-PCR analysis showed much higher fold change in the expression levels of some of the DEGs compared to the RNA-seq results, while a few (~10%) showed completely contradictory results (Fig. S3c).

Identification of key salt tolerance-related genes.
To better understand the relevance of the transcriptome data obtained from A. officinalis roots, the key salt tolerance-related genes were identified by aligning the DEG sequences of A. officinalis roots with published (GEO database) root transcriptome/microarray sequences of Bruguiera gymnorhiza, rice and Arabidopsis obtained upon salt treatment. While 75 genes were obtained by alignment with rice, 21 and 14 genes were identified by alignment with Bruguiera gymnorhiza and Arabidopsis, respectively (Table 3). A total of 93 salt tolerance-related genes were obtained after removal of the repetitive genes and these are listed in Table 4. Based on their GO function, these identified genes were predicted to be involved in metabolic processes, defense and stress, signaling, transport, transcription-related processes, trafficking and cytoskeleton. Among the 93 identified genes, 13 were present in more than one dataset (highlighted in Table 4) which indicates that these could play an important role in rendering salt tolerance to plants. However, the importance of other genes cannot be ignored. The roles of some of these identified genes such as, hexokinase 46 , cationic peroxidase 47 , Trihelix TF 48,49 , NAC domain containing protein 50 , 14-3-3 51 and calmodulin 52 are well studied under salt stress. However, no studies have been carried out on many of the other genes identified. Therefore, further experimental validation would be required to understand the precise roles of these identified genes under salt stress.

Mechanism of salt tolerance in A. officinalis.
Mere identification of candidate salt tolerance-related genes in the roots of A. officinalis is not sufficient to understand the broad regulatory network that involves the functioning of these gene products in rendering salt tolerance. We reasoned that phytohormone signaling, Ca 2+ signaling and specific TFs should play important roles under salt stress to regulate many signaling pathways. Hence, all the identified DEGs that are predicted to be involved/associated with ABA, auxin and ethylene signaling pathways were analyzed in more detail, with the idea that they might reveal important signaling modules for mediating salt tolerance. In total, ~100 unigenes were ABA responsive while 65 and 61 were responsive to auxin and ethylene, respectively. While 11 of these genes were common to ABA and Auxin, 12 were common to auxin and ethylene and 17 were common to ethylene and ABA (Fig. 2a). Finally, 10 genes were found to be common in all the three pathways. In order to understand the potential roles of these genes in salt tolerance of A. officinalis, a broad signaling-network was created using the published information regarding these pathways 15,53 . Upon placing the DEGs in these known pathways, it was evident that most of the DEGs involved in ethylene, auxin and Ca 2+ signaling were up-regulated while those involved in ABA signaling were down-regulated (Fig. 2b). These results imply that several ABA-independent signaling pathways could also play a major role in salt tolerance of A. officinalis. Hence, the expression profiles of most of these genes were validated by temporal gene expression analysis using qRT-PCR (Figs S4 and S5). Role of LRR-RLK and phytohormone signaling. A number of leucine-rich repeat receptor-like kinase (LRR-RLK) genes were up-regulated in the salt-treated roots, suggesting their possible role in perception of the stress signals. Although, the exact function of many LRR-RLK genes in plants have not been understood yet, RLKs are shown to be involved in cell to cell signaling under various environmental stresses by functioning as receptors to various signals 54 . Like other RLKs, LRR-RLK could be involved in phosphorylation of MAPKs, which is supported by the up-regulation of different MAPKs (MAPK3, 8 and 9) in our study (Fig. 2b). Moreover, expression of several genes involved in ethylene biosynthesis such as methionine synthase, S-Adenosyl methionine synthetase (SAM2), SHMT and ACC oxidase were observed to be upregulated in response to salt treatment in both RNA-seq and qRT-PCR experiments (Fig. 2c). In addition, Trihelix TF which is known to interact with AP2/ERFs 49 , hexokinase1 known to be involved in ethylene signaling 55 , cationic peroxidase and glutamate synthase that are induced by ethylene leading to proline synthesis under salt stress [56][57][58] are all identified as key salt tolerance-related genes (Table 4). Therefore, we hypothesize that ethylene signaling could be playing a major role in salt tolerance of mangrove roots. This hypothesis is further supported by the up-regulation of a number of AP2/ERF TFs such as ERF 1B, 14,24,110 and 114 (Fig. 2b). The expression of some of these ERFs was also confirmed by qRT-PCR analysis (Fig. 2d). Auxin [indole-3-acetic acid (IAA)] is essential for plant growth and development. It provides key signal for the formation of lateral roots in many plants. It is produced via tryptophan-dependent and -independent biosynthetic pathways and maintains its homeostasis by processes such as degradation, conjugation to amino acids and directional transport 59,60 . In the current study, a number of auxin responsive genes such as probable indole-3-acetic acid-amido synthetase GH3.1, auxin-induced protein 5NG4, ARFs (6, 25 and 1) auxin-binding protein ABP19a, auxin-responsive protein IAA11, auxin-induced protein AUX22D, Auxin transporter ABCB10, probable auxin efflux carrier PIN6, 14-3-3 and stelar potassium outward rectifying channel SKOR were up-regulated (Figs 2c and 3c), strongly suggesting involvement of auxin signaling in the roots of A. officinalis in response to salt treatment. In further support of this hypothesis, ARF25, PIN6 and 14-3-3 were also identified as the key salt tolerance-related genes by comparative genomic analysis (Table 4).

Continued
Role of transcription factors. In addition to the up-regulation of ARFs and ERFs, a number of TFs such as NAC2, NAC7, WRKYs (9 and 22), bHLHs (130 and 137) and Trihelix were also upregulated in A. officinalis roots ( Fig. 2d and Table 2). Among these TFs, NAC2/NAC7 could be common downstream components of both auxin and ethylene signaling pathways. The induction of these TFs could be acting in an ABA-independent manner to support lateral root development 61 . Overall, our results suggest a potential involvement of a number of TFs and a crosstalk between auxin and ethylene signaling in response to salt treatment. This is further supported by down-regulation of a number of genes involved upstream and downstream of ABA signaling pathway and many ABA responsive TFs such as MYBs, ABFs and bZIPs. Our results also provide possible link between upregulation of WRKY 9 and officinalis. The Trihelix TF, could also play an important role in stress signaling because this TF has been identified as a key salt tolerance-related gene (Table 4) and was highly induced by salt in Arabidopsis 49 . The association of hormones and different TFs was further evident in the interaction network created using ARACNE and CYTOSCAPE to identify important stress-responsive pathways ( Fig. 3a and b). Among different potential genes, the AoARF25 and AoERF114 were preferred as specific nodes due to their significant upregulation in response to salt treatment in the roots of A. officinalis. In addition, AoARF25 was also identified as a salt tolerant gene by comparative genomic analysis ( Table 4). The network analysis with AoARF25 showed interactions with 76 up-regulated genes (Table S3) including PIN6 (an auxin efflux carrier), MAPKs (MAPK8 and MMK1), COP9 signalosome complex subunit 4 (CSN4), bHLH TF, serine/threonine-protein kinase (SOS2) and potassium channel (SKOR). Whereas, AoERF114 was possibly interacting with proline-rich receptor-like protein kinase (PERK2), betaine aldehyde dehydrogenase (BADH), cold-shock proteins (CSP1, CSP3), extensin (HRGP), pyruvate dehydrogenase (PDH) and elongation factor 1-gamma 3. To further validate the expression pattern of these genes, qRT-PCR analysis was performed in both A. officinalis (Fig. 3c) and Arabidopsis by profiling the expression of some of the homologous genes (Fig. 3d). Both the network and the gene expression results provide extensive information regarding the involvement of complex interactions of phytohormones in response to salt treatment in A. officinalis.
Role of Ca 2+ signaling. Calcium is one of the most important second messengers required for plant signaling networks under abiotic stresses. Many external stimuli like salt stress are known to increase Ca 2+ levels in the cytosol within seconds through various Ca 2+ transporters and pumps 64 . In the current study, several genes related to Ca 2+ signaling were up-regulated (Figs 2b and S5b). Calcium-transporting ATPases (ACAs), Ca + /H + exchangers (CAXs) and CNGCs were differentially expressed, which could be leading to the Ca 2+ fluxes during salt stress. The CNGC20 could be involved in Ca 2+ influx across PM, while ACAs (ACA12, ACA2) and CAX2 could be involved in efflux across PM and tonoplast, respectively. The increased Ca 2+ levels are sensed by the calcium sensors such as CaMs and CBLs (SOS3) 65 , both of which were up-regulated in the current study. Further, SOS3 interacts with specific Ser/Thr kinases (CIPKs/SOS2) and this SOS3-SOS2 complex would activate various downstream targets under salt stress 53,64 . This complex would activate Na + /H + antiporters SOS1 and NHX1, leading to Na + efflux across PM and Na + compartmentalization into vacuoles, respectively 53,66 . They are also known to block Na + uptake in the roots by HKT1, leading to salt tolerance in plants. Concomitantly, both NHX1and NHX6 were up-regulated; whereas HKT1 was down-regulated in our study. Although, a number of SOS1 genes were identified, they were not differentially regulated with 24 h salt treatment. Also, the vacuolar ATPase, VHA was up-regulated, which has previously been shown to be activated by this complex. Identification of CAX2 and calmodulin as salt tolerance-related genes (Table 4) makes their role even more significant in salt tolerance. Overall, our results suggest that in addition to phytohormones, Ca 2+ signaling could play an important role in salt tolerance of A. officinalis. However, further experiments are required to confirm the roles of these identified candidate genes in salt tolerance of mangroves.  Table 4. Key salt tolerance-related genes identified from the root transcriptome of A. officinalis: The DEG sequences of A. officinalis were aligned with 4 of the published root transcriptome and microarray data that were obtained from the roots of Bruguiera gymnorhiza, rice and Arabidopsis in response to salt treatment. Column 1 shows the unigene ID, while column 2 represents the ID of the reference gene. Percent similarity between the sequences of A. officinalis and the reference plant is shown in column 3 while e-value is given in the column 4. Column 5 indicates the bit score and the gene name is given in the column 6. The genes that were present in more than one data set are highlighted in bold.
Scientific RePoRTs | 7: 10031 | DOI:10.1038/s41598-017-10730-2 Ethylene response factor (AoERF114) plays an important role in salt tolerance: a case study using Arabidopsis mutant aterf115. Integration of the gene expression data, network analysis and the validation results suggest the importance of ethylene signaling via ERFs in salt tolerance of A. officinalis. In our study, both AoERF114 and AoERF1B were significantly up-regulated (Fig. 2d). Although, ERF1B was identified as a salt tolerance-related gene (Table 4), its expression was suppressed upon salt treatment in Chrysanthemum 67 . Therefore, we chose AoEFR114 for a more elaborate study. Based on the phylogenetic tree generated using the deduced amino acid sequence of AoERF114 and other members of the family from the database, AtERF115 emerged to be one of its homologs in Arabidopsis (Fig. 4a). Sequence alignment of the derived amino acid sequences of AoERF114 and AtERF115 showed 54% identity and 66% similarity between the two. The AP2 domain characteristic of the AP2/ERFs consisting of YRG and RAYD elements 68 was also conserved in both (Fig. 4b). Moreover, eight conserved amino acid residues that are involved in the interaction with the DNA GCC box 69 are present in AoERF114 and AtERF115. The Arabidopsis homozygous insertional mutant aterf115 (AT5G07310.1, SALK_021981 C) was obtained from TAIR to study the effect of salt treatment. The mutant was more sensitive to salt compared to wild-type (WT) seedlings (Fig. 5). Seed germination of aterf115 was significantly reduced (more than half) upon 75 mM and 100 mM NaCl treatment ( Fig. 5a and b). In addition, salt treatment significantly affected the seedling growth of aterf115 on agar plates ( Fig. 5c and d). Moreover, we found that the roots of aterf115 seedlings were shorter than that of WT seedlings when treated with 75 mM and 100 mM NaCl (Fig. 5c). In order to test whether AtERF115 responds to salt treatment, we obtained Arabidopsis lines with GUS expression driven by the promoter of AtERF115 (pAtERF115::GUS line). The GUS expression patterns in the roots showed that AtERF115 gene was induced in response to 3 to 24 h of salt treatment (Fig. 6a). Similarly, the transcript levels of AtERF115 increased upon salt treatment after 3 and 6 h as shown by qRT-PCR analysis (Fig. 6b). In addition, the expression profile of selected, known targets of ERFs were tested and were found to be significantly up-regulated by salt treatment in Arabidopsis roots (Fig. 6b). While NAC2 showed a twofold increase (after 0.5 h of salt treatment), HAK5 and RD29 showed 18-fold (after 24 h) and 80-fold (after 3 h) increase, respectively. To independently verify this data, the expression profiles of these selected target genes were checked in the aterf115 mutant seedling roots, and they were significantly reduced (Fig. 6c). Overall, these findings suggest that the ERF115 TF could be involved in ethylene signaling by regulating some of these genes.

Discussion
In this study, a comprehensive transcriptomic analysis from the roots of A. officinalis in response to salt treatment was carried out in order to identify salt-responsive genes. Despite the advancement in genome sequencing techniques, genomic information for many non-model plants is unavailable. Transcriptome profiling using mRNA-sequencing facilitates rapid generation of large datasets leading to identification and quantification of transcripts even in the absence of a reference genome sequence 70 . While transcriptome studies for a few salt secretor mangrove species have been carried out 32-34 , such information on A. officinalis is missing. In general, roots provide the first line of defense against salinity as they are in direct contact with the saline soil. This necessitates them to exhibit anatomical, physiological and molecular changes in order to adapt to such harsh environments. Therefore, the primary and important mechanisms of salt tolerance may reside in the roots. Salt tolerance is a complex phenomenon which involves the interaction of many genes that brings about tissue tolerance to osmotic stress, ion homeostasis and detoxification 1 . In the current study, several groups of potential salt-responsive genes, which could contribute to the salt tolerance of A. officinalis, were identified. Some of the important salt-related genes are discussed in relation to their known functions from other plant species.
Mangroves have been shown to accumulate high levels of organic solutes such as proline, glycinebetaine, polyols and sugars in order to overcome the salinity-induced osmotic stress 5 . In support of this, among the group of genes that affect metabolites, we observed significant up-regulation of choline monooxygenase (CMO) and BADH involved in glycinebetaine biosynthesis 71 as well as Hexokinase1 (HXK1) and trehalose 6-phosphate phosphatase (TPPA) involved in trehalose biosynthesis 72 in A. officinalis roots treated with salt. BADH and CMO have also been reported from various species like Suaeda, Halogeton, Atriplex, sugar beet etc. 22,25,73,74 . We also observed up-regulation of several key genes encoding enzymes that are involved in reactive oxygen species (ROS) scavenging and detoxification under salt stress 75,76 . These include peroxidases, catalases, glutathione peroxidases and glyoxylases. These ROS scavenging enzymes are well studied in plants and are known to be associated with various abiotic stresses. Additionally, genes related to flavonoid biosynthesis were also up-regulated in our study. Flavonoids have been shown to enhance salt tolerance by mitigating oxidative damage in soybean 77 . Up-regulation of these genes suggests that the oxidative stress is induced by salt treatment, and detoxifying as well as ROS scavenging enzymes are active in the mangrove roots as part of the metabolic adaptation towards salt tolerance as seen in other halophytes 5,74 .
Another major group of DEGs identified in our study comprised of genes involved in regulating ion uptake and transport. Roots have remarkable ability to regulate the plant's Na + and Cl − concentrations. Among the genes known to confer salt tolerance are those that are associated with ion uptake, transport to shoots, root ion homeostasis and water status 78 . Non-selective cation channels such as CNGCs are known to be involved in uptake of Na + , K + and Ca 2+ 79 while NHXs are involved in Na + , K + compartmentalization and pH homeostasis, which function by utilizing the pH gradient generated by V-ATPases 80 . K + transporters are essential to maintain the ionic balance which is altered under salt stress 81 . In A. officinalis, several plasma membrane and tonoplast transporters were up-regulated (Table S2). The CNGC could be involved in ion uptake while NHXs along with V-ATPases and K + transporters could function in ion sequestration and homeostasis. The precise function of ABC transporters is still unknown. However, they are implicated in various functions including transport of heavy metals, osmolyte, fatty acids, auxin and Na + 82-84 . Overall, the data suggest that various genes involved in osmotic adjustment, ion homeostasis, detoxification and metabolic processes are up-regulated and could play important roles in salt tolerance of A. officinalis. However, the functionality of these genes needs further experimental validation.
With the aid of comparative genomic approach, key salt tolerance-related genes in A. officinalis roots were identified, while the temporal expression profiles of a few of them were validated by qRT-PCR. Presence of some of these genes (13) in more than one dataset strengthens their importance in salt tolerance of plants. Among these genes, Glutamate synthase1 was shown to be involved in proline synthesis under salt stress in tomato and its activity was increased by ethylene in Hevea leaves 57,58 . Hexokinase1 plays important role in sugar and ethylene signaling 31,55 . These findings suggest that both metabolism-related genes could be involved in accumulation of osmolytes required for osmotic balance under salt stress in A. officinalis along with the other genes discussed earlier. Interestingly, among the four classes of TFs identified (Trihelix, NAC, ARF, ERF), Trihelix and NAC were present in more than one dataset and hence they can be important candidates for future studies in understanding salt tolerance mechanisms in plants. NACs and Trihelix TFs were also found to be differentially regulated by salt in Suaeda maritima 22 .
Although we have identified 93 genes as salt tolerance-related genes, the relevance and importance of other DEGs cannot be ignored. The changes in expression may not be the same across various species compared, because they may have different mechanisms of response, involving different molecular elements. Even if the same stress treatment is applied, it is expected that two plant species may experience different levels of stress, and accumulate a given transcript at different levels. Considering the broad diversity of salt tolerance mechanisms in plants, diverse gene expression profiles under salt treatment is not unusual. Hence, it will be important to study all the genes that are responsive to salt treatment in our attempts to unravel the tolerance mechanism.
The results based on gene network and signaling network analyses suggest that there is a crosstalk between auxin and ethylene signaling in response to salt treatment, which may operate in an ABA-independent pathway in A. officinalis. This is further supported by the observation that many genes known to be upstream and downstream of ABA signaling pathway were down-regulated (Fig. 2b). Moreover, many ABA dependent TFs such as MYBs, ABFs and bZIPs were also down-regulated (Table 2). In addition to salt and osmotic stresses, mangrove roots regularly experience submergence stress, which requires them to adapt to hypoxic conditions. Ethylene biosynthesis has been shown to be increased in roots under hypoxic condition and involvement of ethylene in inducing aerenchyma formation has also been well studied in other plant systems 85,86 . This explains the up-regulation of several AP2/ERF transcription factors and genes involved in ethylene signaling in A. officinalis. This TF either alone or in combination with other ERFs and TF families could be involved in the control and regulation of ROS accumulation and signaling that is required for adaptation to salt and sub-ambient oxygen concentration. Both higher salt levels and waterlogging inhibit root elongation 87 . While there is ion toxicity to the roots in the former, there is reduced oxygen supply in the latter. Furthermore, studies indicate that the stress hormone ethylene, which accumulates in waterlogged plants, can contribute to the regulation of lateral and adventitious root formation in a complex crosstalk with auxin 88 . Many NAC TFs are found to be responsive to auxin, ABA and abiotic stresses in Arabidopsis 61 . Similarly, we hypothesize that transcription factors NAC2/NAC7 could be promoting lateral root formation under submerged conditions in coordination with auxin and ethylene signaling pathways in the mangrove system reported here. In addition, it was shown that NAC2 in A.thaliana promoted lateral root formation and salt-induced AtNAC2 expression was dependent upon the ethylene and auxin signaling pathways, but not ABA signaling 61 . Under salt stress, ethylene signaling in plants is known to be mediated through ERFs which regulate the stress responsive genes 89 . Some of the salt responsive genes regulated by ethylene are: A dehydration responsive gene RD29 90 , NAC2 61 and high affinity potassium transporter (HAK5) 91,92 . Our results suggest that AoERF114/AtERF115 could be regulating some of these genes to render salt tolerance as all of these genes were suppressed in the salt sensitive aterf115 roots, while their expressions were induced upon salt treatment in the WT (Fig. 6b and c). Therefore, further studies on ERFs (especially ERF115) are necessary to unravel the mechanism of stress response in plants. In addition to phytohormones, Ca 2+ signaling is well known to play important role in salt tolerance of plants 93 . Up-regulation of several genes involved in this signaling pathway, underlines their importance in A. officinalis. The NHXs (NHX2 and NHX6) and vacuolar ATPases could be important for Na + /K + homeostasis in A. officinalis roots, similar to their function in other plants 94,95 . NHX1 and V-ATPases from several halophytes have proven to increase salt tolerance in glycophytes 25 . However, NHXs were not responsive to salt treatment in the leaves of halophytes like Suaeda and Halogeton 22,74 . Overall, identification of a number of up-regulated genes associated with ethylene, auxin as well as Ca 2+ signaling provides critical information regarding the involvement of these signaling pathways in salt tolerance of mangroves.
In conclusion, a comprehensive transcriptome profile of A. officinalis roots is provided in this study. Our data helped to identify numerous salt tolerance-related genes as part of an overall list of DEGs in response to salt treatment, 93 of which would conceivably be playing meaningful roles in conferring salt tolerance in mangroves and other plants. The transcriptome data together with our results from Arabidopsis mutant (aterf115) analysis helped to reveal an important role for this ERF in salt tolerance. Our study also revealed the interplay of various A. officinalis genes involved in ethylene-, auxin-and Ca 2+ -mediated signaling pathways in a salt-responsive manner. This information may be used for future studies on salt tolerance in plants.

Materials and Methods
Plant materials and growth conditions (A. officinalis). The propagules of Avicennia officinalis L. (A. officinalis) were collected during fruiting seasons from the mangrove swamps in Singapore (Berlayer Creek and Sungei Buloh Wetland Reserve). The seedlings were maintained in NaCl-free conditions by growing in potting mixture (Far East Flora, Singapore), until they reached the four-node stage (~2 months) in a greenhouse (25-35 °C, 60-90% relative humidity; 12 h photoperiod), after which they were carefully transferred to pots containing sand and were allowed to adapt for two days by watering with half-strength Hoagland's solution. The plants were then treated with half-strength Hoagland's solution containing 500 mM NaCl for 24 hours. and frozen in liquid nitrogen for total RNA isolation. For histochemical GUS expression analysis, 5-day old seedlings were treated with 50 mM NaCl for various time periods (0 h, 1 h, 3 h, 6 h and 24 h). For seed germination studies, the sterilized and cold stratified seeds were sown on MS Agar plate with and without NaCl and allowed to germinate as mentioned above. The number of germinated seeds was counted from day 1 to day 4 and the root lengths were measured and photographed 7 days after germination.

RNA isolation.
Total RNA was isolated from roots of control and treated (500 mM NaCl for 24 h) greenhouse-grown A. officinalis using Qiagen RNeasy kit (QIAGEN) and DNase treated (RNase-free DNase set, QIAGEN) according to the manufacturer's instructions. The quality of RNA samples was determined using a 2100 Bioanalyzer (Agilent Technologies). For each sample, at least 20 µg of total RNA was sent to Beijing Genomics Institute for Illumina sequencing (commercial service). For qRT-PCR experiments, total RNA was isolated from the roots of control and treated (500 mM NaCl for varying time periods; 0 h, 0.5 h, 1 h, 2 h, 4 h, 8 h, 12 h and 24 h) greenhouse-grown A. officinalis and control and treated (150 mM NaCl for varying time periods; 0 h, 1 h, 3 h, 6 h and 24 h) roots of one-week-old Arabidopsis seedlings as described above. An aliquot of this RNA (1 µg) was used to synthesize cDNA using Maxima first strand cDNA synthesis kit for qRT-PCR (Thermo Scientific) following manufacturer's instructions. cDNA library preparation, sequencing and transcriptome de novo assembly. For each sample, mRNAs were purified using oligo (dT)-attached magnetic beads and fragmented into small pieces (100-400 bp). The cDNA library was prepared by synthesizing the first and second strand cDNAs, using the mRNA fragments as templates primed with random hexamers. The synthesized cDNAs were end repaired, 3' adenylated and ligated with sequencing adaptors. Suitable fragments (~200 bp) were selected by agarose gel electrophoresis and enriched by PCR amplification. Finally, these cDNA libraries were sequenced using Illumina HiSeq ™ 2000 sequencer (Beijing Genomics Institute, BGI, Shenzhen, Guangdong, China). Image data obtained from the sequencing machine was transformed by base calling into sequence data (raw reads) and stored in fastq format. Transcriptome de novo assembly was performed using the short read program Trinity (version release-20121005) 97 . The Trinity software first combined clean reads with a specific length of overlap to form longer fragments without Ns, forming contigs. Next, the contigs were connected to obtain consensus sequences that contained the least Ns and could not be extended on either end. Such sequences were defined as unigenes. Finally, the sequence orientations of the all-unigenes were determined by Blastx against NCBI non-redundant (Nr) protein database, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Clusters of Orthologous Groups (COG) with e-value cut off of <10 −5 . Unigenes that could not be aligned to any of the four databases were scanned using EST Scan 98 , which produced a nucleotide sequence (5′-3′) direction and amino sequence of the predicted coding region. The transcriptome data of this work has been deposited to the NCBI website (GEO GSE73807).

Data analysis.
For functional annotation, unigene sequences were first aligned using Blastx to the Nr, Swiss-Prot, KEGG, and COG protein databases (E-value < 10 −5 ), which retrieved proteins with the highest sequence similarity to A. officinalis unigenes in addition to their protein functional annotations. Sequence searches were performed by querying the NCBI Nr protein database using the Blastx algorithm (E-value < 10 −5 ) 99 . After Nr annotation, the Blast2GO program 100 was used to obtain Gene Ontology (GO) annotations and the WEGO software 101 was used to perform GO functional classification of all unigenes to determine the distribution of gene functions at the macro level. KEGG annotation was carried out to obtain pathway annotations for Scientific RePoRTs | 7: 10031 | DOI:10.1038/s41598-017-10730-2 unigenes. Later, unigenes were aligned to the COG database to predict and classify potential functions based on known orthologous gene products using pathfinder software (version release 63.0). Gene expression analysis was carried out using reads per kilobase per million reads (RPKM) method 102 . For a given unigene, RPKM values were generated using SOAP (version release 2.21). A rigorous algorithm was used to identify differentially expressed genes (DEGs) in salt-treated roots compared to untreated roots. False discovery rate (FDR) ≤ 0.001, the absolute value of log2Ratio ≥2 and P-value ≤ 0.001 was used as the threshold to judge the significance of differential gene expression 103 . For pathway and GO enrichment analysis, all DEGs were mapped to KEGG and GO databases (http://www.geneontology.org/). By using hypergeometric test, significantly enriched GO terms were identified in comparison with the genome background. In addition, the DEGs were classified into various GO categories, based on the published databases and reports on particular genes. To identify important salt tolerance-related genes, the sequences of the DEGs were aligned with the published transcriptome/microarray sequences obtained in response to salt treatment from roots of Arabidopsis, rice and a mangrove species Bruguiera. The main criteria for choosing these species was that the transcriptomic/microarray sequences were obtained in response to salt treatment from the roots of the plants. The commonality among the various datasets (from the different plant species) used was that they were all "salt responsive datasets from roots". From Arabidopsis, two published datasets were used 104, 105 with the GEO IDs; GDS3216 and GSE46208. Similarly, two published datasets were used from rice 106,107 with the corresponding GEO IDs; GSE20746 and GSE14403 and one published data with the GEO ID GSE10942 was used from Bruguiera gymnorhiza 31 . Histochemical GUS staining. Transgenic Arabidopsis seedlings containing pAtERF115::GUS fusion constructs were treated as described above. GUS histochemical staining was performed by vacuum-infiltrating the seedlings immersed in GUS staining solution (0.1 M sodium phosphate buffer pH 7.0, 10 mM EDTA, 0.1% Triton-X, 2 mM 5-bromo-4-chloro-3-indolyl glucuronide (X-Gluc)) for 5 min followed by overnight incubation in the dark at 37 °C without shaking. Staining solution was removed and several washes with 50% ethanol was performed until the chlorophyll was bleached and tissues cleared. The image of blue colored whole seedlings with various salt treatments was recorded using a stereo microscope (Leica DIC 310 FX). GUS-stained tissues and plants shown in this paper represent the typical results of at least six independent plants for each treatment.
Quantitative real-time PCR (qRT-PCR) analysis. The qRT-PCR for differentially expressed genes was performed using the Stepone Real-Time PCR machine (Applied Biosystems) with the following programme: 20 s at 95 °C followed by 40 cycles of 03 s at 95 °C and 30 s at 60 °C. The SYBR Fast ABI Prism PCR kit from KAPA was used for qPCR analysis. The reaction mixture consisted of 5.2 μL master mix (provided in the kit), 0.2 μM FW primer, 0.2 μM RV primer, 3.4 μL nuclease-free water, and 1 μL sample cDNA template for a final volume of 10 μL. All of the data were analyzed using the StepOne TM Software (v2.1, ABI). The primers were designed using the sequences obtained by RNA sequencing and are listed in Supplemental Table S4. Constitutively expressed Ubiquitin 10 was used as internal control.
Network analysis. Differentially expressed up-or down-regulated genes were extracted from the RNA sequencing data of the root samples of A. officinalis. Gene networks for selected genes were constructed using Algorithm for the Reconstruction of Accurate Cellular Networks (ARACNE) algorithm 108 . ARACNE uses the mutual information of the features to determine the connection between genes. The features included in these networks were gene expression (RPKM) and transcript-to-SWISPROT protein alignment score, among others. Based on ARACNE output, the final gene network graphs were created using Cytoscape 109 .