Next Article in Journal
Nationwide Desert Highway Assessment: A Case Study in China
Previous Article in Journal
Work and High-Risk Alcohol Consumption in the Canadian Workforce
Previous Article in Special Issue
Forage as a Primary Source of Mycotoxins in Animal Diets
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Identification of Candidate Genes and Physiological Pathways Involved in Gonad Deformation in Whitefish (Coregonus spp.) from Lake Thun, Switzerland

1
Computational and Molecular Populations Genetics Lab, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland
2
Liverpool Microarray Facility, School of Biological Sciences, University of Liverpool, L69 7ZB Liverpool, UK
3
Centre for Fish and Wildlife Health, University of Bern, Laenggass-Strasse 122, PO-Box 8466, 3001 Bern, Switzerland
4
Institute of Clinical Chemistry, University Hospital, University of Bern, Inselspital, 3010 Bern, Switzerland
*
Author to whom correspondence should be addressed.
Int. J. Environ. Res. Public Health 2011, 8(7), 2706-2733; https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph8072706
Submission received: 28 April 2011 / Revised: 7 June 2011 / Accepted: 15 June 2011 / Published: 30 June 2011
(This article belongs to the Special Issue Advances in Environmental Biology)

Abstract

:
In 2000, fishermen reported the appearance of deformed reproductive organs in whitefish (Coregonus spp.) from Lake Thun, Switzerland. Despite intensive investigations, the causes of these abnormalities remain unknown. Using gene expression profiling, we sought to identify candidate genes and physiological processes possibly associated with the observed gonadal deformations, in order to gain insights into potential causes. Using in situ-synthesized oligonucleotide arrays, we compared the expression levels at 21,492 unique transcript probes in liver and head kidney tissue of male whitefish with deformed and normally developed gonads, respectively. The fish had been collected on spawning sites of two genetically distinct whitefish forms of Lake Thun. We contrasted the gene expression profiles of 56 individuals, i.e., 14 individuals of each phenotype and of each population. Gene-by-gene analysis revealed weak expression differences between normal and deformed fish, and only one gene, ictacalcin, was found to be up-regulated in head kidney tissue of deformed fish from both whitefish forms, However, this difference could not be confirmed with quantitative real-time qPCR. Enrichment analysis on the level of physiological processes revealed (i) the involvement of immune response genes in both tissues, particularly those linked to complement activation in the liver, (ii) proteolysis in the liver and (iii) GTPase activity and Ras protein signal transduction in the head kidney. In comparison with current literature, this gene expression pattern signals a chronic autoimmune disease in the testes. Based on the recent observations that gonad deformations are induced through feeding of zooplankton from Lake Thun we hypothesize that a xenobiotic accumulated in whitefish via the plankton triggering autoimmunity as the likely cause of gonad deformations. We propose several experimental strategies to verify or reject this hypothesis.

1. Introduction

In the year 2000, commercial fishermen reported the sudden appearance of deformed reproductive organs in whitefish (Coregonus spp.) from Lake Thun, Switzerland. A detailed first survey of the variation of gonad morphology in whitefish of this lake suggested that 35% of the examined fish had aberrant gonad morphology [1]. A more detailed study [2] indicated that, although a certain proportion of the observed variation in gonad morphology can be attributed to natural variation, the quantity and quality of gonadal deformations in Lake Thun whitefish goes clearly beyond that level of natural variability [1]. Furthermore, Bittner et al. [2] also found that these deformations occurred predominantly at particular spawning sites of different whitefish forms that coexist in sympatry in Lake Thun.
Many studies were initiated in search of potential causes for the unusual gonad deformations. Most of them focused on chemical contaminants, in particular on endocrine-disrupting compounds since these substances can alter gonad morphology in fish [36]. To date, no unequivocal evidence for chemicals as causative agents have been discovered [79]. It appears that neither inbreeding effects nor recent hybridization between genetically distinct whitefish forms are involved [10]. Inheritance of gonad deformations [11] and improper breeding of whitefish at the local hatchery [1] have also been excluded. Thus, despite the intensive research having been conducted to date, the causes of gonad deformations in whitefish from Lake Thun remain unknown with the exception of some rearing experiments that point to a role of the lake plankton [7,11].
In recent years, gene expression profiling by means of microarray analyses has become an increasingly popular tool to assess the impact of potential stressors on biological pathways at the level of mRNA transcription, especially in fish, which are increasingly used as models for studies of environmental genomics [12]. Analyses of the transcriptome provide insight into the molecular control of biological and physiological processes and specifically into response to disease, toxins, environmental challenge and other stressors [13,14]. Many of the stressor exposures for which there exists data have been accomplished in laboratories under controlled conditions. So far, only few microarray studies have been conducted with field samples because of the risk of high noise effects due to variation of sex, genotype, age and intrinsic genetic variability in gene expression in natural populations [1518] that makes it difficult to sort out specific stressor effects. In these studies, usually a gradient approach was used where local fish populations under stress were compared to less impacted sites or to “reference” sites where the stressor was absent. Thus, European flounder (Platichthys flesus) from the polluted Tyne were compared to the unpolluted Alde estuaries in the UK [19]. In another study, transcriptomic profiles of wild Fundulus from several superfund sites (highly polluted) and nearby reference sites along the Atlantic coast were compared [16,20,21]. However, to our knowledge, arrays have not been used to date to compare fish from the same environment but displaying different manifestations of a phenotypic trait.
Here, we applied gene expression profiling to identify new candidate genes and physiological processes associated with observed gonad deformations in whitefish from Lake Thun, and thus, to obtain insight into new potential causes. By using an in-situ synthesized rainbow trout oligonucleotide array [22], we compared expression levels of 21,492 transcripts in head kidney and liver tissue between male whitefish with normal and deformed gonads. Since we aimed at screening for gene expression patterns pointing to particular xenobiotic substances potentially contained in the lake plankton, we chose to study tissues of the main organs with respect to detoxification processes or excretion of toxic compounds. Liver is an important organ for many detoxification processes and the head kidney is a key excretory organ of toxic compounds in fish. Transcription profiles were analyzed with two complementary approaches. First, in order to identify candidate genes, we used an ANOVA-based gene-by-gene approach (cf. [23]) to infer expression differences between groups of fish in pairwise comparisons associated with the deformation of gonads. Secondly, using the Gene Ontology (GO) database, we categorized annotated transcript probes into functionally meaningfully groupings according to molecular functions, biological processes or cellular locations in order to identify physiological pathways that are associated with deformed gonads.

2. Experimental Section

2.1. Sampling

Lake Thun harbours several distinct sympatric whitefish forms [10]. A large-scale epidemiological survey revealed that gonad deformations occurred predominantly in the summer-spawning “Brienzlig” whitefish form and additionally in one population of the winter-spawning “Albock” form [2]. The two forms display significantly different, though partially overlapping distributions in gill rakers counts, which is an important adaptive and highly heritable morphological trait in whitefish [24].
For microarray analyses, we included only male whitefish in our study in order to avoid additional, sex-related variance in gene expression and because males showed significantly higher frequencies of gonad deformations compared to females [2]. We aimed at comparing the two most distinct gonadal phenotypes: (i) fish with near perfectly symmetrical gonads (normal gonads) and (ii) fish showing adhesions/fusions of testis to the peritoneal wall and the lateral trunk musculature (deformed gonads).
Fish were caught with bottom gill nets and their gonad morphology was macroscopically immediately analyzed. For individuals, which fulfilled our criteria of normal and deformed gonads, liver and head kidney samples were extracted, subsequently transferred into RNA stabilization reagent (RNA Later, Ambion) and stored at −30°C prior to RNA extraction. Two hundred and sixty eight (268) male whitefish were sampled from two Brienzlig spawning sites at depths between 100–200 m in September 2005. Fusions (i.e., deformed gonads) were observed in 52 individuals (19.4%). We pooled samples from the two sites, since the two spawning populations were not significantly different in allele frequencies at 11 microsatellite loci [25]. Furthermore, 231 male whitefish were collected from an Albock spawning site at a depth of about 50m in December 2005 as an independent biological replicate. In this sample, 25 fish (10.82%) with fusions were present.
The age of each sampled fish was determined using the annulus-criteria of [26]. Additionally, the first gill arch was removed in order to determine gill raker counts, which were counted from the anterior left gill arch in the laboratory using a binocular. Furthermore, a small liver sample from each fish was stored in absolute ethanol for subsequent genotyping at 11 microsatellite loci as described in [25].

2.2. Experimental Design

