Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Transcriptomic response to parasite infection in Nile tilapia (Oreochromis niloticus) depends on rearing density

BMC Genomics volume 19, Article number: 723 (2018) Cite this article

Abstract

Background

Captive animal populations, be it for food production or conservation programmes, are often maintained at densities far beyond those in natural environments, which can have profound effects on behaviour, immune and stress levels, and ultimately welfare. How such alterations impact transcriptional responses to pathogen infection is a ‘different kettle of fish’ and remains poorly understood. Here, we assessed survival and gene expression profiles of infected fish reared at two different densities to elucidate potential functional genomic mechanisms for density-related differences in disease susceptibility.

Results

Utilising a whole-transcriptome sequencing (RNAseq) approach, we demonstrate that rearing density in tilapia (Oreochromis niloticus) significantly impacts susceptibility to the oomycete Saprolegnia parasitica, via altered transcriptional infection responses. Tilapia held at low densities have increased expression of genes related to stress, likely due to increased aggressive interactions. When challenged with Saprolegnia, low-density fish exhibit altered expression of inflammatory gene responses and enhanced levels of adaptive immune gene suppression compared to fish reared at higher density, resulting in significantly increased mortality rates. In addition, Saprolegnia infection substantially perturbs expression of circadian clock genes, with fish reared at low-density having higher levels of molecular clock dysregulation.

Conclusions

Our results reveal the wide-scale impact of stocking density on transcriptional responses to infection and highlight the need to incorporate circadian biology into our understanding of disease dynamics in managed animals.

Background

Increasing global human populations are escalating spatial conflicts for food production. As competition for land and water grows, food security and sustainable intensification of food production remain at the forefront of worldwide political agendas [1]. In the past, rearing densities in intensive food animal production were primarily optimised for maximal growth and/or productivity, and hence economic benefit [2]. Today, however, there is increasing consideration of animal rearing density on behaviour [3, 4], stress levels [5, 6], immune levels [5, 7], and ultimately welfare [3, 6].

Intensification of aquaculture increasingly relies on maintaining fish stocks at densities vastly different to those found in natural populations. Stocking/rearing density and disease management are considered pivotal factors determining fish farm productivity and profitability [8,9,10]. In fish, sub-optimal rearing density can increase cortisol levels, likely due to altered social interactions [11, 12]. Wide-scale transcriptional differences due to stocking conditions (potentially due to chronic social stress) appear to include a number of immunologically important pathways [7]. However, critically, it is unknown how these transcriptome-level alterations translate to actual infection responses and thus disease susceptibility. Examining the effects of rearing density on disease susceptibility in fish at a transcriptomic level will not only provide critical insight for effective disease management in aquaculture, but also for other managed animal populations such as terrestrial livestock and captive breeding programmes for conservation.

Infections due to Saprolegnia oomycetes result in significant economic losses in fish farming [13]. S. parasitica, the species most commonly associated with fish saprolegniasis, typically elicits strong host inflammatory responses, yet is capable of actively suppressing several components of the adaptive immune system including T-helper cell cytokines, antigen presentation machinery and immunoglobulins [14]. For tilapia aquaculture, second only to carp in global economic importance [15], devastating outbreaks typically occur in outdoor rearing facilities during cooler months, resulting in high mortality and severe economic losses [16, 17]. Although malachite green is highly effective at controlling saprolegniasis, since its worldwide ban in 2002 there are limited effective treatment options, leading to a global resurgence of the disease [18]. It is therefore essential that conditions influencing susceptibility to Saprolegnia infection are elucidated to effectively manage risks and reduce reliance on pharmaceutical interventions in aquaculture.

Extreme high fish densities increase plasma cortisol in tilapia, which is typically taken as evidence for chronic stress [19]. However, in some fish species, including tilapia, increased levels of aggression between individuals are observed at lower densities [11, 20] which may also act as a potent chronic social stressor. We hypothesise that a sub-optimal, low rearing density in tilapia will significantly impact the expression of immune gene pathways, leading to a decreased resistance to infection challenge. In this study, we experimentally challenged juvenile Nile tilapia (Oreochromis niloticus), reared under two densities, to S. parasitica and compared infection resistance and skin and gills transcriptome-wide gene expression (quantified via RNAseq) between uninfected and infected fish. We use our data to characterise the wider functional genomic effects of stocking density on disease responses and make recommendations for the improvement of approaches used in the management of disease dynamics in captive animals.

Results

Stocking density and disease-challenge survival

Susceptibility to Saprolegnia infection (i.e. mortality due to infection) was examined in tilapia raised under two density conditions (high-density = 6 kg m− 3, low-density = 1.5 kg m− 3). No mortality was observed in sham-challenged (uninfected/control) fish in either high- (n = 80) or low-density (n = 20) treatments for the duration of the experiment. Overall mortality was 10.2 and 52.1% in high- (n = 29/285) and low-density (n = 38/73) Saprolegnia-challenged treatments respectively. Stepwise survival model selection resulted in a final model including density, length and their interaction, excluding body condition and its interactions with other variables (Additional file 1). In the final model, host survival following Saprolegnia challenge was significantly lower in fish maintained at low density (Cox Hazard Ratio HR = 28.631, P < 0.001; Fig. 1a). In addition, smaller fish were significantly less likely to survive infection (HR = 23.355, P < 0.001) and there was a significant interaction between density treatment and fish length (HR = 0.151, P = 0.026; Fig. 1b). Maintaining fish at low density significantly increased average fish length (Welch’s t = − 2.28, P = 0.023).

Fig. 1
figure 1

a Cox’s proportional hazards survival model of Oreochromis niloticus challenged with Saprolegnia parasitica at two stocking densities. Shading indicates 95% CI bands. b Logistic regression predicting probability of survival during Saprolegnia infection by fish standard length. Shading indicates 95% CI bands. c Example of Saprolegnia infection on Oreochromis niloticus

Full size image

Differential gene expression in uninfected fish

