Common Gene Expression Patterns in Environmental Model Organisms Exposed to Engineered Nanomaterials: A Meta-Analysis.

The use of omics is gaining importance in the field of nanoecotoxicology; an increasing number of studies are aiming to investigate the effects and modes of action of engineered nanomaterials (ENMs) in this way. However, a systematic synthesis of the outcome of such studies regarding common responses and toxicity pathways is currently lacking. We developed an R-scripted computational pipeline to perform reanalysis and functional analysis of relevant transcriptomic data sets using a common approach, independent from the ENM type, and across different organisms, including Arabidopsis thaliana, Caenorhabditis elegans, and Danio rerio. Using the pipeline that can semiautomatically process data from different microarray technologies, we were able to determine the most common molecular mechanisms of nanotoxicity across extremely variable data sets. As expected, we found known mechanisms, such as interference with energy generation, oxidative stress, disruption of DNA synthesis, and activation of DNA-repair but also discovered that some less-described molecular responses to ENMs, such as DNA/RNA methylation, protein folding, and interference with neurological functions, are present across the different studies. Results were visualized in radar charts to assess toxicological response patterns allowing the comparison of different organisms and ENM types. This can be helpful to retrieve ENM-related hazard information and thus fill knowledge gaps in a comprehensive way in regard to the molecular underpinnings and mechanistic understanding of nanotoxicity.

Introduction

The exponential increase in production of engineered nanomaterials (ENMs, i.e., particles and fibers in which at least one dimension is <100 nm) in the last decade has raised concerns about their impact on human and environmental health.[1−4] While studies in a number of model species have demonstrated that ENMs induce toxicity at the phenotype level, and the molecular and cellular mechanisms that lead to ENM-induced toxicity are not yet well understood.[5−8]

In an attempt to study the molecular mechanisms in an unbiased way, the so-called omics approaches (e.g., transcriptomics, proteomics, and metabolomics) have become more often used in nanotoxicology in recent years.[9−15] In fact, several reviews have been published, in which authors gathered individual omics studies and summarized their results to discuss the molecular mechanisms associated with specific ENMs and biological species.[7,16,17] A common conclusion among these reviews is that oxidative stress is the most prevalent molecular mechanism associated with ENM toxicity, found in most omics-based ENM studies, followed by metal homeostasis (metal ENMs dominate the reviewed studies), plant hormone homeostasis, and immune response (vertebrates), which were found to be perturbed in fewer studies. These findings fit well to targeted ENM toxicity studies, which often find an increase in reactive oxygen species (ROS), perturbation of metal ion uptake, and infiltration of immune cells in tissues exposed to ENMs.[18−22] While the fit between the omic and targeted approaches adds confidence to the use of omics, thus far, omics studies of ENM-induced toxicity have not lead to the discovery of any mechanisms beyond what was already known.

Importantly, the published reviews did not consider the differences in experimental design and statistical methods used to analyze the omics data sets between the reviewed studies. It has been shown before that performing a meta-analysis based on data sets that have been analyzed by a standard analysis pipeline can reveal new mechanisms of toxicity that the initial studies have missed.[23,24] The aim of this study was to perform such a meta-analysis for all environmental toxicology-relevant ENM transcriptomic toxicity studies in order to find molecular mechanisms of toxicity common for different species exposed to different ENMs. To do this, we built a pipeline for integrative analysis of transcriptomic data from different microarray platforms, which included inference of differentially expressed genes (DEGs) and functional enrichment analysis. We used the pipeline to reanalyze nanotoxicity gene expression studies from three different species, Arabidopsis thaliana, Caenorhabditis elegans, and Danio rerio, which were exposed to different ENMs. In the article, we describe the pipeline, which is available as an R package, and a semiquantitative analysis of the most common ENM mechanisms discovered using our approach for each of the species and ENM types studied.

Methods

Data Acquisition

We searched for all publicly available transcriptomic (microarray and RNA-Seq) data sets, describing gene expression response upon exposure of environmental model organisms to ENMs. This was done by querying commonly used depositories of gene expression data for the following search terms: “nano*”, “particle”, “NP*”, “ENM”, and “quantum”. We found only two public data repositories comprising relevant data sets: Gene Expression Omnibus (GEO)[25] and ArrayExpress.[26] After the initial search, all nonenvironmentally relevant studies (mammals) and studies on bacteria were filtered out; the final search included all plant and aquatic species, including fungi (search date: 12.01.2018). All matching data sets were retrieved, and the experimental information (e.g., assay type, exposure concentration, and time) was summarized (Table S1).

