Introduction

Inorganic sulfate is the most oxidized and dominant form of sulfur. It is an essential macronutrient for all organisms, which can be incorporated into amino acids, vitamins, cofactors, prosthetic groups, and a variety of primary and secondary metabolites1. Sulfur metabolism and the enzymes involved have been extensively studied in higher plants2. However, sulfur assimilation in algae remains largely unexplored but is presumably more similar to that in higher plants than in the prokaryotes, with a few exceptions on anion translocators3 and O-acetylserine (thiol)-lyase4. Briefly, inorganic sulfate transported into the cells is activated by ATP sulfurylase to form adenosine 5′-phosphosulfate (APS)5. Subsequently, APS can be either phosphorylated into 3′-phosphoadenosine 5′-phosphosulfate (PAPS) for the sulfation of metabolites, or converted into cysteine, methionine and glutathione for the synthesis of reduced sulfur compounds6.

Most photosynthetic organisms have limited sulfur storage in cells and depend solely on continuous uptake of sulfate from surrounding environment7. In freshwater and terrestrial environments, sulfate is often present at low concentration (10–50 μM) thus specific sulfate acclimation responses have been developed by different taxa. Transcriptomic responses towards sulfate deprivation had been extensively studied in photosynthetic organisms, such as Arabidopsis thaliana8,9,10,11, freshwater green alga Chlamydomonas reinhardtii7,12,13,14,15, cyanobacteria Synechocystis sp.16, and marine diatom Emiliania huxleyi17. The common responses during sulfate deprivation of these organisms include up-regulation of genes encoding sulfate transporters, increase in sulfate reducing capacity, synthesis of extracellular arylsulfatase and catabolism of sulfur-containing amino acids/secondary metabolites8.

In contrast with freshwater or terrestrial plants, seaweeds grow in marine waters with high sulfate concentrations ranging from 25 to 28 mM17 and hence sulfate may not be a limiting factor for cell growth and productivity in marine organisms. Many sulfated polysaccharides are synthesized by marine algae including agar18. Agar which is produced in the inner cell wall of red seaweeds may contribute to 6–71% of the total dry weight of seaweeds19. The backbone of agar is made up of two carbohydrate monomers D-galactose and L-galactose, with sulfate as the main side chain substitution20. Although the anionic matrix of sulfated polysaccharides provide heavily hydrated surface and higher flexibility enabling seaweeds to cope with dessication, strong ocean waves, osmotic stress, fluctuation in pH and temperature18, sulfation was found to lower agar gel strength21.

Seaweeds belonging to the genus Gracilaria contribute up to 60–80% of worldwide agar production22. Among these agarophytes, G. changii has a good agar yield (12–29% dry weight) and agar gel strength (180–563 g/cm2), while G. salicornia has a lower agar yield (9–14% dry weight) and poor gel strength (90–167 g/cm2)23,24,25,26. The molecular response of these red seaweeds to sulfate deprivation is largely unknown. We had previously shown that G. changii and G. salicornia collected during the rainy and dry seasons had minor increments in agar yields and gel strengths after sulfate deprivation for 5 days, compared to their respective controls, with a significant increase in agar yield for sulfate-deprivated G. salicornia collected during rainy season25. Although the effects of sulfate starvation on agar yield, gel strength, and other physical properties observed were not prominent, transcriptome sequencing of these samples may reveal dynamic changes in gene expression of these seaweeds to a sulfate free environment. In this study, we profiled and compared the transcriptomes of sulfate-deprivated G. changii and G. salicornia collected during rainy season whose yield and quality of agar were more affected by the sulfate deprivation treatment25. In addition, the transcriptomes of G. changii (with a higher agar yield and gel strength) and G. salicornia (with a lower agar yield and gel strength) were compared.

Results and Discussion

Effects of sulfate deprivation on the sulfate contents of seaweed and agar

Marine environment has stable sulfate content since 40–50 million years ago27. The production of sulfated cell wall polysaccharides (such as agar) could be an adaptation strategy that many seaweeds acquired over geological time to avoid sulfate toxicity in an environment with high, unchanging sulfate. We were interested to find out whether the sulfation level of agar decreased when the sulfate content in seawater was removed. Our data (Fig. 1) showed that sulfate deprivation generally reduced the sulfate content in the two Gracilaria species and their agars, but the differences were only statistically significant in G. changii seaweed sample (p = 0.0095) and G. salicornia agar sample (p = 0.0243). The effects of species and treatment on the sulfate content in the seaweed/agar samples were statistically significant (p < 0.001; Supplementary Table S1). These findings suggested that the sulfate distribution into different biological pathways (including the sulfation of agar) may vary according to the sulfate content in the marine environment and in different Gracilaria species.

Figure 1: Sulfate content in the thalli of G. changii and G. salicornia (a) and their agars (b).
figure 1

