Transcriptome analysis identifies genes involved in adventitious branches formation of Gracilaria lichenoides in vitro

Tissue culture could solve the problems associated with Gracilaria cultivation, including the consistent supply of high-quality seed stock, strain improvement, and efficient mass culture of high-yielding commercial strains. However, STC lags behind that of higher plants because of the paucity of genomic information. Transcriptome analysis and the identification of potential unigenes involved in the formation and regeneration of callus or direct induction of ABs are essential. Herein, the CK, EWAB and NPA G. lichenoides transcriptomes were analyzed using the Illumina sequencing platform in first time. A total of 17,922,453,300 nucleotide clean bases were generated and assembled into 21,294 unigenes, providing a total gene space of 400,912,038 nucleotides with an average length of 1,883 and N 50 of 5,055 nucleotides and a G + C content of 52.02%. BLAST analysis resulted in the assignment of 13,724 (97.5%), 3,740 (26.6%), 9,934 (70.6%), 10,611 (75.4%), 9,490 (67.4%), and 7,773 (55.2%) unigenes were annotated to the NR, NT, Swiss-Prot, KEGG, COG, and GO databases, respectively, and the total of annotated unigenes was 14,070. A total of 17,099 transcripts were predicted to possess open reading frames, including 3,238 predicted and 13,861 blasted based on protein databases. In addition, 3,287 SSRs were detected in G.lichenoides, providing further support for genetic variation and marker-assisted selection in the future. Our results suggest that auxin polar transport, auxin signal transduction, crosstalk with other endogenous plant hormones and antioxidant systems, play important roles for ABs formation in G. lichenoides explants in vitro. The present findings will facilitate further studies on gene discovery and on the molecular mechanisms underlying the tissue culture of seaweed.

through the introduction of exogenous plant hormones. Despite the addition of auxins and cytokinins (CTK) alone or in combination, regeneration from callus is difficult. Furthermore, the main morphological response of explants of Gracilariales to in vitro culture is the direct production of adventitious branches (ABs) within a shorter period of time than that required for callus formation as a result of wounding, as reported for G. blodgettii, Gracilariopsis lemaneiformis, G. tenuistipitata, G. asiatica, G. vermiculophylla, G. domingensis, and G. changii [21][22][23][24][25] . Similar results were reported for Gelidiales, Halymeniales, Ceramiales and Gigartinales [26][27][28] , which illustrates cell totipotency and is an alternative morphological response that could be used for seedling production instead of thallus regeneration from callus. However, few studies have investigated the mechanisms of ABs formation, and improving our understanding of the process of ABs formation and development in seaweed species would be of value.
Studies addressing the role of exogenous plant hormones in the formation of ABs in seaweed explants have reported inconsistent results. The differential sensitivity of explants to plant hormones might be due to an inherent regulatory mechanism rather than the concentration of plant growth regulators 29 . Aguirre- Lipperheide et al. showed that the direct regeneration of micropropagules from the explants of red seaweed is regulated by internal factors rather than exogenous plant hormones or other external factors 30 . Endogenous plant hormones, which play a significant role in almost every aspect of plant growth and development, have been detected in green, brown and red seaweeds 31 . Nevertheless, the evidence for a more direct physiological and developmental role of endogenous plant hormones in algae is limited and inconclusive. Furthermore, auxin is the only endogenous hormone involved in polar transport. The establishment of an auxin gradient and polar auxin transport (PAT) mediated by PIN1, a subfamily of ABC transporters (ABCBs), is essential for callus induction and somatic embryogenesis in model land plants such as Arabidopsis [32][33][34] . Understanding the role of endogenous hormones in the regulation of seaweed growth and development is important to improve tissue culture techniques.
Sequencing technology has developed rapidly in recent years, and transcriptome analysis has become increasingly powerful. RNA sequencing (RNA-seq) enables the rapid and comprehensive analysis of plants and other organisms 35,36 . The intensive study of transcriptome sequencing has provided a cost-effective method for the qualitative and quantitative analysis of gene transcripts in many non-model plant species such as Gracilaria 37 . In the present study, we used the Illumina sequencing platform to examine the transcriptional profiles of samples of G. lichenoides exposed to different treatments, including a control(CK), explants with adventitious branches (EWAB), and explants treated with N-1-naphthylphthalamic acid (NPA), a PAT inhibitor.