To examine the transcriptomic effects of stocking density in uninfected fish, we compared gene expression between sham-challenged (uninfected/control) fish maintained at high and low density at two time points (24 and 48 h post-exposure). In the gill, 2672 and 1268 genes were differentially expressed between uninfected high- and low-density fish at 24 and 48 h respectively. Genes with significantly higher expression levels in low-density fish at both timepoints (n = 154) were enriched for biological processes related to vitamin transport (e.g. GO:0051180: vitamin transport, GO:0032218: riboflavin transport) and included folate transporter 1 (folt) and riboflavin transporter 2 (rft2). Biosynthesis/metabolism genes (e.g. GO:0044210: de novo CTP biosynthetic process, GO:0006541: glutamine metabolic process) were also more highly expressed, including CTP synthase 1 (ctps1) and amidophosphoribosyltransferase (ppat). The most significantly enriched gene ontology (GO) terms of genes more highly expressed at low density were found to be related to circadian processes (e.g. GO:0032922: circadian regulation of gene expression, GO:0009648: photoperiodism), including period (per3) and casein kinase I isoform delta (csnk1d). In addition, upregulated genes in low-density fish gills included “steroid hormone signalling” (GO:0003707, e.g. nr1d2, paqr7) and “response to stress” (GO:0006950, e.g. hspa9, hspd1, trap1). Genes with higher expression in the gills of fish maintained at high density (n = 165) were enriched for biological processes related to lipid and nitrogen compound transport (e.g. GO:0015917: aminophospholipid transport, GO:0010876: lipid localisation), including phospholipid-transporting ATPase (atp8a2) and oxysterol binding like protein 2 (osbpl2). A complete list of enriched GO terms from gill tissues is provided in Additional file 2.

In the skin, 6171 and 2810 genes were differentially expressed between uninfected high- and low-density fish at 24 h and 48 h respectively. Genes exhibiting higher expression levels in low-density fish at both timepoints (n = 474) were highly enriched for processes related to metabolism and biosynthesis (e.g. GO:0006629: lipid metabolic process, GO:0016053: organic acid biosynthetic process), including several lipid synthases (fasn, isyna1, pdss1, ptgs1). Similar to the gills, genes with higher expression in low-density fish also included several GO terms related to circadian processes (GO:0009648: photoperiodism, GO:0007623: circadian rhythm; e.g. per3). Genes with higher expression in the skin of fish maintained at high density (n = 644) were enriched for biological processes related to T cells (e.g. GO:0045580 regulation of T cell differentiation; gata3), skin structure (e.g. GO:0098773 skin epidermis development; nrg1) and protein ubiquitination (GO:0016567; rnf19a). A complete list of enriched GO terms from skin tissues is provided in Additional file 3.

Genes previously identified as markers of stress and/or cortisol production in fish and other vertebrates were specifically compared between uninfected fish at high and low density in both skin and gill tissues. In the gills, corticotropin-releasing factor receptor 1 (crhr1) was more highly expressed in low density fish at 24 h. In the skin of low density fish, metallothionein (mt) showed higher expression and glucocorticoid receptor (nr3c1) reduced expression at 24 h (Fig. 2). Furthermore, two skin gene co-expression modules (SK4, SK13; Additional file 4) significantly associated with density (higher expression in low-density fish) included significant enrichment of GO terms related to stress responses such as “cellular response to stress” (SK4; e.g. dmap1) and “cortisol secretion” (SK13; e.g. ptpn11).

Fig. 2
figure 2

Differential expression of stress-related genes in uninfected Oreochromis niloticus at two stocking densities. Asterisks denote significant differences in expression (contrasts between tissues and/or time points not shown)

Full size image

Gene expression responses to Saprolegnia infection shared between densities

The numbers of genes found to be differentially expressed among control (uninfected) and Saprolegnia-challenged fish and their GO term enrichments are summarised in Additional files 5 and 6. In gill tissues, a total of 175 genes were more highly expressed in high and low density infected fish at both time points (24 and 48 h post-exposure). This gene set was highly enriched for metabolic processes including glucose (GO:0006006) and lipid (GO:0006629) metabolism (e.g. hk2, eno1). Genes downregulated in all infected fish groups were predominantly enriched by GO terms related to circadian processes (Additional file 6), but also included “MHC class II protein binding” (GO:0042289; e.g. cd74). Shared early responses (24 h) to Saprolegnia included increased gene expression related to osmotic stress (GO:0006970, e.g. slc31a1), epidermal repair (GO:0048730, e.g. tp63), mucous secretion (GO:0070254, e.g. p2ry2), and response to fungus (GO:0009620, e.g. ncf1). Eight gene co-expression modules were significantly associated with infection, but not density (Additional file 8). Three of the four co-expression modules (G1, G4, G5) negatively correlated with infection status (i.e. lower expression in infected fish) and were enriched for immune response terms including “antigen processing and presentation” (GO:0019882, e.g. cd83, MHC II alpha and beta chains) and “interleukin-12 production” (GO:0032615, e.g. mast2) (Additional file 8). Two of the four modules positively correlated with infection status were enriched for a number of innate and adaptive immune response terms (Additional file 8), including both T cell and mast cell processes (e.g. rac1, rc3h1).

In the skin, genes with decreased expression in infected individuals at both densities were highly enriched for immune functions. The most significant GO term for the 75 genes found to exhibit lower expression in infected fish in both density treatments and time points was “immune response” (GO:0006955) and were also enriched for “interferon-gamma-mediated signalling pathway” (GO:0060333) Additional file 10. This gene set included suppressor of cytokine signalling (socs1), beta-2-microglobulin (b2m), T-cell surface glycoprotein CD4 (cd4), and C-C motif chemokine 17 (ccl17). In addition, genes sharing decreased expression in infected fish in both density treatments, during at least one time point, were enriched for a number of adaptive immune terms including “MHC class II protein complex” (GO:0042613, e.g. HLA class 2 gamma chain, H-2 class 2 alpha chain) and “T cell activation” (GO:0042110, e.g. pag1, gata3) (Additional file 5). The 148 genes found to have higher expression in the skin of infected fish from both time points and density treatments were not significantly enriched for immune related functions. However, shared early responses to infection (24 h post-exposure in both densities) included increased expression of genes associated with inflammation and wound healing (e.g. mast cell migration, epithelial cell migration, interleukin-1 receptor activity; Additional file 5). In the weighted gene co-expression network analysis (WGCNA), no association between co-expression gene modules and density treatment was found, although overall nine gene modules were significantly associated with infection (Additional file 7). Modules positively correlated with infection status (i.e. increased expression in infected individuals) included enrichment for mast cell mediated immune pathways (SK10, SK11, SK14).