DW, dry weight. Error bars denote ± standard deviation (SD), n = 5. Asterisks (*) and (**) show significant difference (p < 0.05) and highly significant difference (p < 0.01), respectively, between control and sulfate-deprivated samples. Different letters indicate significant differences (p < 0.05) analyzed by Duncan Multiple Range Test.

Full size image

De novo transcriptome assembly and functional annotation

RNA-Sequencing (RNA-Seq) produced a total of 51.8 to 53.7 million raw reads for G. changii. The total sequencing throughput for G. salicornia was in the range of 35.9–37.2 million raw reads. A summary of the transcriptome sequencing data and assembly is shown in Supplementary Table S2. After the quality filtering process, a high percentage of clean reads (>99%) was obtained for both species, indicating that the sequencing reads were of good quality. Velvet assembly of paired-end reads and iterative clustering of the contigs produced 15,846 and 20,671 unigenes for G. changii and G. salicornia, respectively (Table 1). Although the sequencing depth for G. salicornia samples was lower compared to that of G. changii, the total number of unigenes obtained from the RNA-Seq data of G. salicornia was about 30% higher.

Table 1 Statistics of clustering and functional annotation of unigenes.
Full size table

The statistics on functional annotation of G. changii and G. salicornia unigenes are summarized in Table 1. A total of 10,341 unigenes (65.26%) from G. changii and 11,736 unigenes (56.78%) from G. salicornia have at least one match to the databases with an E-value < 10−5. However, the numbers of unigenes with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations are smaller compared to the numbers of unigenes that matched to the National Center for Biotechnology Information (NCBI) databases. G. salicornia has a higher number of unknown genes (43.22%) compared to G. changii (34.74%). The unannotated unigenes in these two Gracilaria species could be too short for similarity search or shared low similarities with annotated sequences in the databases.

Many annotated unigenes from G. changii (36.67%) and G. salicornia (38.62%) have high similarities to those in the red seaweed Chondrus crispus (Fig. 2a). However, some unigenes from these two Gracilaria species also share sequence similarity with the ciliated protozoa Tetrahymena thermophila (4.90–5.64%), Paramecium tetraurella (4.72–5.29%), Ichthyophthirius multifiliis (3.10–3.83%), Stylonychia lemnae (0.81–1.63%) and Oxytrichia trifallax (0.90–1.59%), which are taxonomically distant from red algal lineage. Since the transcriptomes were obtained from field samples, we do not exclude the possibility that some of these sequences may have originated from co-existing organisms/eukaryotes that were not removed by the cleaning process.

Figure 2: Functional annotation of unigenes from G. changii and G. salicornia.
figure 2

(a) Distribution of unigenes by species according to the best match of BLASTx results to non-redundant (NR) database at NCBI. (b) Gene Ontology (GO) annotation based on three main categories: cellular component, molecular function and biological process. (c) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways classification of the annotated unigenes.

Full size image

The distribution of GO categories for the two Gracilaria transcriptomes is very similar, with only minor differences observed (Fig. 2b). For example, unigenes assigned in two ontologies “chemoattractant activity” (GO: 0042056) and “protein tag” (GO: 0031386) were present in G. changii but absent in G. salicornia. Comparison of KEGG annotation between the two Gracilaria species showed that G. changii had a higher percentage of unigenes assigned to the category “cell growth and death” while G. salicornia had a higher percentage of unigenes assigned to “translation”, “transcription” and “transport and catabolism” categories (Fig. 2c). Reciprocal BLAST search of G. changii and G. salicornia unigenes revealed 4,018 homologous matches (Supplementary Table S3). We also identified a set of genes that are only expressed in G. changii (e.g., genes encoding galactose-2,6-sulfurylases I, polysaccharide pyruvyltransferase and thiosulfate sulfurtransferase rhodanese-like protein) or G. salicornia (e.g., genes for glutathione synthetase, glutamate dehydrogenase and carbon-sulfur lyase) (Supplementary Table S3). Some of these genes may be responsible for the major morphological and physiological differences (including the yield and quality of agar) between the two Gracilaria species.

Differential gene expression analysis

In this study, RNA pooled from the thalli of three individual plants was sequenced. Although the differences in gene expression of individual seaweeds may not be captured, the gene expression and biological noise from individual biological replicates were averaged by pooling biological samples28. Biological averaging method, in which individual biological replicates were replaced with pooled biological replicates, has been employed in many real time PCR, microarray and RNA-Seq experiments29,30,31. On the other hand, gene expression data derived from high throughput sequencing are highly replicable32 and the pooled samples generally average the technical variation33.

