1 INTRODUCTION

As an economically important kelp, Saccharina japonica (Laminariales) plays an irreplaceable role in marine ecosystems (Tonon et al., 2011). It has become the largest biomass harvested among all seaweeds because of its extensive cultivation in East Asia, Europe and North America for its three major carbohydrates, alginate, laminarian and mannitol (Nyvall et al., 2003; Takeda et al., 2011; Kraan, 2013). Kelp has also become ideal feedstocks for the production of bioethanol, and their cultivation does not require arable land, fertilizer or fresh water (Wargacki et al., 2012; Enquist-Newman et al., 2014; Yang and Li, 2014). According to the 2016 China Fishery Statistics Yearbook, the annual S. japonica production of China in 2015 reached 2.089 2 million tons, accounting for more than 90% of the world total. Despite the long history of cultivation, we currently have a very limited knowledge of its metabolism.

Saccharina japonica is cultivated upside down under a floating rope net with the holdfast and meristem facing the sunlight and with the blade in the shade, making S. japonica sporophyte persistently experience adverse stresses during its development. Environmental shifts including constantly aggravating pollution, continuously changing climate and runaway intensification of aquaculture scale may interact and worsen these stresses (Ye et al., 2015). The life cycle of S. japonica is alternative between heteromorphic sporophyte and gametophyte. This study focused on the sporophyte because it is cultivated and harvested. The sporophyte of S. japonica can be divided into a series of stages different in morphology and physiology (Wu et al., 2015). The first is mushroom stage at which S. japonica is thin and tender and about one meter long; the second is adult stage at which the blade becomes hard and tough, and the content of organic matter increases significantly; the third is mature stage at which S. japonica stops to grow and sporangia develops on blade surface; and the last is aging stage at which the blade dies gradually. The S. japonica sporophyte changes in morphology and physiology at different developmental stages; however, the physiological mechanisms underlining these changes are still unknown, and the metabolism processes of S. japonica during its development are rarely documented.

With the popularization of massively parallel sequencing technologies, RNA sequencing (RNAseq) has been used to decipher the molecular mechanisms that underline the physiological processes of diverse organisms. The aim of this study was to disclose the difference of expressed genes and the metabolic pathways they defined at different developmental stages of S. japonica through RNAseq.

2 MATERIAL AND METHOD 2.1 Content determination of chemicals

Sporophytes were collected from Sanggou Bay, Rongcheng, Shandong Province in China. The sporophytes were cultivated on a rope net floating permanently at the surface of the seawater. The rope net consisted of parallel master ropes anchored to the seafloor and sporophyte hanging ropes knotted between the master ropes at an equal distance. Each sporophyte hanging ropes had the sporophytes of Dongfang No. 7, a variety of S. japonica (Li et al., 2016). At four different developmental stages (mushroom, adult, mature and aging) corresponding to four different months (March, April, June and July, 2015, respectively), sporophytes were collected from three points along a sporophyte hanging rope, each at surface, middle and bottom positions and three sporophyte hanging ropes each stage. In total, 36 sporophytes were collected, which were sun-dried, ground to powders ≥425 μm, mixed, sealed into a plastic bag and stored in dry air and at room temperature. Alginate content was determined with the method described by Shang et al. (2011) while mannitol and iodine contents were assayed following the method described early (Ji, 2004). Simultaneously, the blade length, width, color and flexibility were measured in field and recorded as described early (Li et al., 2016).

2.2 RNA extraction and transcriptome sequencing

Sporophytes were simultaneously collected with those for determining alginate, mannitol and iodine contents, each at surface, middle and bottom position of a hang rope and two hanging ropes each stage. In total, 24 sporophytes were collected, six each stage. Tissues from the growing point each sporophyte of the same developmental stage were cut, grouped, fast frozen in liquid nitrogen and used to extracting total RNA with TransZol Plant, a kit of Invitrogen, following manufacturer's instructions. The concentration and purity of RNA was checked on the NanoDrop 2000c Spectrophotometer (Thermo Fisher Scientific, USA). The quality of RNA was determined through electrophoresis in 1% agarose gel in order to ensure that the extracted RNA was not degraded and contaminated by genomic DNA. The mRNA was purified using poly-T attached magnetic beads, fragmented and used as the templates of the firststrand cDNA synthesis using random hexamer primers and M-MuLV reverse transcriptase (RNase H free). The second-strand cDNA was synthesized, subsequently, using RNase H and DNA polymerase I. The double stranded cDNA was end polished with exonuclease and polymerase. After nucleotide "A" addition at the 3′ ends of cDNA fragments, NEBNext Adaptor was linked on. The cDNA in different lengths were selected using AMPure XP system (Beckman Coulter, Beverly, USA) and amplified through PCR, which were purified using AMPure XP system (Beckman Coulter, Beverly, USA), yielding the cDNA library for sequencing. In total, 8 libraries, 2 each stage, were sequenced on an Illumina HiSeqTM2000/-MiSeqTM platform (Novogene Bioinformatics Technology Co., Ltd.).