Density-specific responses to Saprolegnia infection

To examine the potential transcriptomic underpinnings for the differences in disease susceptibility between stocking densities, we identified differentially expressed genes unique to each density treatment. In the gills, enrichment of inflammatory response regulation was found in genes with significantly increased expression in low-density infected fish at 24 h (e.g. tnip1, sbno2, Additional file 6). In contrast, increased gene expression at 24 h unique to high-density fish was enriched for a number of GO terms related to immunoglobulin-mediated immune responses (e.g. msh2, batf, Additional file 6). Furthermore, we found enrichment of GO terms related to T cell functions in genes with increased expression only in infected fish at high-density at both timepoints (e.g. trpm4, Additional files 6 and 13). Downregulated gill genes in low-density fish, at both timepoints, were enriched for terms related to major histocompatibility complex (MHC) II (e.g. MHC class II beta chain), complement and humoral immune responses (e.g. masp1, Additional file 6). Only a single gene co-expression module was significantly correlated with both infection status (lower expression in infected fish) and density treatment (lower expression in low-density fish) in the gills (G9; Additional files 8 and 11). This module of 127 genes (including per3, csnk1d, ciart) was predominantly enriched for GO terms related to circadian functions (e.g. circadian regulation of gene expression), indicating significant disruption of molecular clock expression due to infection (Fig. 3).

Fig. 3
figure 3

a Differential expression of circadian rhythm genes between uninfected and Saprolegnia-infected Oreochromis niloticus. Positive (green) and negative (red) log2 fold changes indicate higher and lower expression in infected fish respectively. Zero change (white) indicates no significant difference in expression. b Network view of G9 gill module (circadian regulation of gene expression) hub genes. Nodes are labelled with official gene symbols when available. Edge line width represents connection strength (weight); thicker lines denote stronger connections. Node size is proportional to number connections. Node colour represents the log2 fold change between control and infected fish at 48 h post-exposure. Node borders represent differential gene expression results. Grey = nonsignificant; red = decreased expression in the infected; green = increased expression in the infected

Full size image

In the skin, similar to early responses of the gill, responses to infection unique to low density fish at 24 h (n = 1965) were enriched for increased expression of inflammatory pathways (e.g. il17c, Additional files 5 and 9). In contrast, early responses (24 h) of infected high-density fish (n = 905) were enriched for metallopeptidase activity (e.g. mmp28, Additional file 5). In addition, while suppression of immune pathways, particularly adaptive immune responses, were found in the skin at both densities, only the low-density group were found to have decreased collagen biosynthetic activity (e.g. p4ha1, Additional file 5). Gene co-expression modules SK1 and SK6, enriched for adaptive immune responses (SK1; e.g. klhl6, trpm4) and wound healing processes (SK6; e.g. cask, Additional file 7), were negatively correlated with infection status and positively associated with density, indicating a lower level of immunosuppression in the skin of high density fish.

Discussion

Here we demonstrate that sub-optimal stocking levels significantly alter fish disease responses at the transcriptional level, and result in higher disease susceptibility. Specifically, we found that tilapia raised at low densities are likely under chronic social stress [11, 20] (evidenced by increased expression of markers related to cortisol production and other stress responses), have altered expression of inflammatory responses and higher levels of pathogen immunosuppression. In addition to extensive perturbation of immune pathways, rearing density and infection appear to dysregulate the expression of molecular body clock genes, demonstrating potentially profound and complex effects of rearing conditions that must be considered for managed animal disease mitigation.

Saprolegnia parasitica is a generalist pathogen, capable of infecting a broad taxonomic range of freshwater teleosts [21]. To date, transcriptional changes due to infection have only been elucidated in salmonid host species [14, 22, 23]. We found the core responses of tilapia to Saprolegnia challenge (i.e. changes in gene expression common to both stocking densities) are broadly similar to those highlighted in salmonid studies; such as early increases in acute inflammation expression including interleukin-mediated skin inflammatory genes (e.g. il1rl2, Additional files 5 and 6). The pathogenicity of Saprolegnia is likely to be due, at least in part, to its ability to suppress the host’s adaptive immune response [14]. Here, we found evidence for immunosuppression in infected tilapia. In both the skin and gills, there is decreased expression of genes related to adaptive immunity, including T cell and MHC-related genes (Additional files 5, 6, 7 and 8), suggesting the immunosuppressive activities of Saprolegnia are not limited to salmonid infections.

There were substantial differences in susceptibility to Saprolegnia infection between the two stocking densities, with significantly higher mortality in low density fish (Fig. 1a). We propose underlying differences in expression of genes between the two stocking densities significantly contribute to the observed variance in Saprolegnia susceptibility. Although we reveal increased expression of some pro-inflammatory pathways in fish at both densities, early responses unique to low-density fish included increased expression of genes related to downregulation of inflammation (e.g. tnip1, sgk1, mst1) in both tissues and a decrease in pro-inflammatory cytokine (interleukin-12, IL12) production in the skin (Additional files 5 and 6). In mammals, stress is known as a potent disruptor of pro-inflammatory responses, including glucocorticoid inhibition of IL12 production [24]. IL12 is critical to producing Th1 responses and interferon-gamma (INFγ) [25], both considered important for clearing Saprolegnia and other oomycete infections [18]. Taken alone, these results would suggest that raising tilapia at low stocking density dampens their ability to mount an inflammatory response to Saprolegnia challenge, resulting in higher mortality rates. However, low-density fish also exhibited increased expression of interleukin-17 (IL17), a pro-inflammatory cytokine produced by Th17 cells [26]. Therefore, low stocking density (and thus perhaps stress) may not universally decrease inflammatory responses, but instead induce an inflammatory pathway less effective against oomycete infections.