About 89–95% and 87–95% of the short reads from G. changii and G. salicornia, respectively, were mapped to their respective unigenes, indicating that the assemblies were of good quality. In this study, only differentially expressed genes (DEGs) with a reads per kilobase per million reads (RPKM) value equal or above 15 were reported because more than 50% of the DEGs which was below this threshold do not have sequence similarity with red algae (Supplementary Fig. S1). This was performed to reduce the possibility of reporting false positive genes originating from other co-existing organisms. The vast majority of genes in this study do not have high expression fold-change in the samples being compared thus we used 1.5 expression fold-change as a threshold value in defining DEGs. Many interesting and relevant genes may not be included in the list of DEGs by setting a higher threshold value. Furthermore, the same threshold value has been used by other studies for the identification of DEGs16,34,35.

In G. changii, 661 genes were up-regulated in sulfate-deprivated samples while 631 genes were down-regulated in sulfate-deprivated samples (Fig. 3). In G. salicornia, 1,787 genes were up-regulated in sulfate-deprivated samples while 1,652 were down-regulated in sulfate-deprivated samples (Fig. 3). The percentage of DEGs was higher in G. salicornia (16.64%) compared to that in G. changii (8.15%).

Figure 3: Number of up-regulated and down-regulated genes in G. changii and G. salicornia treated under sulfate deprivation for 5 days.
figure 3

Genes with a fold change of ≥1.5-fold in transcript abundance are considered as differentially expressed genes in sulfate-deprivated sample compared to the control sample.

Full size image

Verification of DEGs with quantitative reverse transcriptase polymerase chain reaction (qRT-PCR)

The expression of 15 genes from G. changii and G. salicornia, respectively, were verified by qRT-PCR. These randomly selected unigenes covered a wide range of RPKM values and log2Fold change (FC) values including two genes that were not differentially expressed. In general, the expression profiles of the unigenes in G. changii and G. salicornia, generated by RNA-Seq, corroborated with those obtained from qRT-PCR, with R2 values of 0.92 and 0.95 for G. changii and G. salicornia, respectively (Fig. 4). The results showed a high consistency between the results generated by RNA-Seq and qRT-PCR.

Figure 4
figure 4

Verification of differentially expressed genes from G. changii (a) and G. salicornia (b) using qRT-PCR. A total of 15 DEGs from G. changii and G. salicornia, respectively, were verified. FC, fold change.

Full size image

Responses of G. changii and G. salicornia to sulfate deprivation

To understand the responses of red seaweeds to sulfate deprivation, the expression of DEGs involved in sulfur acquisition and assimilation, metabolism of sulfur-containing amino acids, transport systems, and putative agar biosynthesis (Tables 2 and 3) were analysed and discussed.

Table 2 Expression level of selected differentially expressed genes from G. changii.
Full size table
Table 3 Expression level of selected differentially expressed genes from G. salicornia.
Full size table

Sulfate transport

Genes encoding ATP-binding cassette from G. changii (GC_457), S-type anion channel (GS_8662) and ABC transporters (GS_3784, GS_4151 and GS_660) from G. salicornia were found to be up-regulated in their sulfate-deprivated samples (Tables 2 and 3). It is unclear whether these ABC transporters participated in the transportation of exogenous sulfate into seaweeds. Two unigenes encoding putative sulfate permeases in the DASS/SLC13 permease family were identified in G. changii and G. salicornia, respectively, but were not being differentially expressed in response to sulfate deprivation. Meanwhile, the genes encoding high affinity H+/SO42− symporter (designed as SULTR in plants)2 were not identified in the transcriptomes of G. changii and G. salicornia. Our findings suggested that the sulfate transporters in Gracilaria species could be different from those reported in higher plants.

Sulfur acquisition and assimilation

Although the gene for ATP sulfurylase (ATPS) was identified in the Gracilaria species, it was not being differentially expressed in either species, suggesting that the transcription of this gene was not affected by sulfate-deprivation. APS is a branching point for sulfate assimilation into reductive or phosphorylative pathway2. Gene encoding the key enzyme in the reductive pathway, APS reductase was not found in the transcriptomes of both Gracilaria species (Fig. 5). It is unclear whether the seaweeds could obtain reduced sulfur through other enzymes/pathways, or the demand for the reduced sulfur was low during the experiments. The changes in transcript abundance for APS reductase in sulfate-deprivated E. huxleyi17 and Ch. reinhardtii15 were also non-significant compared to their respective controls. In a separate study36, the transcript encoding APS reductase from Ch. reinhardtii accumulated within 4 hours upon sulfur-deprivation before declining to a level comparable to that in the control, but its enzyme activity increased progressively over 24 hours, suggesting that this enzyme is regulated post-transcriptionally.

Figure 5: Effects of sulfate deprivation on the expression level of genes involved in sulfur metabolism.
figure 5

