RNA sequencing analysis to capture the transcriptome landscape during skin ulceration syndrome progression in sea cucumber Apostichopus japonicus

Background Sea cucumber Apostichopus japonicus is an important economic species in China, which is affected by various diseases; skin ulceration syndrome (SUS) is the most serious. In this study, we characterized the transcriptomes in A. japonicus challenged with Vibrio splendidus to elucidate the changes in gene expression throughout the three stages of SUS progression. Results RNA sequencing of 21 cDNA libraries from various tissues and developmental stages of SUS-affected A. japonicus yielded 553 million raw reads, of which 542 million high-quality reads were generated by deep-sequencing using the Illumina HiSeq™ 2000 platform. The reference transcriptome comprised a combination of the Illumina reads, 454 sequencing data and Sanger sequences obtained from the public database to generate 93,163 unigenes (average length, 1,052 bp; N50 = 1,575 bp); 33,860 were annotated. Transcriptome comparisons between healthy and SUS-affected A. japonicus revealed greater differences in gene expression profiles in the body walls (BW) than in the intestines (Int), respiratory trees (RT) and coelomocytes (C). Clustering of expression models revealed stable up-regulation as the main pattern occurring in the BW throughout the three stages of SUS progression. Significantly affected pathways were associated with signal transduction, immune system, cellular processes, development and metabolism. Ninety-two differentially expressed genes (DEGs) were divided into four functional categories: attachment/pathogen recognition (17), inflammatory reactions (38), oxidative stress response (7) and apoptosis (30). Using quantitative real-time PCR, twenty representative DEGs were selected to validate the sequencing results. The Pearson’s correlation coefficient (R) of the 20 DEGs ranged from 0.811 to 0.999, which confirmed the consistency and accuracy between these two approaches. Conclusions Dynamic changes in global gene expression occur during SUS progression in A. japonicus. Elucidation of these changes is important in clarifying the molecular mechanisms associated with the development of SUS in sea cucumber. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-2810-3) contains supplementary material, which is available to authorized users.