In order to control for possible technical and biological effects on individual transcription profiles, in a first step we homogenized the samples with respect to genetic population association and in a second step we matched the composition of the samples with regard to age and gill raker numbers as good as possible. In more detail, assignment tests [27,28] were carried out with the program GENECLASS 2.0 [29] in order to identify outlier genotypes. Individuals with unlikely multilocus genotypes, i.e., P < 0.01 in their respective parental populations were excluded. Outlier individuals with respect to age and gill raker counts were also excluded. Means and standard deviations for the above parameters as well as weight and length measurements of sampled fish are given in Table S1.
The use of the two whitefish forms as independent biological replicates was verified in this study by assessing the expected morphological and genetic differentiation between the two populations. In detail, the Albock and the Brienzlig displayed significant differences in gill raker counts (Table S1) and were highly significantly genetically differentiated, (FST = 0.136; P < 0.00001; 95% confidence interval: 0.064–0.216). In contrast, no genetic and morphological differentiation between normal and deformed fish within both forms was observed (all P > 0.05; Table S1).
In order to estimate the reproducibility of the array experiments, we replicated the analysis of four head kidney and four liver samples on a different array slide using independently extracted and labelled RNA-samples from the same tissue samples. We observed a high mean squared correlation coefficient (R2) between slide replicates of 0.97 for liver tissue and of 0.94 for head kidney tissue, respectively. Reproducibility of spot replicates on a slide was similarly assessed by comparing the raw intensity measurements of spot replicates within all arrays: R2 was 0.99 for liver (range: 0.983–0.993) and for head kidney (range: 0.989–0.994), respectively, indicating the high quality of the arrays used in this study.
We used the two-colour microarray-based approach and carried out non-competitive hybridization experiments. In total we analyzed transcription profiles of 112 samples (56 liver and 56 head kidney samples). Liver samples were labelled with Cy3 and head kidney samples with Cy5 and were then co-hybridized on the same array. Since we used 4 × 44 k arrays hybridizations were carried out simultaneously for four different whitefish samples incorporating liver and head kidney samples adding to a total of eight array experiments on each slide. One sample per group (i.e., Albock deformed, Albock normal, Brienzlig deformed, and Brienzlig normal) was present on each slide.
To avoid hidden systematic biases we carefully randomized all stages of the microarray experiment (RNA extraction, cDNA labelling and hybridizations) by processing different sets of samples of the four population groups. To avoid further systematic biases hybridizations and washings were simultaneously carried out for two batches of seven slides. The technical replicate slide was hybridized and washed individually.

2.3. Microarrays

Since no microarray designed for Coregonus species was available to date, we used an in situ-synthesized oligonucleotide array chip originally designed against rainbow trout targets [22]. Owing solely to sequence divergence, the number of features that a single-species microarray can detect in targets from another species is expected to decrease with increasing phylogenetic divergence. The use of a microarray generated from one species to measure gene expression in another is termed “heterologous hybridization” [30]. We tested hybridization of Coregonus spp. targets to rainbow trout probes in a preliminary experiment to labelled cDNA from two rainbow trout and two whitefish specimens from liver and head kidney tissues respectively on 22 k custom arrays fabricated by Oxford Gene Technologies (OGT, Oxford, UK) [22]. We found a comparable number of expressed genes in both species in liver and head kidney tissue (Table S2). In addition, log2 transformed raw intensity frequency histograms showed also comparable patterns between the two species as well as for both tissues (Figure S1). These results validate the heterologous use of these probes for Coregonus species.
For our array experiment we used 15 Agilent 4 × 44 k custom arrays with four identical arrays containing two copies of the 21,492 unique 60-mer oligonucleotide gene probes (i.e., spot replicate) as well as negative control probes. The performance of these arrays is excellent in that the dynamic range, spot quality and spot-to-spot consistency exceeds that of typical cDNA arrays by some margin [22].

2.4. cRNA Synthesis, Hybridizations, Imaging and Feature Extraction

Total RNA was extracted automatically using the RNeasy Plus kit on the QiaCube robot (both Qiagen) using 10mg of frozen tissue stabilized in RNA Later (Ambion). Quality and quantity assessment of the extracted RNA were performed using a NanoDrop ND-1000 UV-VIS Spectrophotometer v.3.2.1 and Bioanalyzer (Agilent). All RNA samples had RIN values of > 8 (average mean RIN for liver = 9.36; average mean RIN for head kidney = 9.34) with the exception of one liver sample (RIN = 7).
Sample preparation, labelling, hybridization, washing, scanning and feature extraction was performed using Agilent’s two-colour microarray-based gene expression analysis v.5.5 protocol. We used Agilent’s Low RNA Input Linear Amplification Kit Plus to generate fluorescent cRNA. Quality (i.e., integrity) of the labelled cRNA was checked with Agilent’s Bioanalyzer using the RNA 6000 Nano LabChip kit. To assess the quantity of cRNA we used the NanoDrop ND-1000 UV-VIS Spectrometer. Hybridizations were carried out overnight for 17 h at 65 °C.
Slides were washed and then scanned immediately on an Agilent DNA microarray scanner (5 μm resolution) to minimize the impact of environmental oxidants on signal intensities. Finally, gene expression levels were quantified using Agilent’s Feature Extraction software v.9.5.3. Median local background signal intensities were subtracted from raw median signal intensities of each spot. Non-uniform spots, non-uniform local backgrounds as well as saturated spots were flagged by Feature Extraction software. For liver, out of the 60 array experiments a mean of 0.09% (range: 0.016–1.047%) of the relevant 42,984 spots per array were flagged. For head kidney slightly more spots, i.e., 0.372% (0.184–1.04%) per array were flagged. For subsequent analysis, the technical replicate slide was excluded.

2.5. Normalization and Data Filtering

Supported by raw signal intensity distribution of spot replicates, raw data were log2 transformed. In order to make data comparable across slides, we applied a robust non-linear method for normalization using array signal distribution analysis and cubic splines [31]. We then applied a spot replicate filter, where the absolute distance between pairs of spot replicates was accounted for. In cases where features with absolute distance > 4 times the standard deviation of the population mean between each other were found, the outlying feature was discarded. For liver, a mean of 1.59% features (range 0.61–4.27%) per array and for head kidney, a mean of 0.51% features (range 0.16–1.88%) per array were excluded by spot replicate filtering. Then data were normalized again. As background level of unspecific hybridization, we considered the mean of negative control probes plus two standard deviations in each tissue. Only probes with mean signal intensities above this level in at least one of the groups were considered for further statistical treatment. Additionally, we applied a missing data filter, i.e., features with less than 50% of data in at least one of the groups were excluded from further analyses.

2.6. Gene-by-Gene Analysis

In order to identify genes involved in gonad deformation, pairwise hierarchical analysis of variance (ANOVA) between following combinations of the four groups were performed: (i) Albock normal vs. Albock deformed; (ii) Brienzlig normal vs. Brienzlig deformed; (iii) all normal fish (Albock and Brienzlig pooled) vs. all deformed fish and additionally a “form” comparison (iv) all Albock (normal and deformed pooled) vs. all Brienzlig. ANOVAs were calculated separately for each gene in both tissues separately. The following nested ANOVA model was considered:
y i j k = μ + a i + b i j , + c i j k
where yijk is the log2 gene expression intensity for a given gene probe for the k-th spot replicate in the j-th individual in the i-th experimental group. The effects for gonadal phenotype (a), individual (b) and spot replicate (c) are assumed to be independent and random. Therefore total variance can be calculated:
σ T 2 = σ a 2 + σ b 2 + σ c 2
where σa2, σb2and σc2 are the associated variance components. Note that under this setting, σa2 is also equal to the covariance in gene expression levels between two individuals from the same group [32], and σb2 is the covariance in gene expression levels between two spot replicates from the same individual. Therefore, in a way similar to what is done with genetic data, we can define the following intra-class correlation coefficient:
M S T = σ a 2 / σ T 2
MST is the correlation in gene expression between two individuals from the same group relative to two individuals from two different groups. MST can also be considered as the proportion of total variance in gene expression differences that can be explained due to differences between the groups.
Significance levels of gene-by-gene analysis were estimated using 10,000 permutation tests, as described elsewhere for genetic data [33]. If none of the permutations revealed a mean difference larger than the observed, additional permutations were carried out until at least two permutations were present above criteria (this only accounted for very few probes and was only carried out for comparisons of normal vs. deformed fish). Advantage of our ANOVA based approach was to neglect assumptions regarding normal distributions of normalized gene expression levels as well as equal variances between groups and among individuals (see [23]). In order to adjust for multiple hypotheses testing, false discovery rate (FDR) was applied [34].

2.7. Gene Ontology Analysis

Information about genes is increasingly being captured and organized in the Gene Ontology (GO) nomenclature [35]. Enrichment analysis identifies GO categories that are enriched among the genes that are differently expressed between groups. Gene Ontology terms for the 21,492 gene probes were obtained by downloading corresponding GO terms from the Gene Ontology Annotation (GOA) database using the quickGO tool [36]. Out of the 21,492 unique gene probes on our array, which in some cases were annotated to identical protein identifiers (ID), 8,269 IDs possessed corresponding GO annotations (all of which were IDs annotated to SwissProt database).
Since gene-by-gene analysis did not yield large sets of significantly different expressed genes between normal and deformed fish (see results) all genes were included in the GO analyses by ranking them according to their P-values. Furthermore, in order to focus on the identification of GO categories primarily associated with gonad deformations, we only used co-regulated gene probes (i.e., same direction of regulation of gene expression) of the two within whitefish form comparisons of normal vs. deformed fish. The pooled analysis, i.e., the comparison between all normal fish vs. all deformed fish was discarded (see discussion). For GO analysis, in order to reduce redundancy due to groups of gene probes having the same protein IDs, we included only the gene probe of each such group into the GO analyses showing lowest combined P-values (as calculated by the procedure of Fisher’s Combination [37]).
We used the Wilcoxon rank test procedure implemented in the program FUNC [38] to identify enrichment of particular gene annotations. FUNC uses permutations in order to model the distribution under the null hypothesis that gene associated variables (i.e., P-values) are independent of gene annotation for significance testing. We selected the three root ontologies (i.e., GO taxonomies: “biological process”, “molecular function”, “cellular component”) on which the tests were performed and restricted tested categories to those containing a minimum of 20 genes (default value). Results of GO analysis were further narrowed down with a “refinement analysis”, also implemented in FUNC, using only GO categories with arbitrary threshold of P < 0.01. This procedure excluded GO categories that were significant solely because they contained significant descendant categories. Consequently, the list of categories was limited to the most specific significant GO terms. Only categories with number of genes < 200 were considered from the resulting list of significant GO groups due to otherwise confusing and unspecific large gene lists.