The sulfur metabolism pathway was modified from Takahashi et al.2. Each square and round box represents an unigene from G. changii (GC) and G. salicornia (GS), respectively. Yellow and white rectangles indicate transcripts that are found to be present or absent in the transcriptomes of both Gracilaria species, respectively, whereas grey and pink rectangles indicate transcripts that are only present in G. changii and G. salicornia, respectively. Dotted lines indicate enzymatic reaction encoded by unknown gene(s) from a gene family. All expression values are based on the biological averaging of pooled samples from three individual plants. APS, adenosine 5′-phosphosulfate; APK, APS kinase; APR, APS reductase; ATPS, ATP sulfurylase; BPNT, 3′(2′), 5′-bisphosphate nucleotidase; CBL, cystathionine beta-lyase; CGS, cystathionine gamma-lyase; GGCS, γ-glutamylcysteine synthetase; GSHS, glutathione synthetase; GST, gluthathione-S-transferase; MS, methionine synthase; OAST-TL, O-acetylserine (thiol)-lyase; PAPS, 3′-phosphoadenosine 5′-phosphosulfate; SAH, S-adenosylhomocysteine; SAHH, SAH hydrolase; SAM, S-adenosylmethionine; SAMS, SAM synthethase; SAT, serine acetyltransferase; SH, sulfohydrolase/sulfurylase; SiR, sulfite reductase; SO, sulfite oxidase; SOT, sulfotransferase; ST, sulfate transporter/permease.

Full size image

In this study, a gene encoding sulfite oxidase (GS_3043) was down-regulated in sulfate-deprivated G. salicornia (Table 3; Fig. 5). This gene was reported to be up-regulated in sulfate-deprivated Ch. reinhardtii15, regenerating free sulfate from excess sulfite recycled from cysteine and sulfur-containing amino acids. The down-regulation of this gene may suggest that regeneration of free sulfate from sulfite may reduce in sulfate-deprivated G. salicornia and sulfite is possibly an important intermediate for the biosynthesis of sulfide, cysteine and other sulfur-containing amino acids during sulfate deprivation.

Sulfate deprivation of G. changii decreased the transcript abundance of a gene encoding APS kinase (GC_2611) which is involved in sulfate phosphorylation (Table 2; Fig. 5). This suggested that sulfate deprivation may reduce PAPS production and subsequently the incorporation of sulfate into sulfur-containing metabolites (including agar), as shown by the slight decline in sulfate content in agar (Fig. 1b). However, APS kinase was not found to be differentially expressed in G. salicornia with significantly lower sulfate content in its agar (Fig. 1b) thus suggesting that there are other mechanism(s) that reduced the sulfate content in the agar of this species.

A decrease in the transcript abundance of APS kinase might lower the amount of cellular PAPS, which generally affects the transcriptional activity of sulfotransferases37. True enough, a gene encoding sulfotransferase 1C2A-like (GC_4742) from G. changii was down-regulated following sulfate deprivation (Table 2).

In G. salicornia, genes encoding heparan-sulfate-6-O-sulfotransferase (GS_1740) and carbohydrate sulfotransferase (GS_6282) were down-regulated in sulfate-deprivated sample (Table 3; Fig. 5). The decrease of these two transcripts in sulfate-deprivated G. salicornia coincided with a significant decrease of sulfate content in its agar (Fig. 1b). However, the involvement of these two carbohydrate sulfotransferases in the sulfation of galactose sugar requires further verification. On the other hand, the expressions of three sequences encoding glycolipid sulfotransferases (GS_8537, GS_8540 and GS_8541) were up-regulated in sulfate-deprivated G. salicornia suggesting that the demand for sulfation in glycolipid could be higher than that in the carbohydrate, when the seaweeds were deprived of sulfate sources.

Redistribution and recycling of sulfur

Extracellular arylsulfatase was found to be up-regulated in the green alga Ch. reinhardtii15,38 and marine diatom E. huxleyi17 to recycle sulfate from alternative sources during sulfate deprivation. In this study, gene encoding arylsulfatase was not found in the transcriptomes of Gracilaria species. However, the transcript abundance of a gene encoding galactose-2,6-sulfurylase I (GC_13944) was increased by 1.86 fold in sulfate-deprivated G. changii (Table 2; Fig. 5). The changes in gene expression coincide with the minor reduction of sulfate content in the agar samples of sulfate-deprivated G. changii, suggesting the putative role of galactose-2,6-sulfurylase I in recycling sulfur from sulfated galactans. Galactose-2,6-sulfurylase is a novel enzyme which has only been reported in the red alga C. crispus, and was proposed to catalyse the conversion of D-galactose-6-sulfate into 3,6-anhydro-D-galactose39. However, the transcripts encoding galactose-2,6-sulfurylase I were not found in the transcriptome of G. salicornia. It is unknown whether a higher gel strength in G. changii agar compared to that in the G. salicornia was attributed to the gene expression and activity of the species-specific galactose-2,6-sulfurylase I. Genes encoding rhodanese and cysteine/taurine dioxygenase had been reported to be involved in the recycling of internal sulfur15, but they were not differentially expressed in both Gracilaria species.