Data Analysis Using a Statistical Platform

The content and structure of data sets were manually checked, and only the data sets which fulfilled the following three criteria were selected for analysis

Minimum of three replicates.

Available genomic annotation (gene ontology (GO)).

Complete technical information and raw data needed for data reanalysis.

All data sets which passed these quality criteria (Table S1) were further processed. A pipeline was developed that reanalyzes the raw data of all data sets in order to unify data normalization, significance testing of differential expressed genes (DEG), and functional analysis (Figure ). The pipeline can handle heterogenic data from different microarray manufacturers (e.g., Agilent, Affymetrix, and NimbleGen) and from different environmental relevant model organisms, such as A. thaliana, C. elegans, and D. rerio. It can handle single channel and dual channel arrays and automatically recognizes the format version of the binary files including Agilent, Celera, and Genepix files. The built R package, with detailed instructions on how to install and use the package, a vignette and example microarray files, is available at https://github.com/alxbetz/mira. While we only provide the final results of our analysis in the manuscript, we have also added all the intermediate results, provided as Supporting Information files.

Figure 1

Pipeline work flow and the applied tools to process and analyze different data types. Starting with the raw microarray image reader files, the data are first background corrected[84] and then, one of two methods is applied for between-sample normalization: for Agilent arrays, quantile normalization[85] is used and for Affymetrix arrays, robust multiarray averaging is used. Then, the probe sequences are realigned to the transcriptome (Arabidopsis thaliana: ensemble version 91.10, Caenorhabditis elegans: 93.260, and Dario rerio: 93.11.) using “bowtie”.[86] We allow 1 mismatch per 20 bases to account for variation between genome versions. If one probe matches to multiple targets with the same accuracy, we discard it, and in cases where multiple probes map to the same target, the normalized fluorescent intensities are averaged per gene. To assess differential gene expression, we fit a gene-wise linear model using “limma”, followed by a t-test and FDR multiple-testing correction.[87] Alternatively, a FC ranking can be calculated using FCROS. Finally, either the list of DEGs or the FC rank can be used as input for the functional enrichment analysis using a Kolmogorov–Smirnov-test implemented in “topGO”. Herein, functional analysis based on a linear model (dotted arrow) was not possible because of the low number of DEGs in some studies. The functional analysis was performed with FCROS (solid arrow).

Pipeline Work Flow to Process Microarray Data

The raw data output from the different microarray image analyzers was normalized and analyzed in a three-step work-flow, as detailed in Figure . Most of the studies, which we selected, revealed a low effect size. Application of a standard, uniform p-value and fold-change (FC) cutoff resulted in a very low number of detected DEGs (Section ). Therefore, we decided to assess functional enrichment based on fold-change rank ordering statistics (FCROS)[27−30] instead. Briefly, in the first step, the FC rank of each gene in all pairwise comparisons of treatment and control replicates is computed. Then, the mean of ranks per gene across all comparisons is calculated. The resulting distribution is approximately normal, and the mean and variance of this empirical score distribution are used as parameter estimates for a normal distribution. An f-value is calculated based on this normal distribution and the mean rank value.

Functional Enrichment and Visualization of Response Profiles

To assess the GO term enrichment, we used a Kolmogorov–Smirnov test[31] implemented in the R package “topGO”[32] together with the FC rank. We selected the 100 top-ranked GO terms with the lowest f-value for each contrast and then computed the overlap of GO-terms between contrasts, organisms, and ENM types. Response profiles were determined by assessing the GO-terms which are associated with the commonly observed toxicity mechanisms. The score value indicates the number of GO-terms which were found for each category and was normalized to the number of tested contrasts for each organism or ENM type. The different contrasts are listed in Table .

Table 1

Total Number of DEGs for Each Tested Contrast Using Standard Statistical Thresholds (FDR < 0.05, log FC > 1.2) in Comparison to DEG Numbers Reported in the Respective Publicationsa