2.8. Hierarchical Clustering

Finally, in order to visualize expression patterns of particular genes involved in gonad deformations, we performed a two-way clustering analysis [39], simultaneously classifying genes and groups based on mean gene expression differences. In detail, we included all genes that were in common among certain GO categories between the two whitefish forms for each tissue respectively (see results). For each gene probe, group means were divided by the global mean across all four groups in order to account for heterogeneity in magnitudes of expression between the different genes. As clustering algorithm, we applied unweighted pair group method with the arithmetic mean (UPGMA) [40] and used centered Pearson correlation coefficients as distance measure. Hierarchical clustering analysis was carried out using the MeV v.4.3.01 software [41].

2.9. Real-Time qPCR Analysis

A new real-time PCR assay was established for the Ictacalcin gene. Prior to the design of the assay, sequence analyses were preformed for this gene to obtain detailed sequence information on the regions to be amplified in the real-time PCR reaction. The regions of interest were amplified on a GeneAmp® PCR System 9700 (Applied Biosystems) in a reaction volume of 25 μL, using cDNA from at least one whitefish individual analysed as a template and the following primers sequences: 5′-CACAGGTCCAGCAGGGTATG-3′ and 5′-AAGAACTCATTGCACATCATGG-3′. PCR conditions using Taq polymerase (Qiagen) were the following: Initial activation and denaturation at 95 °C for 10 min, followed by five cycles with denaturation at 95 °C for 1 min, annealing at 55 °C for 30 s, elongation at 72 °C for 1 min, and 25 cycles with denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s and elongation at 72 °C for 1 min with a subsequent final step of 10 min at 72 °C. Terminator Ready Reaction Mix “Big Dye” version 3.1 (Applied Biosystems) and both primers for each gene to sequence both strands in a 10 μL reaction volume. The following thermal profile was used for cycle sequencing: Initial denaturation for 50 s at 90 °C, followed by 25 cycles with denaturation at 90 °C for 10 s, annealing at 50 °C for 10 s and elongation at 60 °C for 4 min. Subsequently, products were cleaned using the DyeEx 96 spin kit (Qiagen) and separated on an ABI Prism 3130XL Genetic Analyser (Applied Biosystems).
For RT-PCR, a custom TaqMan® gene expression assay (Applied Biosystems) was used. The primers for Icatcalcin have been designed as 5′-CCTGAGCAAGGGAGAACTCA-3′ and 5′-CCTTTGCCTGGTCAGTGTTT-3′ and as probe 5′-TCTGCTCAACGCTGAGCTTGGAGAGA-3′. Transcript levels were normalized against 18S, for which primers and TaqMan probes are commercially available (Applied Biosystems). RT-PCR was performed with the 7500 fast real-time PCR system (Applied Biosystems) using 96 wells plates (MicroAmp, fast optical 96-well reaction plate and–MicroAmp optical adhesive film, Applied Biosystems) according to the Custom TaqMan® gene expression assays protocol in duplicate. For the estimation of 18S levels, cDNA was diluted 1:1,000 in RNAse-free water before addition to the 25-μL reaction. As internal standard for every PCR plate, we used a four replicates of of cDNA from a reference individual to determine interassay variability. For each assay, the amplification efficiency was determined based the slopes of standard curves. Relative quantification of the gene expression data was done using the method proposed by Paffl [42].

3. Results and Discussion

3.1. Gene-by-Gene Analysis

In liver, 7,013 and in head kidney 17,834 gene probes showed signal intensities above background, all of which were considered for statistical analyses. In general, we observed only weak differences in gene expression between normal and deformed fish. The subtle differences were more pronounced in head kidney compared to liver. However, no single gene showed significantly different gene expression between normal and deformed fish after P-value adjustments for multiple hypotheses testing with FDR ≤ 5% in both tissues and whitefish forms respectively. In contrast, when comparing transcription profiles between the two whitefish forms 2,764 genes in the liver (39.4%) and 12’820 genes in the head kidney (71.9%) showed significant expression differences (FDR ≤ 5%; Table 1).
At the nominal P-value of 0.01 in the liver 78 gene probes with FDR of 88% (i.e., nine gene probes assigning as true positives) were found, but FDR of 100% was observed for the Brienzlig comparison (Table 1a). In the head kidney, in the Albock, 264 gene probes and in the Brienzlig 528 gene probes were observed at the 1% nominal level with FDRs of 68% and 31% respectively (i.e., 84 gene probes assigning as true positives for the Albock and even 364 for the Brienzlig; Table 1) indicating weak differences in gene expression between normal and deformed fish at the level of individual genes.
Investigating congruencies among the two whitefish forms, in the liver, no gene was found at the nominal P-value of 0.01 in the three comparisons of normal vs. deformed fish and also not in the intersection between the two within-form comparisons (Figure 1a). In contrast, in head kidney, applying the same criteria, one gene, encoding for Icatcalcin, was found in all three comparisons. Another five genes were observed in the intersection of the two within-form comparisons (Figure 1b). Out of these six genes only Ictacalcin was co-regulated in both forms with up-regulation in deformed fish as compared to normal fish. Ictacalcin was also expressed at much higher levels in the head kidney as compared to the liver (Figure 2). In support of this finding two independent gene probes of Ictacalcin accounted for the top and third position by means of P-value rankings in the pooled comparison of all normal vs. all deformed fish (Table S3). No significant difference in Ictacalcin expression between deformed and normal fish was observed when reanalysing the head kidney samples using RT qPCR (Figure 3), although mean and median gene expression levels of deformed fish were higher in deformed fish in both whitefish forms.
When the two within-form comparisons of normal vs. deformed fish were contrasted against the “form comparison”, we observed a general increase in the numbers of genes in the intersections of the Venn diagrams (as indicated by arrows; Figure 1c and 1d). In other words, many of the genes that were present in the comparisons of normal vs. deformed fish, were also present in the comparison between the two whitefish forms. Therefore, when considering the pooled comparison, a substantial number of genes was affected by the strong differences in gene expression between the two whitefish forms. This finding can be described as a “form effect”. The “form effect” is well illustrated by hierarchical clustering analysis (see below). A large number of genes were present, which showed similar expression differences between normal and deformed fish, but which were expressed on different levels with regard to the two whitefish forms (Figure 4a and 4b).
In order to test, whether the strong differences between the two whitefish forms are caused by a technical bias, we carried out the analysis by normalizing the raw data for each form separately. We found that population normalized intensity measures did not change noticeably as compared to a normalization procedure carried out across both whitefish forms (R2 > 0.9999 in both groups). Therefore, based on these results and the above described “form effect” we included only the two within-form comparisons of normal vs. deformed fish in GO analysis.
In summary, while gene expression differences between normal and deformed fish were generally weak, strong differences were found, when the two whitefish forms were compared. These strong differences may be explained by ecological and life-history differences between the Brienzlig and Albock whitefish forms (i.e., spatial and temporal differences in spawning as well as feeding ecology), which require complex physiological adaptations. Additionally, the two whitefish forms were genetically distinct (see Experimental Section) and genetic divergence is well-known to affect global transcription profiles (e.g., [23,43]).
An important recent finding provides a potential explanation for the weak differences in gene expression between normal and deformed fish. In a series of long-term rearing experiments in the laboratory [11] showed, that whitefish fed with zooplankton originating from Lake Thun developed gonad deformations at comparable frequencies to wild fish regardless of genetic background (even whitefish from different lakes additionally to Lake Thun were used) and origin of water (even spring water besides water from Lake Thun was used). These results offer strong evidence that the feeding on zooplankton is directly associated with observed gonad deformations. The whitefish sampled for microarray analyses were caught on their spawning sites, i.e., they have not been feeding for some time, maybe even weeks, which was clearly visible by their entirely empty intestinal tracts during the sampling period (D. Bittner; pers. observation). Therefore, we might have measured gene expression levels at a point in time, where an initial possibly stronger response had already faded away.

3.2. GO Analysis