Metabolism of glutathione and methionine

Overall, the genes involved in the metabolism of glutathione and methionine showed different responses in sulfate-deprivated G. changii and G. salicornia. The gene encoding γ-glutamylcysteine synthetase (GC_2731), which is responsible for the biosynthesis of L-γ-glutamylcysteine from L-cysteine and L-glutamate, was down-regulated in sulfate-deprivated G. changii (Table 2; Fig. 5). Meanwhile, the enzyme which catalyses the biosynthesis of glutathione from L-γ-glutamylcysteine was not found in the transcriptome of G. changii (Fig. 5). Our results suggested that the glutathione synthesis in G. changii could be depressed by sulfate deprivation in G. changii. The genes encoding for γ-glutamylcysteine and glutathionine synthetase were not differentially expressed in sulfate-deprivated G. salicornia (Fig. 5).

A gene encoding peptide methionine sulfoxide reductase (GC_2663) was down-regulated in sulfate-deprivated G. changii (Table 2). This result was in contrast with that in G. salicornia, in which the unigene encoding peptide methionine sulfoxide reductase (GS_3594) was up-regulated in response to sulfate deprivation (Table 3). The increase in the abundance of transcript encoding peptide methionine sulfoxide reductase (which catalyses the regeneration of methionine from methionine sulfoxide) and methionine aminopeptidase (GS_14648) in sulfate-deprivated G. salicornia suggest that methionine might be recycled from proteins during sulfate deprivation.

Meanwhile, a gene encoding methionine synthase (GS_4274) was down-regulated in sulfate-deprivated G. salicornia while the transcript encoding S-adenosylmethionine (SAM) synthethase (GS_4407) was accumulated in its sulfate-deprivated samples (Table 3). This suggested that under sulfate deprivation, methionine was recycled from proteins and converted to SAM which is important for methylation process. The response was similar to that of A. thaliana treated under sulfate deprivation9, but is in contrast with those in sulfate-deprivated G. changii, Ch. reinhardtii15 and E. huxleyi17.

Biosynthesis of dimethylsulfoniopropionate (DMSP)

DMSP is an osmolyte and cryoprotectant accumulated in many algal species in response to abiotic stress and nutrient deprivation40,41,42. DMSP biosynthesis is dependent on APS reductase activity which regulates the synthesis of its precursor, methionine, through the reductive pathway of sulfur assimilation43. Genes involved in the biosynthesis of DMSP have been identified and characterized, including 2-oxoglutarate-dependent aminotransferase, NADPH-linked reductase and SAM-dependent methyltransferase44. These genes were not being found (2-oxoglutarate-dependent aminotransferase and NADPH-linked reductase) or being differentially expressed (SAM-dependent methyltransferase) in G. changii and G. salicornia possibly due to a relatively low content of DMSP in Gracilaria spp. compared to that in green algae and phytoplankton45,46.

Putative agar biosynthesis

The biosynthetic pathway of agar is less explored and the genes involved have yet been elucidated47. In this study, we identified a number of putative genes which could be involved in the biosynthesis of agar precursors from G. changii and G. salicornia (Fig. 6). DEGs related to the biosynthesis of fructose-6-phosphate (a central metabolite used for the synthesis of agar monomer), elongation of polysaccharide chain and modification of its side chain substitutions were found in sulfate-deprivated samples (Fig. 6). The agar yield of sulfate-deprivated G. salicornia increased significantly25, coincides with the increase in the transcript abundance of genes encoding galactose-1-phosphate uridylyltransferase and GDP-mannose-3′,5′-epimerase (Table 3; Fig. 6), suggesting their putative roles in regulating the biosynthesis of agar precursors. The gene expression level of these two genes had been shown to correlate with the agar content in Gracilaria and Gracilariopsis samples24,48,49.

Figure 6: Effects of sulfate deprivation on the expression level of putative genes involved in agar biosynthesis.
figure 6

The putative agar biosynthetic pathway was modified from Lee et al.47. Each square and round box represents an unigene from G. changii (GC) and G. salicornia (GS), respectively. Yellow and white rectangles indicate transcripts that are found to be present or absent in the transcriptomes of both Gracilaria species, respectively, whereas grey rectangle indicates transcripts that are present in G. changii only. Dotted lines indicate enzymatic reaction encoded by unknown gene(s) from a gene family. All expression values are based on the biological averaging of pooled samples from three individual plants. FBPA, fructose-1,6-bisphosphate aldolase; FBA, fructose-1,6- bisphosphatase; GALT, galactose-1-phosphate uridylyltransferase; GME, GDP-mannose-3′,5′-epimerase; GMP, GDP-mannose pyrophosphorylase; GPI; glucose-6-phosphate isomerase; GTs, glycosyltransferases; MPI; mannose-6-phosphate isomerase; MTs, methyltransferases; PFK, phosphofructokinase; PGM, phosphoglucomutase; PMM, phosphomannose mutase; PTs, pyruvyltransferases; STs, sulfotransferases; UGP, UTP-glucose-1-phosphate uridylyltransferase.