study	species	sample source	treatment	DEG (pipeline)	DEG (reported)	statistics (reported)	reference
GSE80461	A. thaliana	leaves	CeO2	0	0–221b	varyingb	Tumburu et al. 2016
GSE80461	A. thaliana	roots	CeO2	44	24–1066b	varyingb
GSE80461	A. thaliana	leaves	TiO2	10	38–2196b	varyingb
GSE80461	A. thaliana	roots	TiO2	779	136–2276b	varyingb
GSE46958	A. thaliana	roots	Au	1187	n.a.	FC > 2; p < 0.05	Taylor et al. 2014
GSE32521	C. elegans	L3 larvae	Au	12	797	FC > 1.5; p < 0.05	Tsyusoko et al. 2012
GSE70509	C. elegans	L1 larvae	Ag (aged)	0	n.a.	n.a.	not published
GSE70509	C. elegans	L1 larvae	AgNO3	0	n.a.	n.a.
GSE70509	C. elegans	L1 larvae	Ag (pristine)	0	n.a.	n.a.
NERC	C. elegans	L1 larvae	AgNO3	0	0–213	varyingc	Schultz et al. in prep.
NERC	C. elegans	L1 larvae	Ag52PVP (uncharged)	0	0–210	varyingc
NERC	C. elegans	L1 larvae	Ag12PVP (uncharged)	0	0–297	varyingc
NERC	C. elegans	L1 larvae	Ag12MUA (negative)	0	0–449	varyingc
NERC	C. elegans	L1 larvae	Ag12AUT (positive)	0	0–441	varyingc
NERC	C. elegans	L1 larvae	Ag (unfunctionalized)	271	604–2459	varyingc
GSE73427	D. rerio	larvae (48 hpf)	Si_BaP	3	n.a.	n.a.	Duan et al. 2016/Hu et al. 2016
GSE73427	D. rerio	larvae (48 hpf)	Si	0	2515	FC > log 2; p < 0.05
GSE61186	D. rerio	larvae (120 hpf)	AgNO3	53	n.a.	n.a.	not published
GSE61186	D. rerio	larvae (120 hpf)	Ag	5	n.a.	n.a.
GSE77148	D. rerio	larvae (96 hpf)	ZnO	0	445	FC > 1.5; p 0.05	Choi et al. 2016
GSE77148	D. rerio	larvae (96 hpf)	ZnSO4	0	653	FC > 1.5; p 0.05
GSE50718	D. rerio	larvae (72 hpf)	Ag (150 nm)	43	7538 (total)	n.a.	not published
GSE50718	D. rerio	larvae (72 hpf)	Ag (50 nm)	156	7538 (total)	n.a.
GSE50718	D. rerio	larvae (72 hpf)	AgNO3	329	7538 (total)	n.a.
GSE41333	D. rerio	larvae (48 hpf)	PAMAM-G3	35	230	FC > log 2; p < 0.001	Oliveira et al. 2014
GSE41333	D. rerio	larvae (48 hpf)	PAMAM-G4	499	220	FC > log 2; p < 0.001

The applied statistics and the reference are reported if available.

Varying statistics: FC > 2 & (p 0.01–0.1); p 0.05 & (FC 2–8).

Varying statistics: p 0.05 & (FC 1.4 and 2); with/without Benjamin Hochberg MSC.

Results & Discussion

Available Nanorelated Transcriptomics Studies

The search for transcriptomic studies in GEO and ArrayExpress returned 46 nanorelated and environmentally relevant publicly available studies (Table S1). The majority of studies were performed with metal-based nanoparticles (silver, titanium, zinc, gold, cerium, and copper), while a few studies were based on carbon (e.g., multiwalled nanotubes), silica, quantum dots, or polyamidoamine (PAMAM) (Figure left). Impact of these ENMs was tested in various environmental organism groups: plants (A. thaliana, Solanumly copersicum), fish (Danio rerio, Pimephales promelas, Oryzias latipes), nematodes (C. elegans), crustaceans (Daphnia magna, Hyalella azteca), clitellate (Eisenia fetida, Enchytraeus albidus), algae (Chlamydomonas reiinhardtii), and molluscs (Mytilus galloprovincialis) (Figure right).

