Sequencing data interpretation
Totally, 30 kelp samples were subjected to transcriptome sequencing and data analysis (Table 1; Additional file 1: Figure S1). All the RNA integrity numbers(RINs) were between 7.2-9.2, which showed that all the samples were qualified for deep sequencing (Additional file 2: Table S1). High-throughput sequencing generated 45.69-84.83 million of 150-bp paired-end reads in each library (Table 2). After data filtering, 1,751,262,386 high quality reads were produced (98.20% of clean reads) with an average Q30% > 96.00%. rRNA removed reads were mapped with our previous completed S. japonica genome (NCBI: MEHQ00000000) with a mapping ratio of c. 80%, except for Ap1-3 sample which had a rather low ratio of 50.46% (Table 2; Additional file 2: Table S1). Assembled transcriptomes were annotated according to 24,419 reference genes. In total, 20,940 known genes (ID starts with “GENE_”) and 4,264 novel genes (ID starts with “XLOC_”) were obtained. There are 1,957 novel genes showing high identities with sequences from E. siliculosus, among which 34.3% were annotated as unknown or hypothetical proteins. KEGG pathway annotation showed two main categories of “Metabolism” and “Genetic information processing” (Additional file 3: Figure S2). Statistics of transcriptomes sequencing output are in Additional file 2: Table S1.
Table 1 The collection background of S. japonica samples
Collection date
|
Sample ID
|
Tissue site
|
Replicates
|
Seawater temperature
|
22nd January
|
JaB
|
Basal
|
JaB-1, JaB-2, JaB-3
|
5.2 ℃
|
4th March
|
MhB
|
Basal
|
MhB-1, MhB-2, MhB-3
|
4.6 ℃
|
10th April
|
ApB
|
Basal
|
ApB-1, ApB-2, ApB-3
|
5.6 ℃
|
Ap1
|
1/3
|
Ap1-1, Ap1-2, Ap1-3
|
Ap2
|
2/3
|
Ap2-1, Ap2-2, Ap2-3
|
ApD
|
Distal
|
ApD-1, ApD-2, ApD-3
|
10th May
|
MyB
|
Basal
|
MyB-1, MyB-2, MyB-3
|
9.0 ℃
|
MyD
|
Distal
|
MyD-1, MyD-2, MyD-3
|
16th June
|
JuB
|
Basal
|
JuB-1, JuB-2, JuB-3
|
13.2 ℃
|
JuD
|
Distal
|
JuD-1, JuD-2, JuD-3
|
Table 2 The output and quality control of the RNA-Seq data
|
Maximum
|
Minimum
|
Average
|
Clean data (bp)
|
12,724,976,400
|
6,852,831,900
|
8,916,922,600
|
HQ clean data (bp)
|
12,214,001,205
|
6,606,166,213
|
8,595,871,647
|
Q30 (%)
|
96.69%
|
95.05%
|
96.00%
|
GC (%)
|
57.26%
|
55.72%
|
56.52%
|
Clean reads No.
|
84,833,176
|
45,685,546
|
594,461,501
|
HQ clean reads No.
|
82,982,428
|
44,857,256
|
58,375,413
|
% HQ clean reads
|
98.47%
|
97.82%
|
98.21%
|
Mapped reads No.
|
65,017,939
|
32,881,766
|
45,563,696
|
Mapping ratio
|
82.76%
|
50.46%
|
80.33%
|
Total genes
|
25,204
|
Trend analysis and functional enrichment of differentially expressed genes (DEGs)
DEGs in different kelp developmental periods and tissue portions were clustered into 29 and 25 profiles, respectively (Additional file 4: Figure. S3). We selected two representative profiles: profile29 (985 genes) with increasing DEGs expression and profile0 (1319 genes) with decreasing trend from January to June. Twenty DEGs encoding ribosomal proteins (RPs) were enriched in the profile29 (Q < 0.05) (Additional file 5: Table S2). “Oxidative phosphorylation”, “photosynthesis-antenna proteins”, “photosynthesis and carbon fixation pathways” were enriched in profile0, with a decreasing trend from juvenile to mature sporophytes (Q < 0.05) (Additional file 5: Table S2). While for DEGs in different tissues, we selected 3 profiles with increase pattern [profiles 25 (367 genes), 22 (322 genes) and 16 (303 genes)] and 2 with decrease pattern [profiles 9 (1619 genes) and 0 (1359 genes)] (Additional file 5: Table S2). DEGs in “ABC transporters and RNA transport pathways” were highly enriched (p < 0.05) and gradually increased from basal to distal blade. “Polyunsaturated fatty acids (e.g. arachidonic acid and linoleic acid) metabolisms” were also enriched in this pattern (Q < 0.05). Profile9 and profile0 enriched “photosynthesis-antenna protein”, “secondary metabolites biosynthesis and thiamine (VB1) metabolism”, and complex metabolic pathways known as “microbial metabolism in diverse environments” under the KEGG nomenclature (Q < 0.05). Notably, “imm upregulated 3” gene family were highly expressed in juvenile sporophytes, with 28 genes enriched in profile0 (Additional file 6: Table S3). In the basal blade, 32 “imm upregulated 3” genes were enriched (Additional file 6: Table S3).We analyzed the expression profiles of genes related with energy-producing metabolisms (Table 3). Although there is no obvious expression tendency from January to June, genes encoding key enzymes in TCA cycle and glyoxylate cycle were highly expressed in the distal blade of kelp, whereas those crucial genes in glycolysis, gluconeogenesis and oxidative phosphorylation were highly expressed in basal blade compared with distal blade.
Table 3 The expression levels of key enzymes involved in energy-producing metabolisms between basal and distal blade of kelp
Pathway
|
Rate-limiting enzyme
|
Gene ID
|
Fold change
|
ApD vs ApB
|
MyD vs MyB
|
JuD vs JuB
|
|
TCA cycle
|
Citrate synthase
|
GENE_007839
|
-2.29
|
-2.11
|
-11.25
|
Oxoglutarate dehydrogense
|
GENE_002415
|
-1.78
|
-2.93
|
-3.95
|
GENE_002283
|
-1.75
|
-1.36
|
-1.86
|
XLOC_013477
|
1.38
|
-1.73
|
-2.86
|
Glyoxylate cycle
|
Isocitrate lyase
|
GENE_025653
|
-7.81
|
-46.48
|
-19.09
|
Malate synthase
|
GENE_024070
|
-2.58
|
-9.05
|
-14.19
|
Glycolysis
|
Glucokinase
|
GENE_019969
|
3.24
|
12.18
|
7.35
|
GENE_027705
|
7.38
|
9.60
|
5.18
|
Pyrophosphate- phosphofructose kinase
|
GENE_017153
|
7.63
|
13.50
|
5.19
|
XLOC_033414
|
3.65
|
6.02
|
4.30
|
GENE_023282
|
1.18
|
1.34
|
0.84
|
Pyruvate kinase
|
GENE_006234
|
1.13
|
1.63
|
1.68
|
Gluconeogenesis
|
Fructose-1,6-bisphosphatase
|
GENE_010339
|
2.90
|
2.40
|
2.85
|
GENE_004642
|
1.60
|
1.14
|
1.27
|
Phosphoenolpyruvate carboxykinase (ATP)
|
XLOC_018117
|
46.29
|
38.48
|
18.73
|
GENE_015129
|
28.30
|
30.78
|
22.01
|
Oxidative phosphorylation
|
ATP synthase
|
GENE_008643
|
1.75
|
2.32
|
0.69
|
GENE_004216
|
4.77
|
9.00
|
3.29
|
GENE_020559
|
1.85
|
2.24
|
0.94
|
GENE_025468
|
2.63
|
2.36
|
1.47
|
GENE_024356
|
3.00
|
1.44
|
0.57
|
Alginate and mannitol content variations during the kelp developmental periods
Mannitol and alginate contents in S. japonica were detected at different periods from January to June (Fig. 1; Additional file 7: Figure S4). Average content of alginate in individual kelp collected from each month was c. 30% without significant difference from month to month. However, mannitol content gradually increased from 1.56% (in March) to 9.63% (in June), a 6.18-fold difference between the mature and the juvenile sporophytes (Fig. 1). The contents of these two metabolites varied in different tissue portions. For instance in April, alginate content was up-regulated from basal blade (22.71%) to distal blade (35.06%) whereas mannitol was down-regulated from 7.04% to 0.58% (Fig. 1; Additional file 7: Figure S4). Alginate content was constantly higher in the distal blade than in the basal blade (1.16-1.54-fold change). On the contrary, the basal blade contained more mannitol than the distal blade, especially in mature sporophytes in June (25.04-fold). In general, mannitol was up-regulated from juvenile to mature sporophytes and the variations of mannitol and alginate showed opposite patterns from basal to distal blade (Fig. 1).
Fig. 1 Contents of alginate and mannitol detected in S. japonica samples collected from different developmental stages and tissue parts. (a) Content variations of alginate from March to June. (b) Content variations of mannitol from March to June.
Identification of alginate/mannitol-related genes and their transcriptional profiles
Based on the annotation of our previously sequenced S. japonica genome, we screened 134 genes encoding enzymes catalyzing alginate biosynthesis, including 3 MPI, 2 PMM, 3 GMD, 1 GT2 and 125 MC5E genes. We identified 6 genes encoding enzymes involved in the mannitol metabolism, of which 3 sequences (M1PDH and M1Pase) were for mannitol biosynthesis and the rest 3 sequences (M2DH and FK) were for degradation (Additional file 8: Table S4).Transcriptional profiles of these genes from RNA-Seq data are shown in Fig. 2. Transcripts of MPI1 (GENE_021848), MPI3 (GENE_013986), PMM1 (GENE_007314) and GMD3 (GENE_022063) were down-regulated from January to June, whereas other genes did not exhibit specific transcriptional pattern, especially for MC5E gene family which contains 125 genes with very diverse transcriptional profiles. Fig. 2a shows the expression levels of 3 representative MC5Es, among which MC5E70 (GENE_007019) and MC5E122 (XLOC_006798) were highly expressed in the distal blade. PMM2 (GENE_006655) and GMD3(GENE_022063) exhibited opposite expression patterns: PMM2 decreased and GMD3 increased from basal blade to distal blade (Fig. 2a). M1PDH1 (GENE_003979) and M1Pase (XLOC_010181) expression were gradually down-regulated and M2DH (GENE_006978 and GENE_006979) and FK (GENE_018623) were remarkably up-regulated from basal to distal blade. Expression levels of all the 6 genes in mannitol cycle were higher in juvenile sporophytes than in later stages (Fig. 2b).
Fig. 2 Transcriptional patterns of the genes involved in alginate and mannitol metabolism. The order of each row is: JaB, MhB, ApB, MyB, JuB, ApB, Ap1, Ap2, ApD. MPI1: GENE_021848; MPI2: GENE_013980; MPI3: GENE_013986; PMM1: GENE_007314; PMM2: GENE_006655; GMD1: GENE_022030; GMD2: GENE_008524; GMD3: GENE_022063; GT2: GENE_006305; MC5E1: GENE_007233; MC5E70: GENE_007019; MC5E122: XLOC_006798; M1PDH1: GENE_011959; M1PDH2: GENE_003979; M1Pase: XLOC_010181; M2DH1: GENE_006978; M2DH2: GENE_006979; FK: GENE_018623.
Identification of alginate-/mannitol- co-expressed genes and pathways via module-trait correlations
WGCNA analysis resulted in 22 distinct modules. Additional file 9: Figure S5 shows the hierarchical cluster tree for modules of co-expressed genes with each branch constituting a module and each leaf as one gene. Additional file 10: Table S5 lists the number of genes clustered in each module. Alginate and mannitol module-trait correlation analysis was conducted based on WGCNA data. Fig. 3a shows correlations from -1 (green) to 1 (red), which revealed that the “brown4” module was closely related to alginate content (r = 0.84, p = 6 × 10-9) and “black” module was highly correlated to mannitol (r = 0.80, p = 1 × 10-7). Fig. 3a shows the opposite module-trait correlation pattern between alginate and mannitol biosynthesis: positive correlated modules for alginate concentrated mannitol-negatively correlated modules and vice versa. Fig. 3b shows the top 2 modules for each correlation analysis: “brown4” and “darkgreen” were correlated with alginate whereas “black” and “darkslateblue” were correlated with mannitol.
Although the “brown4” and “black” modules were highly correlated with alginate and mannitol contents, none of their biosynthetic genes were found in these two modules (Additional file 8: Table S4). Genes involved in alginate biosynthesis: MPI2 (GENE_013980), PMM1 (GENE_007314), GMD3 (GENE_022063), MC5E70 (GENE_007019) and MC5E122 (XLOC_006798) appeared in “darkorange” and “mediumpurple3” modules (Fig. 3a). The two modules clustered genes with higher expression levels in the distal blade than those in the basal blade of kelp, as the “brown4” module (Additional file 11: Figure S6a,b). M1PDH1 and M1Pase (for mannitol synthesis) fell into the module of “greenyellow” (Fig. 3a), in which gene expression levels were relatively higher in basal blade (Figs. 3b; Additional file 11: Figure S6c). However, their transcripts were highly up-regulated in juvenile sporophytes, but mannitol content was higher in adult kelp. M2DH and FK (for mannitol degradation) fell into “darkgreen” and “bown4” modules (Fig. 3a) which showed the opposite expression patterns with M1PDH1 and M1Pase (for mannitol synthesis).
Fig. 3 Module-trait correlations and gene expression patterns of the top 2 modules correlated with alginate and mannitol contents. (a) Module-trait relationships and corresponding p values. The color scale on the right shows correlations from -1 (green) to 1 (red). Panels on the left represent alginate and mannitol content as traits. Other panels show the expression variations of each biosynthetic gene as a trait. (b) Expression patterns of each selected module. “Brown4” and “Darkgreen” were highly correlated with alginate content, whereas “Black” and “Darkslateblue” were highly correlated with mannitol content
The “brown4” and “darkgreen” modules indicated that the expression levels of those genes correlated with alginate biosynthesis were constantly lower in the basal samples from January to June (JaB, MhB, ApB, MyB toJuB), but were higher in the distal blade (ApD, MyD and JuD) (Fig. 3b). The enriched pathways in “brown4” module included TCA cycle, carbon metabolism and amino acid metabolism etc (Additional file 12: Table S6). In the “darkgreen” module, pathways were concentrated on DNA and protein regulation, including DNA replication (7 DEGs), pyrimidine metabolism (12 DEGs), nucleotide excision repair (7 DEGs), and ubiquitin mediated proteolysis (11 DEG) (Additional file 12: Table S6). Seven representative genes were summarized in Table 4.
Table 4 The description of representative genes in DNA and protein regulation pathways which are highly correlated with “darkgreen” module (p < 0.05)
Pathway
|
Gene ID
|
Connectivity
|
Annotation
|
DNA replication
|
GENE_017950
|
216.81
|
Cdc21-like protein
|
GENE_029087
|
94.91
|
DNA polymerase
|
Pyrimidine metabolism
|
GENE_012005
|
112.52
|
CTP synthase
|
XLOC_025500
|
119.87
|
RNA polymerase II
|
Nucleotide excision repair
|
XLOC_012937
|
146.32
|
Transcription factor II H
|
Ubiquitin mediated proteolysis
|
XLOC_014294
|
71.42
|
Ubiquitin-conjugating e2 j1
|
GENE_001248
|
84.28
|
Ubiquitin-conjugating enzyme 1
|
Genes correlated with mannitol content were up-regulated with the kelp growth and developments (“darkslateblue” module). The expression levels of these genes were higher in the kelp basal blade, with an opposite pattern compared with genes correlated with alginate content (Fig. 3b). Interestingly, “greenyellow” module (M1PDH- and M1Pase-correlated module in Fig. 3a) enriched pathways of energy generation, photosynthesis and photomorphogenesis (p< 0.05) (Table 5). We listed all the genes in these pathways (Additional file 13: Table S7) and found that: 1) twenty-two genes were annotated in the most significantly enriched pathway “oxidative phosphorylation”; 2) more than half of the annotated light harvesting complex protein (LHC) genes fell into “greenyellow” module, together with the genes in photosynthesis e.g. PSII, cytochrome b6/f complex, electron transport and ATPase; 3) the complete biosynthetic pathway from ζ-carotene to violaxanthin, the precursor of fucoxanthin, was annotated; 4) one CRY-DASH gene which encodes cryptochrome in response to blue light was annotated; and 5) the carbon fixation pathway from ribose-5-phosphate to F6P was significantly enriched with at least 7 genes annotated.
Table 5 Genes description in enriched pathways from “greenyellow” module (p < 0.05)
Pathway
|
DEGs in pathway
|
Gene ID
|
Annotation
|
Oxidative phosphorylation
|
22 (6.23%)
|
GENE_024356
|
ATP synthase gamma chain
|
GENE_019595
|
NADH dehydrogenase subunit 10
|
Photosynthesis - antenna proteins
|
18 (5.1%)
|
GENE_022269
|
Light harvesting protein lhcf6
|
Photosynthesis
|
9 (2.55%)
|
GENE_014910
|
Photosystem II 12 kDa extrinsic protein
|
Circadian rhythm - plant
|
5 (1.42%)
|
GENE_003847
|
Cryptochrome 3
|
GENE_017333
|
Phytochrome-like protein 3
|
Carotenoid biosynthesis
|
7 (1.98%)
|
GENE_000415
|
Zeta-carotene desaturase
|
GENE_021497
|
Lycopene beta cyclase
|
GENE_027759
|
Cytochrome P450
|
GENE_017807
|
Violaxanthin de-epoxidase
|
GENE_006795
|
Flavoprotein Monooxygenase
|
Carbon fixation in photosynthetic organisms
|
12 (3.4%)
|
GENE_000154
|
Phosphoglycerate kinase
|
GENE_008411
|
Fructose-bisphosphatase
|
GENE_002455
|
Malate dehydrogenase
|
Validation of the expression of representative alginate/mannitol-related genes
Five genes (MPI2, PMM1, GMD3, MC5E70 and MC5E122) correlated with alginate content, and 3 genes (M1PDH1, M1Pase, M2DH) correlated with mannitol were selected for verification with real-time quantitative PCR (RT-qPCR) assay. Their expression patterns detected with RNA-Seq and RT-qPCR were greatly consistent (Fig. 4), indicating the reliability of high-throughput transcriptomes sequencing. M1PDH1 expression levels were much higher in the adult sporophytes in June than that in earlier developmental stages. It was conflicting with the RNA-Seq but consistent with the extensive accumulation of mannitol in late developmental stages.
Fig. 4 Expression levels of the genes involved in alginate and mannitol metabolism by RNA-Seq and real-time qPCR validation. Thehistogram shows the RT-qPCR results and the curve indicates the expression levels from RNA-Seq results
Screening of transcription factors (TFs) from co-expression analysis
Totally, 57 TFs were identified with the online prediction tool PlantTFDB, and their connectivities range from 2.56 to 238.42. All these TF gene sequences were listed in Additional file 14: Table S8. Eight genes encoding heat shock transcription factors (HSFs) were distributed in 7 modules, being the most abundant TFs. GENE_023741 in “darkolivegreen” was the HSF with highest connectivity and expression levels, which is presumed to have central regulating function. Three MYB or MYB DNA binding protein transcription factors and 3 histone-like TFs (NFY2, NFYB3 and NFYC4) were identified, among which NFY2 (GENE_001494 in “darkolivegreen”) and NFYB3 (GENE_001506 in “darkolivegreen”) exhibited high expression levels and high connectivities. We searched TFs in 4 modules correlated with mannitol and alginate biosynthesis (Table 6; Additional file 15: Table S9). In total, 15 TFs were screened in “black” and “greenyellow” modules, of which 3 HSFs (HSF1, HSF3 and HSFA4C), 1 MYB3R and 1 E2F were identified. Only 4 TFs were found in alginate-correlated “brown4” and “darkorange” modules, of which one RNA polymerase and one histone-like TF were identified (Table 6; Additional file 15: Table S9).
Table 6 The identified transcription factors in the mannitol/alginate-correlated modules
Gene ID
|
Name
|
Module
|
Pfam domain
|
Annotation
|
|
GENE_023779
|
HSF1
|
Greenyellow
|
HSF
|
Heat Shock transcription factor (E. siliculosus)
|
|
GENE_025481
|
MYB3R
|
Black
|
SANT
|
myb transcription factor (Nannochloropsisgaditana)
|
|
GENE_008289
|
E2F
|
Black
|
E2F_CC-MB
|
transcription factor E2F (E. siliculosus)
|
|
XLOC_031444
|
Chrac1
|
Darkorange
|
CBFD_NFYB_HMF
|
histone-like transcription factor family (CBF/NF-Y) (E. siliculosus)
|
|
GENE_024306
|
rpa12
|
Darkorange
|
ZnF_C2C2
|
RNA polymerase I transcription factor TFIIS subunit RPA12 (K. flaccidum)
|
|
Protein interaction networks predicted for mannitol and alginate biosynthesis
We screened the co-expressed genes with alginate and mannitol metabolisms by pearson correlation coefficient ≥ 0.6 or ≤ -0.6. The cytoscape representation is shown in Fig. 5. Alginate-biosynthetic genes highly interacted with genes in DNA and protein regulation. RNA polymerase II (RNApol) and oxoglutarate dehydrogense (ODH)were hub genes affecting alginate content (Fig. 5a). Mannitol-biosynthetic genes were highly correlated with light-responsive reactions. Light harvesting complex protein (LHC), malate dehydrogenase (MDH) and HSF1 play the most important role in the mannitol metabolism (Fig. 5b).
Fig. 5 Cytoscape representation of co-expressed genes with alginate and mannitol metabolisms selected by Pearson correlation coefficient ≥ 0.6 or ≤ -0.6. (a) Protein interaction network for alginate biosynthesis. Genes involved in alginate pathway (red), DNA and protein regulation pathways (green), energy-producing metabolisms (purple) and transcription factor (blue) were indicated. (b) Protein interaction network for mannitol metabolism. Genes involved in mannitol biosynthesis (red), mannitol degradation (yellow), light-responsive reactions (green), energy-producing metabolisms (purple), carbon fixation (brown) and transcription factor (blue) were indicated. Full names for all genes were listed in Abbreviations