Full size image

In addition, most of the genes encoding galactosyltransferases (GTs), sulfotransferases (STs), methyltransferases (MTs) and pyruvyltransferase (PTs) were down-regulated in G. changii in response to sulfate deprivation (Fig. 6). In G. salicornia, some of the genes encoding GTs and STs showed different expression patterns from those in G. changii, they were being up-regulated in response to sulfate deprivation (Fig. 6). However, the gene encoding PT was not being identified in sulfate-deprivated G. salicornia. The differential expression of these genes in G. changii and G. salicornia may lead to higher gel strength in the former, as sulfation and methylation on the agar were reported to produce a weaker gel50,51. In most organisms, glycosyltransferases (including glucosyltransferases and galactosyltransferases) represent a large diverse group of enzymes with more than 30 gene families52. To our knowledge, glycosyltransferase which is responsible for agar polymerization has not been identified. In this study, we have identified a few DEGs encoding glycosyltransferase. Two of them (GS_4319 and GS_6663) were up-regulated in sulfate-deprivated G. salicornia sample (Table 3 and Fig. 6), which had an increase in agar yield upon sulfate deprivation25.

Other pathways

The genes involved in photosynthesis, such as those encoding phycobilisome linker polypeptide, photosystem II Q and light harvesting protein, were differentially expressed in sulfate-deprivated G. changii and G. salicornia, with different expression patterns (Tables 2 and 3). In Ch. reinhardtii, the genes encoding different light harvesting proteins have different expression patterns under sulfate deprivation15, suggesting that these proteins may respond differently to sulfate deprivation.

Sulfate deprivation regulated a number of genes involved in carbon metabolism, suggesting a tight interconnection between carbon and sulfur metabolisms. For example, genes encoding transaldolase, fructose-1,6-bisphosphate aldolase, and glycoside hydrolases were differentially expressed in sulfate-deprivated G. changii and G. salicornia, but the expression patterns were species-dependent (some of these genes were up-regulated in one species but down-regulated in another species, and vice versa) (Tables 2 and 3). There are a number of genes that were only differentially expressed in sulfate-deprivated G. changii or G. salicornia, such as genes encoding glucose-6-phosphate isomerase (GC_1234), sedopheptulose-1,7-bisphosphatase (GC_7206), glyceraldehyde-3-phosphate dehydrogenase (GS_5744 and GS_6394), phosphoglycerate kinase (GS_5992) and 6-phosphogluconate dehydrogenase (GS_1698) (Tables 2 and 3).

The genes encoding superoxide dismutase (GC_2437) and glutathione S-transferase (GS_5546) from G. changii and G. salicornia, respectively; were up-regulated in sulfate-deprivated sample (Tables 2 and 3), similar to that observed in sulfate-deprivated A. thaliana9. Meanwhile, a few transcripts encoding haloperoxidase-like proteins, including vanadium-dependent bromoperoxidases, were decreased by sulfate deprivation in both Gracilaria species. Although the biological role of haloperoxidase-like protein remains unclear, these genes have been shown to be differentially expressed in red seaweeds under various abiotic stresses, such as light, salinity and temperature53,54,55. It is noteworthy that the increase of these stress-related transcripts could be a general stress response.

Conclusions

The findings of this study offer a snapshot of the global transcriptional responses upon sulfate deprivation in two Gracilaria species, i.e. G. changii and G. salicornia. Among the DEGs were genes involved in putative agar biosynthesis, sulfur acquisition and assimilation, metabolism of sulfur-containing amino acids, transport systems and oxidative stress. These two species responded differently to sulfate deprivation, possibly due to differences in their cellular sulfate storage and demand as well as synthesis of sulfated polysaccharides, thus affecting the expression of genes encoding enzymes related to sulfate phosphorylation, removal of sulfate, methionine and SAM cycle, and recycling of sulfur. In addition, sulfate deprivation responses of Gracilaria species were found to be different from those in Ch. reinhardtii, E. huxylei and A. thaliana, suggesting a unique acclimation response in marine macroalgae. However, these differences could also be contributed by gene expression profiling methods used, physiological states of the sample, experimental design, severity of sulfate deprivation, and duration of treatment.

Materials and Methods

Seaweed materials and sulfate deprivation treatments