In addition to differences in inflammatory responses to Saprolegnia between fish held at different stocking densities, we revealed important variation in adaptive immune gene expression. Genes with increased expression in the gills of only high-density (less susceptible) fish were enriched for immunoglobulin production and B cell activation at 24 h. Increased expression of genes such as batf and msh2 suggest immunoglobulin class-switch recombination [27, 28] plays a role in effective defences against Saprolegnia infections. T cell mediated immunity was enriched in genes with increased expression at both time points only in high-density fish (Additional file 6). These genes included several transient receptor potential cation channel subfamily M 4 members (e.g. trpm4), key to regulating Th responses [29], providing further support for the impact of rearing density on T helper cell activity. While there was evidence of the immunosuppressive activity of Saprolegnia in infected fish from both density treatments, these results suggest that fish at high density (with likely lower levels of stress) are less susceptible to adaptive immune pathways disruption, resulting in lower mortality rates. Future work to examine how the gene groups and pathways we identify here are perturbed under various chronic stressors (e.g. temperature, pollutants etc) and other immunosuppressive pathogen species will be highly informative for understanding the impacts of chronic stressors on disease dynamics in both wild and captive animal populations.

Perhaps the most striking finding of this study was the extent to which expression of genes regulating circadian rhythms was affected by stocking density and infection (Fig. 3a). The vast majority of genes central to vertebrate circadian rhythm regulation [30] were differentially expressed between uninfected and infected fish. Genes of both the negative (e.g. arntl and clock) and positive (e.g. per and cry) arms of circadian transcriptional loops, plus associated regulators (e.g. rora), exhibited decreased expression in infected fish. However, the extent of variation in circadian expression profiles due to infection differed by stocking density, with greater changes observed in fish at low density (Fig. 3b). Moreover, some of the most significant differences in gene expression between stocking densities in unchallenged fish are associated with circadian rhythm functions. In humans and other mammalian models, there is increasing recognition of the intricate links between stress, disease, and circadian rhythms [31]. Despite this, little or no consideration has been given to the implications for managed animal populations. Here, we provide the first indication that, in fish, management practices may substantially impact circadian functioning and thus immunity to pathogens. Time-series of expression during infection (i.e. to determine if circadian gene expression is continually suppressed during infection or expression oscillations are shifted) and further studies under different levels of stress will be highly valuable to our understanding of the impacts of rearing conditions on animal immunity. For aquaculture, where photoperiod manipulation is a common management tool [32, 33], we posit it is essential optimal husbandry conditions are implemented to minimise stress, thus decreasing disease susceptibility and ultimately increasing not only welfare, but economic productivity.

Conclusions

This study demonstrates the complex functional genomic effects of rearing density on fish immune responses and disease susceptibility. We find that densities that likely induce increased social stress in tilapia result in altered inflammatory responses to Saprolegnia, increased immunosuppression of adaptive immune pathways by the pathogen, and thus higher mortality due to disease. By employing a transcriptome-wide approach, we revealed a high level of dysregulation of molecular clock genes due to stress and during infection. Our results provide not only an important in-depth transcriptomic resource for future studies of rearing stress and immunity in fish, but also highlights the need to incorporate circadian infection biology into our understanding of disease dynamics in managed animals.

Methods

Tilapia rearing conditions

Four hundred and sixty 3-week-old mixed-sex Nile tilapia fry (standard length 11 to 29 mm, mean = 19.2 mm, weight 0.31 to 4.29 g, mean = 1.54 g) were obtained from a commercial facility (FishGen Ltd.). Fry were visually and microscopically determined free of parasitic infections upon arrival. Fry were randomly allocated to glass aquaria (90 × 30 × 38 cm) and reared for 9 weeks under two stocking conditions. In the low density (LD) treatment, fry (n = 94) were maintained at 1.5 kg m− 3. In the high density (HD) treatment, fry (n = 366) were maintained at 6 kg m− 3, thus maintaining 4 times higher density than LD treatment throughout rearing. Both density treatments are within typical fry stocking density regimes in commercial tilapia aquaculture [34]. Fish were fed ad libitum twice daily with trout pellets (Skretting Nutra Parr) and maintained under 12:12 h light: dark conditions at 24 °C (± 0.5 °C) with external filtration. Water quality measurements (ammonia, nitrite, nitrate) and 50% water changes were performed every other day in both treatments to ensure equivalent abiotic conditions (ammonia, nitrite = 0 ppm, nitrate = 20–30 ppm).

Experimental infections

Saprolegnia parasitica (species confirmed with rDNA ITS sequencing according to [35], isolated in 2015 from a naturally infected Salmo salar), were used for experimental infections. Mycelial cultures were maintained at 24 °C on potato dextrose agar plates. Prior to zoospore production, mycelial plugs were transferred to glucose yeast extract broth and incubated for 3 days at 24 °C. Mycelia masses were washed in sterile dechlorinated water and then maintained for 3 days in 50% dechlorinated water, 50% water from tilapia stock tanks, at 15 °C to induce zoospore production. Heat-killed zoospores were produced by incubating at 60 °C for 1 h. A haemocytometer was used to quantify zoospore concentrations and zoospore suspensions were equilibrated to 24 °C before use [36].