Background
The transcriptome is a set of all RNA transcripts including rRNA, tRNA, mRNA and non-coding RNA produced in one type of cell or a population of certain types of cells at a particular stage in an organism [1]. Unlike the genome, which is roughly fixed for a certain cell type, the transcriptome is considered to be highly dynamic [2][3][4][5][6][7]. These transcriptomic changes are the prelude to the impact of protein translation on the phenotype of the organism. Transcriptome analysis is thus essential for elucidating the underlying molecular constituents of cells and tissues in various biological progresses. With the development of RNA-sequencing (RNA-seq) technology, it provides a far more precise measurement of the levels of transcripts and their isoforms than other methods [8]. Many biologically related issues, such as the expression levels of specific genes, the presence of novel transcripts and fusion transcripts and the events of differential splicing, allelespecific expression and RNA editing can be determined accurately by RNA-seq technology [1,[9][10][11]. Among several RNA-seq technologies, the Illumina sequencing platform is a particularly attractive approach to transcriptomic studies, because it is relatively low cost, and its coverage and depth are far greater than other currently available sequencing technologies [12][13][14][15]. Extensive studies have taken the advantage of the RNA-seq strategy to interpret the dynamic transcriptome in many species, such as maize [16][17][18], Drosophila [2], locust [4], Xenopus [5], and sea cucumber [6]. Other studies have also focused on whole transcriptome analysis to provide an improved understanding of the processes of tumorigenesis and cancer metastasis in human [19,20].
The sea cucumber A. japonicus is one of the most valuable sea foods in East Asian countries because it has a high nutritional value. With increasing market demands for related products, sea cucumber aquaculture has developed rapidly in recent years, and has now become a vigorous industry along the northern coasts of China. However, various diseases in sea cucumber caused by bacteria, viruses or protozoa have severely limited the stable development of this important industry. Among these diseases, SUS is the most common and serious, usually resulting in high mortality. The SUS is characterized by skin ulceration, evisceration, general atrophy, swollen mouth and anorexia [15,21,22]. Some studies have focused mainly on the investigation, isolation and identification of the pathogens responsible for SUS epidemics [22][23][24][25][26][27]. Recently, investigations of the microRNAs (miRNAs) [15,28] and genes involved in immune pathways [29] associated with SUS in A. japonicus have been reported. However, these molecular studies do not provide a comprehensive understanding of the immune response against SUS in A. japonicus and further studies of the molecular mechanisms involved in the SUS epidemics of A. japonicus are required. SUS infection is a dynamic progress. After infection, skin ulceration begins with the appearance of one or more small white ulcerative specks, followed by deep and enlarged ulcerative lesions, leading to exposure of the underlying muscle and spicules. Finally, the ulcerative specks develop into extensive lesions and in severe case, many of the infected sea cucumbers lose the ability to attach to the tanks [21,22,26]. During the onset and progression of SUS, besides the appearance of obvious white ulceration on the skin, the molecular changes triggered by the disease comprise a process that is both dynamic and complicated. Systematic investigations of the multiple pathways involved in the multi-stage SUS progression and its polygenic regulation are required to clarify the origin and development of the disease. In the present study, we screened 21 cDNA libraries from populations of SUS-affected and healthy A. japonicus by high-throughput sequencing using the Illumina HiSeq™ 2000 platform. We compared the differentially expressed genes (DEGs) in four different tissue types (body walls, BW; intestines, Int; respiratory trees, RT; and coelomocytes, C) in SUS-affected and healthy samples, and also investigated the dynamic changes in the transcriptome of BW during the three-stage SUS progression. Furthermore, GO and KEGG enrichment analyses were conducted to elucidate the main functions of DEGs and to provide an overview of the regulation of genes during SUS progression. These results expand our understanding of the complex molecular mechanism of SUS, provide a framework for further investigations of SUS progression, and will be useful in determining strategies aimed at the prevention of SUS caused by bacterial pathogens in sea cucumber.

Transcriptome sequencing and analysis
To obtain a comprehensive gene expression profile for SUS progression in A. japonicus, we first constructed 13 cDNA libraries (from D1 to D13) for RNA-seq that represented the ulcerative body wall and distal normal tissue at all three SUS stages (Fig. 1), as well as four tissues (BW, Int, RT and C) obtained from SUS-affected and healthy individuals. Based on the analysis of the gene expression profiles in the first 13 cDNA libraries, we constructed eight cDNA libraries (from T1 to T8) for RNA-seq that represented the ulcerative BW of the three SUS stages and the BW of healthy individuals using the Illumina HiSeq™ 2000 platform. A total of 553 million (553,212,906) raw reads were generated from 21 libraries; these were deposited in the NCBI SRA database (accession number: SRP050068). After trimming the raw reads, a total of 542 million (541,944,042) high-quality clean reads were obtained ( Table 1). The map reference transcriptome containing 161,174,656 Illumina reads of

Analysis of the gene expression profiles during SUS progression
One of the primary goals of RNA-seq is to compare the gene expression levels from different samples. Analysis of the first 13 cDNA libraries revealed the sequential large-scale gene expression profiles of healthy and SUSaffected A. japonicus individuals, as well as the DEGs of ulcerative and normal BW in the same individuals. Transcripts with at least two-fold differences (log 2 Ratio ≥1, FDR ≤0.001) were classified as DEGs. As shown in Table 3, the largest number of DEGs was observed when the profiles of the Int, RT and C were compared to those of the BW in the healthy group (16,

Stage-series gene expression modeling during SUS progression
The patterns of DEGs for SUS stages I, II and III were assigned as three serial points in the progression of SUS.
To classify the dynamics of the SUS development transcriptome on a global scale, we performed gene expression profile clustering. The seven representative expression patterns of SUS progression are presented in Fig. 2. Visual inspection of these expression models suggested the existence of a diverse and complex pattern of regulation during SUS progression. The first model that represented the main expression pattern comprised Cluster 1, Cluster 4 and Cluster 7, containing 4,245 DEGs. These clusters showed higher levels of DEGs compared to healthy samples and were sustained stably with advancing SUS progression. In contrast, Cluster 2 showed sustained lower levels of DEG expression compared to healthy samples at the three SUS stages. Three expression models based on Cluster 3, Cluster 6, and Cluster 9, revealed a similar expression pattern at SUS stages I and III, while a distinct expression pattern was observed during the transition period at SUS stage II. Moreover, Cluster 5, Cluster 6 and Cluster 8 presented significantly higher levels of expression specific for stages I, II and III, respectively. GO and KEGG enrichment analysis of the stably upregulated expression pattern (4,245 DEGs) during SUS progression showed 175 GO terms (25 GO terms in the category "cellular component", 13 GO terms in the category "molecular function", and 137 GO terms in the category "biological process") (Additional file 1). GO analysis demonstrated that these DEGs with positive regulation function were involved in a broad range of physiological processes. The most notably abundant GO terms in the stably up-regulated expression model were "protein binding", "signal transduction", "cell communication", "signaling" and "single organism signaling". This main stage-series gene expression model of SUS progression associated with biological pathways was evaluated by enrichment analysis of DEGs. Significantly enriched pathways were identified and are presented in Fig. 3 (Additional files 2, 3, 4, 5 and 6; Q value ≤ 0.05). For the 35 significantly enriched pathways identified in this expression model, 14 were related to signal transduction. Eight significantly enriched pathways consisting of "natural killer cell mediated cytotoxicity", "leukocyte transendothelial migration", "Fc epsilon RI signaling pathway", "B cell receptor signaling pathway", "T cell receptor signaling pathway", "Toll-like receptor signaling pathway (TLR pathway)", "Fc gamma R-mediated phagocytosis" and "chemokine signaling pathway" were related to the immune system. Seven significantly enriched pathways were related to cellular processes, with apoptosis and regulation of autophagy pathways being noteworthy. Three significantly enriched pathways, "dorsoventral axis formation", "osteoclast differentiation" and "axon guidance", were related to development.

Investigation of DEGs during SUS progression
Based on the gene expression profiles in the first 13 cDNA libraries, we focused on the changes in the transcriptome in the BW of sea cucumber during the SUS progression from stage I to stage III. Thereafter, we resequenced samples from ulcerative BW of the three SUS stages and healthy BW. Compared to healthy A. japonicus, the number of DEGs in the BW of SUS stages I, II and III were 5,478, 5,589 and 4,437, respectively (Additional files 7 and 8). Combined with the results of the first 13 libraries, we further investigated the DEGs involved in the immune responses from extracellular interaction with bacteria to the activities within the nucleus. Based on GO and KEGG pathway annotations, manual blast and literature searches, 92 DEGs with Nr annotations were selected and divided into four main functional categories: (1) attachment/pathogen recognition (17 genes); (2) inflammatory reactions (38 genes); (3) oxidative stress response (7 genes); and (4) apoptosis (30 genes). A subset of these candidates is listed in Table 4. Among them, the vast majority of these DEGs were up- regulated in the process of SUS progression in A. japonicus. The comparison resulted in a trend concordance between the two sets of sequencing data.

Expression validation using qRT-PCR
To validate the reliability of the RNA-seq results, twenty DEGs from Table 4 were detected by qRT-PCR at three stages of SUS progression (SUS stage I, II and III) and H samples. The cytochrome b (Cyt b) gene was chosen as the reference gene [32]. The 20 DEGs were selected for their clear background information in the function of immune responses, and some of them were involved in the complement pathway, TLR signaling pathway and Apoptosis. The Pearson's correlation coefficient (R) was used to assess the consistency of DEGs expression profiles between these two methods. All the DEGs at three stages of SUS progression showed a consistent expression pattern with R value ranging from 0.811 to 0.999 (Fig. 4).

Changes in gene expression profiles during SUS progression
Using RNA-seq we have generated an extensive transcriptome profile for SUS-affected A. japonicus. The transcriptome profile allowed us to look at the changes in gene expression associated with the progression of the disease. Through construction of the first 13 cDNA libraries and comparison of the DEGs in four different tissues taken from SUS-affected A. japonicus against those of healthy A. japonicus samples, as well as investigation of the transcriptome associated with the three disease stages, we acquired a broad understanding of the dynamic transcriptome during the progression of SUS. Recently, 4,858 DEGs were found in A. japonicus coelomocytes by comparisons  [29]. However, our data showed that between the SUS and healthy individuals, the biggest difference in expression profiles was found in BW. White skin ulceration is the main symptom of SUS; therefore, we further analyzed the changes in the transcriptome in the BW of sea cucumber during SUS progression from stage I to stage III.

Pathways and DEGs during SUS progression
During the progression of SUS caused by V. splendidus challenge, immune regulation was an important event in this host-pathogen interaction process. In fact, most of the significantly enriched pathways involved in the immune system were observed in the up-regulated expression model during SUS progression. Of the eight significantly enriched pathways that were related to immune system, the "TLR pathway" has been reported in A. japonicus [32][33][34][35] response to bacterial challenge and in sea star Pycnopodia helianthoides response to sea star wasting disease [36]. Several signal transduction pathways in the up-regulated expression model during SUS progression, including the insulin [37], TGF-beta [38], MAPK [29], NF-kappa B [39,40], and Notch [41] signaling pathways, have been reported to involve in innate immunity in invertebrate. In the studies, we further investigated DEGs involved in attachment/pathogen recognition, inflammatory reactions and apoptosis.

Attachment/pathogen recognition
Immune responses against invasive pathogens depend primarily on the recognition of pathogen components by innate pattern recognition receptors (PRRs). C-type lectins, which are a type of PRRs involved in innate immunity, specifically recognize sugars on the surface of pathogens in the presence of Ca 2+ , and cause a series of immune responses to resist pathogen invasion [42]. It is worth noting that four C-type lectins, including GalNAc-specific lectin [43], showed downregulated expression and presented a downward trend during SUS progression (Table 4). Mannan-binding Ctype lectin expression was also down-regulated at 48 h and 72 h after LPS challenge in A. japonicus [32]. Further studies are required to clarify the correlation between negative-regulation of this group of receptors and the occurrence of SUS in A. japonicus. Cell adhesion molecules play important roles in many facets of the immune system. Integrins are the   largest family of cell adhesion molecules that mediate cell-to-cell and cell-to-matrix interactions in a broad range of physiological and pathological processes [44].
The essential roles of integrins in developmental and innate immunity in invertebrates have been well studied, especially in phagocytosis [45,46], encapsulation [47] and degranulation [48]. Integrin β1 is upregulated in hemocytes in response to various microbes in Spodoptera exigua [49] and integrin β is also involved as a cell adhesion receptor in the immune responses against microbe challenge in the shrimp  Litopenaeus vannamei [50,51]. Intriguingly, six integrin family genes were up-regulated during SUS progression in A. japonicus. The functions of integrins in A. japonicus response to SUS need to be further studied.

Inflammatory reactions
The innate immune system is the host's first line of defense against infection. It triggers diverse humoral and cellular activities via signal transduction pathways which are conserved both in invertebrates and vertebrates [52]. In our studies, the DEGs included in inflammatory reactions were mainly involved in the complement pathway and TLR signaling pathway. The complement system is a major component of the innate immune system involved in defending against all the foreign pathogens [53]. The key genes of the complement system participating in this defense include complement components C3, C3-2, and factor B, all of which were up-regulated during SUS progression. In our earlier study, the expression of two complement component C3 genes in A. japonicus coelomocytes was also clearly up-regulated after the animal was challenged with LPS [54]. Zhong et al. reported that LPS challenge of A. japonicus induced significant up-regulation of the expression of a homologue of complement B factor (AjBf-2) in the coelomocyte and body wall [55]. C3 plays a pivotal role in the activation of both the classical and alternative complement pathways as well as the lectin pathways [56]. As the second complement component in the alternative pathway, Bf binds to the activated C3b and C3(H 2 O) [57]. In the sea urchin, these two complement components may be part of a simple complement system that is homologous to the alternative pathway in vertebrates [58]. The essential role of TLRs in the activation of innate immunity by recognizing specific patterns of microbial components is well established [59]. While individual TLR genes were not captured, downstream signaling components traditionally associated with TLR pathway were found. For up-regulated DEGs, the important molecules of this pathway, such as MyD88, IRAK4-like, TOLLIP, TAK, NF-kB p105 subunit, NF-kB Rel and stress-activated protein kinase JNK-like, were screened during SUS progression. It is noteworthy that TRAF1 and TRAF5-like, and serine/threonine-protein kinase TBK1-like were associated with the MyD88-independent TLR pathway. Two adaptor molecules in the A. japonicus Toll signaling cascade, MyD88 and TRAF6, have been isolated and characterized. The expression levels of these two genes have been shown to increase significantly during V. splendidus challenge [33]. The NF-kB homologues, Aj-rel and Aj-p105, have also been identified, and shown to be involved in LPS-induced immunity in A. japonicus [34]. The important molecules in this pathway such as MyD88, NF-kB Rel, TRAF6 and LPS-induced TNFα were all shown to be differentially expressed in A. japonicus coelomocytes after LPS challenge by RNA-seq analysis in our early studies [32]. As to the MyD88-independent pathway, two Toll-like receptor genes (AjTLR3 and AjToll) were identified, which were functionally involved in the immune responses against Gram-positive and Gram-negative bacteria, fungi and double-stranded RNA (dsRNA) viruses [35].

Apoptosis
Apoptosis is an essential process in metazoans and is critical for the formation and function of tissues and organs [60]. Dysregulation of the apoptotic processes often leads to serious consequences in humans, such as neurodegenerative diseases [61], cancer [62], and autoimmunity [63]. Some apoptosis-related proteins, such as sarcoma oncogene (Src), vitronectin and vinculin, have been identified and displayed time-dependent depressed expression in the coelomocytes of V. splendidus-challenged A. japonicus [28]. In the present study, we identified the DEGs involved in the apoptosis pathway and determined their expression profiles during SUS progression. There are two major pathways leading to apoptosis: an extrinsic pathway initiated by death receptors and an intrinsic pathway that occurs through the mitochondria. In the extrinsic pathway, procaspase 8 is activated by receptors for FasL and TNF through the recruitment of intracellular death domaincontaining proteins, such as FADD. In the intrinsic pathway, procaspase 9 is activated by cytochrome C released from the damaged mitochondria [64]. Our results showed down-regulated cytochrome C expression during SUS progression. Apoptosis regulator Bcl-2 is a family of evolutionarily related proteins that govern mitochondrial outer membrane permeabilization and perform either proapoptotic (such as BAX) or anti-apoptotic (such as Bcl-2) functions. Interestingly, our studies showed that these two opposing regulators were up-regulated expression during SUS progression. Another important apoptosis regulator, baculoviral inhibitor of apoptosis repeat-containing protein 2 (BIRC2), was identified with up-regulated expression during SUS progression. In addition, expressions of two transcription factors NF-kB p105 subunit and NF-kB Rel that suppress apoptosis were also up-regulated during SUS progression. In the study of human cancer, molecular targeting therapies have been focused on the regulation of apoptosis by Bcl-2 family proteins [65], IAPs [66] and NF-kB [67]. Since apoptotic signals are complicated and regulated at several levels, the mechanism underlying the regulation of apoptosis in SUS of sea cucumber is worthy of further exploration.

Conclusion
The development of SUS in sea cucumber is a complex process in which tens of thousands of genes showed significantly different expression during the progression of the disease. Systematic investigation of the polygenic regulation and multiple pathways involved in the multistage progression of SUS is required to elucidate the dynamic mechanism of SUS progression. The findings reported in this study will be useful for further studies into the origins and development of SUS in sea cucumber. In further investigations, we will focus on the network biology approach to comprehensively depict the miRNA-mRNA and mRNA-protein networks relevant to SUS in A. japonicus.

Animals
Animals used in this research were obtained from commercial sea cucumber catches, therefore approval from any ethics committee or institutional review board was not necessary. Healthy A. japonicus (10-12 g) were collected from Zhuanghe, Liaoning Province (China). The animals were acclimated in the laboratory for 1 week, and then subjected to artificial infection. The sea water used in the experiment was filtered through sand, and then through 300-μm nylon sieves. Twenty-five percent of the seawater in the tank was exchanged daily. The animals were maintained at 12°C, pH 8.1, with salinity of 32 and continuous aeration.

Artificial infection and sample collection
V. splendidus used to infect A. japonicus was previously isolated and identified in our laboratory from SUSaffected A. japonicus according to the method reported by Deng et al. [22]. The bacteria were cultivated in 2216E medium at 28°C for 24 h. The bacterial cells were then harvested by centrifugation at 1,000 × g for 5 min, and then re-suspended in 0.22-μm-filtered seawater. For the wounded immersion infection, 200 healthy A. japonicus individuals were cultured in five rectangular tanks (80 cm in length, 45 cm in width, and 45 cm in height). A cut (0.5 cm × 0.5 cm) was made in the body wall of the sea cucumbers using a scalpel and V. splendidus was added to the tank (final concentration, 5 × 10 9 CFU mL −1 ), with 25 % seawater replenished daily; V. splendidus was supplemented and the bacterial concentrations were maintained. In the sample collection, white skin ulceration was considered to be the most important mark to distinguish diseased and healthy individuals [15,28]. The progression of SUS was divided into three stages. In stage I, the animals showed one small white speck of skin ulceration (diameter <0.2 cm). The animals retained the ability to attach to the surface of the tank and did not eviscerate. In stage II, the animals exhibited 2 to 3 larger white specks (diameter >0.2 cm). The animals continued to exhibit the ability to attach to the surface of the tank and did not eviscerate. In stage III, the individuals showed several deep and extensive ulcerations, lost the ability to adhere to the tank and displayed evisceration. A. japonicus that cultured under normal conditions without the treatments of cut and bacterial challenge served as healthy controls. Thirty individuals were selected from each of the SUS stages (stage I, stage II, and stage III) and the ulcerative body wall and its distal normal tissue were sampled for RNA analysis (Fig. 1). In addition, samples of the Int, RT and C were also obtained at SUS stage II. The C were collected by centrifugation at 1,000 × g for 5 min. Samples of BW, Int, RT and C were separated from 30 healthy A. japonicus individuals. For each sampled individual, the size of BW tissues is about 0.2 cm × 0.2 cm × 0.2 cm. The quality of sampled Int and RT tissues is less than 100 ng. The volume of C is about 1 mL. All samples were frozen immediately in liquid nitrogen and then stored at −80°C before RNA isolation.

RNA extraction, cDNA library construction and sequencing
Total RNA was extracted from the specimens with Trizol reagent (Invitrogen, Carlsbad,CA,USA) according to the manufacturer's instructions. The quality and concentration of total RNA were measured by Nanodrop 1000 (Thermo). The initial amount of high-quality total RNA was 1 μg per library (generated from 30 individuals). Subsequently, the mRNA in the total RNA was enriched using Oligo (dT) magnetic beads and sheared into short fragments by adding fragmentation buffer, followed by first-and second-strand cDNA synthesis using random hexamer primers. The cDNA fragments were subjected to an end-repair process, addition of "A" base, and ligation of sequencing adapters. After agarose gel electrophoresis, suitable fragments were selected and used as templates for the PCR amplification to create the final cDNA libraries. We first constructed 13 cDNA libraries (from D1 to D13), representing the ulcerative BW and their distal normal tissues at all three stages of SUS development (Fig. 1) and four tissues (BW, Int, RT and C) of SUS-affected and healthy samples. High-throughput sequencing was conducted using the Illumina HiSeq™ 2000 platform to generate 50-bp reads.

Data processing, assembly and annotation
The original image data were converted to sequence data by base-calling and saved as fastq files. The raw reads were then cleaned by discarding adaptors, low-quality reads (quality scores <20), reads with unknown bases greater than 5 %, and reads of less than 20 nt. De novo transcriptome assembly combined with Illumina reads obtained from GenBank [28,30] was carried out with the shortread assembly program, Trinity [68]. The assembled sequences from the Illumina reads, 454 sequences and ESTs were used in further processing of the assembly according to the method reported by Zhou et al. [30]. Blastx analysis Table 5 Primers used for qRT-PCR validation   Gene  Primer Sequence (5′-3′ of unigenes longer than 200 bp was conducted against the Non-Redundant (Nr) database, Swiss-Prot, COG and KE GG database with an e-value cutoff of 1E-5. GO classification was analyzed by Blast2GO software (v.2.5.0) based on best hits of Nr annotation. In addition, transcripts were also annotated with the NCBI non-redundant nucleotide (Nt) database using Blastn with an e-value cutoff of 1E-5.

Screening and analysis of DEGs
The unigene expression was calculated using the RPKM (Reads Per kb per million reads) method [9]. To identify the DEGs in the first 13 A. japonicus libraries, a rigorous algorithm was developed for statistical analysis according to "the significance of digital gene expression profiles" [69]. The FDR (false discovery rate) is used in multiple hypothesis testing to correct for P-values [70]. In our analysis, the FDR ≤ 0.001 and RPKM ratio larger than 2 (|log 2 Ratio| ≥ 1) were used as the threshold to judge the significance of differences in gene expression. The dynamic expression profiles of DEGs during the threestage SUS progression were visualized by clustering analysis, based on the K-means method using the Euclidean distance algorithm. The enrichment analysis was carried out based on an algorithm presented by KOBAS [71], with the entire transcriptome set as the background. The P-value was approximated by a hypergeometric distribution test and FDR multipletesting correctionwas used to corrected-P-value [70]. The GO enriched cutoff was a corrected-P-value ≤ 0.05, and the KEGG enriched cutoff was a Q value ≤ 0.05.
Transcriptome sample preparation, re-sequencing and DEG identification After analysis of the gene expression profiles in the first 13 cDNA libraries, we re-sequenced samples of the ulcerative BW of the three SUS stages and healthy BW. The three stages of SUS-affected samples were also infected with V. splendidus. Subsequently, ulcerative BW of individuals during the three SUS stages and healthy BW were separated for cDNA library construction (from T1 to T8) and RNA-seq. Two biological replicates of each sample were set up. High-throughput sequencing was conducted using the Illumina HiSeq™ 2000 platform to generate 100-bp paired-end reads. After quality control for raw reads, clean reads were mapped to the reference transcriptome and their abundance estimated using RSEM (RNA-seq by expectation maximization). DEGs were detected by using the DESEQ program from the R-Bioconductor package.

Expression validation using qRT-PCR
To validate the RNA-seq results, 20 DEGs was employed for qRT-PCR (Mx3005p™ detection system, Agilent Stratagene, Santa Clara, CA, USA). Total RNA from samples used in the first RNA-seq was reverse-transcribed into cDNA templates with the PrimeScript™ RT reagent Kit (TaKaRa, Otsu, Japan) according to the manufacturer's instruction. Reactions were incubated at 37°C for 15 min, and then 85°C for 5 s. According to the sequence information in the transcriptome database, primers were designed using Primer 5.0 software according to rigorous criteria. The primer information was provided in Table 5.
The cytochrome b (Cytb) gene was chosen as the reference gene [32]. Optimal primer pairs were examined by checking the melting curve at the end of each PCR reaction to confirm the specificity of PCR product. The qRT-PCR amplification was conducted in a volume of 20 μL containing 10 μL of 2× SYBR Premix Ex Taq™ [72] was used to calculate the expression differences between control and SUS groups (SUS stage I, II and III). The correlation of DEGs expression profiles between RNA-seq and qRT-PCR was assessed by Pearson correlation coefficient (R). The nearer the scatter of points is to the straight line, the higher the strength of association between the variables. The value of R ranges from −1 to +1 [73,74].