Collection of seaweed materials and sulfate deprivation treatment were performed as described by Lee et al.25. Briefly, Gracilaria changii (Xia et Abbott) Abbott, Zhang et Xia and G. salicornia (C. Agardh) Dawson without any cystocarpic structures were collected from the mangrove swamp at Morib, Selangor, Malaysia (02°45.878′ N; 101°25.976′ E) in October 2012 (rainy season). The seaweeds were washed from visible mud, epiphytes and epibionts, before acclimating in aerated aquaria containing 30 part per thousand (ppt) artificial seawater (ASW: Instant Ocean Synthetic Sea Salt, Aquarium Systems Inc., Mentor, OH, USA) at pH 7 and 25 °C.

Sulfate deprivation treatment was conducted on G. changii (CC, control sample; CT, treated sample) and G. salicornia (SC, control sample; ST, treated sample). The thalli of G. changii and G. salicornia were treated for 5 days in aerated aquaria containing sulfate-free ASW (450 mM NaCl, 370 mM KCl, 9 mM CaCl2.2H2O, 49 mM MgCl2.6H2O and 2 mM NaHCO3), with a control in normal ASW (450 mM NaCl, 370 mM KCl, 9 mM CaCl2.2H2O, 23 mM MgCl2.6H2O, 26 mM MgSO4.7H2O and 2 mM NaHCO3). The seaweeds were treated for five days because we could not maintain the seaweeds in a healthy condition (without necrosis) in the laboratory for a longer period. Growing in a marine environment with high and stable sulfate content, the seaweeds may have accumulated sufficient sulfate that may last for sometime. Due to that, we have chosen a sulfate-free treatment for a short duration (5 days). The treatment was performed under laboratory condition with natural photoperiod (12:12 h light:dark cycle, 8–10 μmol photons m−2s−1 during the light cycle), while the temperature (25 °C), pH (7.8) and salinity of seawater (35 ppt) were kept constant throughout the experiment. All the samples were collected from the field or harvested during day time (1200–1400 local time) to minimize the diurnal patterns among the samples.

Measurement of sulfate content

The total sulfate content in seaweed and agar (that were pooled from thalli of a few individual plants, respectively) was determined using a modified BaCl2-gelatin turbidimetric method56. Approximately 15 and 40 mg of dried seaweed and agar samples, respectively, were hydrolysed in 0.5 M HCl (in a ratio of 1:150 w/v) at 110 °C for 10 hours. The solution was cooled to room temperature and filtered through Whatman No. 1 filter paper (Whatman, Hillsboro, OR). The seaweed or agar hydrolysate (0.2 ml each) was added to 3.8 ml of 3% (w/v) trichloroacetic acid and 1 ml of BaCl2-gelatin reagent (40 mM BaCl2, 0.3% w/v gelatin), mixed vigorously and incubated at room temperature for 30 min. The sample was replaced with 0.5 M HCl in the blank, while a “control” which consisted of gelatin solution without BaCl2 was prepared for each sample to eliminate artifacts contributed by ultraviolet light-absorbing materials produced during acid hydrolysis. Spectrophotometric measurement for each sample and its “control” was performed against the blank solutions at an absorbance of 360 nm. The optical reading obtained from the “control” was subtracted from that of the seaweed or agar hydrolysate before the sulfate content was calculated from a standard curve prepared with 1.15–5.74 mM K2SO4, per mg of seaweed or agar.

Statistical analyses

All statistical analyses were performed using SAS statistical software version 9.3 (SAS Institute Inc., Cary, North Carolina, USA). Two-way ANOVA test was used to determine the significant effects of species and treatment, and their interactions on the sulfate content in seaweed and agar. Statistical difference in sulfate content between the control and sulfate-deprivated samples (from five replicates, respectively) was determined with Student’s t-test while Duncan Multiple Range Test (DMRT) was used to determine the significant differences among the experimental groups (e.g. CC, CT, SC and ST).

Isolation and quality assessment of RNA

Total RNA was isolated from 10 g of each Gracilaria sample (i.e., CC, CT, SC and ST) which was pooled from the thalli of three individual plants each, as described by Chan et al.57. Genomic DNA was removed using RNase-free DNase I (New England Biolabs, Beverly, MA) following the manufacturer’s instructions. The purity and quality of the RNA were evaluated using Nanodrop spectrophotometer (Implen UK Ltd, Essex, UK) and Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The RNA samples which met the following criteria were processed for RNA-Seq and qRT-PCR analyses: absorbance ratios A260/280 and A260/230 that were between 1.9–2.1 and more than 2, respectively, 25 S rRNA:18 S rRNA ratio near to 2:1 and RNA Integrity Number (RIN) more than 6.5.

RNA-Seq, de novo assembly and annotation of the transcriptomes