Among the gene probes showing signal intensities above background, 3,449 in liver and 7,226 in head kidney had GO annotations and were consistently either up- or down-regulated in deformed compared to normal fish. In an enrichment analysis of GO categories performed on these gene sets revealed an elevated number of co-regulated genes in specific GO categories, i.e., even though single genes were not significantly differently expressed, a combination or group of genes contained within GO categories showed significant difference in gonadal phenotype correlated gene expression. Significant GO categories were identified in all four analyses (i.e., normal vs. deformed fish in Albock and Brienzlig for both tissues respectively), after P-value adjustment for multiple hypotheses testing (Table S4). After refinement analysis (see Material and Methods section), in both tissues significant GO terms related to the immune system were most prominent (Table 2). These were present in all four analyses (i.e., in both whitefish forms and tissues respectively) with particular representation in the comparison of normal vs. deformed Albock in the liver where five out of seven significant GO categories were involved in immunity (Table 2a). Strong overlap of corresponding genes among the immune system categories was found. Additionally, we found strong congruencies of genes between GO categories “proteolysis” of the Albock and “extracellular region” of the Brienzlig in the liver (Table 2a, b) and between GO terms related to GTPase activation and Ras protein signal transduction in both forms in the head kidney (Table 2c, d).
Due to the structure of the gene ontology, a particular gene can be annotated in multiple GO categories simultaneously. In the liver, 252 unique genes were present among the seven significant GO terms in the Albock and 326 unique genes among the four significant GO terms in the Brienzlig out of which 102 genes were found in common between the two whitefish forms (see Figure 4a for gene identities). In the head kidney, 318 unique genes were observed among the twelve significant GO terms in the Albock and 172 unique genes among the six significant GO groups in the Brienzlig out of which again a relatively large number of genes, (71 in total, see Figure 4b for gene identities) were found in common between the two forms. These findings indicated substantial overlap in GO analysis between the two whitefish forms and both tissues, respectively.
Hierarchical clustering analysis of the four groups across the 102 common genes in liver, and the 71 shared genes in head kidney showed a separation of whitefish form rather than gonadal phenotype, i.e., Albock normal clustered with Albock deformed fish and vice versa. This reflects the strong global gene expression differences between the two whitefish forms (Figure 4). In liver, most of the particular 102 shared genes were involved in immune system processes and a smaller number in the biological process “proteolysis” and the cellular compartment “extracellular region” (Figure 4a). In head kidney, again many of the 71 shared genes were involved in immunity, but also a large number of genes were related to the molecular function “GTPase activation” and the biological process “Ras protein signal transduction (Figure 4b).
In the liver 70 genes were shared among the immune system related GO categories “immune response”, “humoral immune response”, “innate immune response”, “activation of immune response” and “immune effector process” in the Albock (Table 2a), and “immune system process” in the Brienzlig (Table 2b). Among these 70 common genes, 29 were down- and 41 genes were up-regulated in deformed fish. A total of 19 genes were found in common among all six GO categories comprising both whitefish forms. These genes were primarily involved in activation of the complement system (Figure 5a). A continuous up-regulation of these genes was observed with the exception of two genes.
Another major group of immune system related genes encoded for histocompatibility antigens. Relatively large mean expression differences in both forms between normal and deformed fish were further observed for the “proteasome subunit beta type 8 precursor” gene, which plays an important role in antigen processing. All these genes were down-regulated in deformed fish of both forms. In contrast, a few immunoglobulin genes of the heavy or kappa chains were up-regulated (Figure 5a). Furthermore, in the liver, 44 genes were in common between GO terms “proteolysis” in the Albock and “extracellular region” in the Brienzlig. It is interesting that 18 of these 44 genes were in common with the 19 specific genes involved in complement activation, which is indicated in Figure 5b. Additionally, a small group of genes encoding for coagulation factors and collagenase precursors was down-regulated in deformed fish; however, 33 out of the 44 shared genes were up-regulated (Figure 5b).
Head kidney, also displayed a remarkable overlap between the two whitefish forms of genes in GO categories related to the immune system (Table 2c, d). The GO terms concerned were “regulation of immune system process” in the Albock and “leukocyte differentiation” and “T cell activation” in the Brienzlig (Table 2c, d) with 28 shared genes between the two whitefish forms (eleven genes were present in all three GO terms; Figure 5c). In contrast to genes involved in the up-regulation of the complement system in the liver, no clear pattern in the direction of gene expression regulation was observed in head kidney. However, as in the liver, two genes encoding for histocompatibility antigens were down-regulated (Figure 5c). Furthermore, in the head kidney, a second overlap of genes among significant GO categories between the Albock and Brienzlig concerned categories related to “small GTPase regulator activity” and “regulation of Ras protein signal transduction” (Table 2c, d). In more detail, a total of 22 genes were in common among the GO categories “small GTPase regulator activity” and “regulation of Ras protein signal transduction” in the Albock and “GTPase activator activity” in the Brienzlig (Figure 4d). Nine of these 22 shared genes were up- and 13 were down-regulated in deformed fish. Many of these genes encoded for GTPase activating proteins and a tendency in down-regulation of GTPase activating proteins was observed (Figure 5d).

3.3. Deviation in Regulation of Immune System Processes