Materials and Methods
G. lichenoides culture. G. lichenoides thalli were collected from the Dongshan coast, Fujian Province. The fresh samples were cleared of mud and epiphytes using a soft brush and filtered seawater, and cultured in tanks with f/2 medium 38 at 25 ± 1 °C under cool-white fluorescent illumination (30 μ mol photons m −2 s −1 ) in a12:12 light:dark (L:D) cycle. The water inside the tanks was agitated with constant aeration and the culture medium was replaced every 3 days.
Tissue culture and sample collection. Healthy thalli of G. lichenoides were examined under a microscope, and those with the least contamination of epiphytes were selected and sterilized according to the methodology described by Shen et al 39 . Briefly, thalli were washed three times in autoclaved seawater with a writing brush, three times in an ultrasonic cleaner (Branson, Beijing, China) and three times in 0.7% potassium iodide (KI).
The sterilized thalli were cut into short explants (5-7 mm) using a sterile surgical blade and inoculated in a 200-mL culture vessel for liquid media with f/2 medium (21 ± 1 °C, 30 μ mol photons m −2 s −1 irradiance, 14-h: 10-h, light-dark cycle), which was replaced every 3 days. The effect of auxins (indole-3-acetic acid, IAA) and NPA (1 and 100 μ M, respectively) and their combination on callus formation or AB induction in G. lichenoides explants was determined. A treatment control (without IAA or NPA) was prepared simultaneously. After exposure to the different conditions, the samples were subjected to extraction of RNA. RNA extraction,de novo transcriptome assembly and annotation. The transcriptional profiles of the differently treated samples of G. lichenoides (CK, EWAB and NPA) were examined using the Illumina sequencing platform. Total RNA was extracted according to Method 3 as described by Chan et al 40 . The extracted RNA was qualified and quantified using a Nanodrop ND-1000 Spectrophotometer (Nanodrop Technologies,Wilmington, DE, USA) and all samples showed a 260/280 nm ratio between 1.9 and 2.1 with no sign of degradation. The construction and sequencing of cDNA libraries were performed by BGI (Shenzhen, China) on an Illumina HiSeq 2000 platform (San Diego, CA, USA) in accordance with the manufacturer's instructions.
A filter was used to remove low-quality reads, which were defined as those with adaptor or ambiguous bases ('N' < 5%), as well as all reads with nucleotides that had Phred 41 quality scores < 20, where a Phred score of 20 corresponds to a 1% expected error rate. Ultimately, clean reads were obtained by this filtering process and assembled into unigenes with the Trinity program 42 (release-20130225, http://trinityrnaseq.sourceforge.net/). Since there is no reference genome of G. lichenoides, a k-mer value cutoff of 25 was used after removing redundant nucleotide sequences by Tgicl (v2.1, http://sourceforge.net/projects/tgicl/files/tgicl%20v2.1/).
All the unigenes were used for functional annotation by BLASTX against various protein databases, such as non-redundant protein (Nr) databases, the Swiss-Prot database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and the Cluster of Orthologous Groups (COG) database, as well as the NCBI non-redundant nucleotide sequence (Nt) database with an E value cutoff of 10 −5 43 , and the best aligning results were selected to annotate the unigenes. If the aligning results from different databases were in conflict with each other, the results from the Nr database were preferentially selected, followed by the Swissprot, KEGG and COG databases. CDS (coding sequence) was obtained by blasting against these four databases, as well as by using ESTScan 44 . Gene ontology (GO) annotation for the unigenes was based on the best alignment of Swissprot obtained using the Blast2GO program 45 , and WEGO software 46 was then used to plot the GO annotation results. Pathway assignments were performed according to the KEGG pathway database 47 . Differential expression analysis of unigenes. The specifically expressed genes were screened from the distinct G. lichenoides explants by differential expression analysis. The expression levels of unigenes were calculated using reads per kilobase transcriptome per million mapped reads (RPKMs) to eliminate the effect of gene length and sequencing level on the calculation of gene expression 48 . The calculated results were directly used to analyze differences in gene expression levels between the different samples using the statistical comparison method 49 . The results of multiple hypothesis testing were corrected using the false discovery rate (FDR) control method, and the ratio of RPKMs was converted to fold-changes in the expression of each gene in two samples simultaneously. Significant differentially expressed genes were screened according to the FDR threshold of ≤ 0.001 and the absolute value of |log 2 Ratio| ≥ 1. When differentially expressed genes were identified, all significantly differentially expressed genes (DEGs) were mapped to the KEGG database to conduct pathway enrichment analysis to identify significantly enriched metabolic pathways or signal transduction pathways with differentially expressed genes using hypergeometric tests. The calculated q-values were subjected to Bonferroni correction, with a corrected q-value ≤ 0.05 as a threshold. Unigenes were annotated with the NR, NT, Swiss-Prot, KEGG, COG and GO databases, and the number of unigenes annotated with each database was counted.

Results
Tissue culture of G. lichenoides. In all experiments, the main morphogenetic response of the explants of G. lichenoides in vitro was the production of direct ABs from the cortex within 8 days rather than callus formation as a result of wounding. The rate of ABs formation was 94.9%. ABs became mature with rhizoids 4 months after leaving cut segments (Fig. 1). NPA inhibited ABs formation in the explants of G. lichenoides, and this effect was not reversed by the addition of IAA. However, ABs formed after transfer of the explants to normal medium.
Illumina sequencing, assembly and sequence analysis. A total of 17,922,453,300 nucleotide (nt) paired-end reads with a read length of 100 bp were obtained conveniently and efficiently using the three cDNA libraries, which were sequenced with Illumina HiSeq ™ 2000. The G + C percentages were all approximately 52% (Table 1). The percentage of Q20 bases (those with a base quality > 20 and an error rate < 0.01) was higher than 99%,   Table 1. Output statistics of sequencing. * Total reads and total nucleotides are actually clean reads and clean nucleotides. The total nucleotides should be greater than contract provision. The Q20percentage is the proportion of nucleotides with a quality value higher than 20. The N percentage is the proportion of unknown nucleotides in clean reads. The GC percentage is the proportion of guanidine and cytosine nucleotides among total nucleotides. indicating the high quality of the raw sequencing reads of the sample (Table 1). After stringent quality assessment and data filtering performed using the Trinity de novo assembly program, the next-generation short-read sequences were assembled into 21,294 unigenes, which represented putative protein-coding transcripts. These consisted of the partial non-overlapping coding regions of G. lichenoides, with an N50 [length of the smallest transcripts in the set containing the fewest (largest)] length of 5,055 nt and a mean length of 1,883 nt. An overview of the assembly results is provided in Table 2.
Functional annotation of the transcriptome. Regarding  In the Nr annotation, 6,709 sequences had perfect matches with E values < 10 −45 and 5.2% of the matches had a similarity over 95% (Fig. 2a,b). The species distribution showed that 30.0% had top matches to the Galdieria sulphuraria genes (Fig. 2c). The next two top matching species were the Cyanidioschyzon merolae strain 10D and  Table 2. Statistics of assembly quality. * N50: length of the smallest transcripts in the set that contain the fewest (largest). Guillardia theta CCMP2712, accounting for 4.5 and 3.7%, respectively. Only 1.16% of the sequences showed matches to Gracilaria. GO analysis was conducted on the annotated 58,214 unique sequences using the Blast2GO program. The unique sequences were successfully assigned to biological processes, cellular components, and a molecular function GO category, respectively (Fig. 3). The most abundant were "cellular processes" (4,856 members) and metabolic processes (4,606 members) in the biological process category, and cells (5,651 members), cell parts (5,650 members) and organelles (4,677 members) for the cellular component category. Most of the sequences were classified into catalytic activity (3,744 members) and binding (3,211 members) functions in the molecular function category. Additional information derived from the GO annotation was as follows: of the unigenes analyzed, 2,123 were annotated into the response to stimulus category, 580 into reproduction, 496 into reproductive processes, 336 into transporter activity, and 40 into antioxidant activity, which may provide an explanation for the in vitro regeneration mechanism of macroalgae explants.
To further evaluate the completeness of our transcriptome library and the effectiveness of our annotation process, all annotated unigenes were aligned to the COG database to predict and classify possible functions. In total, 9,490 sequences had a COG classification. Among the 25 COG categories, "general function prediction only" represented the largest group (2,882 members), followed by "translation, ribosomal structure and biogenesis" (2,443 members) and "post-translational modification, protein turnover, chaperones" (1,568 members) (Fig. 4).
The KEGG pathway database records the networks of molecular interactions in the cell and variants specific to particular organisms. Pathway-based analyses are important to further understand the biological functions and interactions of genes. In total, 10,611 sequences were assigned to 126 KEGG pathways. The pathways with the most representation among the unique sequences were the metabolic pathways (2,654 unigenes), followed by Ribosome (1,233 unigenes) and biosynthesis of secondary metabolites (1,089 unigenes).A hallmark of dedifferentiation and the subsequent formation of ABs is the ability to control cell division and cell wall accumulation associated with polysaccharide metabolism. This is associated with the expression of the corresponding hydrolytic enzymes and their accumulation of storage reserves and secondary metabolites. Furthermore, some important pathways related to plant growth, development, morphogenesis and response to stress stimuli, including ABC transporters (222 unigenes), ubiquitin-mediated proteolysis (152 unigenes), glutathione metabolism (119 unigenes), peroxisome(106 unigenes), carbon fixation in photosynthetic organisms (92 unigenes), phosphatidylinositol signaling system (67 unigenes), plant hormone signal transduction(61 unigenes), zeatin biosynthesis(38 unigenes) and others, have also been annotated successfully. Global analysis of differentially expressed genes. RSEM software was applied to map the clean reads in different libraries to the corresponding assembled consensus sequences. A strict binomial distribution algorithm was used to identify DEGs, considering the differences in library size for the differential selection of DEGs among libraries. Pairwise comparisons between the above three groups (CK, EWAB and NPA) were performed for differential expression analysis (Fig. 5). To gain an understanding of the genes expressed in this time course study, we selected genes that were differentially expressed between two groups by log 2 (fold change) > 10.
In the CK/EWAB pairwise comparison, 9,525 genes showed differential expression, including 2,800 up-regulated genes and 6,725 down-regulated genes in EWAB. In the CK/NPA pairwise comparison, there were 9,891 DEGs, including 3,651 up-regulated genes and 6,240 down-regulated genes in NPA (Fig. 5). GO term enrichment analysis was performed to evaluate significantly over-represented GO terms in the differentially expressed genes, which resulted in a total of 49 GO terms with a p-value < 0.05 and FDR < 0.01. These GO terms included functions and processes such as reproduction, reproductive processes (Fig. 3), auxin polar transport, and response to auxin stimulation (Table 3). Furthermore, 4,087 (42.9%) and 4, 217 (42.6%) unigenes were annotated to the KEGG pathway in CK/EWAB and CK/NPA, respectively. The majority of those annotated genes between CK/EWAB and CK/NPA had similar expression patterns and were related to plant growth and development, cell vegetation and differentiation, and response to stimulus, mainly including ABC transporters, ubiquitin-mediated proteolysis, plant hormone signal transduction and other pathways ( Table 4). The process of in vitro regeneration of explants from macroalgae was complex.
Between EWAB and NPA, 4,264 DEGs were discovered, of which 2,688 genes were up-regulated and 1,576 were down-regulated in NPA (Fig. 5). GO term enrichment analysis included functions and processes such as auxin polar transport (one up-regulated unigene), auxin-mediated signaling pathways (two up-and two down-regulated unigenes), response to the auxin stimulus(six up-and two down-regulated unigenes), and auxin biosynthetic processes(one up-regulated unigene) ( Table 3). In addition, 2,108 (49.4%) unigenes were annotated to the KEGG pathway, mainly including the ABCB1 family of ABC transporters (five up-and one down-regulated unigenes), ubiquitin-mediated proteolysis (10 up-and one down-regulated unigenes), plant hormone signal transduction (10 up-regulated unigenes) and antioxidant systems (18 up-andthree down-regulated unigenes) ( Table 4). Given the known role of auxin in regulating plant embryogenesis, the competence for AB induction might be the result of regulated auxin responses of these cells. These annotations provide a valuable resource for investigating specific processes, functions, and pathways and allow for the identification of novel genes involved in these pathways.
Heterogeneous SSR detection. SSR detection performed using MISA software with unigenes as a reference resulted in the identification of 3,287 high-quality SSRs. This dataset included 1,143 SSRs with two nucleotides repeated more than six times, 1,558 SSRs with three nucleotides repeated more than five times, and 228 SSRs with five or more nucleotides repeated more than four times. The SSR discovery revealed that tri-nucleotide repeats are the most abundant because any addition or reduction in repeat numbers would not cause frame shift mutations.   Among TriNR, GCG was the most common motif, followed by GGC and CGC. The most abundant motifs were GC in DiNR, GCGT in TetraNR, CCAAA in PentaNR, and CCGCGT in HexaNR (Fig. 6).

Discussion
In the present study, although callus induction was not observed, ABs were successfully formed on the cut surfaces of G. lichenoides explants in response to most treatments, showing the capacity for regeneration after wounding. Similar results were obtained in Gelidiales, Halymeniales, Ceramiales and Gigartinales, and other species of Gracilaria 25-28 . Huang and Fujita reported that callus was not induced from explants of G. textorii using IAA 15  These studies suggest that the differential sensitivity of macroalgae explants to plant hormones could be attributed to an inherent regulatory mechanism rather than the concentration or combination of plant hormones 29 . In our study, NPA, an inhibitor of PAT, effectively prevented the formation of ABs in G. lichenoides, suggesting that PAT plays an important role in the process of regeneration of macroalgae. This was consistent with results obtained in higher plants [32][33][34] . Moreover, the detection of increases in DNA, RNA and protein content prior to or concomitant with morphological responses indicates changes in gene expression 51 . The development of STC techniques in the near future and their combination with molecular genetics could provide support for the same biotechnological applications used in higher plants in the genomic age 26,52,53 . Transcriptome sequencing is one of the most important tools for expression pattern identification and gene discovery 54,55 . Here, we used transcriptome sequencing with the Illumina platform to regenerate macroalgae explants in vitro. Functional enrichment showed the presence of many regeneration-related genes in the transcriptome (Tables 3 and 4   Gracilaria sp 31 . In the present study, 49 genes related to auxin biosynthesis (four transcripts), transport (11 transcripts), metabolism (nine transcripts), response to the auxin stimulus (27 transcripts), auxin response factor (ARF) (one transcript), and ABCB1 (26 transcripts), were differentially expressed during ABs formation (Table 3 and 4).
In general, plant development and physiology depend on a unique, plant-specific process, namely cell-to-cell or polar auxin transport. PAT is controlled by efflux proteins, including PIN and ATP binding cassette type B/P-glycoprotein/multidrug resistance (ABCB/PGP/MDR) proteins. ABCB-and PIN-mediated auxin efflux can function independently and play identical cellular but separate developmental roles 56 . E.siliculosus, for instance, has several ABCB efflux IAA transporters that participate in the polarized transport of IAA and stabilize PIN1 57 . ABCB1 and ABCB19 were recently identified as binding proteins of the synthetic auxin-efflux inhibitor, 1-N-Naphtylphthalamic acid (NPA) [58][59][60][61][62] . TWISTED DWARF1 (TWD1) binds to NPA, which disrupts the TWD1-ABCB1 interaction 58,63 . This disrupts ABCB1 activity, suggesting that TWD1 and ABCB1 represent essential components of the NPA-sensitive-efflux complex 63 . Several lines of evidence suggest that PIN proteins do not themselves act as direct targets of NPA 34,62,64,65 . Here, 26 auxin transport-associated ABCB1 transcripts were found to be differentially expressed. Nineteen ABCB1 transcripts were down-regulated and one transcript was up-regulated in the CK/EWAB group, whereas 22 (five up-and 17 down-regulated) and six (five up-and one down-regulated) ABCB1 transcripts were differentially expressed in the CK/NPA and EWAB/NPA groups, respectively ( Table 4). The annotated ABCB1 genes mainly down-regulated, which are direct binding proteins of NPA, during the formation of ABs indicating that they might negatively regulate PAT. By contrast, five ABCB1 unigenes were up-regulated by NPA, indicating that those five genes might be induced by NPA and prevent ABs formation by inhibiting PAT. Unigene11327_All was obviously up-regulated in the CK/EWAB group and down-regulated in the EWAB/NPA group, indicating that this gene might promote PAT. Unexpectedly, influx (AUX1) and efflux (PIN) proteins were not annotated in our study. This could be attributed to NPA may not be strictly specific for PIN proteins. It is thus possible that the auxin efflux transporter of G. lichenoidesis different from that of higher plants.
ARF, which is a type of transcription factor with B3 DNA binding domains, mediates auxin responses and regulates the expression of auxin-induced genes by binding to auxin response elements (AuxRE) in specific gene promoter regions 66 . Because Aux/IAA proteins can function as active repressors, dimerization between an ARF transcriptional activator bound to an AuxRE and an Aux/IAA protein would result in the repression of auxin response genes in higher plants. When auxin concentrations increase, Aux/IAA repressors are degraded rapidly by the proteasome through the action of SKP1-Cullin-F-box (SCF) E3 ligases, namely SCF TIR1/AFB , leading to the dissociation of Aux/IAA proteins from the ARF activators and the activation or repression of auxin response-related genes [67][68][69] . In the present study, only unigene8184_All, a ARF transcript, was down-regulated in group CK/EWAB. This suggests that unigene8184_All negatively regulates genes related to auxin. Moreover, the ubiquitin-mediated proteolysis pathway plays a central role in most hormone signaling pathways 70 , and we detected certain unigenes associated with the SCF-complex, including two up-and five down-regulated unigenes in group CK/ EWAB, three up-and five down-regulated unigenes in group CK/NPA, and three up-regulated unigenes in group EWAB/NPA (Table 4). However, we did not identify a transcriptional inhibitor similar to AUX/IAA and auxin specific targeting factor TIR1. Similarly, an AUX-IAA-ARF-mediated auxin signaling cascade is not present in microalgae, suggesting that these organisms have alternative pathways mediating auxin responses 71 . Therefore, although the mechanism itself may be conserved, it may rely on different transcriptional inhibitors recognized by an IAA-protein complex, such as AUX/IAA, that can target them for degradation. This could explain the inability of auxin to induce callus or promote AB formation, which differs from its effect on higher plants, during the process of macroalgae tissue culture.
CTK, which was detected previously in Gracilaria sp 72 , mediates cell division and shoot apical meristem development 73,74 . Therefore, CTK might play an important role in the regeneration process of G. lichenoides explants, as indicated by the down-regulation of nine transcripts associated with the zeatin biosynthesis pathway in group CK/ EWAB, and the up-regulation of seven and six zeatin transcripts, IPT, and cytokinin dehydrogenase [EC: 1.5.99.12] in groups CK/NPA and EWAB/NPA, respectively. Recently, Cheng et al. showed that AtIPT5 expression is negatively regulated by ARF3 through direct binding to its promoter, which affects stem cell initiation and meristem formation in Arabidopsis thaliana 75 . This suggests that the down-regulation of ARF genes may have led to the up-regulation of IPT genes.
Abscisic acid (ABA) is an important plant hormone that plays a role in several regulatory processes in response to abiotic stresses, and in seed dormancy, germination, and the regeneration of somatic embryos [76][77][78][79][80][81] . In the absence of ABA, PP2Cs (protein phosphatase 2C) directly inactivates SnRK2 (SNF1-related protein kinase) by dephosphorylating multiple residues in the kinase activation loop 82,83 . In the presence of ABA, PYR/PYL/RCAR (pyrabactin resistant/pyrabactin resistant-like/regulatory component of ABA Receptor), an ABA receptor, binds to and inhibits PP2Cs, resulting in the accumulation of phosphorylated SnRK2s and the subsequent phosphorylation of ABFs (ABA response element binding factors), activating downstream genes in the ABA signal transduction pathway 84 . Further, Singh et al. proposed that ABA, which is considered an auxin inhibitor, may reduce auxin to an optimum level required for regeneration and suppress excessive callus formation in cultures 85 . In the present study, in group CK/ EWAB, the annotated transcripts involved in the ABA signal transduction pathway included PP2C (four down-), SnRK2s (one down-), ABF (two down-regulated) transcripts, whereas no PYR/PYL/RCARs were observed. However, alternative ABA receptors may exist in G. lichenoides. Our results suggest that the ABA content of G. lichenoides explants increases during ABs formation, which reduces auxin to an optimum level to allow the formation of ABs. By contrast, in groups CK/NPA and EWAB/NPA, PP2C was up-regulated, accounting for the decrease in the ABA content of G. lichenoides explants when the auxin content increased after PAT was inhibited by NPA, which suppressed the induction of ABs.
Ethylene (ET) is related to the differentiation and growth of plant tissues cultured in vitro 86 . In the present study, CTR1, an important gene for ET signal transduction, was annotated successfully. CTR1 plays a negative regulatory role in the ET response pathway of Arabidopsis thaliana 87 . The ETR family of receptors directly activates the kinase activity of CTR1 in the absence of ET, which could activate CTR1 to phosphorylate downstream targets and turn off the ET response pathway. Binding of ET to the receptors inhibits their activity, resulting in the inactivation of CTR1 88 . Moreover, the formation of embryos in Coffea canephorais inhibited in response to inhibitors of ET synthesis (Co 2+ ions) and ET activity (Ag + ions). By contrast, exogenous ET gas partially alleviates the effect of Co 2+ ions 86 . Furthermore, ET stimulates auxin biosynthesis and basipetal auxin transport toward the elongation zone, where it activates a local auxin response leading to inhibition of cell elongation 89 . Muday et al. suggested that ET modulates auxin synthesis, transport, and signaling through unique targets and responses in a range of tissues to fine tune seedling growth and development 90 . We found two annotated CTR1 unigenes that were down-regulated in group CK/EWAB, leading to the activation of ET signal transduction and enhancing the formation of ABs. However, six CTR1 unigenes were up-regulated in group CK/NPA or EWAB/NPA, indicating that ET signal transduction was suppressed because PAT was inhibited. Consequently, our results suggested that auxin and ethylene largely function synergistically to control ABs formation.
Brassinosteroids (BRs) are a group of naturally occurring novel steroidal plant hormones that regulate plant growth and development by inducing an array of physiological changes [91][92][93] . They are involved in a wide range of activities including seed germination, vascular formation, cell growth, reproductive growth, and production of flowers and fruit. BR treatment results in the elongation of excised segments of hypocotyls (seedling stems), principally because of directional cell expansion 94 . However, BRs can also promote leaf senescence 95 and negatively regulate phototropism 96 . In the present study, seven BR biosynthesis enzymes such as det2 and DWF1were down-regulated in group CK/EWAB, indicating that BRs plays a negative role in the regeneration process of G. lichenoidesin vitro. One BRs biosynthesis enzyme, SMT1, was up-regulated in response to NPA, which could be attributed to auxin-mediated biosynthesis of BR. Recently, Nemhauser et al. proposed the transcriptional regulation of ARFs by BRs 94 . The effects of auxin and BRs on a number of Aux/IAA genes could favor the formation of specific transcriptional complexes, promoting growth. Our results indicated that auxin could suppress BRs biosynthesis, whereas BRs might bind to inhibitors of auxin and then prevent PAT, which leaded to fail to induce ABs.
The role of salicylic acid (SA) is evident in seed germination, glycolysis, flowering, fruit yield, ion uptake and transport, photosynthetic rate, stomatal conductance, and transpiration 97 . Systemic-acquired resistance (SAR) is a broad-spectrum resistance in plants that involves the up-regulation of a battery of pathogenesis-related (PR) genes 98,99 and TGA, basic leucine zipper (bZIP) transcription factors, are required for the establishment of the SA-dependent SAR 100 . Several reports have described the effect of exogenous SA on enhancing somatic embryogenesis in plants 101,102 . In the present study, both TGA and PR1 transcripts were up-regulated to enhance the activity of defense programs in explants in response to wounding, osmotic stress, and inhibition of PAT.
The initiation of regeneration from callus may be related to the balance between different endogenous plant hormones 103,104 . Thus, we presumed that the mechanism of direct adventitious induction rather than callus formation of the G. lichenoides explants in vitro might be mediated by the cross-talk among different plant hormones, especially auxin. Antioxidant systems. Stress plays an important role during somatic embryogenesis similar to plant hormones 105 . Jin et al. revealed that somatic embryogenesis involves the activation of stress responses, which may regulate the balance between cell proliferation and differentiation 106 . Excised explants, which are used for callus or ABs induction in macroalgae, are subjected to sterilization procedures, cut into pieces, and cultured in an artificial nutritional environment. They therefore undergo considerable stress during this process, as indicated by the formation of reactive oxygen species (ROS). Recently, the morphological development of callus-like cells in S.japonica was closely linked to the activity of NADPH oxidase in the plasma membrane following ROS production 107 . Furthermore, ROS-induced signaling is entwined with plant hormonal responses. ET biosynthesis is an early O 3 response, and later, SA and ABA are produced 108 . In contrast to the increase in signaling via the stress hormones SA, ABA, and ET, ROS treatment causes auxin signaling to be transiently suppressed. The auxin transport inhibitor TIBA causes a similar expression profile to O 3 treatment regarding auxin-responsive genes, supporting a role for auxin transport in the regulation of apoplastic ROS responses 109 . In the present study, response to ROS was observed successfully in the GO enrichment analysis (Table 5).Genes related to ROS were mainly up-regulated by NPA, whereas SA, ABA, and BR-related genes were up-regulated. Our results suggest that ROS in combination with plant hormones, especially auxin, may provide a convenient reference to clarify the mechanism of direct formation of ABs instead of callus in the regeneration process of explants of macroalgae.
However, an inability of the explants to effectively scavenge ROS could lead to severe oxidative stress, causing damage to intracellular organelles and impairing the physiological processes related to cell regeneration. Thus, the antioxidant system, which mainly includes glutathione S-transferase (GST), NADPH, superoxide dismutase (SOD), and ascorbate peroxidase, is important for scavenging ROS during ABs formation.   110 . Similarly, we observed 10, 17, and 11 GST genes with marked expression changes in groups CK/EWAB, CK/ NPA, and EWAB/NPA, respectively (Table 3). More interestingly, 10 GST genes were obviously up-regulated in the EWAB/NPA group. This is consistent with the function of certain GSTs in tissue dedifferentiation, which is mediated by alterations in the redox status of the cell induced by changes in the endogenous levels of important plant growth hormones such as auxin 111 .

GO ID GO NAMES
We found that five, seven, and four SOD unigenes showed differential expression in groups CK/EWAB, CK/ NPA, and EWAB/NPA, respectively. Increased SOD activity enhances embryogenic cell differentiation during somatic embryogenesis 112 . Lin and Lai also demonstrated that a gradual increase in SOD activity during early somatic embryogenesis and late zygotic embryo development may play a key role in cell differentiation and the removal of ROS in longan somatic embryogenesis 113 . Therefore, variations in SOD activity might play an important role in blocking ROS production during the regeneration of G. lichenoides explants in vitro. ABs formation might be determined partially by the balance between the production and scavenging of ROS.
ROS is considered a factor controlling cell elongation in higher plants 114 . For example, rhizoid elongation in F. serratus may be caused by ROS produced by NADPH oxidase 115 . Moreover, NADPH oxidase shows Ca 2+ sensitivity 116 , suggesting that high Ca 2+ levels activate NADPH oxidase in the plasma membrane of callus-like cells and promote ROS production, resulting in cell wall loosening and cell elongation. Accumulation of ROS also occurs during multiple stages of root development, initially in response to transient cytosolic calcium increases as a result of a mechanical stimulus applied to a growing root tip 117 , and later in response to auxin during cell elongation and lateral root initiation 118 . In the present study, four, four, and two NADPH-related genes showed marked expression changes in the comparison groups CK/EWAB, CK/NPA, and EWAB/NPA, respectively. Meanwhile, one up-and fourdown-regulated genes were detected in the CK/EWAB group, while two transcripts, Unigene12429_All (11.5) and Unigene12914_All (11.3), were obviously up-regulated in response to NPA. In addition, seven calm transcripts (four up-and three down-regulated) were observed in the CK/EWAB group, whereas 14 transcripts (11 up-and three down-regulated) and 14 transcripts (nine up-and five down-regulated) were differentially expressed in CK/ Figure 7. IAA-protein similar as AUX/IAA was able to suppress ARFs expression. Thus, ARFs would become active after IAA-protein was identified and degraded by SCF-complex in the presence of auxin [67][68][69]71 , and then suppress downstream auxin genes related to the induction of ABs. In the presence of ABA, ABA response factor such as PYR/PYL/RCARs bind to and inhibit PP2Cs, which allows the accumulation of phosphorylated SnRK2s and the subsequent phosphorylation of ABFs, leading to the activation of downstream genes of the ABA signal transduction pathway 85 . ABA is considered an auxin inhibitor, and may reduce auxin to the optimum level required for regeneration of G. lichenoides explants in vitro; CTR1 encodes a negative regulator of the ET response pathway 89 , which might modulate auxin synthesis, transport, and signaling with unique targets and responses resulting in AB formation. Auxin could suppress BR biosynthesis, whereas BR binds to inhibitor proteins of auxin, resulting inthe prevention of PAT. ARFs could directly bind to the promoter of IPT and negatively regulate IPT expression 75 , affecting AB formation. The type-B ARRs act as transcription factors and their phosphorylation activates the transcription of the cytokinin-regulated genes, which negatively regulate auxin signaling 126 . ROS-induced signaling is entwined with plant hormonal responses. Systemicacquired resistance (SAR) is a broad-spectrum resistance in plants that involves the up-regulation of a battery of pathogenesis-related (PR) genes 100,101 and basic leucine zipper (bZIP) transcription factors TGA are required for the establishment of the SA-dependent plant defense response systemic-acquired resistance 102 . ET biosynthesis is an early O 3 response, and later, SA and ABA are produced 110 . In contrast to the increase in signaling via the stress hormones SA, ABA, and ET, ROS treatment caused auxin signaling to be transiently suppressed. The auxin transport inhibitor TIBA gave a similar expression profile toO 3 treatment for auxin-responsive genes, supporting a role for auxin transport in the regulation of apoplastic ROS responses 111 . The antioxidant system mainly includes GST, NADPH, and SOD, which are very important for scavenging the ROS production, which effectively promotes AB formation; auxin and brassinosteroid signaling and biosynthesis and auxin transport might be linked by an emerging upstream connection involving calcium-calmodulin and phosphoinositide signaling 96 during explant regeneration of G. lichenoides in vitro.
NPA and EWAB/NPA, respectively. It is possible that ROS production regulated by changes in calcium concentration and the activity of NADPH oxidase plays an important role in the regeneration process of explants of G. lichenoides.
Therefore, high levels of antioxidant activities relating to scavenging the ROS production might also determine effectively AB formation of explantsin G. lichenoides. Phosphatidylinositol (PI) signaling system. The PI metabolic pathway plays significant roles in many developmental processes, such as gene expression, hormone function, cell differentiation, cell responses to environmental factors, and interaction with other signaling pathways such as the calcium-related signaling network [119][120][121][122] . Connett et al. investigated the breakdown of membrane-bound phosphatidylinositol (PI) in homogenates of soybean (Glycine max) callus, which could be stimulated by Ca 2+ 123 . Similarly, calcium is essential for the regeneration of many plants 124 . In the present study, 13, 29, and 22 PI genes showed marked expression changes in the comparison groups CK/EWAB, CK/NPA, and EWAB/NPA, respectively, suggesting that PI signal transduction plays an important role during the formation of ABs. The addition of NPA up-regulatedcalm, PI5K, PIPK, PI(3)P5K, PLC and other PI transcripts. Further, phosphoinositide signaling has been reported to be triggered by auxin 125 . Hardtke et al. showed that auxin and BR signaling and biosynthesis and auxin transport might be linked by an emerging upstream connection involving calcium-calmodulin and phosphoinositide signaling 96 . Their findings together with our results suggest that PI signal transduction mediated by auxin plays an important role during formation of ABs.
In conclusion, the mechanism of direct AB induction rather than callus formation in G. lichenoides explants in vitro is complex and may be mediated by cross-talk among multiple elements, including plant hormones, antioxidant systems and other signal transduction pathways (Fig. 7). We investigated the mechanism of macroalgae tissue culture using the Illumina sequencing platform for the first time and identified important genomic factors that could help clarify the process of regeneration of macroalgae in the future.