Figure 2

Coverage of nanomaterials (left) and organisms (right) of publicly available transcriptomic studies found in GEO and ArrayExpress. ENMs marked as “others” were copper (4%), quantum dots (2%), PAMAM (2%), and polystyrene (2%). Group of plants included Arabidopsis thaliana (8) and Solanum lycopersicum (3), fish were Dario rerio (8), Pimephales promelas (2), and Oryzias latipes (1), the group of bacteria included Escherichia coli (3), Nitrosomonas europaea (3), Bacillus cereus (1), Pseudomonas aeruginosa (2), and Shewanella oneidensis (1), all nematode studies were performed with Caenorhabditis elegans, crustaceans included Daphnia magna (6) and Hyalella Azteca (1), clitellata were Eisenia fetida (2) and Enchytraeus albidus (2), and the group of “others” included Chlamydomonas reinhardtii (2) and Mytilus galloprovincialis (1).

We performed a similar search for scientific articles in Web of Science and Scopus, using the same keywords as for the data sets and then filtering using the term transcript*. In total, over 600 different publications were found. While the majority of the publications were based on mammalian data, approximately 150 studies were environmentally relevant. The number of scientific studies is therefore over three-fold greater than the number of data sets in public depositories, the conclusion being that the scientific community is still not doing enough to openly share its data.

Study Selection Based on Three Defined Quality Criteria

Three quality criteria were established to select appropriate studies for the meta-analysis: a minimum of three replicates, established GO, and availability of the complete raw data in a processable format (Figure ). From the initial list, 10 studies had less than three replicates, however, a minimum of three biological replicates should be standard in order to have the statistical power to calculate DEGs and account for the high variances of gene expression levels. Thirteen studies were performed with organisms that have no GO annotation available; this includes, for example, D. magna, E. albidus, O. latipes, or P. promelas. Whereas functional analysis with nonmodel organisms using related sets of annotation is possible, misinterpretation may be a consequence, for example, of overrepresentation of highly-conserved marker genes.[33] Thus, we excluded these. For 12 studies, the technical information or raw data were insufficient: two studies provided no cDNA information, four studies had no probe sequences, three studies were without raw data, and for three studies, the data organization did not allow data processing (e.g., duplicates within the data frames) (Table S1).

Figure 3

Defined selection criteria for the selection of studies suitable for the meta-analysis: in total, 46 data sets were identified, 10 studies had less than three replicates, for 13 studies the genomic or GO annotation was insufficient, and for 12 studies, the reported data were insufficient. Eleven studies were selected for the meta-analysis.

Studies with incomplete reporting were removed as proper reporting is crucial to reproduce experiments and also to perform high-quality meta-analysis. In the future, it is vital that authors report data compliant with the existing reporting guidelines, such as minimum information about a microarray experiment (MIAME).[34] The checklists provided by MIAME ensure that the requested information is adequate to process and reanalyze the data.[35]

Overall, we selected 11 complete microarray data sets for reanalysis, comprising three environmentally relevant model organisms (A. thaliana, C. Elegans, and D. rerio) exposed to different ENMs (silver, titanium, cerium, gold, silica, zinc, and PAMAM) (Table S2).

Pipeline

The statistical analysis, reported in the publications which utilized the complete microarray data sets, differed in terms of data processing (background correction, normalization, and alignment) and the choice of statistical thresholds to calculate DEGs.[36−43] We established an R-scripted pipeline that can semiautomatically process heterogenous microarray data sets of different model organisms. The pipeline handles data from three microarray manufacturers (Agilent, Affymetrix, and Nimblegen) including single- and dual-channel arrays (Agilent) (Figure ). Functional analysis can be either performed with a list of DEGs or the FC rank based on FCROS. When analyzing data from different technologies, the use of a pipeline ensures unified and coherent statistical analysis and reliable output. No RNA-Seq data set fulfilled the selection criteria and was included in the work flow of the pipeline.

Principal Component Analysis of Gene Expression