2.3 Sequence data processing

The sequencing records were translated into the raw reads in FASTA format through CASAVA Base Calling (Zhang et al., 2015; Miao et al., 2016). Clean reads were obtained after removing those either containing adaptor or ploy-N or at low quality. Clean reads were assembled using Trinity (Grabherr et al., 2011) to generate unigenes (here after genes). Differential expression analysis was performed using the DEGseqR package (1.12.0) (Anders and Huber, 2010). The abundance of gene transcripts were estimated based on FPKM (fragment per kilobase of exon length per million reads) values and the abundance of gene transcripts was considered to be significantly differential if q-value < 0.05 (Trapnell et al., 2010). Differentially expressed genes (DEGs) were selected by comparing two adjacent stages of S. japonica development. GO enrichment of DEGs was carried out using the GOseqR Package based on a hypergeometric test (Young et al., 2010). GO functional classifications were performed with WeGO software (Ye et al., 2006). DEGs were also enriched into KEGG pathways (http://www.genome.jp/kegg/) with KOBAS software (Mao et al., 2005). The enrichment P-values were adjusted with Benjamin and Hochberg method. In order to appropriately describe the function of expressed genes, their corresponding models in the S. japonica genome (Ye et al., 2015) were Blast screened against SwissProt for their functional homologs. In addition, the genes involved in sugar biosynthesis were screened out through blasting the expressed gene assemblage with the enzyme amino acid sequences retrieved from NCBI as queries with BLAST v 2.2.9 (Altschul et al., 1997).

3 RESULT AND DISCUSSION 3.1 Morphological and metabolite shifts at different development stages

The morphological changes of S. japonica at different sporophyte developmental stages mainly manifested in the length, color and flexibility of blade (Fig. 1). The length of blade increased from mushroom to adult, and simultaneously the color became deeper and deeper and the flexibility became tougher and tougher. The most protrusive morphological changes included the blade thickening from mushroom to adult stages and the appearance of sporangia from adult to mature stages. A large number of sporangia protruded from surface of the entire blade at mature stage and were observable with naked eyes at both mature and aging stages. The percentile content of alginate was found to be similar among developmental stages while that of mannitol showed a gradual increase from early to latter developmental stages and that of iodine showed a difference before adult and after mature stages (Table 1).

3.2 Transcriptome profiling

In total, 204 992 741 raw reads were obtained for 4 developmental stages, of them, 196 238 591 were clean, accounting for 95.73% of the total. Of the clean reads, 79.96% were mapped onto the S. japonica genome (Ye et al., 2015), and 49.55% to the gene models of the genome. From the clean reads, 7 800 gene models were found to be transcribed, accounting for 58.53% (7 800/13 327). In total, 915 genes were found to be differentially expressed between stages, of them, 64 were between mushroom and adult stages (27 up-regulated and 37 down-regulated), 539 between adult and mature stages (348 up-regulated and 191 down-regulated), and 302 between mature and aging stages (99 up-regulated and 203 downregulated) (Table 2, Fig. 2).

Against SwissProt, only 2 697 genes of 7 800 expressed in total were annotatable. Against more data bases, the number of annotatable genes should be increasable; however, an obvious boost of such number was not expected. We found that the total abundance of the annotatable genes accounted for 50.44% of the total. Of the expressed, we found 778 genes differentially expressed between at least two stages. Of these differentially expressed genes (DEGs), only 372 were annotatable, which accounted 19.23% of the total abundance of DEGs (Table 3).

It should be noticed that S. japonica was not a model organism. It is not close to the plant lineage consisting of red and green algae and high plants but rather phylogenetically close to diatoms, dinoflagellates, ciliates among others (Keeling et al., 2005; Baldauf, 2008; Lane and Archibald, 2008). Till present, the genes isolated and characterized in S. japonica were very limited. As a brown algal model, the genome of Ectocarpus siliculosus has been sequenced (Cock et al., 2010), however, a large percentage of genes remained non-annotated. In addition, the genetic and physiological difference between S. japonica and E. siliculosus should be large although they are all brown large. Such scenario may have restricted the annotation of S. japonica genes and also description of its physiological processes, thus it was understandable that less genes were annotatable, especially when the annotation is against only one data base (here SwissProt).

3.3 Weird transcripts

Very interestingly, we found that 1 208 of the 7 800 expressed and 2 697 annotated were virus associating genes and the total abundance of their transcripts accounted for about 10.21% of the total expressed. Such weirdness was observed also in the DEGs. Of 778 DEGs, 372 were annotatable and 209 were virus associating genes with total abundance of their transcripts accounted for 8.90%. It was true that the expressed genes were annotated against SwissProt only, and some of virus associating genes may be functionally assignable against other data bases. However, such a large portion of virus associating genes indicated that the S. japonica genome contained a large portion of virus originating genes, and these genes were active indeed in the current evolved genome. This observation implied that the S. japonica genome was weird. Delaroque and Boland (2008) found that the genome of the brown alga E. siliculosus contains a series of viral DNA pieces, suggesting an ancient association with large dsDNA viruses. It seemed that S. japonica genome did contain also a large portion of virus originating genes, thus evolving in a specific route, which is worthy of further study.

We identified these DESs by referring to the published genome sequence (Ye et al., 2015). Our findings generated multiple, magic and mystery stories to be unveiled. It was obvious that the S. japonica genome should be resequenced with newly available tools and strategies (Eid et al., 2009; Lieberman-Aiden et al., 2009; Burton et al., 2013; Vanburen et al., 2015) to meet the standards of the reference in order to solve the questions such as less annotated genes (Table 3) and possibly missed gene models. Similarly the genes implying S. japonica genome evolution included a large set of functionally repetitive genes. We found that 439 of 2 697 annotated genes were multiple copy genes, which either were or contained the domain of known proteins. Again these genes were worthy of further inspection of their structures and functions (Table 3).

3.4 Gene expression pattern associating with structural and storage sugar changes

In order to describe the physiological mechanisms underlining the morphological and biochemical changes at different developmental stages, we tried to enrich the genes differentially expressed between at least 2 stages, unfortunately, we failed to find any metabolism and signaling pathway. As was pointed out above, the annotated transcripts were less than the expected in this non-model brown alga, and those annotated against SwissProt covered a large portion of virus genes or virus protein domain fusing genes. As an alternative strategy, we turned to find out the genes controlling the biosynthesis of structural and storage sugars, which were identified in recent years from the genome of E. siliculosus (Michel et al., 2010a, b; Tonon et al., 2017) using a local BLAST approach and the previously identified genes as queries. In total, 19 genes were found, which controlled the biosynthesis of 5 types of sugars including sulfated fucan, alginate, mannitol, trehalose and laminarin (Table 4; Fig. 3).

A FPKM value represents the normalized abundance of a gene transcript. In order to describe the biosynthesis intensity of a sugar, the transcript abundances of the genes involved in its biosynthesis were summed up first in the direction of developmental stages and then in the direction of the list of all genes involved in its biosynthesis with the total sum averaged among the genes involved in its biosynthesis as the intensity of its biosynthesis. Following the same philosophy, the abundances of all genes were summed up in both directions and averaged as the whole genome transcription (background) and used as the reference to the biosynthesis intensity of a sugar. It was found that the transcripts abundance associating with sugar biosynthesis was about 2.13 folds of all expressed genes (10.434/4.893), indicating that the biosyntheses of structural and storage polysaccharides were very important cellular processes of S. japonica. Interestingly, the abundance of genes responsible for structural (alginate) and storage (mannitol and laminarin) polysaccharides was similar to each other (11.538/12.741). From the averaged sum of the abundance of genes involved in the biosynthesis of each polysaccharide, we found that the biosynthesis of laminarin was the strongest process (17.781), which was followed by that of alginate (11.538), mannitol (7.701), trehalose (7.331) and sulfated fucan (4.510). At different developmental stages, the total abundance of genes involved in the biosynthesis of alginate and laminarin was similar among stages, however, that of genes involved in the biosynthesis of mannitol increased about 2 folds from mushroom and adult stages to mature and aging stages. Such trend explained our observations that the content of alginate was almost constant at different developmental stages while that of mannitol increased gradually (Table 1).

3.5 Persistent expression of defense associating genes

Of 2 697 genes annotated against SwissProt, we found that a set of defense and cell recurring genes highly expressed in S. japonica sporophyte development, and many of them expressed differentially among stages, which included those encoding vanadium-dependent bromoperoxidase, heat shock 70 protein, and peroxiredoxin and peptide methionine sulfoxide reductase (Table 5). On average, the sum abundance of the transcripts of these genes at 4 developmental stages was 3.40 (24.258/7.138) and 4.96 (24.258/4.893) folds of all annotated and all expressed genes, respectively, indicating that S. japonica sporophytes persistently respond possible pathogen and environmental stresses.

Such physiological status is possible as strong UV irradiation, worse pollution, changing climate, intensive cultivation among others may have been making S. japonica sporophytes face to continuous stresses.

The pathways of phytohormone biosynthesis in higher plants have been well described, and the role of phytohormones in the regulation of metabolic processes in algae is no longer in doubt (Kiseleva et al., 2012); however, the hormone metabolism at gene level in different groups of algae remains largely unknown although all known phytohormones are found in various algae, and auxin metabolism and function seem to be characterized in E. siliculosus (Le Bail et al., 2010). Similarly, oxylipins, a group of oxygenized polyunsaturated fatty acids (PUFAs) and derivatives, have been described in diverse algae, especially in brown algae like S. japonica. The mechanism of polyunsaturated fatty acids (PUFAs) oxidation at least via a lipoxygenase step and then alterative and diverse subsequent reactions collectively is referred to oxylipin pathway (Andreou et al., 2009). The biosynthesis of diverse oxylipin is the response of organisms to both abiotic and biotic stresses, which forms a few rings along a chain from the recognition and reception of stress factors and pathogenic initiators and effectors to the transcriptional reprogramming of genome, namely oxylipin signaling, a widely observed response mechanism of organism also referred to the innate immunity (Ritter et al., 2014). As the intensively cultivated kelp species, S. japonica is facing a set of biotic and abiotic stresses, including continuously raising seawater temperature, content of carbon dioxide, strong fluctuation of sunshine, strengthening UV irradiation, territorial sources of metal pollutants, intensifying pathogens due to continuous and high density cultivating among others. The genome sequence S. japonica is available currently (Ye et al., 2015), however the gene assemblage encoding the enzymes functioning in oxylipin synthesis and signaling have not been well characterized yet. The analyses in a model brown alga, E. siliculosus, have shown that this brown alga genome contains neither jasmonate (a group of oxylipins) receptor nor its synthesizing enzyme encoding genes (Cock et al., 2010; Ritter et al., 2014). We speculated that S. japonica may do so because they are closely related each other. It has been hypothesized that S. japonica responds to both abiotic and biotic stresses by up-taking and metabolizing iodine, bursting H₂O₂, synthesizing both octadecanoids and eicosanoids derived oxylipins, and reprogramming gene transcription of its genome with unknown transcriptional networks (Küpper et al., 1998; Cosse et al., 2009; Ritter et al., 2014; Thomas et al., 2014; Tsuda and Somssich, 2015).

Vanadium-dependent bromoperoxidase is a radiative oxygen species (ROS) detoxifying enzyme specific for brown algae. Heat shock 70 may facilitate the refolding of damaged proteins and aid to protecting protein aggregation under oxidative stress. The other two proteins associate also with the stress response. In fact, these genes were picked up from those annotated against only SwissProt. We are sure that homologous searching should find more stress responding genes. However, our findings should have implied that cultivated S. japonica persistently responds to abiotic and biotic stresses. In recent years, environmental changes may have placed the cultivated S. japonica at a threshold, a little lower making it develop healthily, and a little higher causing the outbreak of diseases. Breeding high disease resistant varieties is appreciated while reforming S. japonica cultivating system may realize such purpose fast.

4 CONCLUSION

Transcripts of kelp (Saccharina japonica) at different sporophyte developmental stages were a little weird. The transcription of associating genes and the metabolite shift of sugars at were consistent with each other. The content of alginate was almost constant at different developmental stages while that of mannitol increased sharply. S. japonica sporophytes persistently respond possible pathogen and environmental stresses by expressing defense genes. Further research is required to reveal the mechanism of S. japonica defense.

5 DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author on request.