In both density conditions, fish were divided into challenged (HD; n = 285, LD; n = 73) and sham-challenged control (HD; n = 80, LD; n = 20) groups, maintaining 1.5 kg m− 3 and 6 kg m− 3 in low- and high-density treatments respectively. The sample size for Saprolegnia-challenged groups was based on previous Saprolegnia differential susceptibility experiments using other stressors [37]. To avoid the overriding effects of acute stress due to confinement and individual isolation [38], fish were simultaneously challenged with S. parasitica in these treatment groups. Immediately prior to exposure, all fish were shaken in a hand net in groups of five for 1 min to facilitate infection [39]. In the high-density treatment, live zoospores were added directly to aquaria to achieve a concentration of 5 × 106 zoospores L− 1. To ensure low density fish were exposed to equivalent 1) number of infective zoospores per individual and 2) concentration of organic matter, a mixture of 1:3 live: heat-killed zoospores was added directly to each aquarium to achieve a concentration of 5 × 106 zoospores L− 1. Filtration was turned off during exposure and zoospore solutions were completely changed in all tanks every 6 h to minimise waste accumulation. For uninfected controls, sham-challenged fish were also net-shaken, returned to holding aquaria, and water changed every 6 h to mimic the Saprolegnia-exposed groups. After 24 h, water was completely changed in all tanks and filtration turned back on.

For the 24-h exposure period and duration of infection experiment, fish were visually inspected hourly (under red light during dark periods) for signs of lethargy, reduced swimming ability and mycelial growth. Fish showing severe signs of infection were immediately euthanised with an overdose of Tricaine Methanesulfonate (MS222, 500 mg L− 1) according to Home Office Schedule 1. All fish were weighed, standard length measured, and body condition calculated (Fulton’s Condition Index; K = weight/standard-length3). At 24 and 48 h post-exposure, five fish per treatment combination (high/low density, challenged/sham-challenged) were euthanised with an overdose of MS222 (500 mg L− 1) and samples of gill (all arches) and skin (standardised right and left flank area immediately posterior to opercula) were immediately preserved in RNAlater and stored at − 80 °C until RNA extraction. Skin and gill tissues were chosen as they are primary sites of infection of Saprolegnia. Individuals selected for gene expression analyses were size-matched among density and infection treatments. For infected fish, to ensure gene expression signals were not simply due severe physiological stress near death, individuals were selected showing mycelial growth but no other signs of pathology (e.g. loss of swimming/righting ability, lethargy, hypo/hyperventilation). All tissue samples used for gene expression analyses were taken at the same time of day (11.00 am). In the absence of reliable methods to quantify Saprolegnia infection burdens, here we refer to disease “susceptibility/resistance” as mortality/survival post-challenge. Statistical analyses of survival were performed in R. The effect of stocking density, fish size, and body condition on Saprolegnia challenge survival was examined using Cox proportional hazards (R package survival [40]). Continuous variables (length, condition) were scaled and centred on the mean. The initial model predicted survival using all variables (density treatment, length, condition) and their two-way interactions. Stepwise model selection on Akaike information criterion (AIC) was then performed. The interaction of density and length on probability of survival was further examined using plots of binary logistic regressions.

Transcriptome sequencing

Total RNA was extracted from each tissue sample separately using AllPrep DNA/RNA Micro kit (Qiagen). RNA was quantified using Qubit High-Sensitivity RNA assays (ThermoFisher Scientific) and integrity was determined using Agilent 4200 Tapestation RNA assays (Agilent Technologies). All samples had RNA integrity values greater than 8.0. Libraries were generated using the Illumina Stranded TruSeq mRNA sample preparation kits (high-throughput protocol), as per the manufacturer’s instructions (Illumina, San Diego, CA). Libraries’ quality was quantified and assessed using Agilent 4200 Tapestation prior equimolar pooling. The library pool was run four times on Illumina NextSeq500 (2 × 75 PE) to achieve a minimum of 20 million read pairs per sample. Raw reads are available at the NCBI Short Read Archive under Accession Number PRJNA413738.

Gene expression analyses

Reads were filtered and trimmed using Trimmomatic version 0.32 [41] to: (1) trim any Illumina adapter sequence, (2) trim the 5′ and/or 3′ end of reads where quality score dropped below Q20, (3) trim anywhere within each read where a 5-bp window drops below Q20, and (4) discard any trimmed reads less than 36 bp long. This ensured only the highest quality reads were used for subsequent analyses. Trimmed reads were then mapped to the O. niloticus genome (v1.1 [42]) using HISAT2 version 2.0.5 [43] and quantified using RSEM version 1.2.30 [44]. A summary of sequence read data and mapping is presented in Table 1. Gene expression profiles were visualised using principal component analysis (PCA) plots of z-scaled and centred read counts to ensure no individual outliers driving expression differences [45] (Additional file 12). Differential expression tests were performed using the R package DESeq2 [46], comparing 1) high- and low-density control (uninfected) fish at both time points (24 and 48 h post-exposure), and 2) infected and uninfected (control) fish at both time points (24 and 48 h post-exposure). Transcripts were filtered to include only those with at least two counts per million mapped reads (TPM) in at least two individuals. Although fish were size matched between groups, as fish size was significant in models of survival, length was included as a covariate in expression testing and tissues were analysed separately. False discovery rate (FDR) corrected P values of less than 0.05 were considered for significant differential expression. Overlap of differentially expressed genes between treatment groups were determined using Venny version 2.1 [47]. Gene ontology (GO) functional enrichment tests (with FDR correction) were carried out via Blast2GO [48] to detect significantly overrepresented biological processes, molecular functions, and cellular components of groups of differentially expressed genes shared/unique to particular treatment groups.

Table 1 Summary of Oreochromis niloticus transcriptome sequence data
Full size table

Weighted gene co-expression network analyses

Differential gene expression (DGE) is usually identified using exact tests carried out on each gene separately; however, due to the need of correcting p-values with stringent multiple testing methods, only genes with the largest differences in expression are typically identified. An alternative for quantifying transcriptional responses is weighted gene co-expression network analysis (WGCNA [49]). This method was used to quantify transcriptional responses of fish to stocking density and infection challenge by Saprolegnia, enabling the identification of networks (modules) of co-expressed genes (i.e. genes that show consistent expression profiles across samples), and thus potentially identifying functionally important genes with only subtle changes in expression that may otherwise (i.e. via DGE) not have be detected. Trimmed mean of M-values (TMM) normalised fragments per kilobase per million mapped expression values were analysed using the R package WGCNA [49]. Our gene modules were defined using the dynamicCutTree function and TOMType “signed” with a minimum module size of 100. A module eigengene distance threshold of 0.25 was used to merge highly similar modules. Modules were considered robust if the average module adjacencies were significantly higher than observed in a set of 1000 random module permutations. Modules were then correlated with time, stocking density, infection and fish length to identify gene networks significantly associated with factors of interest. GO term enrichment tests of each significant gene module were performed using Blast2Go as described above. Each gene within a module was ranked by its module membership (kME), calculated by WGCNA. Network hub genes were defined as those ranked in the top 100 module membership values and with the highest 150 network connection weights. Hub gene network connections were exported to Cytoscape [50] for visualization.