Principal component analysis (PCA) is a commonly used tool to identify relationships and capture patterns in microarray datasets with multiple features.[44] We found high variability across all treatments for all three organisms and no clear separation of ENMs compared to the control treatments (Figure S1A–H). Separation was only found according to the experimental design, for example, when comparing different tissue types of A. thaliana, samples of leaves clearly separated from root samples (Figure S1A). This may be explained by the exposure path and experimental design of A. thaliana experiments, where normally roots are the first and major target resulting in different toxicological responses.[36] Further, separation by PCA was found when analyzing the single studies, for example, samples of silica exposed D. rerio larvae separated from respective control samples (Figure S1E).

Gene Expression Analysis of Selected Studies

Using the pipeline, data sets were processed and completely reanalyzed in order to identify statistically DEG for each treatment (contrast) as first output. When applying standard statistical threshold values (pFDR < 0.05, log FC > 1.2) and using a linear model (“limma”) in order to assess DEGs, out of 28 contrasts (Table ), in 14 contrasts no DEGs were found (Figure ). The highest number of DEGs was found for A. thaliana exposed to gold ENMs (1187) and titanium ENMs (779), C. elegans exposed to unfunctionalized silver ENM (271), and D. rerio exposed to PAMAM (499) or 50 nm-sized silver ENM (156). No DEGs were found for A. thaliana leaves exposed to cerium ENM, C. elegans treated with functionalized (PVP, MUA, and AUT) ENMs or aged silver, and in D. rerio exposed to silica and zinc.

Figure 4

Fraction of DEGs for A. thaliana, C. elegans, and D. rerio upon treatment with different nanomaterials. Differential gene expression was assessed using standard statistical threshold values (FDR < 0.05, abs(log FC) > 1.2).

For most data sets, our DEG numbers (pipeline) did not match the DEG numbers reported in the publications that originally analyzed the data sets (Table ; “DEG (reported)”). The reason is that less stringent statistics were used in those publications, in particular, in most studies no multiple sample correction (MSC) was performed.[36−38,40−42] This is normally done by implementing a false discovery rate (FDR) or adjusting the p-value. Only two studies accounted for multiple testing,[39,43] and the DEG numbers reported in the two are in the range of our pipeline results (Table ). For example, PAMAM G3 and G4 resulted in 230 and 220 DEGs[43] compared to the output of our pipeline with 35 and 499 DEGs, respectively. In one study the authors used both analysis with and without MSC;[37] however, only the results without MSC were considered for functional analysis. Omission of MSC can be justified if one is not worried about false positives;[45,46] however, the use of MSC is considered as standard in microarray analysis.[47] Further, the data of three studies are not published and no statistical thresholds were reported (GSE70509, GSE61186, and GSE50718).

Overall, the selection of appropriate statistical methods, including the choice of threshold values, is crucial to provide comparable output. The aim of this study was to analyze all the data sets in the same manner in order to better compare between the studies and synthesize the results. However, because the output in terms of DEGs was low in many studies, it was not possible to perform functional analysis with linear models. Therefore, we applied fold-change rank ordering statistics (FCROS),[27] which uses FC-based rank calculations instead of classical statistical testing. It was proposed to be more favorable for data sets with a high biological variability and also eludes the issue of MSC.[27]

In order to estimate the uncertainty in our data analysis pipeline, we repeated the analysis with varying inputs: we deleted 10% of all samples from the data set and repeated this process until all combinations of samples had been tested. The resulting confidence intervals (Figures S2 and S3) show that our results are robust with respect to small variations in the input.

Common Mechanisms of Nanotoxicity Found across Select Studies

Based on FCROS, functional enrichment analysis of gene ontologies (GO) was performed. The output in the form of GO terms was clustered into different functional categories (Table S3), revealing mechanistic functions specific to the organisms and ENM types (Figure ). The findings were visualized in radar charts, and by this, commonly expressed patterns across all treatments (Figure A), individual organisms (Figure B–D), and ENM types (Figures E,F and S4) were found. We identified several mechanisms which are commonly related to ENM toxicity but also found mechanisms with less-described responses.

Figure 5

Overlapping GO terms across organisms and ENMs were determined by calculation of the top ranked genes per contrast using FCROS followed by GO enrichment and classification of GO terms. Overlap of GO-term categories across organisms and ENMs are listed, whereby the color indicates the percentage of all observed GO-terms found in this category. Gray color indicates that no GO term was present for this category.

Figure 6