Massive paired-end RNA Seq was performed on each Gracilaria sample using Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA) at a read length of 90 bp. The RNA-Seq data were deposited in European Nucleotide Archive (ENA) under the accession number PRJEB13899. The raw reads obtained from RNA-Seq were analysed with FastX toolkit v0.0.13.2 (http://hannonlab.cshl.edu/fastx_toolkit/) to remove low-quality reads with Phred quality scores <20, adaptor and ambiguous bases ‘N’. De novo assembly of clean reads obtained from the four samples was performed with Velvet (v1.2.08)58 with default parameters and optimized k-mers at 67, 69, 65, and 65 for CC, CT, SC, and ST, respectively. Assembled sequences that were less than 100 nucleotides were discarded. Iterative clustering of contigs from control and treated samples was performed with The Gene Index Clustering Tool (TGICL)59 to generate two sets of non-redundant contigs (unigenes) for G. changii and G. salicornia, respectively.

Sequence based alignments were performed against the following public databases: non-redundant protein sequences (NR), UniProtKB/SwissProt (Swiss-Prot), GO and KEGG. BLASTx program60 with default settings and a cut-off E-value of 10−5 was used to match the nucleotide sequences to those in the NR and SwissProt databases at the NCBI. The GO annotation was performed with BLAST2GO61 based on the most significant SwissProt matches and the results were visualized using Web Gene Ontology Annotation Plot (WEGO)62. KEGG annotation was performed with the KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/kegg/kaas/)63. The most significant and informative matches were used to assign putative function to each unigene. Reciprocal BLASTn64 was conducted on the two sets of unigene with a cut-off E-value of 10−30.

Identification of differentially expressed genes (DEGs)

The paired-end RNA-Seq reads from each sample were mapped to respective transcriptomes with Bowtie265. The number of reads mapped to each unigene was converted into RPKM66, to normalise the differences in transcript length and sequencing depth. DEGs between control and sulfate-deprivated samples of G. changii and G. salicornia, respectively, were identified by DEGseq67, an R package which uses the MA-plot-based method with Random Sampling model (MARS). In the absence of replicates, DEGseq treated the counts of reads that were mapped to a gene from two samples (control and sulfate-deprivated samples in this study) as independent values and tested them with the following hypotheses, H0: p1 = p2 = p versus H1: p1≠p2; p1 and p2 denote the probability that a read obtained from sample 1 (e.g. control) and 2 (e.g. sulfate-deprivated), respectively. The two estimates generated from above were then used for the calculation of Z-score which can be converted to two-sided p-value67. Transcripts with RPKM read ≥15 in at least one of the samples, p-value less than 0.001 and fold change ≥1.5-fold were defined as DEGs (Supplementary Table S4).

Verification of DEGs using qRT-PCR

Verification of DEGs was performed on the same RNA samples for RNA-Seq experiment. Affinity Script QPCR cDNA Synthesis Kit (Agilent Technologies, USA) was used to reverse transcribe 2.5 μg of total RNA from each sample into the first strand cDNA. The DEGs and the primers used for qRT-PCR analysis are listed in Supplementary Table S5. Primer3 software ver. 0.4.0 (http://frodo.wi.mit.edu)68 was used to design primers that are specific to the 3′-untranslated region flanking the open reading frame, to avoid amplification from possible isoforms/gene family members. Melting curve and standard curve (with serial dilutions of cDNA samples) were generated to check for amplification of single PCR product, absence of primer dimers and PCR efficiency (90–105%).

All qRT-PCR analyses were performed using Brilliant III Ultra-Fast SYBR Green QPCR Master Mix (Agilent Technologies, USA) and iQTM5 real time detection system (Bio-Rad, USA). The reaction mixture (20 μl) consisted of 100 ng of first-strand cDNA template, 1X Brilliant III Ultra-Fast SYBR Green QPCR master mix, 200–400 nM of each respective primers and molecular grade H2O (bioWORLD, USA). The PCR amplification was performed with the following conditions: initial denaturation at 95 °C for 3 min, followed by 40 cycles of 95 °C for 5 sec and 60 °C for 10 sec; and a melt curve analysis with temperature increasing from 55 °C to 95 °C with an increment of 0.5 °C/10 sec. A control sample without template was included in all experiments to rule out the presence of contaminants and primer-dimer.

Two internal reference genes for G. changii (encoding histone and cytosolic fructose-1,6-biphosphatase, respectively) and G. salicornia (encoding 40S ribosomal protein S15a and tRNA-dihydrouridine synthase-1-like, respectively) were selected based on their stability ranking and pairwise variation as analysed by geNorm version 3.569. All samples in the qRT-PCR analysis were performed in three technical replicates. The raw Ct values obtained were analysed using 2−ΔΔCt method70,71, and normalized against the expression of internal reference genes. The normalized expression level in all samples was expressed as relative fold change to that in the control samples.

Additional Information

How to cite this article: Lee, W.-K. et al. Transcriptome profiling of sulfate deprivation responses in two agarophytes Gracilaria changii and Gracilaria salicornia (Rhodophyta). Sci. Rep. 7, 46563; doi: 10.1038/srep46563 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.