RNA-seq analysis reveals significant transcriptome changes in turbot (Scophthalmus maximus) suffering severe enteromyxosis

RESEARCH ARTICLE Open Access RNA-seq analysis reveals significant transcriptome changes in turbot (Scophthalmus maximus) suffering severe enteromyxosis Diego Robledo 1† , Paolo Ronza 2† , Peter W Harrison 3 , Ana Paula Losada 2 , Roberto Bermúdez 4 , Belén G Pardo 5 , María José Redondo 6 , Ariadna Sitjà-Bobadilla 6 , María Isabel Quiroga 2* and Paulino Martínez 3,5 Abstract Background: Enteromyxosis caused by the intestinal myxozoan parasite Enteromyxum scophthalmi is a serious threat for turbot (Scophthalmus maximus, L.) aquaculture, causing severe catarrhal enteritis leading to a cachectic syndrome, with no therapeutic options available. There are still many aspects of host-parasite interaction and disease pathogenesis that are yet to be elucidated, and to date, no analysis of the transcriptomic changes induced by E. scophthalmi in turbot organs has been conducted. In this study, RNA-seq technology was applied to head kidney, spleen and pyloric caeca of severely infected turbot with the aim of furthering our understanding of the pathogenetic mechanisms and turbot immune response against enteromyxosis. Results: A huge amount of information was generated with more than 23,000 identified genes in the three organs, amongst which 4,762 were differently expressed (DE) between infected and control fish. Associate gene functions were studied based on gene ontology terms and available literature, and the most interesting DE genes were classified into five categories: 1) immune and defence response; 2) apoptosis and cell proliferation; 3) iron metabolism and erythropoiesis; 4) cytoskeleton and extracellular matrix and 5) metabolism and digestive function. The analysis of down-regulated genes of the first category revealed evidences of a connexion failure between innate and adaptive immune response, especially represented by a high number of DE interferon-related genes in the three organs. Furthermore, we found an intense activation of local immune response at intestinal level that appeared exacerbated, whereas in kidney and spleen genes involved in adaptive immune response were mainly down-regulated. The apoptotic machinery was only clearly activated in pyloric caeca, while kidney and spleen showed a marked depression of genes related to erythropoiesis, probably related to disorders in iron homeostasis. The genetic signature of the causes and consequences of cachexia was also demonstrated by the down-regulation of the genes encoding structural proteins and those involved in the digestive metabolism. Conclusions: This transcriptomic study has enabled us to gain a better understanding of the pathogenesis of enteromyxosis and identify a large number of DE target genes that bring us closer to the development of strategies designed to effectively combat this pathogen. Keywords: RNA-seq, Transcriptome, Turbot, Enteromyxum scophthalmi, Enteromyxosis, Immune response, Apoptosis, Erythropoiesis, Cytoskeleton, Digestive function * Correspondence: [email protected] † Equal contributors 2 Departamento de Ciencias Clínicas Veterinarias, Facultad de Veterinaria, Universidad de Santiago de Compostela, Lugo 27002, Spain Full list of author information is available at the end of the article © 2014 Robledo et al.; licensee BioMed Central. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Robledo et al. BMC Genomics 2014, 15:1149 http://www.biomedcentral.com/1471-2164/15/1149 Background Turbot (Scophthalmus maximus, L.) is a marine flatfish which has been intensively cultured in Europe for more than 25 years. Turbot aquaculture production currently accounts for over 10,000 tons/year in Europe, and has rapidly increased in China during the last decade, reaching more than 60,000 tons in 2011 [1]. Although its entire life cycle is routinely carried out at farm facilities, major challenges related to growth rate and disease outbreaks are the main concerns for the turbot industry [2-5]. The myxozoan genus Enteromyxum includes three intestinal species, which all cause serious problems for seawater aquaculture. Enteromyxum’s virulence varies depending on the specific parasite-host interaction, rangingfromhighlypathogenicinsomespeciespairs, such as E. scophthalmi infection of turbot, to chronic (E. leei in gilthead sea bream, Sparus aurata) or subclinical (E. fugu in tiger puffer, Takifugu rubripes)diseasesignsin others. E. scophthalmi is a serious threat for cultured turbot, that spreads rapidly in farm facilities due to direct fish-to-fish transmission, and causes a cachectic syndrome which eventually leads to death [6]. Diseased fish present anorexia, anaemia, weight loss, poor conversion rates and delayed growth, with mortality rates reaching up to 100% in many cases [6,7]. The main histological changes are catarrhal enteritis and lymphohaematopoietic depletion of spleen and kidney, the severity of which increases with the progression of the infection [8]. Although some drugs have been able to lower the mortality rates [9], currently there is no effective treatment for enteromyxosis. Understanding the disease pathogenesis through the study of host-parasite interaction and turbot immune response is critical in order to develop effective treatments and apply preventive measures. Numerous recent studies have been focused on elucidating these processes, mainly through histopathological, immunoenzymatic and serological assays [10-16]. Nevertheless, while PCR-array and microarray based molecular profiling of gilthead sea bream response to E. leei has been recently published [17,18], gene expression characterization of E. scophthalmi-infected turbot is lacking. Transcriptome analysis is an invaluable tool for the elucidation of the biological processes behind host-parasite interactions, and in the last decade this approach, mainly based upon microarrays, has been extensively used in fish immunology and pathology [19], and specifically in turbot for analyzing furunculosis and scuticociliatosis [20,21]. Additionally, the identification of relevant immune gene variants conferring tolerance to parasites is essential in order to develop marker assisted selection programmes that can lead to increased resistance [17,22]. RNA-seq is a powerful technique for the analysis of gene expression due to its higher sensitivity and specificity in comparison to microarrays, along with its ability to detect new genes, rare transcripts, alternative splice isoforms, and novel SNPs which can be used for association studies [23-25]. For these reasons it is rapidly becoming the technology of choice for transcriptomic studies [19,26]. In this study, a gene expression analysis of E. scophthalmiinfected turbot was carried out using RNA-seq on the two major lymphohaematopoietic organs, head kidney and spleen, and also on the pyloric caeca, the target intestinal region where the parasite infection starts in this species [27]. Fish categorised with a severe infection were used with the aim of capturing the gene expression signatures associated with advanced stages of the disease as a first reference to investigate the genetic mechanisms underlying the pathogenesis of enteromyxosis. Our findings constitute the basis of future studies aimed at investigating resistance-related genes and associated genetic variants that could be applied in breeding programmes. This is the first study to tackle the molecular basis of lesion development and the immune response underlying enteromyxosis in turbot. Results Histopathology Experimentally-infected (recipient, RCPT) fish selected for RNA-seq analysis presented catarrhal enteritis characterized by severe parasitic load along the entire gastrointestinal tract associated with moderate to severe inflammatory infiltrates and lining epithelium detachment.Inmostgutsegments,apoptoticfigures in both epithelium and lamina propria were observed, while signs of epithelium regeneration were sporadically annotated. Spleen and head kidney showed moderate to severe cellular depletion, with occasional observation of indicators of apoptosis. There were no significant histological changes in the other organs of recipient fish (RCPT) neither in any samples from control fish (CTRL). RNA-seq A total of ~170 million 100 bp pair-end reads were sequenced, accounting for on average 15 million reads per sample post-filtering. Filtered reads were mapped to the turbot genome identifying a total of 54,864 transcripts and 23,063 genes in the three organs. The average number of raw, filtered and mapped reads for the samples of each organareshowninTable1.Foreachorganthecontroland infected samples were hierarchically clustered according to Table 1 RNA-Seq sample statistics Sample Raw reads Trimmed reads Aligned reads Head kidney 17,431,026 15,929,258 14,321,389 Spleen 16,393,159 14,705,171 13,124,477 Pyloric caeca 17,094,324 15,341,340 9,963,567 Average raw reads obtained by Illumina sequencing, average trimmed reads remaining after filtering and average reads aligned to the turbot genome per sample are shown for each organ. Robledo et al. BMC Genomics 2014, 15:1149 Page 2 of 17 http://www.biomedcentral.com/1471-2164/15/1149 their transcript expression (Figure 1) confirming the correct classification of infected and control samples. Differential expression analysis A total of 4,762 differentially expressed (DE) genes were identified across the three organs when comparing RCPT and CTRL fish. The number of DE genes (upand down-regulated) for each organ and those shared between organs is shown in Figure 2. A high number of downregulated genes (3,062) were detected, 68.5% more than up-regulated (1,817). Pyloric caeca showed the highest amount of DE genes, almost double that of the other two organs. The percentage of DE annotated transcripts was similar in the three organs: 44.1% in head kidney, 181 up-regulated and 400 down-regulated; 42.3% in spleen, 229 up-regulated and 353 down-regulated; and 46.9% in pyloric caeca, 562 up-regulated and 851 down-regulated. Log 2 fold change (FC) values ranged from 11.26 to −11.18 in head kidney, from 13.29 to −12.83 in spleen, and from 12.13 to −15.18 in pyloric caeca. Common DE genes between the three organs A group of 117 DE genes were shared between the three organs. Among them, 48 were successfully annotated: 11 upand 26 down-regulated, and 11 either upor downregulated depending on the organ (Additional file 1: Table S1, ST1). Amongst the shared up-regulated genes, some were involved in innate immune response and antigen presentation such as IL4I1, involved in the lysosomal processing and presentation of antigens; ALOXE3,anenzymethat participates in leukotrienes metabolism; and MASP3b, which plays a prominent role in the activation of the lectin complement pathway. Amongst the down-regulated genes, two sets were of particular interest, one directly involved in immune response and the other related to cell and tissue structure disruption. The first group included several interferon (IFN)-related genes involved in antiviral immune response, such as IRF7, Gig2, IFIT-1, cHERC5 and IFI44. Other two interesting genes were related to major histocompatibility complex class I molecules (MHC-I), involved in the presentation of intracellular-derived antigens, and CAT, which encodes an essential antioxidant enzyme for cell protection against oxidative damage. The second set included several genes related to cytoskeleton: FMNL1 (cytoskeletal organization and cell morphology and motility); TMOD4 (geometry of themembraneskeleton);GRXCR1 (architecture of actin filament-rich structures); and SNPH (microtubule-associated protein). COL1A2, involved in extracellular matrix structure and organization, was also found in this group. A group of 12 up-regulated genes in pyloric caeca, but down-regulated in head kidney and spleen, included three genes which promote apoptosis (KCNN3,EGR and TNFRSF10b) and the immunoglobulin light chain, which plays a key role in the adaptive immune response. The only up-regulated gene in pyloric caeca and head kidney, but down-regulated in spleen was the complement component C3, essential for the activation of the complement pathway. Finally, the only gene up-regulated in spleen and head kidney and down-regulated in pyloric caeca was aminopeptidase n, which encodes an enzyme that participates in the final digestion of peptides, but also in processing other peptide molecules such as hormones, neuropeptides and MHC class II-bound antigen peptides. DE genes classification and organ-specific expression We inspected the list of organ-specific DE transcripts and made a selection of the most interesting genes, which were grouped in five key broad functional categories based on GO term characterisation and on the current literature in the field. The chosen categories were 1) immune and defence response, 2) apoptosis and cell proliferation, 3) iron metabolism and erythropoiseis, 4) metabolism and digestive function and 5) cytoskeleton and extracellular matrix. The selected genes and the group to which they belong to are listed in Additional file 2: Table S2 (ST2) for head kidney, Additional file 3: Table S3 (ST3) for spleen and Additional file 4: Table S4 (ST4) for pyloric caeca. A selection of the 50 most relevant DE genes from these five categories is presented in Table 2 and their expression shown in a heatmap (Figure 3). 100 Spleen 1 A Spleen C Spleen 2 Spleen 3 100 100 Kidney 2 B Kidney C Kidney 1 Kidney 3 100 100 Pyloric 1 C Pyloric C Pyloric 2 Pyloric 3 100 Figure 1 Samples hierarchical clustering by organ. Hierarchical clustering of all diseased and control samples for ASpleen, BHead kidney and CPyloric caeca. Approximately unbiased P-values, computed by multi-scale bootstrap resampling, are displayed on branch nodes and clusters of samples with an approximately unbiased P-value > 0.95 are indicated with a dashed red box, indicating strong support. Robledo et al. BMC Genomics 2014, 15:1149 Page 3 of 17 http://www.biomedcentral.com/1471-2164/15/1149 Immune and defence response A large number of genes related to this function were found in the three organs, but, while the number of down-regulated genes was comparable, pyloric caeca showed the most numerous group of up-regulated genes. This included chemokines, chemokine receptors, immunoglobulin chains, interleukins and several other genes involved in both innate and adaptive immune response. In particular, a broad representation of components of the inflammatory reaction pathway like G-CSFR, Figure 2 DE genes Venn diagrams. Venn diagrams of A) all DE genes, B) up-regulated DE genes, C) down-regulated DE genes in the three organs (head kidney, spleen and pyloric caeca) are shown. The total number of DE expressed genes in each tissue and the number of unique and common genes between them is displayed. Robledo et al. BMC Genomics 2014, 15:1149 Page 4 of 17 http://www.biomedcentral.com/1471-2164/15/1149 Table 2 Selection of the 50 most relevant DE genes Gene Short name Log-FC head kidney Log-FC spleen Log-FC pyloric caeca Category Associated function MHC class I antigen MHC I −6.3 −4.9 −5.4 1 Antigen processing and presentation Interferon regulatory factor 7 IFR7 −1.8 −1.6 −3.6 1 Positive regulation of interferons production Interferon-induced protein 44 IFI44 −3.2 −3.6 −2.2 1 Interferon-inducible protein Interferon-induced protein with tetratricopeptide repeats 1 IFIT-1 −4.7 −5.0 −2.1 1 Interferon-inducible protein Gig2-like protein Gig2 −10.4 −12.8 −7.4 1 Interferon-inducible protein Interferon gamma receptor alpha chain IFNGR1 −2.7 −4.4 -- 1 Regulation of interferon-gamma-mediated signalling pathway Interferon regulatory factor 3 IFR3 −1.8 −2.1 -- 1 Positive regulation of interferons production Toll-like receptor 9 TLR9 −3.0 −3.0 -- 1 Innate immune response; Positive regulation of interferons production Mannose-binding lectin-associated serine protease-3b MASP3b 7.6 4.9 4.5 1 Complement activation Epidermis-type lipoxygenase 3-like ALOXE3 5,7 6.8 10.3 1 Leukotriene metabolic process L-amino-acid oxidase-like IL4I1 4.8 5.9 9.5 1 Innate immune response Interleukin-17a f-1 IL-17A/F-1 -- -- −3.8 1 Inflammatory response Interleukin-22 IL22 −4.5 -- −1.7 1 Inflammatory response CD83 antigen CD83 -- −1.8 2.5 1 Defence response CD209 antigen-like CD209 -- −4.3 2.2 1 Innate immune response Cytochrome b-245 heavy chain CYBB −4.0 -- 3.7 1 Inflammatory response CC chemokine CCL -- 1.9 4.0 1 Inflammatory response Interleukin-1 receptor accessory protein IL1RAP -- −4.2 2.6 1 Inflammatory response Lipopolysaccharide-induced tumor necrosis factor-alpha factor LITAF -- −1.6 6.2 1,2 Regulation of cytokine production-Apoptotic process Immunoglobulin light chain IGlc −1.8 −2.0 2.7 1 Antigen binding T-cell surface glycoprotein cd4 CD4 -- −1.6 -- 1 T cell receptor signalling pathway T-cell receptor beta chain TCRB -- −2.0 -- 1 T cell receptor signalling pathway Perforin-1-like PRF1 −3.2 -- -- 1 T cell-mediated cytolysis Granzyme A/K GZM-A/K −1.7 -- -- 1 T cell-mediated cytolysis Catalase CAT −2.0 −1.8 −2.5 1 Hydrogen peroxide catabolic process Superoxide dismutase SOD -- -- −2.1 1 Removal of superoxide radicals Glutathione s-transferase theta-1 GSTT1 -- -- −2.0 1 Oxidation-reduction process Caspase-3-like CASP3 -- -- 5.5 2 Positive regulation of apoptotic process Cytochrome c CYTC 2.5 -- 4.3 2 Apoptotic DNA fragmentation TNF receptor-associated factor 2-like TRAF2 -- -- 3.6 2 Regulation of apoptotic process Tumor necrosis factor receptor superfamily member 10b-like TNFRSF10B −3.2 −2.2 2.6 2 Regulation of apoptotic process Hemoglobin subunit alpha-d HBAD −5.2 −3.5 -- 3 Oxygen transport Hemoglobin subunit beta-2 HBB2 −4.7 −3.8 -- 3 Oxygen transport Hemoglobin subunit beta-1 HBB1 −4.6 −3.2 -- 3 Oxygen transport Band 3 anion exchange protein SLC4A1 −6.6 −4.5 -- 3 Erythrocytes differentiation Gata-binding factor 2-like GATA2 −3.5 −2.2 -- 3 Erythrocytes differentiation Mitoferrin-1 SLC25A37 −3.6 −3.5 -- 3 Erythrocytes maturation TAL1 −2,8 −1,6 -- 3 Robledo et al. BMC Genomics 2014, 15:1149 Page 5 of 17 http://www.biomedcentral.com/1471-2164/15/1149 Table 2 Selection of the 50 most relevant DE genes (Continued) T-cell acute lymphocytic leukemia protein 1 Erythrocytes differentiationErythrocytes maturation Hepcidin HEPC 2.7 2.0 -- 3 Iron metabolism Aminopeptidase n ANPEP 2.2 7.3 −9.2 4 Protein metabolic process Intestinal-type alkaline phosphatase 1 ALPI -- -- −11.7 1,4 Metabolic process Acidic mammalian chitinase CHIA -- -- −14.1 4 Carbohydrate metabolic process Apolipoprotein a-iv precursor APOA4 -- -- −6.5 4 Lipoprotein metabolic process Gastric inhibitory polypeptide GIP -- -- −3.1 4 Response to nutrient levels Cocaineand amphetamine-regulated transcript CART -- -- 6.7 4 Negative regulation of appetite Gastrin cholecystokinin-like peptide GAST-CCK -- -- −4.4 4 Digestion Collagen alpha-2 chain COL1A2 −2.0 −1.5 −2.4 5 Extracellular matrix structural constituent Tropomodulin 4 TMOD4 −5.6 −3.9 −2.3 5 Tropomyosin binding Formin-like protein 1-like FMNL1 −4.0 −3.0 −2.8 5 Actin cytoskeleton organization Alpha actin ACTA −3.9 -- −2.7 5 Skeletal muscle fiber development Statistically significant fold changes are shown for each organ for 50 relevant genes associated with enteromyxosis. Categories: 1) Immune and defence response; 2) Apoptosis and cell proliferation; 3) Iron metabolism and erythropoiesis; 4) Metabolism and digestive function; 5) Cytoskeleton and extracellular matrix. Non significant differences have been marked as “--”. Figure 3 Heatmap of 50 selected genes. Heatmap of the fifty selected genes presented in Table 2. Displayed are EdgeR [101] normalized counts for each sample and gene. Sample names are displayed at the bottom of the figure whilst gene symbols are shown to the right and have been hierarchically clustered according to their pearson correlation. The category assigned to each gene is also shown with a color code. Robledo et al. BMC Genomics 2014, 15:1149 Page 6 of 17 http://www.biomedcentral.com/1471-2164/15/1149 immune-responsive gene-1, pand e-selectin, as well as the transcription factors AP-1 and CEBPB, was detected. Moreover, genes such as LITAF, which promotes the expression of the pro-inflammatory cytokine TNF-α,and CYBB, a superoxide-generating enzyme of phagocytes, were up-regulated in pyloric caeca but down-regulated in spleen and head kidney, respectively. Lymphoid organs, spleen and head kidney, showed a similar number of upand down-regulated genes involved in inflammation and acute-phase response, including shared up-regulated genes like hepcidin, heat shock proteins and prostaglandin synthases. On the other hand, several genes related to B and T cells (such as those encoding immunoglobulin molecules and the T cell-related proteins CD4,TCRB, granzyme and perforin) were down-regulated in these organs. Moreover, in spleen there was a depression of two genes considered to be markers for dendritic cells, CD83 and CD209, which were, conversely, up-regulated in pyloric caeca. Spleen and pyloric caeca, in turn, showed a common up-regulation of the C-type lectin MRC1, while other two members of this family were up-regulated, but only in spleen (MASP1) or in pyloric caeca (CLEC4M). Noticeably, more genes related to IFNs were identified among down-regulated genes, like Gig1,IFNGR1 and IFR3 in head kidney and spleen or GVINP1 and IRF4 in spleen and pyloric caeca. Moreover, the TLR9, also involved in defence against intracellular pathogens, was down-regulated in spleen and head kidney. Also, some regulated genes involved in Th17 cells response were detected, particularly the down-regulation of Th17 cytokines: IL22 in head kidney and pyloric caeca and IL17 and its receptor in pyloric caeca. In this organ, it was also remarkable the down-regulation of four genes related to anti-oxidant defence (MRSA,SOD,GSTT1 and TXNDC17) and the up-regulation of several genes involved in tissue remodelling and repair (e.g. MMP13, PLAT,FGF10,VNT). Spleen and head kidney, as well, showed up-regulation of HGF, a cytokine acting in tissue repair but also in modulation of immune response, and of two HGF-related genes (HGFR in spleen and HGFAC in head kidney). Finally, a few genes known to be induced under hypoxic and/or oxidative stress conditions were found to be up-regulated, including the angiopoietin related proteins showing an increase between 4 and 7.4 FC in the three organs, and HIGD1A, the adrenomedullin genes and a cytochrome c oxidase mitochondrial subunit in pyloric caeca and head kidney. Apoptosis and cell proliferation Apoptosis and cell proliferation DE genes were found mainly in pyloric caeca. Several genes participating in the apoptotic process, especially the caspase-3 (FC = 5.5) and cytochrome c (FC = 4.3), which are essential players in the execution phase of apoptosis, were up-regulated. In general, in this organ we found more prothan anti-apoptotic genes, but also other genes involved in cell proliferation, such as the PCNA,FGF10 and cyclins b1 and a2, were up-regulated. In head kidney and spleen a few genes belonging to this group, like the pro-apoptotic cytochrome c in head kidney and clusterininspleenwereup-regulated. Iron metabolism and erythropoiesis In head kidney and spleen, the main lymphohaematopoietic organs in teleosts, we observed down-regulation of several genes related to haematopoiesis. The expression of genes involved in erythrocyte maturation and differentiation, like TAL1,GATA2 and mitoferrin-1, was depressed in both organs, and we also observed a dramatic decrease in the expression of genes related to oxygen transport. For example, haemoglobin subunit alpha-d and haemoglobin subunit beta-2 showed a −5.2 and a −4.7 FC, respectively, in head kidney, and a −3.5 and −3.8 in spleen. The band 3 anion transport protein gene, the major glycoprotein of the erythrocyte membrane, also suffered a notable down-regulation in head kidney (FC = −6.6) and spleen (FC = −4.6). On the other hand, two genes related to iron homeostasis, hepcidin (up-) and ferritin heavy subunit (down-), were regulated in both organs. Metabolism and digestive function This group of genes was analyzed in pyloric caeca to evaluate intestinal function during enteromyxosis, and were mainly down-regulated. Most of these genes showed high expression in the control sample, while its expression was practically undetectable in infected individuals. That was the case of the digestive enzymes CHIA (−14.1), ALPI (−11.7), CYP7A1 (−9.1) or CPO (−6.0). Also the FABP2 and APOA4 genes, involved in lipid metabolism, showed very highly expression in control samples, but −4.1 and −6.5 FCs, respectively, in infected samples. Moreover, there was a depression of genes induced by food intake (GIP,CCK2, gastrin-cholecystokinin-like peptide) and of the gene coding for the galanin type I receptor, an orexigenic petide. On the other hand, two anorexigenic genes (CARTand CGRP) were up-regulated. Cytoskeleton and extracellular matrix Several myosin, collagen, actin, tubulin, coronin and spectrin genes, were down-regulated in the three organs. Pyloric caeca and head kidney exhibited the highest number of down-regulated genes. Of particular interest were collagen alpha-1, alpha actin and the different myosin genes, that were abundantly expressed in control samples, showing FCs ranging from −6.1 (myosin heavy chain) to −1.8 (collagen alpha-1). Likewise, the TPM4, Robledo et al. BMC Genomics 2014, 15:1149 Page 7 of 17 http://www.biomedcentral.com/1471-2164/15/1149 which was highly expressed in spleen and pyloric caeca of CTRL fish, was down-regulated in RCPT samples. GO enrichment analysis The full transcriptome of the three organs was annotated and GO terms for each sequence were obtained. A Fisher exact test (FDR corrected p-value = 0.05) was used to compare DE sequences with the background transcriptome to obtain the enriched GO terms for each organ (Figure 4). Oxygen binding was clearly overrepresented in both spleen and head kidney of RCPT fish, likely indicating alterations in the erythrocyte machinery, as mentioned earlier. Lipid metabolism and catalytic activity were enriched categories in pyloric caeca, which might evidence problems in the digestive function. Extracellular space or extracellular region GO terms were present in all three organs. GO enrichment was also performed for up-regulated and down-regulated genes separately, obtaining an additional up-regulated GO category in spleen, peptidase activity, and in pyloric caeca, cell cycle. Discussion This is, to our knowledge, the first report of a RNA-seq transcriptomic analysis applied to the study of a fish-parasite model. We investigated turbot at an advanced enteromyxosis stage, selected on the basis of histopathological evaluation. This approach allows the analysis of fish with a more uniform health status, minimizing interindividual variation, and consequently, enabling a more consistent identification of regulated genes on a reduced number of animals. Figure 4 GO terms enrichment. GO enrichment (p < 0.01 FDR corrected) for DE genes in A) head kidney, B) spleen and C) pyloric caeca. The percentage of sequences with the associated GO term present in the full organ transcriptome is shown in blue, while the percentage of sequences with the GO term in the DE gene group is shown in red. Robledo et al. BMC Genomics 2014, 15:1149 Page 8 of 17 http://www.biomedcentral.com/1471-2164/15/1149 This study advances our understanding of how the E. scophthalmi infection develops and the determination of the clinical signs and lesions characteristic of infection in turbot. Figure 5 depicts the cascade of events leading to severe enteromyxosis in turbot considering in particular the transcriptomic changes found in the current study. We found far more DE genes in pyloric caeca (3022) than in either kidney or spleen (1316 and 1377, respectively). This is perhaps not that surprising since the intestine is the target tissue of Enteromyxum spp. infection and the lesions caused by the disease are mostly restricted to the gastrointestinal tract. Additionally, the most characteristic clinical signs of the disease, such as weight loss and anorexia, are attributable to the alteration of the normal physiology of the digestive system. Furthermore, spleen and kidney shared 321 down-regulated genes, an interesting result which can be attributed to the common lymphohaematopoietic functions and the cellular depletion observed in both organs in late stages of the disease [8]. Figure 5 Enteromyxosis flowchart. Flowchart showing the main processes involved in severe turbot enteromyxosis. The flowchart has been divided in five sections according to the most representative processes occurring during Enteromyxum infection: I, blue, immune and defence response; II, purple, apoptosis and cell proliferation; III, red, iron metabolism and erythropoiesis; IV, green, metabolism and digestive function; V, orange, cytoskeleton and extracellular matrix. Robledo et al. BMC Genomics 2014, 15:1149 Page 9 of 17 http://www.biomedcentral.com/1471-2164/15/1149 maximus) along the infection process with Aeromonas salmonicida using an immune-enriched oligo-microarray. Mar Biotechnol 2011, 13(6):1099–114. 21. Pardo BG, Millán A, Gómez-Tato A, Fernández C, Bouza C, Álvarez-Dios JA, Cabaleiro S, Lamas J, Leiro JM, Martínez P: Gene expression profiles of spleen, liver, and head kidney in turbot (Scophthalmus maximus) along the infection process with Philasterides dicentrarchi using an immune-enriched oligo-microarray. Mar Biotechnol 2012, 14(5):570–82. 22. Rodríguez-Ramilo ST, Fernández J, Toro MA, Bouza C, Hermida M, Fernández C, Pardo BG, Cabaleiro S, Martínez P: Uncovering QTL for resistance and survival time to Philasterides dicentrarchi in turbot (Scophthalmus maximus). Anim Genet 2013, 44(2):149–57. 23. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res 2008, 18(9):1509–17. 24. Morozova O, Hirst M, Marra MA: Applications of new sequencing technologies for transcriptome analysis. Annu Revi Genom Hum G 2009, 10(1):135–51. 25. Nielsen R, Paul JS, Albrechten A, Song YS: Genotype and SNP calling from next-generation sequencing data. Nat Rev Genet 2011, 12(6):443–51. 26. Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 2009, 10(1):57–63. 27. Redondo MJ, Palenzuela O, Álvarez-Pellitero P: Studies on transmission and life cycle of Enteromyxum scophthalmi (Myxozoa), an enteric parasite of turbot Scophthalmus maximus.Folia Parasit 2004, 51(2–3):188–98. 28. Dahl MR, Thiel S, Matsushita M, Fujita T, Willis AC, Christensen T, Vorup-Jensen T, Jensenius JC: MASP-3 and its association with distinct complexes of the mannan-binding lectin complement activation pathway. Immunity 2001, 15:127–35. 29. Hajela K, Kojima M, Ambrus G, Wong KH, Moffatt BE, Ferluga J, Hajela S, Gal P, Sim RB: The biological functions of MBL-associated serine proteases (MASPs). Immunobiology 2002, 205:467–75. 30. Gadjeva M, Takahashi K, Thiel S: Mannan-binding lectin–a soluble pattern recognition molecule. Mol Immunol 2004, 41:113–21. 31. Fujita T, Matsushita M, Endo Y: The lectin-complement pathway–its role in innate immunity and evolution. Immunol Rev 2004, 198:185–202. 32. Kania PW, Sorensen RR, Koch C, Brandt J, Kliem A, Vitved L, Hansen S, Skjodt K: Evolutionary conservation of mannan-binding lectin (MBL) in bony fish: identification, characterization and expression analysis of three bona fide collectin homologues of MBL in the rainbow trout (Onchorhynchus mykiss). Fish Shellfish Immunol 2010, 29:910–20. 33. Gruden-Movsesijan A, Petrovic M, Sofronic-Milosavljevic L: Interaction of mannan-binding lectin with Trichinella spiralis glycoproteins, a possible innate immune mechanism. Parasite Immunol 2003, 25:545–52. 34. Ambrosio AR, De Messias-Reason IJ: Leishmania (Viannia) braziliensis: interaction of mannose-binding lectin with surface glycoconjugates and complement activation. An antibody-independent defence mechanism. Parasite Immunol 2005, 27:333–40. 35. Evans-Osses I, Ansa-Addo EA, Inal JM, Ramirez MI: Involvement of lectin pathway activation in the complement killing of Giardia intestinalis. Biochem Bioph Res Co 2010, 395:382–6. 36. Cestari I, Evans-Osses I, Schlapbach LJ, de Messias-Reason I, Ramirez MI: Mechanisms of complement lectin pathway activation and resistance by trypanosomatid parasites. Mol Immunol 2013, 53:328–34. 37. Sahu A, Lambris JD: Structure and biology of complement protein C3, a connecting link between innate and acquired immunity. Immunol Rev 2001, 180:35–48. 38. Cuesta A, Munóz P, Rodríguez A, Salinas I, Sitjà-Bobadilla A, Álvarez-Pellitero P, Esteban MA, Meseguer J: Gilthead seabream ( Sparus aurata L.) innate defence against the parasite Enteromyxum leei (Myxozoa). Parasitology 2006, 132:95–104. 39. van Vliet SJ, Garcia-Vallejo JJ, van Kooyk Y: Dendritic cells and C-type lectin receptors: coupling innate to adaptive immune responses. Immunol Cell Biol 2008, 86:580–7. 40. Vázquez-Mendoza A, Carrero JC, Rodríguez-Sosa M: Parasitic infections: a role for C-type lectins receptors. Biomed Res Int 2013, 2013:456352. 41. Geijtenbeek TB, Torensma R, van Vliet SJ, van Duijnhoven GC, Adema GJ, van Kooyk Y, Figdor CG: Identification of DC-SIGN, a novel dendritic cell-specific ICAM-3 receptor that supports primary immune responses. Cell 2000, 100:575–85. 42. Lugo-Villarino G, Balla KM, Stachura DL, Bañuelos K, Werneck MBF, Traver D: Identification of dendritic antigen-presenting cells in the zebrafish. Proc Natl Acad Sci U S A 2010, 107(36):15850–5. 43. Y-h H, Zhang M, Sun L: Expression of Scophthalmus maximus CD83 correlates with bacterial infection and antigen stimulation. Fish Shellfish Immunol 2010, 29(4):608–14. 44. Rauta PR, Nayak B, Das S: Immune system and immune responses in fish and their role in comparative immunity study: a model for higher organisms. Immunol Lett 2012, 148:23–33. 45. Sitjà-Bobadilla A, Palenzuela O, Riaza A, Macias MA, Álvarez-Pellitero P: Protective acquired immunity to Enteromyxum scophthalmi (Myxozoa) is related to specific antibodies in Psetta maxima (L.) (Teleostei). Scand J Immunol 2007, 66(1):26–34. 46. Keller AD, Maniatis T: Identification and characterization of a novel repressor of beta-interferon gene expression. Gene Dev 1991, 5(5):868–79. 47. Le Bon A, Tough DF: Links between innate and adaptive immunity via type I interferon. Curr Opin Immunol 2002, 14(4):432–6. 48. Schroder K, Hertzog PJ, Ravasi T, Hume DA: Interferon-gamma: an overview of signals, mechanisms and functions. J Leukocyte Biol 2004, 75(2):163–89. 49. McCall MBB, Sauerwein RW: Interferon-γ—central mediator of protective immune responses against the pre-erythrocytic and blood stage of malaria. J Leukocyte Biol 2010, 88(6):1131–43. 50. Zou J, Secombes CJ: Teleost fish interferons and their role in immunity. Dev Comp Immunol 2011, 35(12):1376–87. 51. Rappocciolo G, Birch J, Ellis SA: Down-regulation of MHC class I expression by equine herpesvirus-1. J Gen Virol 2003, 84(2):293–300. 52. Vincent IE, Zannetti C, Lucifora J, Norder H, Protzer U, Hainaut P, Zoulim F, Tommasino M, Trépo C, Hasan U, Chemin I: Hepatitis B Virus Impairs TLR9 Expression and Function in Plasmacytoid Dendritic Cells. PLoS One 2011, 6(10):e26315. 53. Gainey MD, Rivenbark JG, Cho H, Yang L, Yokoyama WM: Viral MHC class I inhibition evades CD8+ T-cell effector responses in vivo but not CD8+ T-cell priming. Proc Natl Acad Sci U S A 2012, 109(47):E3260–7. 54. Song W, Kao W, Zhai A, Qian J, Li Y, Zhang Q, Zhao H, Hu Y, Li H, Zhang F: Borna disease virus nucleoprotein inhibits type I interferon induction through the interferon regulatory factor 7 pathway. Biochem Bioph Res Co 2013, 438(4):619–23. 55. Babu S, Blauvelt CP, Kumaraswami V, Nutman TB: Diminished expression and function of TLR in lymphatic filariasis: a novel mechanism of immune dysregulation. J Immunol 2005, 175(2):1170–6. 56. Asteal F: A study of Immune responses to L. mexicana Antigens and Immunogenicity of L. donovani Centrin-3. In PhD thesis: Nottingham Trent University, School of Science and Technology. Nottingham; 2011. 57. Srivastava S, Pandey SP, Jha MK, Chandel HS, Saha B: Leishmania expressed lipophosphoglycan interacts with Toll-like receptor (TLR)-2 to decrease TLR-9 expression and reduce anti-leishmanial responses. Clin Exp Immunol 2013, 172(3):403–9. 58. Young ND, Cooper GA, Nowak BF, Koop BF, Morrison RN: Coordinated down-regulation of the antigen processing machinery in the gills of amoebic gill disease-affected Atlantic salmon (Salmo salar L.). Mol Immunol 2008, 45(9):2581–97. 59. Peck A, Mellins ED: Precarious balance: Th17 cells in host defense. Infect Immun 2010, 78(1):32–8. 60. van de Veerdonk FL, Gresnigt MS, Kullberg BJ, van der Meer JW, Joosten LA, Netea MG: Th17 responses and host defense against microorganisms: an overview. BMB Rep 2009, 42(12):776–87. 61. Fujino S, Andoh A, Bamba S, Ogawa A, Hata K, Araki Y, Bamba T, Fujiyama Y: Increased expression of interleukin 17 in inflammatory bowel disease. Gut 2003, 52(1):65–70. 62. Kelly MN, Kolls JK, Happel K, Schwartzman JD, Schwarzenberger P, Combe C, Moretto M, Khan IA: Interleukin-17/interleukin-17 receptor-mediated signaling is important for generation of an optimal polymorphonuclear response against Toxoplasma gondii infection. Infect Immun 2005, 73(1):617–21. 63. Zenewicz LA, Yancopoulos GD, Valenzuela DM, Murphy AJ, Stevens S, Flavell RA: Innate and adaptive interleukin-22 protects mice from inflammatory bowel disease. Immunity 2008, 29(6):947–57. 64. Glass EJ: The molecular pathways underlying host resistance and tolerance to pathogens. Front Genet 2012, 3:263. 65. Wynne JW, O’Sullivan MG, Cook MT, Stone G, Nowak BF, Lovell DR, Elliott NG: Transcriptome analyses of amoebic gill disease-affected Atlantic salmon (Salmo salar) tissues reveal localized host gene suppression. Mar Biotechnol 2008, 10(4):388–403. Robledo et al. BMC Genomics 2014, 15:1149 Page 16 of 17 http://www.biomedcentral.com/1471-2164/15/1149 66. Bates JM, Akerlund J, Mittge E, Guillemin K: Intestinal alkaline phosphatase detoxifies lipopolysaccharide and prevents inflammation in zebrafish in response to the gut microbiota. Cell Host Microbe 2007, 2(6):371–82. 67. Bol-Schoenmakers M, Fiechter D, Raaben W, Hassing I, Bleumink R, Kruijswijk D, Maijoor K, Tersteeg-Zijderveld M, Brands R, Pieters R: Intestinal alkaline phosphatase contributes to the reduction of severe intestinal epithelial damage. Eur J Pharmacol 2010, 633(1–3):71–7. 68. Lee C, Chun J, Hwang SW, Kang SJ, Im JP, Kim JS: The effect of intestinal alkaline phosphatase on intestinal epithelial cells, macrophages and chronic colitis in mice. Life Sci 2014, 100(2):118–24. 69. Keklikoglu N, Koray M, Kocaelli H, Akinci S: iNOS expression in oral and gastrointestinal tract mucosa. Digest Dis Sci 2008, 53(6):1437–42. 70. Hansen R, Thomson JM, El-Omar EM, Hold GL: The role of infection in the aetiology of inflammatory bowel disease. J Gastroenterol 2010, 45(3):266–76. 71. Marchiando AM, Shen L, Graham WV, Edelblum KL, Duckworth CA, Guan Y, Montrose MH, Turner JR, Watson AJ: The epithelial barrier is maintained by in vivo tight junction expansion during pathologic intestinal epithelial shedding. Gastroenterology 2011, 140(4):1208–18. e1201-1202. 72. Watson AJ, Hughes KR: TNF-alpha-induced intestinal epithelial cell shedding: implications for intestinal barrier function. Ann NY Acad Sci 2012, 1258:1–8. 73. Panaro MA, Cianciulli A, Mitolo V, Mitolo CI, Acquafredda A, Brandonisio O, Cavallo P: Caspase-dependent apoptosis of the HCT-8 epithelial cell line induced by the parasite Giardia intestinalis.FEMS Immunol Med Mic 2007, 51(2):302–9. 74. Lovegrove FE, Gharib SA, Patel SN, Hawkes CA, Kain KC, Liles WC: Expression microarray analysis implicates apoptosis and interferon-responsive mechanisms in susceptibility to experimental cerebral malaria. Am J Pathol 2007, 171(6):1894–903. 75. Wyllie AH: “Where, O Death, Is Thy Sting?”A Brief Review of Apoptosis Biology. Mol Neurobiol 2010, 42(1):4–9. 76. Bienvenu A-L, González-Rey E, Picot S: Apoptosis induced by parasitic diseases. Parasite Vector 2010, 3(1):106. 77. James ER, Green DR: Manipulation of apoptosis in the host parasite interaction. Trends Parasitol 2004, 20(6):280–7. 78. Castellanos-Gonzalez A, Yancey LS, Wang HC, Pantenburg B, Liscum KR, Lewis DE, White AC Jr: Cryptosporidium infection of human intestinal epithelial cells increases expression of osteoprotegerin: a novel mechanism for evasion of host defenses. J Infect Dis 2008, 197(6):916–23. 79. Okunishi K, Dohi M, Nakagome K, Tanaka R, Mizuno S, Matsumoto K, Miyazaki J, Nakamura T, Yamamoto K: A novel role of hepatocyte growth factor as an immune regulator through suppressing dendritic cell function. J Immunol 2005, 175(7):4745–53. 80. Benkhoucha M, Santiago-Raber M-L, Schneiter G, Chofflon M, Funakoshi H, Nakamura T, Lalive PH: Hepatocyte growth factor inhibits CNS autoimmunity by inducing tolerogenic dendritic cells and CD25 + Foxp3+ regulatory T cells. Proc Natl Acad Sci U S A 2010, 107(14):6424–9. 81. Benkhoucha M, Molnarfi N, Schneiter G, Walker P, Lalive P: The neurotrophic hepatocyte growth factor attenuates CD8+ cytotoxic T-lymphocyte activity. J Neuroinflammation 2013, 10(1):154. 82. Ganz T: Hepcidin and iron regulation, 10 years later. Blood 2011, 117(17):4425–33. 83. Neves JV, Caldas C, Wilson JM, Rodrigues PNS: Molecular mechanisms of hepcidin regulation in sea bass (Dicentrarchus labrax). Fish Shellfish Immunol 2011, 31(6):1154–61. 84. Pereiro P, Figueras A, Novoa B: A novel hepcidin-like in turbot (Scophthalmus maximus L.) highly expressed after pathogen challenge but not after iron overload. Fish Shellfish Immunol 2012, 32(5):879–89. 85. Ganz T: The role of hepcidin in iron sequestration during infections and in the pathogenesis of anemia of chronic disease. Isr Med Assoc J 2002, 4(11):1043–5. 86. Tizard IR: Veterinary Immunology: An Introduction, 8th Edition. India: CBS Publishers & Distributors; 2009. 87. Estensoro I, Benedito-Palos L, Palenzuela O, Kaushik S, Sitjà-Bobadilla A, Pérez-Sánchez J: The nutritional background of the host alters the disease course in a fish–myxosporean system. Vet Parasitol 2011, 175(1–2):141–50. 88. Salem M, Silverstein J, Rexroad CE 3rd, Yao J: Effect of starvation on global gene expression and proteolysis in rainbow trout (Oncorhynchus mykiss). BMC Genomics 2007, 8:328. 89. Pérez-Jiménez A, Cardenete G, Hidalgo MC, García-Alcázar A, Abellán E, Morales AE: Metabolic adjustments of Dentex dentex to prolonged starvation and refeeding. Fish Physiol Biochem 2012, 38(4):1145–57. 90. Antonopoulou E, Kentepozidou E, Feidantsis K, Roufidou C, Despoti S, Chatzifotis S: Starvation and re-feeding affect Hsp expression, MAPK activation and antioxidant enzymes activity of European Sea Bass (Dicentrarchus labrax). Comp Biochem Phys A 2013, 165(1):79–88. 91. Hayes J, Volkoff H: Characterization of the endocrine, digestive and morphological adjustments of the intestine in response to food deprivation and torpor in cunner, Tautogolabrus adspersus. Comp Biochem Phys A 2014, 170:46–59. 92. Kyriazakis I, Tolkamp BJ, Hutchings MR: Towards a functional explanation for the occurrence of anorexia during parasitic infections. Anim Behav 1998, 56(2):265–74. 93. Barber I: Parasites, behaviour and welfare in fish. Appl Anim Behav Sci 2007, 104(3–4):251–64. 94. Colditz IG: Six costs of immunity to gastrointestinal nematode infections. Parasite Immunol 2008, 30(2):63–70. 95. Wykes LJ, Fiorotto M, Burrin DG, Del Rosario M, Frazer ME, Pond WG, Jahoor F: Chronic low protein intake reduces tissue protein synthesis in a pig model of protein malnutrition. JNutr1996, 126(5):1481–8. 96. Lenaerts K, Sokolovic M, Bouwman FG, Lamers WH, Mariman EC, Renes J: Starvation induces phase-specific changes in the proteome of mouse small intestine. J Proteome Res 2006, 5(9):2113–22. 97. Redondo MJ, Palenzuela O, Riaza A, Macias A, Álvarez-Pellitero P: Experimental transmission of Enteromyxum scophthalmi (Myxozoa), an enteric parasite of turbot Scophthalmus maximus.J Parasitol 2002, 88(3):482–8. 98. Lohse M, Bolger AM, Nagel A, Fernie AR, Lunn JE, Stitt M, Usadel B: RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Res 2012, 40(W1):W622–7. 99. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg S: TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol 2013, 14(4):R36. 100. Langmead B, Salzberg SL: Fast gapped-read alignment with Bowtie 2. Nat Meth 2012, 9(4):357–9. 101. Robinson M, Oshlack A: A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol 2010, 11:R25. 102. Suzuki R, Shimodaira H: Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics 2006, 22:1540–2. 103. Edgar R, Domrachev M, Lash AE: Gene Expression Omnibus: NCBI gene expression and hybridization data repository. Nucleic Acids Res 2002, 30(1):207–10. doi:10.1186/1471-2164-15-1149 Cite this article as: Robledo et al.:RNA-seq analysis reveals significant transcriptome changes in turbot (Scophthalmus maximus) suffering severe enteromyxosis. BMC Genomics 2014 15:1149. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit Robledo et al. BMC Genomics 2014, 15:1149 Page 17 of 17 http://www.biomedcentral.com/1471-2164/15/1149