This is the first network level analysis of an environmental study integrating multilevel omic datasets. We discovered that the overall molecular state of the flounder liver (transcriptomics and metabolomics) is representative of the chemical contaminant burden of the sediments. Network reconstruction showed that the interface between transcriptional and metabolic network domains is linked to fish morphometric indices and is predictive of environmental exposure. In-depth analyses of predictive networks have identified putative novel pathways representative of responses to exposure. This approach provides a framework both for prediction of chemical pollutants in complex mixtures and for prediction of the health outcomes for exposed animals.
The molecular state of flounder liver samples reflects chemical exposure
The chemical exposures predicted from CTD interactions were partly confirmed by chemical data () despite the complexity of the environment, potentially including mixture effects, bioaccumulation and non-chemical stressors. Additional stressors were indicated that had not been chemically measured. Taking the Brunsbuttel site as an example, ethinyl-estradiol was a predicted contaminant and serum VTG protein, a canonical marker of endocrine disruption, was induced relative to the Alde (). Perfluorooctane sulfonic acid 
and other persistent organic pollutants including PCBs, dieldrin and endosulfan 
have been detected at elevated concentrations in the Elbe estuary and floodplain, and were all identified by our approach. Additional chemicals highlighted included systhane and vinclozolin fungicides, the halogenated aromatic hydrocarbon pesticide lindane, chlorine and tetradecanoylphorboyl acetate (TPA). It is uncertain whether these compounds are in fact present at this site as, for example, the presence of TPA appears unlikely. However TPA is a well-known tumour promoter 
, so detection of its associated gene expression changes might be viewed as a biomarker of effect, not necessarily of a specific exposure. At a number of sites flavonoids and flavonols, such as epicatechin gallate, were predicted, potentially indicating plant-derived exposures not of anthropogenic origin. At Morecambe Bay and the Tyne the prediction of paraquat perhaps reflected an oxidative stress response rather than the presence of this particular compound.
These results support the use of a knowledge-based approach to infer chemical exposure profiles from molecular responses and validate the underlying assumptions in the study. Predictions from interrogation of the CTD database (; Table S3
) differed between sites suggesting that the approach can be sufficiently sensitive to specific differences in the exposure profiles. However, we do not propose that these associations necessarily indicate the presence of each specific contaminant at each site, for example ‘tobacco smoke pollution’ in the Mersey, we instead hypothesise that these represent the effects of related stressors, for example, AhR inducers at the Mersey site.
Interpretation of the molecular interaction network
The development of a modular network, representing the integration between molecular and physiological readouts, provided us with an interpretive framework to analyse the complex molecular signatures linked to exposure. One of the most interesting findings is that the modules that predict environmental exposure with greatest accuracy represent the interface between metabolite and transcriptional networks and link to higher level indicators of fish health, such as condition factor and hepatosomatic index ().
Consistent with this observation, network modules at the interface between metabolite and transcriptional networks were also differentially regulated in response to single chemical laboratory exposures. It should be borne in mind that the environmentally sampled fish have been chronically exposed to pollutants, and that chronic exposure can result in different responses than acute exposure 
. In addition bioavailability, mixture effects, metabolism and bioaccumulation affect compound-specific responses within the livers of these fish. This is illustrated by the modules containing genes that responded to 16-day treatments of flounder with individual toxicants (, Figure S2
). While all toxicants induced changes in the metabolic-interface genes, they also affected the secondary area of the network that related more to acute stress and immune response (, area B), in contrast to the differences between environmental sites, where only one module (40) in this area was affected.
The characterisation of transcripts and metabolites that differed between sites was undertaken to provide insights into the molecular mechanisms that they describe, and to inform on the potential health outcomes for the fish. Canonical pathways that contributed to these differences included those relevant to metabolism of toxicants; AhR signalling, metabolism of xenobiotics by cytochromes P450, the NRF2-mediated oxidative stress response, glutathione metabolism and bile acid bioysnthesis (). Together these describe phase I and II metabolism of xenobiotics, such as aromatic hydrocarbons, and their excretion via the bile. Additional endobiotic metabolic pathways were affected. Changes in glycolysis, pyruvate metabolism, the citric acid cycle and oxidative phosphorylation implied disturbances to the energy pathways of the liver that could reflect the energetic requirements of xenobiotic metabolism and lead to further metabolic disruption. Changes in amino acid synthesis and proteasomal protein degradation also indicated reorganisation of metabolism.
This change in metabolic state and gene expression could be viewed as a successful compensatory response to toxicants and thus of little concern for the health of individual fish and these fish populations. Further examination of the annotation of transcripts and metabolites differing between sites implied that this hypothesis was false. As illustrated in , and shown in and , there is a remarkable overlap between site-predictive modules and modules associated with hepatocellular carcinoma (HCC). Additionally, liver cholestasis -annotated modules overlapped with HCC and site predictive modules and this area of the network was highly associated with bile acid biosynthesis. Apart from this metabolic interface group only one other module (module 40) was predictive of site. This was also associated with hepatocellular carcinoma, and additionally with liver fibrosis, indicative of chronic liver damage, and occurred in an area of the network associated with inflammation. Therefore flounders inhabiting differentially contaminated sites show transcript and metabolite changes that have been associated with liver carcinogenesis in mammals. A question remains as to whether this simply represents the detection of HCC in the liver samples, as histopathology data were unavailable for the fish sampled off Germany. By comparison with studies of tumours from the closely related flatfish dab (Limanda limanda
) this does not appear to be the case. In dab tumours the metabolites choline, phospocholine and glycine were reduced in concentration and lactate increased, an indication of the switch to anaerobic metabolism in the bulk tumours 
. In, for example, the Brunsbuttel samples compared with non-tumour bearing Alde fish, choline, phosphocholine and glycine increased, and lactate decreased. Additionally, transcripts for ribosomal proteins showed co-ordinated induction in bulk tumours from dab, indicative of proliferation 
, but no such induction was apparent from the present samples. The changes in gene expression and metabolites detected in this study do not recapitulate those found in bulk tumours, and may be viewed as indicating either an earlier stage of tumourigenesis or a permissive micro-environment in which hyperplastic tissue may form and lead to tumour formation.
Network analysis reveals potentially novel response pathways
Ingenuity networks, based on mammalian interaction data, permitted more detailed biological characterisation of the site-associated modules. Complete pathways were not recapitulated by these analyses, as only a minority of the transcripts and metabolites from flounder liver were examined. Nevertheless, the analyses highlighted important processes and inferred key regulators. Here the most significant network derived from site-predictive modules is discussed in detail and additional networks are discussed in terms of their key inferred regulators.
The most striking finding from the Ingenuity analyses was the co-ordinated repression of proteasomal subunit genes at the Brunsbuttel site (; Figure S3A
1). This was not so marked at other sites (Figure S1
), indeed at Morecambe Bay these genes were induced in comparison with the Alde fish. Proteasome maturation protein (POMP) has been found to be a critical regulator of proteasomal activity 
and has been shown to be repressed by the halogenated aromatic hydrocarbon 2,3,7,8-tetrachlorodibenzodioxin (TCDD) in an AhR-independent manner 
. Although TCDD concentration was not measured, the mean expression of proteasomal genes was inversely correlated (r
−0.79) with fish liver PCB concentrations but did not correlate well with sediment PAH or PCB concentrations. Tyne fish, for example, displayed relatively high proteasomal gene expression and had low liver PCB but high PAH concentrations (, Figure S1
). Therefore the repression of proteasomal genes may represent a halogenated aromatic hydrocarbon-related response (). In trout Oncorhynchus mykiss
, a proteasome inhibitor reduced PAH-dependent CYP1A induction 
, in contrast to mammalian studies 
. This difference may contribute to the lower inducibility of CYP1A in flounder in comparison with many mammals. Ingenuity analysis also predicted an interaction between the proteasome and NF kappa B, a key regulator of mammalian hepatocarcinogensis 
. The proteasome represses NK kappa B activation, and potentially disruption of proteasomal activity could have extensive additional effects on intracellular protein levels due to its role in the degradation of numerous proteins. We found no significant changes in NF kappa B gene expression between sites, and the consequences of putative activation at the Brunsbuttel site, and repression at the Morecambe site, due to changes in the proteasome, are difficult to predict, as in the early stages of carcinogenesis NF kappa B can have a protective effect, whereas in later stages it can promote tumourigenesis 
Relationship between the AhR pathway and proteasome.
From the Ingenuity networks a number of key regulatory molecules were inferred. These included insulin (, Figure S3A2, S3A4, S3B1
), estrone, luteinizing hormone (LH) and follicle stimulating hormone (FSH) (Figure S3A3 and S3A6
), platelet derived growth factor beta (PDGFBB) (, Figure S3A2, S3B1
), transforming growth factor beta (TGF-beta) (Figure S3A
8), vascular endothelial growth factor (VEGF) (, Figure S3A7, S3B1
), tumour necrosis factor (TNF) (Figure S3A
5), and angiotensinogen (, Figure S3B
Insulin, in fish as in mammals, is a key hormonal regulator of energy, glucose and lipid metabolism, all pathways that were identified as affected by sampling site. By the Ingenuity networks it was linked to protein kinases, metabolites (including glucose and lactate) and the glucose transporter SLC2A4. The most obvious explanation for changes in insulin and related parameters would be differences in diet between fish from different sites. Amino-acid levels are more important regulators of insulin in carnivorous fish such as the flounder than sugars 
. Dietary parameters would be expected to be highly variable depending upon recent feeding history of the fish, which was unknown for these individuals. However, insulin can also be modulated by exposure to toxicants including organophosphates 
that was suggested to lead to an increase in lipogenesis, in agreement with our observations of phospholipidosis in fish from polluted sites (). Mild estrogenic endocrine disruption was suggested by VTG induction in Brunsbuttel fish, and networks shown in Figures S3A6 and S3A3
inferred that estrogen receptor alpha (ESR1), FSH and LH target genes were modulators of the different responses between sampling sites. ESR1 and HNF alpha were linked in Figure S3A
6 and both are involved in hepatic cholestasis, indeed EE2-induced hepatotoxicity has been linked to alterations in bile acid biosynthesis in mice 
PDGFBB is the dimeric form of platelet derived growth factor beta (PDGF-B). Notably, PDGF-B over-expressing mice spontaneously developed liver fibrosis 
, and PDGF-BB was inferred as part of the network deriving from the liver fibrosis-annotated module 40 in our analysis. Additionally PDGF-B over-expressing mice developed hepatocellular carcinoma in response to phenobarbital and diethylnitrosamine treatment and induced TGF-beta and VEGF expression. TGF-beta was inferred to be an important regulator in site-specific responses (Figure S3A
8) and is a well-known mediator of cancer initiation, progression and metastasis, via interaction with the inflammatory response 
. Furthermore, the pro-inflammatory cytokine TNF-alpha, an initiating signal for the innate immune response in fish as well as mammals 
, was also identified by Ingenuity analysis (Figure S3A
5). Release of TNF alpha from Kupffer cells leads to hepatocyte cell death, regeneration and fibrosis that can lead to hepatocellular carcinoma 
. VEGF, best known as a stimulator of angiogenesis, was also highlighted in both the fibrosis-related and carcinoma-related sections of the network, and was linked with cell cycle, oncogenes and tumour suppressor genes (CDKN1A, TP53, MYC). Angiogenesis is a key requirement for the transition from fibrosis to hepatocellular carcinoma 
Angiotensinogen (AGT) is the precursor of angiotensin and was found to be repressed at all sites in comparison to the Alde reference site (Table S1
). Angiotensin is a signal for vasoconstriction in mammals and in fish its expression is related to osmoregulation 
with repression in liver in response to higher salinity. As the sampling sites differed in salinity, alteration of AGT transcription was not a surprising finding. As shown in , AGT was a member of the fibrosis-related module 40 and was predicted to form part of a complex network with VEGF, PDGF and intracellular kinases. Angiotensin has indeed been linked to stimulation of inflammatory liver fibrosis 
, via fibroblast proliferation and production of inflammatory cytokines and growth factors, including TGF-beta. Inhibition of the angiotensin system by antagonism of its receptor 
or inhibition of angiotensin-converting enzyme 
has been shown to reduce hepatic fibrosis.
VEGF, TGF beta, TNF alpha, PDGF and AGT are all intimately related to the progression of fibrosis to cirrhosis and hepatocellular carcinoma in mammals. These molecules were all highlighted as important regulators of the differences between molecular profiles of flounder livers from different sampling sites using an unbiased approach combining network inference and predictive algorithms.
A combination of omics, multiple biomarkers and bioinformatics were used to identify and characterise hepatic molecular changes between fish sampled from several environmental sites. Based on these data, parasite infection, fish morphology and genetics do contribute to the differences between sites, but do not explain the majority of changes seen. For example, within-site tests showed that morphometric parameters and parasite infections could be significantly associated only with a small proportion (<3%) of the gene expression differences between sites (Table S1
, Table S4
). Taken as a whole with our previous studies 
, we find that anthropogenic chemical contamination of the marine environment is a major factor in explaining the molecular differences between fish sampled from these sites.
The different methodologies employed displayed different strengths and weaknesses. Histopathology was a good guide to broad levels of pollution effect, but provided little information upon the nature of the contaminant profile. Protein biomarkers and enzyme activities were useful for categorising sites by major classes of toxicant, but gave little information on the potential health outcomes. 1H NMR metabolomics showed low technical variability, and metabolite profiles alone were more predictive of sampling site than gene expression profiles alone, however the annotation of metabolites is not yet well advanced, limiting the functional information currently available. Transcriptomics exhibited higher variability than metabolomics, but was more informative due to better annotation. Overall the methodologies were highly complementary, allowing analyses that would be impossible if one were limited to a single technique.
The gene expression signatures associated with fish from each sampling site were used to predict the presence of chemical contaminants using the CTD gene expression-chemical interaction database. Mixture effects, other environmental influences and the similarity of certain stressors, such as the metals, might be expected to confound this approach. Additionally the incomplete nature of the flounder microarray and the CTD database and the limited numbers of samples for certain sites, which is a common issue in field studies, reduce the potential of this analysis. Therefore we did not expect to predict all environmental contaminants by this method. While this approach was useful with the current dataset, it may be expected to improve in future as both the CTD database and transcriptomic data become more comprehensive.
Data integration and network analyses were essential; both to predicting health outcomes and to identifying and examining affected biological pathways. They allowed visualisation of the highly complex dataset and facilitated comparison of the effects of different stimuli upon the model system. Modules associated with specific parameters could then be examined in detail, utilising interaction databases (Ingenuity) for further characterisation. Detailed examination of these networks illustrated the changes detected by broader classification of modules by annotation terms. In addition to potential interactions with diet and salinity, the majority of networks contained key regulators of inflammation, hepatic fibrosis and hepatocellular carcinoma. Therefore we propose that network biology approaches can lead to the identification of health impacts of environmental pollutants upon non-model organisms.
The molecular differences between reference and contaminated sampling sites were associated with carcinogenesis, and this outcome is supported by previous histopathology 
. Flatfish hepatic histopathology has long been associated with chemical contamination 
and our results demonstrate the linkages between toxicants and histopathology via alterations in molecular signalling pathways and metabolism.