ENM toxicity response profiles illustrate the general response across all contrasts (A), the different organisms A. thaliana (B), C. elegans (C), D. rerio (D), and the ENMs silver (E) and silica (F). For this, the GO terms were clustered into different functional categories (Table S3). The score value is assessed by the number of GO terms which are present in each functional category and normalizing this number to the number of contrasts that were available for each organism or ENM group. This allows the comparison between species (B vs C vs D) or nanoparticles (E vs F). For example, for A. thaliana, 22 GO terms were related to “energy generation” resulting in a score value of 4.4 considering that five contrasts of A. thaliana were included. In comparison, there were eight different contrasts for C. elegans and 30 GO terms related to “energy generation” resulting in a score of 3.75. The general (average) response was assessed by normalizing the total number of GO terms of each category to the total number of contrasts.

The functional categories present across all treatments were the category of energy generation, general signaling, and DNA metabolism (Figures and 6A). Disturbed energy balance is most likely a consequence of mitochondrial-related ENM effects (e.g. ATP, NAD/NADP).[48−50] Once in the cell, intracellular ENMs can directly impact intracellular transport processes or damage cell organelles including mitochondria. As such, silver and titanium are known to physically impact mitochondria, for example, by changing the permeability of membranes.[51,52] For both, silver and titanium, we have found significant perturbations in the energy generation category (Figures E and S4). Often, mitochondrial damage is related to oxidative stress and the formation of intracellular ROS,[53] which can interfere with calcium uptake, resulting in structural damage of mitochondria.[54] In general, oxidative stress is often discussed as the prevalent mechanism of nanotoxicity;[55] the reactive surface characteristics of ENMs can promote the generation of intracellular reactive hydroxyl radicals.[56] Dramatic increase of these free radicals can result in lipid peroxidation, interference with proteins (e.g., posttranslational modifications), and DNA damage (e.g. histone binding).[8] In our analysis, oxidative stress was indicative by upregulation of GO terms referring to redox cell homeostasis, response to oxygen levels, and metabolic ROS processes. The presence of oxidative stress was found for all ENMs, except for PAMAM particles, but only a few GO terms relating to oxidative stress have been found in our analysis (Figure ), compared to a higher number for, for example, energy generation. This is partly because fewer GO terms are related to oxidative stress (Table S3) but also potentially caused by dynamics of oxidative stress response. It has been shown that the effects of oxidative stress can occur in a time-dependent manner, for example, ROS expression is followed by expression of p38 and p53 which results in DNA damage.[57] The magnitude of ROS formation at a cellular level varies between different ENMs types and its specific surface properties. In some studies, exposure to titanium dioxide even did not result in oxidative stress, although TiO2 exposure is commonly associated with it.[52]

Another commonly reported consequence of exposure to ENMs is genotoxicity.[58,59] We identified upregulation of several DNA-related processes, such as DNA metabolism, DNA repair, and DNA strand break repair mechanisms for all three organisms and all ENM types. Intracellular ENM exposure can lead to single- and double- DNA strand breaks.[8,60] This can be either be due to direct interaction, for example, interference with histone[61] or a consequence of ROS formation, which can also result in mutations.[62,63] Further, we found GO terms describing interference with transcriptional and translational processes and amino acid metabolism. The responses were present across all data sets (Figures and 6). This finding strengthens the paradigm of the cell nucleus as one of the main targets of ENMs, be it through direct or indirect interaction such as via dissolved metal ions.

Further, interference with membrane transport and cytoskeletal components were found in most treatments (Figures and S4). Plasma membranes are the first target of ENMs, which can disturb cell membrane function either by physical interaction or direct permeation.[64−67] Damage of cytoskeletal components and proteins was, for example, indicated by impaired actin filament organization and microtubule polymerization. Microtubules and actin are the major cytoskeletal constituents and are pivotal for mitotic processes. Interference of ENMs can lead to chromosomal aberrations such as polyploidy.[62,68,69] General signaling, which is indicative for impaired signal transduction, was also present for all ENMs, whereby only few GO terms were referred to this category (Figures and S4). These GO terms were inferred with general signaling pathways such as Wnt (Table S2), however, not with immune and inflammatory responses, a commonly reported effect of ENM exposure. ENMs are known to impact immune and inflammatory responses by, for example, affecting secretion of cytokines such as TNF-α, and IL-6.[70,71] However, most inflammation effects have been reported in studies associated with pulmonary exposure to ENMs,[72] which was not relevant for the data sets we used in this study.