Abbreviations

AIC:

Akaike information criterion

bp:

Base pair

DGE:

Differential gene expression

FDR:

False discovery rate

GO:

Gene ontology

HD:

High-density

ITS:

Internal transcribed spacer

LD:

Low-density

MHC:

Major histocompatibility complex

MS222:

Tricaine Methanesulfonate

PCA:

Principal component analysis

RNAseq:

Whole-transcriptome RNA sequencing

Th:

T helper

TMM:

Trimmed mean of M-values

WGCNA:

Weight gene co-expression network analysis

References

  1. Godfray HCJ, Beddington JR, Crute IR, Haddad L, Lawrence D, Muir JF, et al. Food security: the challenge of feeding 9 billion people. Science. 2010;327(5967):812–8.

    Article  CAS  Google Scholar 

  2. Thornton PK. Livestock production: recent trends, future prospects. Philos Trans R Soc Lond Ser B Biol Sci. 2010;365(1554):2853–67.

    Article  Google Scholar 

  3. Campbell, D.L.M., , G. N. Hinch, J. A. Downing and C. Lee. Outdoor stocking density in free-range laying hens: effects on behaviour and welfare. Animal 2017;11(6):1036–1045.

    Article  CAS  Google Scholar 

  4. Cornale P, Macchi E, Miretti S, Renna M, Lussiana C, Perona G, et al. Effects of stocking density and environmental enrichment on behavior and fecal corticosteroid levels of pigs under commercial farm conditions. J Vet Behav. 2015;10(6):569–76.

    Article  Google Scholar 

  5. Yarahmadi P, Miandare HK, Fayaz S, Caipang CMA. Increased stocking density causes changes in expression of selected stress-and immune-related genes, humoral innate immune parameters and stress responses of rainbow trout (Oncorhynchus mykiss). Fish Shellfish Immunol. 2016;48:43–53.

    Article  CAS  Google Scholar 

  6. Calabrese S, Nilsen TO, Kolarevic J, Ebbesson LOE, Pedrosa C, Fivelstad S, et al. Stocking density limits for post-smolt Atlantic salmon (Salmo salar L.) with emphasis on production performance and welfare. Aquaculture. 2017;468:363–70.

    Article  Google Scholar 

  7. Sun P, Bao P, Tang B. Transcriptome analysis and discovery of genes involved in immune pathways in large yellow croaker (Larimichthys crocea) under high stocking density stress. Fish Shellfish Immunol. 2017;68:332–40.

    Article  CAS  Google Scholar 

  8. Engle CR, McNevin A, Racine P, Boyd CE, Paungkaew D, Viriyatum R, et al. Economics of sustainable intensification of aquaculture: evidence from shrimp farms in Vietnam and Thailand. J World Aquac Soc. 2017;48(2):227–39.

    Article  Google Scholar 

  9. Johnson K, Engle C, Wagner B. Comparative economics of US catfish production strategies: evidence from a cross-sectional survey. J World Aquac Soc. 2014;45(3):279–89.

    Article  Google Scholar 

  10. Bondad-Reantaso MG, Subasinghe RP, Arthur JR, Ogawa K, Chinabut S, Adlard R, et al. Disease and health management in Asian aquaculture. Vet Parasitol. 2005;132(3):249–72.

    Article  Google Scholar 

  11. Ellis T, North B, Scott AP, Bromage NR, Porter M, Gadd D. The relationships between stocking density and welfare in farmed rainbow trout. J Fish Biol. 2002;61(3):493–531.

    Article  Google Scholar 

  12. Turnbull J, Bell A, Adams C, Bron J, Huntingford F. Stocking density and welfare of cage farmed Atlantic salmon: application of a multivariate analysis. Aquaculture. 2005;243(1):121–32.

    Article  Google Scholar 

  13. van West P. Saprolegnia parasitica, an oomycete pathogen with a fishy appetite: new challenges for an old problem. Mycologist. 2006;20(3):99–104.

    Article  Google Scholar 

  14. Belmonte R, Wang T, Duncan GJ, Skaar I, Mélida H, Bulone V, et al. Role of pathogen-derived cell wall carbohydrates and prostaglandin E2 in immune response and suppression of fish immunity by the oomycete Saprolegnia parasitica. Infect Immun. 2014;82(11):4518–29.

    Article  Google Scholar 

  15. Naylor RL, Goldburg RJ, Primavera JH, Kautsky N. Effect of aquaculture on world fish supplies. Nature. 2000;405(6790):1017.

    Article  CAS  Google Scholar 

  16. Saad TT, Atallah ST, El-Bana SA. Fish diseases and its economic effect on Egyptian fish farms. J Agri Food Tech. 2014;4(5):1–6.

    Google Scholar 

  17. Chauhan R. Fungal attack on Tilapia mossambicus in culture pond, leading to mass mortality of fishes. Int J Phram Sci Rev Res. 2014;7.

  18. Phillips AJ, Anderson VL, Robertson EJ, Secombes CJ, van West P. New insights into animal pathogenic oomycetes. Trends Microbiol. 2008;16(1):13–9.

    Article  CAS  Google Scholar 

  19. Barcellos LJ, Nicolaiewsky S, De Souza SM, Lulhier F. The effects of stocking density and social interaction on acute stress response in Nile tilapia Oreochromis niloticus (L.) fingerlings. Aquac Res. 1999;30(11–12):887–92.

    Article  Google Scholar 

  20. Adams CE, Turnbull JF, Bell A, Bron JE, Huntingford FA. Multiple determinants of welfare in farmed fish: stocking density, disturbance, and aggression in Atlantic salmon (Salmo salar). Can J Fish Aquat Sci. 2007;64(2):336–44.

    Article  Google Scholar 

  21. Gozlan RE, Marshall WL, Lilje O, Jessop CN, Gleason FH, Andreou D. Current ecological understanding of fungal-like pathogens of fish: what lies beneath? Front Microbiol. 2014;5.

  22. Roberge C, Páez DJ, Rossignol O, Guderley H, Dodson J, Bernatchez L. Genome-wide survey of the gene expression response to saprolegniasis in Atlantic salmon. Mol Immunol. 2007;44(6):1374–83.

    Article  CAS  Google Scholar 

  23. de Bruijn I, Belmonte R, Anderson VL, Saraiva M, Wang T, van West P, et al. Immune gene expression in trout cell lines infected with the fish pathogenic oomycete Saprolegnia parasitica. Dev Comp Immunol. 2012;38(1):44–54.

    Article  CAS  Google Scholar 

  24. Glaser R, Kiecolt-Glaser JK. Stress-induced immune dysfunction: implications for health. Nat Rev Immunol. 2005;5(3):243.

    Article  CAS  Google Scholar 

  25. Hsieh CS, Macatonia SE, Tripp CS, Wolf SF, O'Garra A, Murphy KM. Development of TH1 CD4+ T cells through IL-12 produced by listeria-induced macrophages. Science. 1993;260(5107):547–9.

    Article  CAS  Google Scholar 

  26. Jin W, Dong C. IL-17 cytokines in immunity and inflammation. Emerg Microbes Infect. 2013;2(9):e60.

    Article  CAS  Google Scholar 

  27. Ehrenstein MR, Neuberger MS. Deficiency in Msh2 affects the efficiency and local sequence specificity of immunoglobulin class-switch recombination: parallels with somatic hypermutation. EMBO. 1999;18(12):3484–90.

    Article  CAS  Google Scholar 

  28. Morman RE, Schweickert PG, Konieczny SF, Taparowsky EJ. BATF regulates the expression of Nfil3, Wnt10a and miR155hg for efficient induction of antibody class switch recombination in mice. Eur J Immunol. 2018. https://doi.org/10.1002/eji.201747360.

    Article  CAS  Google Scholar 

  29. Weber KS, Hildner K, Murphy KM, Allen PM. Trpm4 differentially regulates Th1 and Th2 function by altering calcium signaling and NFAT localization. J Immunol. 2010;185(5):2836–46.

    Article  CAS  Google Scholar 

  30. Bell-Pedersen D, Cassone VM, Earnest DJ, Golden SS, Hardin PE, Thomas TL, et al. Circadian rhythms from multiple oscillators: lessons from diverse organisms. Nat Rev Genet. 2005;6(7):544–56.

    Article  CAS  Google Scholar 

  31. Dumbell R, Matveeva O, Oster H. Circadian clocks, stress, and immunity. Front Endocrinol. 2016;7.

  32. Boeuf G, Le Bail P-Y. Does light have an influence on fish growth? Aquaculture. 1999;177(1):129–52.

    Article  Google Scholar 

  33. Taranger GL, Carrillo M, Schulz RW, Fontaine P, Zanuy S, Felip A, et al. Control of puberty in farmed fish. Gen Comp Endocrinol. 2010;165(3):483–515.

    Article  CAS  Google Scholar 

  34. Ronald, N., B. Gladys, and E. Gasper. The effects of stocking density on the growth and survival of Nile tilapia (Oreochromis niloticus) fry at son fish farm, Uganda. J Aquac Res Dev 2014;5(2).

  35. Diéguez-Uribeondo J, Fregeneda-Grandes JM, Cerenius L, Pérez-Iniesta E, Aller-Gancedo JM, Tellería MT, et al. Re-evaluation of the enigmatic species complex Saprolegnia diclina–Saprolegnia parasitica based on morphological, physiological and molecular data. Fung Genet Biol. 2007;44(7):585–601.

    Article  Google Scholar 

  36. Stewart A, Jackson J, Barber I, Eizaguirre C, Paterson R, van West P, et al. Hook, line and infection: a guide to culturing parasites, establishing infections and assessing immune responses in the three-spined stickleback. Adv Parasitol. 2017;98:39.

    Article  Google Scholar 

  37. Carballo M, Munoz MJ, Cuellar M, Tarazona JV. Effects of waterborne copper, cyanide, ammonia, and nitrite on stress parameters and changes in susceptibility to saprolegniosis in rainbow trout (Oncorhynchus mykiss). Appl Environ Microbiol. 1995:2108–12.

  38. Auperin B, Baroiller J-F, Ricordel M-J, Fostier A, Prunet P. Effect of confinement stress on circulating levels of growth hormone and two prolactins in freshwater-adapted tilapia (Oreochromis niloticus). Gen Comp Endocrinol. 1997;108(1):35–44.

    Article  CAS  Google Scholar 

  39. Hoshiai G. Studies on saprolegniasis in cultured coho salmon, Oncorhynchus kisutch Walbaum. Bull Nippon Vet Zoo Coll. 1990;39:154–7.

    Google Scholar 

  40. Therneau T. Package ‘survival’. 2017. https://github.com/therneau/survival. Accessed 15 May 2018.

  41. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  Google Scholar 

  42. Brawand D, Wagner CE, Li YI, Malinsky M, Keller I, Fan S, et al. The genomic substrate for adaptive radiation in African cichlid fish. Nature. 2014;513(7518):375–81.

    Article  CAS  Google Scholar 

  43. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  Google Scholar 

  44. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12(1):323.

    Article  CAS  Google Scholar 

  45. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protocols. 2013;8(8):1494–512.

    Article  CAS  Google Scholar 

  46. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  Google Scholar 

  47. Oliveros JC. VENNY. An interactive tool for comparing lists with Venn Diagrams. 2007. http://bioinfogp.cnb.csic.es/tools/venny/index.html. Accessed 15 May 2018.

  48. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

    Article  CAS  Google Scholar 

  49. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9(1):559.

    Article  Google Scholar 

  50. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  Google Scholar 