GO analysis revealed subtle but congruent changes in gene expression patterns between normal and deformed fish of both, genetically and ecologically differentiated whitefish forms. This suggests that the observed differences in gene expression are of biological relevance and actually associated with the gonadal deformations in these fish. Significant GO categories associated primarily with the immune system and the results were consistent across both whitefish forms and tissues, respectively.
A dysregulation of immune function is commonly observed in situations where fish are stressed (e.g., [4446]), and it thus follows that deformed fish might be more stressed than fish with normally developed gonads. However, since gene expression patterns between normal and deformed fish were rather weak and since we do not know the threshold when to consider differences as dysregulation, here, we use the neutral term “deviation in regulation”. Deviation in regulation of gene expression of immune system processes can be classified into three different major immune responses: (i) primary response of the immune system (i.e., first response when organisms are exposed to exogenous factors such as pathogens or toxicants, that involve inflammatory and acute phase proteins as well as chemokines); (ii) secondary response of the immune system (i.e., a response that follows at a later stage, which is an adapted and specifically directed immune response) and iii) self-directed secondary response of the immune system, i.e., autoimmune response.
Microarrays have been used extensively to describe responses to pathogenic agents, i.e., viral, bacterial and parasitic disease processes (as reviewed by [47]) and exposure to environmental contaminants (e.g., [19,4850]), where deviations in regulation of gene expression patterns of immune system processes were also observed.
The most prominent finding of our analyses concerned the up-regulation of the complement system in deformed fish in the liver. The complement system is a very ancient mechanism of defense that is conserved across many species. Complement is primarily involved in host-pathogen interactions, i.e., defense against bacteria, viruses, fungi and parasites. As part of the innate immune response, complement plays an important role in enhancing the uptake and processing of antigens by antigen-presenting cells. The complement system is composed of a range of soluble plasma proteins that play key roles in innate and adaptive immunity [51]. It is interesting that teleosts, unlike mammals, contain multiple isoforms of complement components, which implies novel roles in fish innate immunity [52]. Furthermore C3/C4 fragments bound to antigen or immune complexes enhance uptake and processing of antigen by antigen-presenting cells, which leads to more effective primary and secondary antibody responses to the antigen [5355]).
Even though involvement of the complement system points towards potential roles of pathogens, thorough histological investigations have failed to demonstrate any signs of inflammation and/or infection in the gonads of deformed fish [1]. Furthermore, potential involvement of pathogens is unlikely, since the two whitefish forms were sampled independently of each other and at different locations. Such a case would be more than coincidental. Additionally, various exposure studies have always reported much stronger differences in gene expression compared to our data, where pathogens were involved (e.g., 1.8–30 fold differences in transcription of complement genes were observed by [49]. However, contrary to all exposure studies that compared exposed fish to non-exposed fish, in the present study, normal and deformed fish were sampled in the same environment, and thus had likely been equally exposed.

3.4. Autoimmune Disease?

Alternatively, the gene expression patterns found in immune system processes in both tissues and in particular the up-regulation of the complement system in the liver could be explained by a chronic autoimmune disease in the testis. In humans, microarrays have been widely used for the study of autoimmune diseases and specific gene expression patterns have been detected [56]. However, data is lacking as to specific hallmarks of gene expression patterns that characterize autoimmune disease, even more so in fish. But autoimmune disease often involves multiple organs [56] and specific GO categories involved in the immune system as well as additional indications in our data fit the idea of an autoimmune activation. In particular following significant GO categories were concerned: “T cell activation”, “GTPase activator activity”, “small GTPase regulator activity”, “regulation of Ras protein signal transduction”, “cell adhesion”, “cell migration”, “extracellular region” and “proteinaceous extracellular matrix” (Table 2).
In comparison with the current literature we found many interesting links of particularities of our data in the context with autoimmunity. Effective immune responses require the appropriate activation and differentiation of peripheral T cells. Defects in appropriate regulation of T cell activation underlie the pathogenesis of many autoimmune disorders in humans [57]. Even though the molecular machinery to control these processes and to prevent autoimmunity is not fully understood, activation of Rho GTPases, is recognized as a key event in the coordination of immune responses and, particularly, in the activation of T cells [57]. These recent findings are very striking in regard to our own results, since “GTPase activation”, which involved Rho GTPase activating proteins, and “T cell activation” were among the identified significant GO categories. However, Rho GTPase activating proteins are molecular switches that also control many other biological processes including complex signaling pathways as well as cytoskeletal reorganization [57].
Cell-cell and cell-matrix adhesion molecules as well as extracellular matrix components are target structures of antibody-mediated autoimmunity. Pathogenic auto-antibodies against these molecules are causally related to disturbances of cell and tissue adhesion [58]. Again our data are consistent with these findings. In liver, we identified related GO terms “extracellular region”, “proteolysis”, “cell adhesion” and “integral to plasma membrane”, whereas in head kidney GO terms “cell migration” and “cell proliferation” were present among the significant categories. Interestingly, proteins of the extracellular matrix, which is part of the extracellular region, have also been shown to play an important role in testis development in male mice [59] and specifically in rainbow trout [60]. In trout, microarray data indicated that the testicular extracellular matrix participates in processes such as growth, adhesion, differentiation, cell migration and patterning, i.e., many processes, which were also evident in our data. In this regard, it is also interesting to note that in mammals the complement system has been shown to participate in gamete development and provides potential bridges through which gametes can interact [61,62]. Developmental regulation of various complement components in reproductive processes was also demonstrated in rainbow trout [60,63,64]. This is particularly interesting when considering that gonad deformations develop right with onset of gonad differentiation [11]. Therefore the gonad deformations could result from dysfunctional cell-cell interactions, i.e., failure in normal cell-cell communication, during development, induced by autoimmunity. Additionally, collagen is well known to play an important role in tissue development. Our data provide a small detail further supporting the autoimmunity hypothesis, since a group of collagen genes were present in GO categories “proteolysis” (Albock) and also in “extracellular region” (Brienzlig).
By contrast, autoimmunity is commonly associated with increased expression of major histocompatibility complex (MHC) antigens and immunoglobulins [65]. This is partly incoherent with our findings, where genes encoding for histocompatibility antigens were consistently down-regulated. However, we analyzed “only” liver and head kidney tissue and not the gonads themselves, the targeted tissue of the possible autoimmune disease, where we would expect up-regulation of histocompatibility antigens according to Herskowitz et al. [65].
With regard to the up-regulation of genes of the complement system in the liver, it is important to mention that it is well known that complement pathway is also one of the major means by which the body recognizes tissue injury [66]. Complement activation in the course of systemic autoimmunity leads to tissue damage in a number of ways including direct lysis of cell, modification of cell function and by contributing to the formation of immune complexes [67] and thus is thought to contribute significantly to end organ damage [68]. In mammals, it has been shown that auto-antibodies in combination with intracellular antigens can trigger cell injury at cell surface by activating the complement system [67]. On the other hand, there are many mechanisms by which the complement system is activated and it is also involved in many other processes (see above).
Although rare, especially in fish, autoimmunity has previously been demonstrated, e.g., in vaccinated farmed Atlantic salmon (Salmo salar), where the presence of auto-antibodies was confirmed [69]. Adhesion of internal organs and tissue damage was found in all vaccinated fish in this study. Interestingly, lesions of testis characteristic of autoimmune attack were also shown in rainbow trout passively immunized with anti-sperm antibodies, which involved the presence of complement components [70]. Mature sperm (and their antigens) develop through the process of spermatogenesis, which arise only in adult organisms, i.e., long after tolerance of non-specific immune system has been established. In the testis, Sertoli cells establish intercellular junctions that are essential for spermatogenesis and the blood-testis barrier is formed by tight junctions between these Sertoli cells. The presence of a blood-testis barrier is thought to provide an environment in which the appearance of specific antigens (possibly also auto-antigens) is facilitated [71]. As mentioned above, it has been demonstrated that autoimmune lesions in the testis could be induced experimentally, which, besides rainbow trout, was also the case in Atlantic salmon by injecting autologous or allogenic testis material into adult fishes [72]. These findings suggest that the blood-testis barrier may have broken down during autoimmune response to the testis, which was confirmed in the study by Lou and Takahashi [71] on Nile tilapia. The Authors provided morphological evidence for the presence of a blood-testis barrier and its break down during immune response. This indicates that the function of the blood-testis barrier may indeed be to prevent autoimmune reaction. Interestingly, it has been shown that impaired Sertoli cell junctional protein expression were the consequence of exposure to different classes of reproductive toxicants [73].

3.5. Indications of Ictacalcin

Studies are accumulating that autoimmune diseases are also stimulated by chronic exposure to various chemicals, i.e., xenobiotic substances, in humans as well as in animals [74,75]. Our data, however, provide only a weak indication pointing to particular chemical substances. We found one gene, Ictacalcin, at the nominal 1% level that showed up-regulation in deformed fish in head kidney of both whitefish forms. Reanalysis of the samples with RT qPCR revealed qualitatively the same but a statistically not significant difference between deformed and normal fish in both forms. Ictacalcin is a Ca2+ binding protein abundant in epithelial cells of the olfactory barbells, skin, and gills of fish suggesting an important role in calcium homeostasis [76,77]. Interestingly, expression of Ictacalcin was significantly increased in juvenile rainbow trout exposed to TCDD dioxin in the head kidney, but not in liver [78], which is in line with our results. The mechanisms of the TCDD-induced increase of Ictacalcin in rainbow trout, however, were not clear. Calcium is the most common intracellular messenger in the signal transduction pathway and is required for cell growth and survival and thus plays an important role in various physiological processes. Therefore, the isolated induction of Ictacalcin alone is not conclusive that whitefish in Lake Thun are exposed in any way to dioxin-like substances. Furthermore, the GO analyses did not reveal any GO categories that are indicative for the involvement of dioxins. In this context, it also is worth noting that the array contains several probes for CYP1A genes, which are known to show the most sensitive response to exposure to dioxins like TCDD and which did not show any sign of up-regulation in deformed fish.

4. Conclusions

Our study showed that the application of gene expression profiling to study natural populations represents a promising approach to address complex ecotoxicological problems. By using such a transcriptomic approach, we identified physiological processes and one candidate gene potentially involved in gonad deformations. A deviation in regulation of gene expression of immune system processes, the most prominent finding of our data, is commonly observed in situations where fish are stressed, thus indicating, that deformed fish are under increased influence of stress. In comparison with current literature, the gene expression patterns found here best fit to a signature of an autoimmune disease in the testes. Based on the recent observations that gonad deformations are induced through feeding of zooplankton from Lake Thun, we hypothesize that a xenobiotic accumulated in whitefish via the plankton triggering autoimmunity as the likely cause of gonad deformations in Lake Thun. As a follow-up study, in order to test the autoimmunity hypothesis, we propose screen normal and deformed whitefish for auto-antibodies against gonad tissue.

Supplementary Materials

ijerph-08-02706-s001.doc

Acknowledgements

We thank the Fisheries Inspectorate of Bern and the professional fishermen of Lake Thun for their co-operation and assistance in the field work. Special thanks go to L. Olohan, L. Rainbow U. Amstutz technical help and helpful discussion. We also thank S. Pashk, A. Casanova and G. Andrey for their help with the qPCR analyses. This work was funded by a joint grant (MICROCORE) from the Swiss Federal Office for the Environment (FOEN) and Fisheries Inspectorate of Bern to C.R. Largiadèr.

References and Notes

  1. Bernet, D; Wahli, T; Kueng, C; Segner, H. Frequent and unexplained gonadal abnormalities in whitefish (central alpine Coregonus sp.) from an alpine oligotrophic lake in Switzerland. Dis. Aquat. Org 2004, 61, 137–148. [Google Scholar]
  2. Bittner, D; Bernet, D; Wahli, T; Segner, H; Kung, C; Largiader, CR. How normal is abnormal? Discrimination between deformations and natural variation in gonad morphology of European whitefish Coregonus lavaretus. J. Fish Biol 2009, 74, 1594–1614. [Google Scholar]
  3. Allen, Y; Scott, AP; Matthiessen, P; Haworth, S; Thain, JE; Feist, S. Survey of estrogenic activity in United Kingdom estuarine and coastal waters and its effects on gonadal development of the flounder Platichthys flesus. Environ. Toxicol. Chem 1999, 18, 1791–1800. [Google Scholar]
  4. Van der Ven, LTM; Wester, PW; Vos, JG. Histopathology as a tool for the evaluation of endocrine disruption in zebrafish (Danio rerio). Environ. Toxicol. Chem 2003, 22, 908–913. [Google Scholar]
  5. Vethaak, AD; Lahr, J; Schrap, SM; Belfroid, AC; Rijs, GBJ; Gerritsen, A; de Boer, J; Bulder, AS; Grinwis, GCM; Kuiper, RV; et al. An integrated assessment of estrogenic contamination and biological effects in the aquatic environment of The Netherlands. Chemosphere 2005, 59, 511–524. [Google Scholar]
  6. Schafers, C; Teigeler, M; Wenzel, A; Maack, G; Fenske, M; Segner, H. Concentration- and time-dependent effects of the synthetic estrogen, 17 alpha-ethinylestradiol, on reproductive capabilities of the zebrafish, Danio rerio. J. Toxicol. Environ. Health, Part A 2007, 70, 768–779. [Google Scholar]
  7. Bernet, D; Liedtke, A; Bittner, D; Eggen, RIL; Kipfer, S; Küng, C; Largiader, CR; Suter, MJ-F; Wahli, T; Segner, H. Gonadal malformations in whitefish from Lake Thun: defining the case and evaluating the role of EDCs. Chimia 2008, 62, 383–388. [Google Scholar]
  8. Bogdal, C; Naef, M; Schmid, P; Kohler, M; Zennegg, M; Bernet, D; Scheringer, M; Hungerbuhler, K. Unexplained gonad alterations in whitefish (Coregonus spp.) from Lake Thun, Switzerland: Levels of persistent organic pollutants in different morphs. Chemosphere 2009, 74, 434–440. [Google Scholar]
  9. Kipfer, S; Segner, H; Wenger, M; Wahli, T; Bernet, D. Long-term estrogen exposure of whitefish Coregonus lavaretus induces intersex but not Lake Thun-typical gonad malformations. Dis. Aquat. Org 2009, 84, 43–56. [Google Scholar]
  10. Bittner, D. Gonad Deformations in Whitefish (Coregonus spp.) from Lake Thun, Switzerland-A Population Genetic and Transcriptomic Approach. PhD Thesis, University of Bern, Bern, Switzerland, 2009. [Google Scholar]
  11. Bernet, D; Wahli, T; Küng, C; Zieri, HR; Segner, H. Zooplankton as the key factor for the development of macroscopical gonad deformation in whitefish (Coregonus lavaretus) from Lake Thun, Switzerland. Aquatic Biology 2011. under revision. [Google Scholar]
  12. Cossins, AR; Crawford, DL. Opinion-Fish as models for environmental genomics. Nat. Rev. Genet 2005, 6, 324–333. [Google Scholar]
  13. Denslow, ND; Garcia-Reyero, N; Barber, DS. Fish ‘n’ chips: the use of microarrays for aquatic toxicology. Mol. BioSyst 2007, 3, 172–177. [Google Scholar]
  14. Miller, KM; Maclean, N. Teleost microarrays: development in a broad phylogenetic range reflecting diverse applications. J. Fish Biol 2008, 72, 2039–2050. [Google Scholar]
  15. Jin, W; Riley, RM; Wolfinger, RD; White, KP; Passador-Gurgel, G; Gibson, G. The contributions of sex, genotype and age to transcriptional variance in Drosophila melanogaster. Nat. Genet 2001, 29, 389–395. [Google Scholar]
  16. Oleksiak, MF; Churchill, GA; Crawford, DL. Variation in gene expression within and among natural populations. Nat. Genet 2002, 32, 261–266. [Google Scholar]
  17. Ranz, JM; Castillo-Davis, CI; Meiklejohn, CD; Hartl, DL. Sex-dependent gene expression and evolution of the Drosophila transcriptome. Science 2003, 300, 1742–1745. [Google Scholar]
  18. Townsend, JP; Cavalieri, D; Hartl, DL. Population genetic variation in genome-wide gene expression. Mol. Biol. Evol 2003, 20, 955–963. [Google Scholar]
  19. Williams, TD; Gensberg, K; Minchin, SD; Chipman, JK. A DNA expression array to detect toxic stress response in European flounder (Platichthys flesus). Aquat. Toxicol 2003, 65, 141–157. [Google Scholar]
  20. Oleksiak, MF; Roach, JL; Crawford, DL. Natural variation in cardiac metabolism and gene expression in Fundulus heteroclitus. Nat. Genet 2005, 37, 67–72. [Google Scholar]
  21. Fisher, MA; Oleksiak, MF. Convergence and divergence in gene expression among natural populations exposed to pollution. BMC Genomics 2007, 8, 108. [Google Scholar]
  22. Olohan, LA; Li, W; Wulff, T; Jarmer, H; Gracey, AY; Cossins, AR. Detection of anoxia-responsive genes in cultured cells of the rainbow trout Oncorhynchus mykiss (Walbaum), using an optimized, genome-wide oligoarray. J. Fish Biol 2008, 72, 2170–2186. [Google Scholar]
  23. Giger, T; Excoffier, L; Amstutz, U; Day, PJR; Champigneulle, A; Hansen, MM; Kelso, J; Largiader, CR. Population transcriptomics of life-history variation in the genus Salmo. Mol. Ecol 2008, 17, 3095–3108. [Google Scholar]
  24. Bernatchez, L. Ecological theory of adaptive radiation An empirical assessment from Coregonine fishes (Salmoniformes). In Evolution Illuminated, Salmon and Their Relatives; Hendry, AP, Stearns, SC, Eds.; Oxford University Press: New York NY, USA, 2004; pp. 175–207. [Google Scholar]
  25. Patterns of morphological changes and hybridization between sympatric whitefish morphs (Coregonus spp.) in a Swiss lake: a role for eutrophication? Mol. Ecol 2010, 19, 2152–2167.
  26. Berg, A; Grimaldi, E. A critical interpretation of the scale structures used for determination of annuli in fish growth studies. Mem. Ist. Ital. Idrobiol 1967, 21, 225–239. [Google Scholar]
  27. Rannala, B; Mountain, JL. Detecting immigration by using multilocus genotypes. Proc. Natl. Acad. Sci. USA 1997, 94, 9197–9201. [Google Scholar]
  28. Paetkau, D; Slade, R; Burden, M; Estoup, A. Genetic assignment methods for the direct, real-time estimation of migration rate: A simulation-based exploration of accuracy and power. Mol. Ecol 2004, 13, 55–65. [Google Scholar]
  29. Piry, S; Alapetite, A; Cornuet, JM; Paetkau, D; Baudouin, L; Estoup, A. GENECLASS2: A software for genetic assignment and first-generation migrant detection. J. Hered 2004, 95, 536–539. [Google Scholar]
  30. Renn, SCP; Aubin-Horth, N; Hofmann, HA. Biologically meaningful expression profiling across species using heterologous hybridization to a cDNA microarray. BMC Genomics 2004, 5, 42. [Google Scholar]
  31. Workman, C; Jensen, LJ; Jarmer, H; Berka, R; Gautier, L; Nielsen, HB; Saxild, HH; Nielsen, C; Brunak, S; Knudsen, S. A new non-linear normalization method for reducing variability in DNA microarray experiments. Genome Biol 2002, 16, 0041–0048. [Google Scholar]
  32. Sahai, H; Ageel, MI. The Analysis of Variance: Fixed, Random and Mixed Models; Birkhäuser: Boston, MA, USA, 2000. [Google Scholar]
  33. Excoffier, L; Smouse, PE; Quattro, JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes-application to human mitochondrial-DNA restriction data. Genetics 1992, 131, 479–491. [Google Scholar]
  34. Benjamini, Y; Hochberg, Y. Controlling the false discovery rate-a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B 1995, 57, 289–300. [Google Scholar]
  35. Harris, MA; Clark, JI; Ireland, A; Lomax, J; Ashburner, M; Collins, R; Eilbeck, K; Lewis, S; Mungall, C; Richter, J; et al. The Gene Ontology (GO) project in 2006. Nucleic Acids Res 2006, 34, D322–D326. [Google Scholar]
  36. Camon, E; Barrell, D; Lee, V; Dimmer, E; Apweiler, R. The gene ontology annotation (GOA) database-An integrated resource of GO annotations to the UniProt knowledgebase. In Silico Biol 2004, 4, 5–6. [Google Scholar]
  37. Sokal, RR; Rohlf, FJ. Biometry: The Principles and Practices of Statistics in Biological Research, 3rd ed; W.H. Freeman: San Francisco, CA, USA, 1994. [Google Scholar]
  38. Prufer, K; Muetzel, B; Do, HH; Weiss, G; Khaitovich, P; Rahm, E; Paabo, S; Lachmann, M; Enard, W. FUNC: A package for detecting significant associations between gene sets and ontological annotations. BMC Bioinf 2007, 8, 41. [Google Scholar]
  39. Hartigan, J. Clustering Algorithms; Wiley: New York, NY, USA, 1975. [Google Scholar]
  40. Sneath, PHA; Sokal, RR. Numerical Taxonomy; Freeman: San Francisco, CA, USA, 1973. [Google Scholar]
  41. Saeed, AI; Sharov, V; White, J; Li, J; Liang, W; Bhagabati, N; Braisted, J; Klapa, M; Currier, T; Thiagarajan, M; et al. TM4: A free, open-source system for microarray data management and analysis. Biotechniques 2003, 34, 374–390. [Google Scholar]
  42. Paffl, W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res 2001, 29, e9. [Google Scholar]
  43. St-Cyr, J; Derome, N; Bernatchez, L. The transcriptomics of life-history trade-offs in whitefish species pairs (Coregonus sp.). Mol. Ecol 2008, 17, 1850–1870. [Google Scholar]
  44. Cuesta, A; Esteban, MA; Meseguer, J. Effects of different stressor agents on gilthead seabream natural cytotoxic activity. Fish Shellfish Immunol 2003, 15, 433–441. [Google Scholar]
  45. Yin, Z; Lam, TJ; Sin, YM. The Effects of crowding stress on the nonspecific immune-response in fancy carp (Cyprinus-Carpio L). Fish Shellfish Immunol 1995, 5, 519–529. [Google Scholar]
  46. Ortuno, J; Esteban, MA; Meseguer, J. Effects of short-term crowding stress on the gilthead seabream (Sparus aurata L.) innate immune response. Fish Shellfish Immunol 2001, 11, 187–197. [Google Scholar]
  47. Douglas, SE. Microarray studies of gene expression in fish. OMICS: J. Integr. Biol 2006, 10, 474–489. [Google Scholar]
  48. Hook, SE; Skillman, AD; Gopalan, B; Small, JA; Schultz, IR. Gene expression profiles in rainbow trout, Onchorynchus mykiss, exposed to a simple chemical mixture. Toxicol. Sci 2008, 102, 42–60. [Google Scholar]
  49. Roberge, C; Paez, 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, 1374–1383. [Google Scholar]
  50. Williams, TD; Diab, A; Ortega, F; Sabine, VS; Godfrey, RE; Falciani, F; Chipman, JK; George, SG. Transcriptomic responses of European flounder (Platichthys flesus) to model toxicants. Aquat. Toxicol 2008, 90, 83–91. [Google Scholar]
  51. Gasque, P. Complement: a unique innate immune sensor for danger signals. Mol. Immunol 2004, 41, 1089–1098. [Google Scholar]
  52. Boshra, H; Li, J; Sunyer, JO. Recent advances on the complement system of teleost fish. Fish Shellfish Immunol 2006, 20, 239–262. [Google Scholar]
  53. Jacquiersarlin, MR; Gabert, FM; Villiers, MB; Colomb, MG. Modulation of antigen-processing and presentation by covalently-linked complement C3b fragment. Immunology 1995, 84, 164–170. [Google Scholar]
  54. Villiers, MB; Villiers, CL; JacquierSarlin, MR; Gabert, FM; Journet, AM; Colomb, MG. Covalent binding of C3b to tetanus toxin: Influence on uptake/internalization of antigen by antigen-specific and non-specific B cells. Immunology 1996, 89, 348–355. [Google Scholar]
  55. Boackle, SA; Morris, MA; Holers, VM; Karp, DR. Complement opsonization is required for presentation of immune complexes by resting peripheral blood B cells. J. Immunol 1998, 161, 6537–6543. [Google Scholar]
  56. Oertelt, S; Selmi, C; Invernizzi, P; Podda, M; Gershwin, ME. Genes and goals: An approach to microarray analysis in autoimmunity. Autoimmun. Rev 2005, 4, 414–422. [Google Scholar]
  57. Pernis, AB. Rho GTPase-mediated pathways in mature CD4* T cells. Autoimmun. Rev 2009, 8, 199–203. [Google Scholar]
  58. Moll, R; Bahn, H; Bayerl, C; Moll, I. Cellular adhesion molecules and components of the extracellular matrix as target structures of autoimmunity. Pathologe 1996, 17, 254–61. [Google Scholar]
  59. Clement, TM; Anway, MD; Uzumcu, M; Skinner, MK. Regulation of the gonadal transcriptome during sex determination and testis morphogenesis: Comparative candidate genes. Reproduction 2007, 134, 455–472. [Google Scholar]
  60. von Schalburg, KR; Cooper, GA; Yazawa, R; Davidson, WS; Koop, BF. Microarray analysis reveals differences in expression of cell surface and extracellular matrix components during development of the trout ovary and testis. Comp. Biochem. Physiol., Part D: Genomics Proteomics 2008, 3, 78–90. [Google Scholar]
  61. Mastellos, D; Lambris, JD. Complement: More than a ‘guard’ against invading pathogens? Trends Immunol 2002, 23, 485–491. [Google Scholar]
  62. Riley-Vargas, RC; Lanzendorf, S; Atkinson, JP. Targeted and restricted complement activation on acrosome-reacted spermatozoa. J. Clin. Invest 2005, 115, 1241–1249. [Google Scholar]
  63. von Schalburg, KR; Rise, ML; Brown, GD; Davidson, WS; Koop, BF. A comprehensive survey of the genes involved in maturation and development of the rainbow trout ovary. Biol. Reprod 2005, 72, 687–699. [Google Scholar]
  64. Von Schalburg, KR; McCarthy, SP; Rise, ML; Hutson, JC; Davidson, WS; Koop, BF. Expression of morphogenic genes in mature ovarian and testicular tissues: Potential stem-cell niche markers and patterning factors. Mol. Reprod. Dev 2006, 73, 142–152. [Google Scholar]
  65. Herskowitz, A; Ahmed-Ansari, A; Neumann, DA; Beschorner, WE; Rose, NR; Soule, LM; Burek, CL; Sell, KW; Baughman, KL. Induction of major histocompatibility complex antigens within the myocardium of patients with active myocarditis: A nonhistologic marker of myocarditis. J. Am. Coll. Cardiol. Found 1990, 15, 624–632. [Google Scholar]
  66. Holers, VM. Complement receptors and the shaping of the natural antibody repertoire. Springer Semin. Immunopathol 2005, 26, 405–423. [Google Scholar]
  67. Toskos, GC; Fleming, SD. Complement in Autoimmunity; Karger: Basel, Switzerland, 2004. [Google Scholar]
  68. Song, W-C. Complement in Autoimmunity; Karger: Basel, Switzerland, 2004. [Google Scholar]
  69. Koppang, EO; Bjerkas, I; Haugarvoll, E; Chan, EKL; Szabo, NJ; Ono, N; Akikusa, B; Jirillo, E; Poppe, TT; Sveier, H; et al. Vaccination-induced systemic autoimmunity in farmed Atlantic salmon. J. Immunol 2008, 181, 4807–4814. [Google Scholar]
  70. Secombes, CJ; Lewis, AE; Laird, LM; Needham, EA; Priede, IG. Role of autoantibodies in the autoimmune-response to testis in rainbow-trout (Salmo-Gairdneri). Immunology 1985, 56, 409–415. [Google Scholar]
  71. Lou, YH; Takahashi, H. The blood-testis barrier and its breakdown following immunization to testis material in the Nile tilapia, Oreochromis-Niloticus. Cell Tissue Res 1989, 258, 491–498. [Google Scholar]
  72. Laird, LM; Wilson, AR; Holliday, FGT. Lesions observed in the testis of precociously maturing male atlantic salmon, Salmo-Salar L. J. Fish Biol 1980, 17, 343–348. [Google Scholar]
  73. Fiorini, C; Tilloy-Ellul, A; Chevalier, S; Charuel, C; Pointis, G. Sertoli cell junctional proteins as early targets for different classes of reproductive toxicants. Reprod. Toxicol 2004, 18, 413–421. [Google Scholar]
  74. Banerjee, BD; Chakraborti, A; Suke, SG; Ahmed, RS; Tripathi, AK. Xenobiotic-induced immune alterations: Implications in health and disease. Indian J. Biochem. Biophys 2008, 45, 7–15. [Google Scholar]
  75. Bigazzi, PE. Autoimmunity caused by xenobiotics. Toxicology 1997, 119, 1–21. [Google Scholar]
  76. Bettini, E; Porta, AR; Dahmen, N; Wang, H; Margolis, FL. Expressed sequence tags (Est) identify genes preferentially expressed in catfish chemosensory tissues. Mol. Brain Res 1994, 23, 285–291. [Google Scholar]
  77. Porta, AR; Bettini, E; Buiakova, OI; Baker, H; Danho, W; Margolis, FL. Molecular cloning of ictacalcin: A novel calcium-binding protein from the channel catfish, Ictalurus punctatus. Mol. Brain. Res 1996, 41, 81–89. [Google Scholar]
  78. Cao, ZJ; Tanguay, RL; McKenzie, D; Peterson, RE; Aiken, JM. Identification of a putative calcium-binding protein as a dioxin-responsive gene in zebrafish and rainbow trout. Aquat. Toxicol 2003, 63, 271–282. [Google Scholar]
Figure 1. Venn diagrams showing the number of genes differently expressed between whitefish with normal (nor) and deformed (def) gonads. (a) and (b) show within-form comparisons contrasted to the pooled normal and deformed fish comparison and (c) and (d) the within-form comparisons contrasted to the comparison between all Albock and all Brienzlig individuals. Only genes with P-value < 0.01 for liver (a & c) and head kidney (b & d) are considered.
Figure 1. Venn diagrams showing the number of genes differently expressed between whitefish with normal (nor) and deformed (def) gonads. (a) and (b) show within-form comparisons contrasted to the pooled normal and deformed fish comparison and (c) and (d) the within-form comparisons contrasted to the comparison between all Albock and all Brienzlig individuals. Only genes with P-value < 0.01 for liver (a & c) and head kidney (b & d) are considered.
Ijerph 08 02706f1aIjerph 08 02706f1b
Figure 2. Box plots of log2 normalized gene expression measurements of Ictacalcin in head kidney (a) and liver tissue (b) of normal (nor) and deformed (def) Albock (AL) and Brienzlig (BR) whitefish. Shown are median (vertical line), 25% and 75% quartiles (box), and minimum and maximum value (error bars). Outliers depicted as small squares represent data points that are more than 1.5 times the interquartile range away from the box.
Figure 2. Box plots of log2 normalized gene expression measurements of Ictacalcin in head kidney (a) and liver tissue (b) of normal (nor) and deformed (def) Albock (AL) and Brienzlig (BR) whitefish. Shown are median (vertical line), 25% and 75% quartiles (box), and minimum and maximum value (error bars). Outliers depicted as small squares represent data points that are more than 1.5 times the interquartile range away from the box.
Ijerph 08 02706f2
Figure 3. Box plots of relative RT qPCR gene expression measurements of Ictacalcin in head kidney tissue of normal (nor) and deformed (def) Albock (AL) and Brienzlig (BR) whitefish. Shown are median (vertical line), mean (cross), 25% and 75% quartiles (box), and minimum and maximum values (error bars).
Figure 3. Box plots of relative RT qPCR gene expression measurements of Ictacalcin in head kidney tissue of normal (nor) and deformed (def) Albock (AL) and Brienzlig (BR) whitefish. Shown are median (vertical line), mean (cross), 25% and 75% quartiles (box), and minimum and maximum values (error bars).
Ijerph 08 02706f3
Figure 4. Transcript abundance patterns for 102 genes in the liver (a) and 71 genes in the head kidney (b) found in common across the significant GO categories in comparisons of deformed and normal whitefish of the Albock and the Brienzlig form, respectively. As clustering algorithm, average linkage (UPGMA) was applied. Centered Pearson correlation coefficients were used as distance measures. The dendrogram on the left groups genes with similar expression and the dendrogram on the top separates populations. The colour bar indicates arbitrary expression units. For each gene probe, the mean gene expression levels of each group divided by the global mean across all four groups are shown. Yellow colour (light) designates low expression and purple colour (dark) represents high expression levels respectively.
Figure 4. Transcript abundance patterns for 102 genes in the liver (a) and 71 genes in the head kidney (b) found in common across the significant GO categories in comparisons of deformed and normal whitefish of the Albock and the Brienzlig form, respectively. As clustering algorithm, average linkage (UPGMA) was applied. Centered Pearson correlation coefficients were used as distance measures. The dendrogram on the left groups genes with similar expression and the dendrogram on the top separates populations. The colour bar indicates arbitrary expression units. For each gene probe, the mean gene expression levels of each group divided by the global mean across all four groups are shown. Yellow colour (light) designates low expression and purple colour (dark) represents high expression levels respectively.
Ijerph 08 02706f4aIjerph 08 02706f4b
Figure 5. Log2 normalized gene expression differences between normal and deformed fish in Albock (AL; black) and Brienzlig (BR; grey) for shared genes among significant GO categories. (a) genes present in immune system related GO categories in the liver; (b) genes present in GO terms encoding for “proteolysis” and “extracellular region” also in the liver; (c) genes present in immune system related GO terms in the head kidney; and (d) genes present in GO categories related to GTPase activation and Ras protein signal transduction also in the head kidney. The presence or absence of a particular gene in the respective GO terms is indicated by black lines. Note that for (a) and (b) only genes present simultaneously in both whitefish forms are shown (i.e., not all genes assigned to specific GO terms). Finally, in (b) the explicit marking of the first 18 genes points out that these genes are also present among the shared genes shown in (a).
Figure 5. Log2 normalized gene expression differences between normal and deformed fish in Albock (AL; black) and Brienzlig (BR; grey) for shared genes among significant GO categories. (a) genes present in immune system related GO categories in the liver; (b) genes present in GO terms encoding for “proteolysis” and “extracellular region” also in the liver; (c) genes present in immune system related GO terms in the head kidney; and (d) genes present in GO categories related to GTPase activation and Ras protein signal transduction also in the head kidney. The presence or absence of a particular gene in the respective GO terms is indicated by black lines. Note that for (a) and (b) only genes present simultaneously in both whitefish forms are shown (i.e., not all genes assigned to specific GO terms). Finally, in (b) the explicit marking of the first 18 genes points out that these genes are also present among the shared genes shown in (a).
Ijerph 08 02706f5aIjerph 08 02706f5bIjerph 08 02706f5c
Table 1. Number and frequencies of genes showing difference in expression levels in head kidney and liver tissue observed in pairwise gene-by-gene analyses at three different nominal P-value thresholds, (i) 5% level; (ii) 1% level and (iii) after correcting for multiple hypothesis testing with FDR ≤ 5%. Comparisons were carried out between the following pairs of samples: (A) Albock normal vs. Albock deformed; (B) Brienzlig normal vs. Brienzlig deformed; (C) all normal (pooled Albock and Brienzlig) vs. all deformed and (D) all Albock (pooled normal and deformed) vs. all Brienzlig. The total number of expressed genes (GEXP) per tissue is indicated. Additionally, for the 5% and 1% levels the corresponding FDR at the respective thresholds are given in parentheses.
Table 1. Number and frequencies of genes showing difference in expression levels in head kidney and liver tissue observed in pairwise gene-by-gene analyses at three different nominal P-value thresholds, (i) 5% level; (ii) 1% level and (iii) after correcting for multiple hypothesis testing with FDR ≤ 5%. Comparisons were carried out between the following pairs of samples: (A) Albock normal vs. Albock deformed; (B) Brienzlig normal vs. Brienzlig deformed; (C) all normal (pooled Albock and Brienzlig) vs. all deformed and (D) all Albock (pooled normal and deformed) vs. all Brienzlig. The total number of expressed genes (GEXP) per tissue is indicated. Additionally, for the 5% and 1% levels the corresponding FDR at the respective thresholds are given in parentheses.
Liver: GEXP = 7,013Comparison5%FDR1%FDRFDR ≤ 5%
A346 (0.049)178 (0.011)0.880
B74 (0.011)110 (0.001)10
C177 (0.025)133 (0.005)10
D3,407 (0.486)0.12,407 (0.343)0.032,764 (0.394)

Head kidney: GEXP = 17,834Comparison5%FDR1%FDRFDR ≤ 5%

A1,162 (0.065)0.77264 (0.015)0.680
B2,548 (0.143)0.35528 (0.030)0.310
C651 (0.037)0.95102 (0.006)0.950
D12,416 (0.696)0.0110,382 (0.582)0.0212,820 (0.719)
Table 2. Gene Ontology (GO) categories that comprised genes that showed significant over-representation in head kidney and liver tissues at the 1% level of low P-value ranking among all genes after refinement analysis (see Material and Methods section) in comparisons of whitefish with deformed and normal gonads. Only GO categories with number of genes < 200 are shown.
Table 2. Gene Ontology (GO) categories that comprised genes that showed significant over-representation in head kidney and liver tissues at the 1% level of low P-value ranking among all genes after refinement analysis (see Material and Methods section) in comparisons of whitefish with deformed and normal gonads. Only GO categories with number of genes < 200 are shown.
(a)Liver: Albock normal vs. Albock deformed
GO root nodeNameGO IDGenesP
biological_processimmune responseGO:0006955700.000
biological_processhumoral immune responseGO:0006959220.001
biological_processinnate immune responseGO:0045087250.001
biological_processactivation of immune responseGO:0002253220.002
biological_processProteolysisGO:00065081600.004
cellular_componentintegral to plasma membraneGO:0005887600.009
biological_processimmune effector processGO:0002252320.009

(b)Liver: Brienzlig normal vs. deformed
GO root nodeNameGO IDGenesP

cellular_componentextracellular regionGO:00055761970.002
cellular_componentproteinaceous extracellular matrixGO:0005578310.002
biological_processimmune system processGO:00023761070.006
biological_processcell adhesionGO:0007155850.009

(c)Head kidney: Albock normal vs. deformed
GO root nodeNameGO IDGenesP

biological_processregulation of cell proliferationGO:00421271140.000
biological_processpositive regulation of cell proliferationGO:0008284540.000
biological_processembryonic developmentGO:0009790860.000
biological_processprotein import into nucleusGO:0006606230.001
molecular_functionATP-dependent helicase activityGO:0008026200.001
biological_processcell migrationGO:0016477520.001
biological_processembryonic development ending in birth or egg hatchingGO:0009792580.002
molecular_functionsmall GTPase regulator activityGO:0005083490.004
biological_processregulation of Ras protein signal transductionGO:0046578270.004
biological_processcellular protein complex assemblyGO:0043623410.009
biological_processregulation of immune system processGO:0002682450.009
biological_processcentral nervous system developmentGO:0007417450.010

(d)Head kidney: Brienzlig normal vs. deformed
GO root nodeNameGO IDGenesP

molecular_functionGTPase activator activityGO:0005096360.003
biological_processleukocyte differentiationGO:0002521320.006
biological_processgamete generationGO:0007276400.006
biological_processT cell activationGO:0042110310.008
cellular_componentsoluble fractionGO:0005625320.009
biological_processmuscle contractionGO:0006936280.010

Share and Cite

MDPI and ACS Style

Bittner, D.; Cossins, A.R.; Segner, H.; Excoffier, L.; Largiadèr, C.R. Identification of Candidate Genes and Physiological Pathways Involved in Gonad Deformation in Whitefish (Coregonus spp.) from Lake Thun, Switzerland. Int. J. Environ. Res. Public Health 2011, 8, 2706-2733. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph8072706

AMA Style

Bittner D, Cossins AR, Segner H, Excoffier L, Largiadèr CR. Identification of Candidate Genes and Physiological Pathways Involved in Gonad Deformation in Whitefish (Coregonus spp.) from Lake Thun, Switzerland. International Journal of Environmental Research and Public Health. 2011; 8(7):2706-2733. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph8072706

Chicago/Turabian Style

Bittner, David, Andrew R. Cossins, Helmut Segner, Laurent Excoffier, and Carlo R. Largiadèr. 2011. "Identification of Candidate Genes and Physiological Pathways Involved in Gonad Deformation in Whitefish (Coregonus spp.) from Lake Thun, Switzerland" International Journal of Environmental Research and Public Health 8, no. 7: 2706-2733. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph8072706

Article Metrics

Back to TopTop