In addition to the common mechanisms discussed above, our analysis revealed several responses and processes which are beyond the commonly described paradigms of nanotoxicity. Response to misfolded proteins and protein folding was found for silica, silver and zinc (Figure ). Structural damage of proteins can lead to adverse effects such as the bundling of actin.[73] However, it has been shown that ENMs can also have chaperone-like characteristics, and thus, can promote protein refolding.[74] In our analysis, chaperone mediated folding was one of the affected GO terms within the category of protein folding. It is unclear if ENMs will result in positive or adverse effects.

Further, we found perturbed neuronal activity in samples of C. elegans and D. rerio with the highest response for silica followed by silver (Figure C–F). Several GO terms related to neuron development, generation and differentiation were impaired as was synaptic transmission. Impairment of neurological functions by silica and silver has been shown before in both D. rerio(75,76) and C. elegans.[77−79]

Only few studies report on ENM-related effects on DNA and RNA methylation and epigenetics and many mechanisms are still unraveled.[80] We found methylation related effects in all organisms and only absent for silica ENMs (Figure ). Methylation is a prominent effect of ENM exposure that has been only reported before in ENM studies that specifically measured it. The detection of methylation GO terms across most studies suggests that the perturbation of cellular methylation is a common consequence of ENM exposure. Since dysfunctional methylation of DNA, RNA or histones can impact cellular functions and also lead to inheritable epigenetic changes,[81−83] we suggest it should be looked at in more detail in future studies.

When comparing our cross-species results with the results of the original studies (which are summarized in Table S4), we see that the main biological pathways implicated in nanotoxicity in those studies (e.g., RNA metabolism in GSE50718 or oxidative stress in GSE41333) have also been found in the meta-analysis, therefore at the level of the individual study we can say that our conclusions coincide well with the original studies. In contrast, since these studies focused mostly on their strongest results, they did not detect and/or discuss other biological pathways, such as DNA/RNA methylation, illustrating the value of the meta-analysis approach.

Applicability for Environmental Risk Assessment

The molecular and cellular mechanisms of ENM induced toxicity are complex, since one ENM often affects more than one target. However, standardized testing in order to compare or group different ENMs is still extremely difficult, because of the large variety of ENMs, biological species and experimental designs used in toxicological studies. Herein, we performed a meta-analysis of existing transcriptomics data sets in order to identify common mechanisms of nanotoxicity across extremely heterogenous studies. Using this approach, we have found known mechanisms of ENM induced toxicity which were described by previous, less quantitative reviews: oxidative stress, mitochondria- or DNA-related toxicity, and translational repression. However, we also found that DNA/RNA methylation is perturbed in most studies, which was not seen outside of specific DNA methylation studies. This demonstrates that such a meta-analysis can also be used to find less-described toxicity mechanisms and potentially even new ones. The package we developed and made available to the community can be used to perform similar meta-analyses in other fields of toxicology.

The toxicological profiles we present are based on a simple scoring method and visualization of the common mechanisms in radar charts. This can be useful for comparison between organism groups or ENM classes but also has the potential to provide input into ENM environmental risk assessment. Available risk assessment tools, such as the Swiss precautionary matrix, require simple input about the common toxicity mechanisms, which our study provides at a qualitative level (low, medium, and high).

Because the high-quality data sets that we were able to use in our meta-analysis are dominated by metal-based nanoparticles, it is difficult to assess how general our findings are for the whole nanotoxicity field. Because we only used nanotoxicity data obtained from three different species with heterogeneous experimental design, it is also difficult to assess how general the findings are across the tree of life. Our results would have been more robust if all the nanotoxicity gene expression studies undertaken thus far would have been annotated according to the field standards and openly shared in the public space. This is the responsibility of the whole scientific community; therefore, we here appeal to all its members, starting with the researchers who need to share the data but also funders, editors, and reviewers who need to demand that the data are shared before funding and the articles are published. It is only through a joined effort that we will be able to make use of the entirety of the information that toxicological science produces.