Download references

Acknowledgments

We thank Emily Matthews for assisting in tissue collection and members of the Cable lab group for fish husbandry. Transcriptomic data analyses were made possible using High Performance Computing Wales, a collaboration between Welsh Universities, the Welsh Government and Fujitsu.

Funding

Funding for this study was provided by the Welsh Government and Higher Education Funding Council for Wales through the Sêr Cymru National Research Network for Low Carbon, Energy and the Environment (NRN-LCEE) AquaWales project.

Availability of data and materials

All sequence data have been submitted to the NCBI Sequence Read Archive (Accession: PRJNA413738) and will be made publicly available upon acceptance of this manuscript for publishing. All other data are available as additional files with this article.

Author information

Authors and Affiliations

  1. School of Biosciences, Cardiff University, Cardiff, CF10 3AX, UK

    Amy R Ellison, Pablo Orozco-terWengel & Jo Cable

  2. Biosciences Department, Swansea University, Swansea, SA2 8PP, UK

    Tamsyn M Uren Webster, Olivier Rey, Carlos Garcia de Leaniz & Sofia Consuegra

  3. Present address: Université Perpignan Via Domitia, IHPE UMR 5244, CNRS, IFREMER, University Montpellier, F-66860, Perpignan, France

    Olivier Rey

Authors
  1. Amy R Ellison
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tamsyn M Uren Webster
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Olivier Rey
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Carlos Garcia de Leaniz
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Sofia Consuegra
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Pablo Orozco-terWengel
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jo Cable
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

All authors designed the study. SC, CG, JC and POW organized funding. AE, OR and TUW. collected data. AE performed analyses. AE wrote the manuscript with contributions and edits from all authors. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Amy R Ellison.

Ethics declarations

Ethics approval and consent to participate

All procedures were performed under Home Office project license PPL 303424 with full approval of Cardiff University Animal Ethics committee. Consent to participate: not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Cox Proportional Hazards model selection: Results of Cox Proportional Hazards models for evaluating the effect of density on survival during Saprolegnia infection, showing the model selection process for the best fit of the data according to AIC values. (DOCX 14 kb)

Additional file 2:

Gene Ontology enrichment of uninfected gill tissues: Full GO enrichment results of gill genes differentially expressed between uninfected fish at high and low density. (XLSX 6567 kb)

Additional file 3:

Gene Ontology enrichment of uninfected skin tissues: Full GO enrichment results of skin genes differentially expressed between uninfected fish at high and low density. (XLSX 3920 kb)

Additional file 4:

Gene Ontology enrichment of skin gene modules: Full GO enrichment results of gill gene co-expression modules. (PDF 80 kb)

Additional file 5:

Summary of skin gene ontology enrichment: Summary of gene ontology (GO) term enrichment of skin genes differentially expressed between control (uninfected) and Saprolegnia-infected Oreochromis niloticus density treatment groups (LD; low density, HD; high density), including total number of differentially expressed genes, most significant GO term, and major biological process clusters determined using ReViGO. ↑ denotes increased expression in infected fish and ↓ denotes decreased expression. (DOCX 38 kb)

Additional file 6:

Summary of gill gene ontology enrichment: Summary of gene ontology (GO) term enrichment of gill genes differentially expressed between control (uninfected) and Saprolegnia-infected Oreochromis niloticus density treatment groups (LD; low density, HD; high density), including total number of differentially expressed genes, most significant GO term, major biological process clusters determined using ReViGO. ↑ denotes increased expression in infected fish and ↓ denotes decreased expression. (XLSX 2288 kb)

Additional file 7:

Summary of skin gene module enrichment: Summary of gene ontology (GO) term enrichment of skin gene co-expression modules defined by WGCNA, including total number of genes per module, most significant GO term, and major biological process clusters determined using ReViGO. Significant correlations with infection status, density treatment, time point, and fish standard length indicated (+; positive correlation, −; negative correlation). (XLSX 2803 kb)

Additional file 8:

Summary of gill gene module enrichment: Summary of gene ontology (GO) term enrichment of gill gene co-expression modules defined by WGCNA, including total number of genes per module, most significant GO term, and major biological process clusters determined using ReViGO. Significant correlations with infection status, density treatment, time point, and fish standard length indicated (+; positive correlation, −; negative correlation). (XLSX 5834 kb)

Additional file 9:

Gene Ontology enrichment of infected gill tissues: Full GO enrichment results of gill genes differentially expressed between infected and uninfected fish at high and low density. (DOCX 21 kb)

Additional file 10:

Gene Ontology enrichment of infected skin tissues: Full GO enrichment results of skin genes differentially expressed between infected and uninfected fish at high and low density. (DOCX 21 kb)

Additional file 11:

Gene Ontology enrichment of gill gene modules: Full GO enrichment results of gill gene co-expression modules. (DOCX 17 kb)

Additional file 12:

Principal Component Analysis plots of gene expression: Principal components analysis (PCA) plots of total gene expression profiles in A) gill, and B) skin of uninfected (CTL) and Saprolegnia-challenged (INF) Oreochormis niloticus at two stocking densities (H; high, L; low) and time points (24; 24 h post-exposure, 48; 48 h). (DOCX 16 kb)

Additional file 13:

Dot plots of GO term enrichment: Summary of GO term enrichment in differentially expressed genes comparing uninfected and Saprolegnia-challenged fish. Dot colour intensity indicates FDR-corrected P value, dot size indicates number of genes. (XLSX 7015 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ellison, A.R., Uren Webster, T.M., Rey, O. et al. Transcriptomic response to parasite infection in Nile tilapia (Oreochromis niloticus) depends on rearing density. BMC Genomics 19, 723 (2018). https://doi.org/10.1186/s12864-018-5098-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5098-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us