Biomarkers of nanomaterials hazard from multi-layer data

doi:10.21203/rs.3.rs-1555233/v1

PDF

Article

Biomarkers of nanomaterials hazard from multi-layer data

https://doi.org/10.21203/rs.3.rs-1555233/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 22 Apr, 2022

You are reading this latest preprint version

There is an urgent need to apply effective, data-driven approaches to reliably predict engineered nanomaterial (ENM) toxicity. Here we introduce a predictive computational framework based on the molecular and phenotypic effects of a large panel of ENMs across multiple in vitro and in vivo models. Our methodology allows for the grouping of ENMs based on multi-omics approaches combined with robust toxicity tests. Importantly, we identified novel mRNA-based toxicity markers and extensively replicated them in multiple independent datasets. We found that models based on combinations of omics-derived features and material intrinsic properties displayed significantly improved predictive accuracy as compared to physicochemical properties alone.

engineered nanomaterials

grouping

multi-omics

predictive toxicology

systems biology

Advances in molecular and cellular biology, along with technological improvements in assay development, are changing the paradigm of chemical safety assessment¹. These developments are also being increasingly adopted for the safety assessment of engineered nanomaterials (ENMs)^2,3. In particular, there is a shift from descriptive toxicology towards a mechanism-based predictive assessment of chemicals and ENMs. Computational approaches in which intrinsic properties of the toxicant are used as predictors of toxicity have been traditionally exploited in the context of Quantitative Structure-Activity Relationship (QSAR) modelling. However, these methods, although valid in principle, do not provide information concerning the biological mechanism of action of chemicals. An important change of paradigm in chemical safety assessment concerns the quantitative analysis of molecular (mechanistic) and functional (phenotypic) effects occurring at multiple levels of the biological organization. To this end, systems biology approaches are applied in order to discover biomarkers of exposure and uncover the mechanism of action of chemical substances⁴. Omics technologies, including transcriptomics and proteomics, are predicted to have a key role in future risk assessment strategies and these methodologies have been extensively used in recent years for the assessment of the biological responses of ENMs in vitro and in vivo^5,6. The combination of different omics-based methods provides a powerful way to evaluate the biological impact of low dose exposure to toxicants⁷, which are necessary to replicate real-life situations for humans and the environment.⁸ Hence, Toxicogenomics provides new insight on the chemical-biological interactions, which goes beyond the traditional structure-activity relationships. Recently, significant efforts have been initiated in order to explain intermediate mechanistic aspects of chemical exposure through the adoption of adverse outcome pathways (AOP) and the definition of biomarkers of the mechanism of action of toxicants. Such efforts are coordinated in Europe in the form of multiple EU-funded projects and initiatives coordinated by the Joint Research Center (JRC) of the European Commission, in the USA by the activities of the Environmental Protection Agency (EPA), as well as globally under the Environmental, Health and Safety (EHS) program of the Organisation for Economic Co-operation and Development (OECD) expert groups such as the Advisory Group on Molecular Screening and Toxicogenomics (EAGMST).

While biomarkers selected based on a priori knowledge have already been used as proxies of apical events in the context of large screening programs including recent US inter-agency collaborations^9,10, toxicogenomics can help to identify new mechanistic biomarkers in a hypothesis-free fashion^2,11.

One of the main challenges in nanosafety concerns the lack of understanding of the rules that govern the biological effects of different ENMs¹². Likewise, specific biomarkers of ENM exposure remain to be identified. Previous work demonstrated the utility of a multipronged assessment of ENMs, using a panel of cell types and cytotoxicity assays reflective of different endpoints, and has shown that it is possible to identify ENMs with similar patterns of biological activity across different cell types¹³. However, the latter study did not provide information on (predictive) biomarkers of ENM toxicity. During the EU-funded FP7-NANOSOLUTIONS project (http://nanosolutionsfp7.com), we implemented a multi-layered, omics-driven systems toxicology approach for ENM grouping and prioritization. The main objective was to derive a set of omics-based biomarkers that, in combination with the intrinsic properties of ENMs, could predict their hazard potential. To this end, data obtained using different in vitro and in vivo assays were homogenized and integrated in order to group the ENMs into three hazard classes (i.e., low or no hazard, intermediate hazard and high hazard). Then, feature selection algorithms and Random Forest (RF) classifiers were used to identify sets of composite (or synergistic) biomarkers able to distinguish the three hazard classes (Fig. 1). The models selected from different data layers, including omics-based data as well as data on intrinsic properties of ENM, were evaluated alone and in combination to test whether the integration of heterogeneous data can improve the predictions.

ENMs, apical toxicity and mechanistic data overview

The present study focuses on a comprehensive toxicity assessment of 31 industrially relevant ENMs. The set of ENMs covers eight different chemistries, spanning both metal and metal oxide nanomaterials (including semiconductor crystals) as well as carbonaceous materials, i.e., multi-walled carbon nanotubes and nanodiamonds (NDs) with three different surface modifications, i.e., amino/ammonium groups, carboxyl/carboxylate groups and poly(ethylene glycol) (PEG)-terminated surfaces, and unmodified ENMs (designated as “core” ENMs), along with variations in size (e.g., 5 nm versus 20 nm Au particles) and shape (i.e., spherical versus rod-shaped TiO₂ particles). We focused on industrially relevant ENMs and the role of surface modifications. Although it is impossible to cover all available ENMs in the context of a single study due to limited resources, those selected here are representative of a range of different physicochemical properties and hazard potentials for robust predictions. Detailed information on the synthesis can be found in ref. ¹⁴. Moreover, the Expanded Methods in the Suppl. Information file provides details on the 31 ENMs and the experimental protocols. A rigorous grid design was applied in order to allow direct evaluation of the effects of the chemistry and surface properties, as well as their combination. The multi-omics datasets included global expression levels of mRNA, miRNA and proteins from two human cell lines (monocyte-like THP-1 cells and bronchial epithelial BEAS-2B cells), mRNA expression from mouse lung tissues, protein corona profiles, and comprehensive characterization of all the ENMs (Suppl. Table S1). THP-1 cells are commonly used as a model for evaluating the cyto- and immunotoxicity of ENMs^15,16, while BEAS-2B cells are a preferred model for the in vitro assessment of pulmonary toxicity including the potential genotoxicity of ENMs^17,18. The multi-omics-based analysis was based on an equipotent, low-concentration exposure (EC₁₀). The EC₁₀ was identified by means of concentration-response studies carried in the THP-1 and BEAS-2B cell lines, as the dose (for each of the 31 EMNs) that elicited 10% cell death. At this subtoxic concentration, the mechanism of action measurements do not mirror the toxicity, since toxic phenotypes are not yet exacerbated. In fact, while the molecular alterations (e.g., at the transcriptomics level) can be intuitively used as a relatively direct measure of toxicity when the exposure doses are high enough to observe toxic phenotypes, they become the indication of important intermediate underlying mechanisms when observed at low, subtoxic, concentrations. Moreover, high doses that allow direct observation of acute toxicity in vitro and/or in vivo are often unable to mimic real life exposure scenarios. For these reasons, we decided to measure the molecular mechanism of action (transcriptomics, proteomics, etc.) in experimental conditions (EC₁₀) in which they would not bluntly mirror toxicity phenotypes, but at which they would instead allow the identification of early mechanistic biomarkers of toxicity, as previously described^19–21.

Grouping of ENMs into hazard categories

The data collected from the in vitro assays (cytotoxicity, genotoxicity and immunotoxicity) was homogenized by applying a point-based conversion system that resulted in a general toxicity score between 1 (no toxicity) and 6 (highest toxicity). Suppl. Tables S2-S5 describe the criteria to apply the point-based conversion system, while Tables S6-S8 illustrate some examples of the homogenization results. The Bayesian information criterion (BIC) was used to determine the number of clusters for each toxicity grouping (Suppl. Figure S1). The BIC was computed varying the number of clusters between 2 and 6. The analysis showed that 2, 3 and 4 are reasonable numbers of groups since they represent a good compromise between the BIC values and the model stability. Based on the result of the BIC analysis, three hazard categories were defined: 1) no-to-low (NoL), 2) medium (M), and 3) high (H) hazard. For the sake of clarity and brevity, only the results using these three hazard categories are shown here. However, results using 2 and 4 hazard categories are also reported in the Suppl. Figure S7.

We then defined three different classification tasks (Fig. 2) based on cytotoxicity (CYT), integration of the in vitro assays (INT), and the in vivo neutrophil infiltration (NEU). In the CYT task, the 31 ENMs were grouped into the hazard categories according to their ability to cause similar patterns of cytotoxicity, while the INT task considers the combination of the homogenized assays by the means of the multi-view clustering algorithm Similar Network Fusion (SNF)²². In parallel, we profiled the immunotoxicity of the 31 ENMs in C57BL/6 mice, based on the total cell counts as well as the number of macrophages, neutrophils, eosinophils, and lymphocytes identified in the bronchoalveolar lavage (BAL) fluid. Neutrophil counts in BAL fluid after 4 days exposure are a well-documented marker of acute inflammation in mice, as we and others have shown^23–26. The neutrophils (NEU) counts were then used to divide the ENMs in the hazard categories: NoL (count < 1), M (1 ≤ count < 10) and H (count ≥ 10).

Exploratory analysis based on PCA and univariate strategies

Next, we performed explorative analysis to assess the distributions of each data layer. For each dataset, the first two principal components were computed and used to plot the projected data points (Suppl. Figures S2-S4). Overall, this analysis highlights that different combinations of omics data layers and cell types convey distinct information about the mechanism of action of ENMs since the distributions of projected data points (ENMs) are remarkably different between each other. While the density of data points in the mRNA of BEAS-2B dataset is the most homogeneous, the protein corona and all the datasets collected from THP-1 are heavily influenced by a few samples. Altogether, the scatterplots highlight the absence of linear separability of the toxicity classes by using all the available molecular features. This is not a surprising phenomenon in omics datasets, in which it is expected that the majority of the total variance of the data is associated with more than two dimensions²⁷. Moreover, to define a baseline classification capability based on the sole predictive power of each independent variable in every dataset, the level of association between every variable in each dataset and the toxicity class of the ENMs according to the different labeling schemes (CYT, INT and NEU) was evaluated by fitting the univariate logistic regression model with repeated random splits (Fig. 3a). We found that single variables are only weak predictors of the ENMs toxicity.

Selecting multivariate biomarkers for ENM safety assessment

We further hypothesized that if any relationship between the intrinsic properties/molecular data and the toxicity grouping is present, it may be encoded in higher order interaction terms. Therefore, we applied multivariate modelling strategies with the objective of finding compact sets of synergistic biomarkers.

We used toxicity-based groupings of the 31 ENMs as target variables, and the omics datasets and the intrinsic properties of the ENMs as predictors. Alternative feature selection and classification methods were used to determine the best composite biomarkers for each classification task (INT, CYT and NEU). We applied standard machine learning (ML) techniques such as Logistic Regression with PCA (LR-PCA), LASSO (a regularization-based method) and varSelRF²⁸, and a recently described method, GARBO. GARBO is a specialized genetic algorithm that uses Fuzzy Logic and RF-based classifiers, to select biomarker sets that optimize the trade-off between classification accuracy and number of biomarkers²⁹. GARBO, LR-PCA and LASSO were applied to select alternative biomarker models from the large-scale omics datasets (mRNA, miRNA and proteomics). Compared to LR-PCA and LASSO methods, GARBO generated composite biomarker models achieving the highest prediction accuracy in most of the used data layers (see Suppl. Figures S5-S6). We systematically evaluated the classification performances of the selected biomarkers when using a different number of groups for each defined classification task (see Suppl. Figure S7). Our evaluation highlights that the selected biomarker models achieve classification accuracies similar to those obtained when defining groupings of three classes. Figure 3 reports the classification performances of univariate models (Fig. 3a) and multivariate-based models selected by using the GARBO algorithm (Fig. 3b).

In the present paradigm, mRNA-, proteomics-, and protein corona-based models reached high accuracy (Fig. 3b). On the other hand, the models selected from the physicochemical properties performed less accurately as compared to the other data layers. As expected, the best models for the prediction of mouse lung neutrophil infiltration (NEU) were those built from the mRNA in vivo data. However, in the BEAS-2B cell line, the mRNA-based models also satisfactorily predict in vivo immunotoxicity (NEU), while proteomics-derived biomarkers have limited predictability both in respect to cytotoxicity (CYT) and immunotoxicity (NEU). On the other hand, protein-based feature sets in THP-1 cells exhibited high accuracy and stability scores for the cytotoxicity (CYT) classification task.

We also investigated whether the integration of different omics data types and physicochemical properties could improve the classification performance of the ENM safety classifier (Suppl. Figure S8). We observed that ensemble classifiers derived from single or multiple data layers in THP-1 outperformed those derived from other biological models, when considering the cytotoxicity classification task (CYT). However, when considering the integrated classification of toxicity (INT), the ensemble classifiers derived from single or multiple data layers in the BEAS-2B cell model systematically achieved the highest classification performances (Suppl. Figure S8).

Biomarkers to predict neutrophil infiltration

Bronchoalveolar lavage (BAL) immune cell identification and counting is a commonly accepted non-invasive procedure for the accurate and confident diagnosis of specific lung diseases³⁰. Furthermore, neutrophil infiltration is a well-known marker of inflammation induced by ENM exposure. Here we asked whether the sets of specific biomarkers previously identified by means of linear regression would predict neutrophil infiltration. As shown in Suppl. Figures S9-S10, all ten models based on the in vivo mRNA data displayed satisfactory performances (R² > 0.6) with respect to neutrophil BAL counts in mice. Predicting in vivo endpoints from in vitro data to facilitate the implementation of the 3R principles in nanosafety is a relevant topic. To this end, we analysed in vitro transcriptomics data to identify biomarkers that could serve as predictors of neutrophil infiltration in vivo (Suppl. Figure S11, Suppl Table S12). Overall, all the models derived from BEAS2B have good predictive performances (R2 > 0.6) and outperform the models generated from THP1 data.

External validation of mRNA-based classifiers

Next, external transcriptome datasets of in vitro and in vivo exposures to ENMs were retrieved from the NCBI GEO database (Suppl. Table S9) in order to validate the top 10 mRNA-based biomarker sets selected for each exposure system (THP-1, BEAS-2B and mouse lung) and classification task (CYT, INT and NEU). Suppl. Tables S10-S12 report the gene sets representing the best mRNA-based biomarker models. Some of the ENMs in the selected external datasets correspond to the same class of ENMs as the ones included in the training set (for instance, TiO₂ and CNTs), while other types of materials were not represented in the training set (e.g., graphene oxide and crocidolite asbestos). For each biomarker set, an RF-based classifier was used to generate class probabilities. The scores from the top 10 RF-based models were averaged to yield one set of class probabilities for each test (NoL%, M%, H%), and the class associated with the highest score was chosen as the predicted class. Figures 4a-c report the level of toxicity predicted by each model (cell type/classification task) on each external dataset. In order to improve the legibility of the evaluation results, the predictions derived from different doses of the same exposure were averaged. The top mRNA models indicated a high hazard priority for MWCNTs (Mitsui-7, or MWCNT-mits7 in Fig. 4). MWCNT-7 (Mitsui-7), a nanomaterial known to cause damage to the lungs³¹ and classified as a potential human carcinogen by IARC³². Crocidolite asbestos, a known carcinogen, was also predicted as hazardous by the models trained on the toxicity classes derived from the integration of different toxicity endpoints. Anatase TiO₂ was predicted to be hazardous. Indeed, the anatase form of TiO₂ is known to be chemically more reactive leading to greater toxicity in vitro and in vivo as compared to the rutile form^33,34. TiO₂-nanobelts were also predicted as highly hazardous and this is in line with the findings of the corresponding original study³⁵, where the authors characterized patterns of gene expression in THP-1 cells and primary small airway epithelial cells exposed to high doses of TiO₂ nanobelts. Figures 4d-f indicate how close the predictions are made by different cell type/mRNA models. It is interesting to observe that THP-1-derived predictive models are closer to in vivo models than those derived using BEAS-2B when focusing on cytotoxicity. This supports the hypothesis that different in vitro models capture different aspects of the chemical exposures, hence collecting complementary data in multiple cell systems aids the in vitro-in vivo extrapolation of predictive biomarkers³⁶. This is an important conclusion that accords well with previous studies aimed at assessing the capacity of in vitro assays to predict relevant in vivo outcomes²³. Thus, it is unlikely that a single cell-based assay (focusing on a single endpoint) will accurately predict the more complex and concerted biological outcomes in vivo. Overall, the ENM safety classifier also yields robust and accurate results for external datasets and demonstrates the feasibility of a toxicogenomic based safety classification of ENMs. Moreover, this is the first study in which predictive models of nanotoxicity are validated in a large collection of manually curated public datasets. Our analysis shows that, despite the profound differences in experimental design, material characterization, and omics technologies used, published data can be of considerable practical utility when properly curated and made available to the community.

External validation of selected molecular markers

We observed that APOE, encoding apolipoprotein E, is consistently found in the THP-1 based models predictive of the cytotoxicity classification task (CYT) (Suppl. Table S10). To determine the validity of this molecular marker, we exposed THP-1 cells to a panel of amorphous SiO₂ ENMs shown to display varying degrees of cytotoxicity (Fig. 5a). We then performed RT-PCR and found that the upregulation of APOE correlated with cytotoxicity when cells were exposed to 10 different SiO₂ ENMs (Fig. 5b). This is thus fully in accordance with the CYT classifier for THP-1 cells (Suppl. Figure S12). However, SPNS2 was not validated as a biomarker of cytotoxicity.

Biological roles of the selected biomarkers

Focusing on the mRNA-based classifiers, several genes of interest were identified from the Suppl. Tables S10-S12. First, it is noteworthy that the BEAS-2B and THP-1 based models do not encompass the same genes nor do the in vitro models display similarities to the in vivo models at the level of the individual mRNAs. Second, for each model, it is noted that certain genes are more prevalent than others. For instance, LDLR, encoding the low-density lipoprotein receptor, is prominently featured in the case of the BEAS-2B based models selected for prediction of in vivo endpoints (NEU). Similarly, as noted above, APOE, encoding a lipid-binding protein, is consistently found in the THP-1 based models predictive of the cytotoxicity classification task (CYT) while CHIL3, encoding the chitinase-like 3 protein, is linked to the mouse lung-based models of in vivo classification. Despite the absence of endogenous chitin, a number of chitinases and chitinase-like proteins that bind but do not degrade chitin have been identified. These proteins play important roles in lung injury, and are also known to play key roles in Th2-dominated disorders such as asthma^37,38. In addition, a previous study using a mouse model of ovalbumin-induced asthma revealed that exposure to graphene oxide increased macrophage production of chitinases, CHI3L1 and AMCase³⁹.

As stated above, all the THP-1 models derived from single or multiple data layers outperformed other biological models for integrated (INT) classification tasks. Interestingly, the THP-1 based model featured three genes, namely AHRR, TMOD1, and NEAT1 (Suppl. Table S10, Suppl. Figure S13). AHRR encodes the aryl hydrocarbon receptor repressor (AhRR), which acts as a tumor suppressor gene in multiple human cancers.⁴⁰ Moreover, AhRR has a major impact on regulating inflammation.⁴¹ TMOD1 encodes the tropomodulin 1 protein (Tmod1), whose role in nanotoxicity is unexplored, but it could be linked to actin cytoskeleton-related responses to ENMs. In fact, SWCNTs were previously shown to reorganize cellular actin structures.⁴² NEAT1 (nuclear paraspeckle assembly transcript 1), in turn, is a long non-coding RNA (lncRNA) known to be upregulated in multiple malignancies.⁴³ Recently, NEAT1, a target gene of the tumour suppressor gene p53, was shown to enable tumorigenesis in vivo by promoting survival of oncogene-targeted cells.⁴⁴ Interestingly, extracellular vesicles enriched in lncRNAs, such as NEAT1, drive fibrosis in a mouse model of ischemic heart disease.⁴⁵ NEAT1 has also been suggested to drive the progression of liver fibrosis⁴⁶ and, more recently, its role in the promotion of pulmonary fibrosis has been shown^47,48. Indeed, fibrosis is a commonly observed adverse outcome related to ENM exposures, especially in the lung^49,50. This supports our hypotheses that integrated methods may be applied to identify early biomarkers predictive of long-term outcomes of exposure. Altogether, these results not only provide a means to establish innovative tools for the prediction of the toxicity of ENMs, but also clarify important aspects related to the intermediate mechanisms of the exposure. This information can be used in further dedicated studies to draft new adverse outcome pathways and refine existing ones.

The chemical industry is undergoing a profound reorientation of research and development towards “safe-and-sustainable-by-design” according to the principles of green chemistry. To fully unleash the power of this new paradigm, the nanotechnology industry needs to gain full access to integrated models that take into consideration both intrinsic and biological properties of the ENMs. These models go beyond the traditional predictive models in which the exposure is linked to an apical endpoint, by adding key elements of the mechanisms underpinning the biological responses. This is the first study in which such an approach has been used to analyze a large collection of industrially relevant ENMs with unprecedented depth. It should still be pointed out that these integrated approaches do not replace traditional toxicity testing or QSAR modelling. Instead, they provide a unique complementary view on the chemical-biological interactions. The work presented here is novel and ground-breaking for the present and future nanosafety for the following reasons:

i) it presents the largest in-depth characterization of intrinsic and biological properties for a selection of 31 industrially relevant nanomaterials. However, certain common ENMs were not included in the panel, due to limited resources, and for reasons having to do with the amenability of some ENMs to surface modification. Nonetheless, the dataset was designed to be as representative as possible, and the PCR validation of the identified biomarkers on an external set of silica ENMs (not included in the original panel) was found to correlate with cytotoxicity. This suggests that the applicability of the identified biomarkers robustly extends beyond the set of ENMs used for their discovery; ii) it showcases a novel integrated modelling approach to define hybrid predictive models of toxicity, comprising both intrinsic and mechanistic properties of ENMs. The best classification performances were obtained with a multiple classifier system integrating the model predictions of RF-based classifiers trained on biomarker sets selected from different omics data types. In our paradigm, predictive models relying exclusively on physicochemical properties of the ENMs achieved lower accuracy. However, their accuracy was improved by building hybrid models in which intrinsic ENM properties and omics-driven mechanism-of-action information was combined. Although carefully selected, the ENMs investigated here do not cover the full spectrum of nanomaterials presently available on the market. However, our validation results and the full availability and reusability of our data/source code in accordance with the FAIR principles⁵¹ will allow expanding and refining of our models as new data will become available; iii) it provides validation of the most accurate models in a large selection of manually curated toxicogenomics datasets as well as newly generated molecular data; iv) it highlights new biomarkers of toxicity that anchor the toxicity potential of ENMs to specific molecular and cellular functions, thus facilitating the generation and refinement of ENM-specific AOPs; v) it proves that profiling the molecular alterations of biological systems after ENM exposure at subtoxic doses, provides not only molecular proxies of toxicity, but novel and not-redundant knowledge on the mechanism of action of the exposure. In sum, the integrated models presented here predict the hazard potential of ENMs and may guide the prioritization of ENMs.

Synthesis and characterization of ENMs

Details on the synthesis and characterization of the 31 ENMs studied herein along with further details on the experimental protocols used to generate the in vitro and in vivo toxicity and omics data, and methods deployed for the validation of selected biomarkers, are found in the Supporting Information file (and refer to ref. ¹⁴ for ENMs).

Computational infrastructure

Data on ENM physicochemical properties, in vitro and in vivo toxicity, and different omics datasets, as described above, were derived from exposures to 31 different ENMs. This large set of data was exploited to identify biomarkers of ENM toxicity, to build predictive classifiers, and to validate the RF-based models on other ENM exposures obtained from external datasets. This was achieved through the implementation of a computational infrastructure consisting of the following steps: 1) identification of targeted classification tasks of toxicity, 2) selection of relevant intrinsic and biological properties enabling toxicity assessment of ENMs, 3) comparison of selected models with simpler computational methods (GARBO vs LASSO and LR-PCA), and 4) testing in vivo based mRNA biomarkers for the prediction of neutrophil count by using a regression-based model, 5) evaluation of ensemble models integrating more than one single-view classifier, and 6) validation by using external datasets (focusing on transcriptomics data due to their greater availability as compared to other omics data).

Grouping of the studied ENMs

Toxicity data generated from high-throughput screening methods were employed to identify shared patterns of toxicity across the selected 31 ENMs. We identified three different groupings of the 31 ENMs, one based on an integration of cyto-, geno- and immunotoxicity data using in vitro models (INT - in vitro), one based solely on cytotoxicity data (CYT - in vitro) and one based on neutrophil counts in BAL fluid (NEU - in vivo). Specifically, we used toxicity assays to measure DNA and chromosome damage in the BEAS-2B cell line (genotoxicity), cell viability (cytotoxicity) in multiple different cell types (BEAS-2B, Jurkat, THP-1), and cytokine profiles for THP-1 cells (immunotoxicity). We decided to conduct a separate classification analysis for the cytotoxicity data since they rely on different cell models and show more consistent toxicity profiles across the 31 ENMs. The in vitro cytotoxicity and genotoxicity assay data were homogenized by applying a point-based classification system. Given a toxicity endpoint, this system assigns a number, ranging between 1 (no toxicity) and 6 (high toxicity), indicating a general toxicity score. The categorization-based system along with the conversion results are reported in the Suppl. Information file. After homogenizing the toxicity assay data, i) the k-means algorithm implemented in R was employed to create groups of ENMs based on their cytotoxicity profiles, here indicated as CYT (more details in next section). Then, a multi-view clustering algorithm, namely Similar Network Fusion (SNF), was used to create a grouping of ENMs based on the integration of their cyto-, geno- and immunotoxicity profiles, here indicated as INT. Finally, we defined a third grouping of ENMs based on the BAL neutrophil levels, here indicated as “NEU” (for neutrophils). In particular, the neutrophil counts were exploited to define a third grouping of the ENMs which consists of the following three categories: NoL (value < 1), M (1 ≤ value < 10) and H (value ≥ 10).

Assessing the number of clusters for grouping

We assessed how data would group according to the toxicity assay data. To this end, for each toxicity score, we performed a cluster analysis by fitting a gaussian mixture model (GMM) with a spherical covariance matrix for each cluster. We compared the results varying the number of clusters, ranging from 2 to 6 using the BIC. The BIC can be interpreted as a goodness-of-fit score that penalizes overly complex models; the lower the score, the better the fit. For each number of clusters, we fitted the GMM 250 times to evaluate the influence of the random initialization on the algorithm and compared the respective distributions of the BICs. A larger variance of BICs indicates that the model depends heavily on the random initialization (this is evident for the case of 5 and 6 clusters, where the chance of overfitting is higher and represented by the many individual models with BIC values very distant from the average). On the other hand, when the number of clusters is kept between 2 and 4 the values of the BIC are compactly distributed around their respective mean values, implying that the clustering algorithm is not so influenced by the random initialization. Here, the results are reported dividing the ENMs into 3 groups. This value was chosen as a trade-off between data variability and stability of the results. For comparison, the same results are reported for different, equally reasonable, groupings k = 2, 4 in Suppl. Figure S1.

Training, validation and testing of RF-based classifiers

In order to address the feature and classification tasks, we partitioned the datasets into training (70%) and testing (30%). Because of the relatively limited number of studied ENMs, the training test was utilized for both model selection and training of RF-based classifiers on the final selected models. The implemented ML-based evaluation strategy defined training, validation, and the testing sets. The classification accuracy was systematically compiled on the test sets (30%), while the training sets (70%) were used to select the biomarker models (or feature subsets) and to train RF-based classifiers. The initial splitting in the training and testing set is repeated 5 times. Then, GARBO utilizes k-fold cross-validation for model evaluation and selection. Note that the univariate analysis is also based on the defined training and test sets.

Exploratory analysis based on PCA and univariate strategies

Prior to PCA analysis, each layer of omics data has been standardized in order for each variable to have an average value of 0 and a standard deviation of 1. After standardization, each dataset was projected on its corresponding first two principal components and plotted for graphical inspection. On each training set (70%), each variable was first standardized, then a univariate logistic regression was fitted and evaluated on the corresponding test set. For each variable in each dataset, the test set (30%) performances were collected at each split. The best-ranking features were chosen according to the best mean test set accuracies achieved and reported in Fig. 3a together with their corresponding 95% confidence interval as a measure of the spread of the estimate.

Marker selection for the prediction of toxicity

While toxicity assay data were used to define different categorizations of the ENMs, omics profiling and the intrinsic properties of ENMs were employed to identify predictive markers of ENM induced toxicity. Given one categorization of ENMs, different feature selection and classification algorithms were used. In particular, we tested two standard methods such as Logistic Regression combined with Principal Component Analysis (LR-PCA), and LASSO. LASSO is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy of the final trained models. Both models LR-PCA and LASSO, were trained using the defined 5-repeated train/test splits, with proportions 70%/30%. In both settings, each variable in the training set has been standardized to have mean 0 and standard deviation 1 (the learned standardization is equally applied to the corresponding test split). For the LR-PCA model, due to data availability, only the first two principal components were systematically used to fit the logistic regression model. For the LASSO model, a further model selection step on the training split has been performed to find the best regularization parameter. In particular, the training set of each split is further divided into a training set and a validation set according to a 5-fold cross-validation scheme. The nested cross-validation scheme was used to estimate the performance of different models indexed by the regularization parameter which varied in the range [\({2}^{-3}\), \({2}^{-1}\)]. The best parsimonious model was selected using the “one standard deviation” empirical rule⁵². After identification of the best parameter, a new model is again fitted on the whole training set this time and the generalization capabilities were evaluated on the corresponding test set. For both LR-PCA and LASSO models, the test set performances across the 5 splits were collected and aggregated into a mean value and a 95% confidence interval and reported in Suppl. Figures S5-S6. Then, we applied a novel feature selection algorithm, GARBO, which enables the selection of very small and highly accurate biomarker models from large-scale genomics data. More detailed information about GARBO can be found in ref.²⁹. LR-PCA, LASSO and GARBO were used to select multiple informative marker sets from the mRNA, miRNA and proteomics data layers. Backward variable elimination algorithm, namely varSelRF²⁸, was used to select feature sets from the intrinsic properties and protein corona profiles. Classification performances were calculated on the testing sets and by using the overall classification accuracy. We used the overall classification accuracy as the main metric for the evaluation of biomarker models, because of the limited number of samples in the testing sets (average of 6). Test accuracy and stability-based metrics were used to select the top 10 mRNA-based biomarker models for the external validation tests and the discussion on the most interesting biomarkers. The stability was calculated with the Dice-Sorensen’s index, which is always in the range of [0, 1]. Stability aims at measuring the capability of the feature selection process in reproducing (more or less) the same feature subsets with different training sample sets. High stability is often correlated to high reproducibility. The classification accuracy and the Dice score were finally used to rank the multiple marker sets (or markers of toxicity) generated by the feature selection and classification algorithms. Specifically, the rank was compiled as the weighted sum of the classification accuracy (w = 0.5) and the Dice stability score (w = 0.5).

Regression-based models for predicting BAL cell counts

The mRNA markers identified by means of the genetic algorithm on the in vivo and in vitro mRNA datasets for the in vivo classification tasks were used to build a linear regression model to predict the neutrophil BAL cell counts. The linear regression modelling was performed by means of the lm function from the stats R library. Before modelling, the neutrophil cell counts were log transformed.

External datasets focusing on mRNA-based features

The RF-based classifiers built upon the top 10 mRNA marker sets selected for each classification task were validated with mRNA profiles derived from external ENM related studies. The mRNA-based models were selected based on a weighted sum between test accuracy and the stability calculated with the Dice coefficient. Suppl. Table S9 reports Gene Expression Omnibus (GEO) ID and a brief description of the selected external (publicly available) ENM studies. The prediction of the class labels for each test was made upon the class probabilities (NoL%, M%, H%) generated by the top 10 RF-based models. The scores from the top 10 RF-based models were finally averaged to yield one set of class probabilities for each test, and the class associated with the highest score was chosen as predicted class label.

Acknowledgements

This work was supported by the European Commission through the Seventh Framework Programme (FP7-NANOSOLUTIONS; grant agreement no. 309329). D.G. was also supported by the Academy of Finland (grant agreement no. 322761) and the EU H2020 project NanoSolveIT (grant agreement no. 814572). The authors wish to thank the following members of the FP7-NANOSOLUTIONS consortium for the provision and characterization of ENMs: Jaimé Ruiz, Didier Astruc, Alexej Antipov, Yirij Fedutik, Carsten Jost, Alexey Kalachev, Alexandros Besinis, Guocheng Wang, Nicky Ehrlich, Zeljka Krpetic, Francesco Muraca, Alejandro Vilchez, Vicenç Pomar Portillo, Jose Luiz Muñoz, Julie Muller, Nathalie Luizie, Zahraa Al-ahmady, Cyrill Bussy, and Kostas Kostarelos. We also thank Michael Persson, Nouryon (Sweden), for the provision of colloidal silica particles, and Vesa Hongisto, Ivica Letunic, and Roberto Tagliaferri for assistance with data management during the course of the project. Finally, we wish to thank Troy Faithfull, Tampere University, for language editing.

Author contributions

D.G. conceived and supervised the project, interpreted the results, and wrote the paper. V.F. developed the algorithms, pre-processed the data, performed the analyses, interpreted the results, and wrote the paper. B.F. interpreted the results and wrote the paper. K.S. coordinated the FP7-NANOSOLUTIONS project with H.A. P.A.S.K., H.A. and D.G. conducted the in vivo studies and generated the transcriptomics (mRNA) data based on the in vivo and in vitro samples. M.F and A.S. performed the analyses, interpreted the results, and wrote the paper. L.A.S. collected, homogenized, preprocessed the public transcriptomics data, and wrote the paper. A.G., G.G. and G.V., R.G., H.N. and B.F. generated the in vitro toxicological data. M.C., K.L., S.M., E.H.L., S.V.C., R.H., F.K. characterized the ENMs and elaborated the dispersion protocols. O.R. and R.L. generated the miRNA data. J.Y. and R.Z. generated the proteomics data. M.M and K.D. generated the protein corona data. T.S. and J.K. coordinated the sample management for the omics data generation. P.R. and L.T. coordinated the data management in the FP7-NANOSOLUTIONS project. The final manuscript was approved by all the authors.

Competing interests

The authors declare no competing interests.

Data availability

The pre-processed datasets described in this paper, along with the scripts used to perform the described analyses, are available from the online repository Zenodo (10.5281/zenodo.4247173)⁵³. The in vitro (GSE148705) and in vivo (GSE157266) omics data have been deposited in the NCBI Gene Expression Omnibus (GEO) data repository.

Nel, A. E. & Malloy, T. F. Policy reforms to update chemical safety testing. Science 355, 1016–1018 (2017).
Fadeel, B. et al. Advanced tools for the safety assessment of nanomaterials. Nat. Nanotechnol. 13, 537–543 (2018).
Pietroiusti, A., Stockmann-Juvala, H., Lucaroni, F. & Savolainen, K. Nanomaterial exposure, toxicity, and impact on human health. Wiley Interdiscip. Rev. Nanomed. Nanobiotechnol. 10, e1513 (2018).
Riebeling, C., Jungnickel, H., Luch, A. & Haase, A. Systems biology to support nanomaterial grouping. Adv. Exp. Med. Biol. 947, 143–171 (2017).
Duan, J. et al. Quantitative Profiling of Protein S-Glutathionylation Reveals Redox-Dependent Regulation of Macrophage Function during Nanoparticle-Induced Oxidative Stress. ACS Nano 10, 524–538 (2016).
Kinaret, P. et al. Network analysis reveals similar transcriptomic responses to intrinsic properties of carbon nanomaterials in vitro and in vivo. ACS Nano 11, 3786–3796 (2017).
Pillai, S. et al. Linking toxicity and adaptive responses across the transcriptome, proteome, and phenotype of Chlamydomonas reinhardtii exposed to silver. Proc Natl Acad Sci USA 111, 3490–3495 (2014).
Krug, H. F. & Wick, P. Nanotoxicology: an interdisciplinary challenge. Angew Chem Int Ed Engl 50, 1260–1278 (2011).
Dix, D. J. et al. The ToxCast program for prioritizing toxicity testing of environmental chemicals. Toxicol. Sci. 95, 5–12 (2007).
Attene-Ramos, M. S. et al. The Tox21 robotic platform for the assessment of environmental chemicals–from vision to reality. Drug Discov. Today 18, 716–723 (2013).
Serra, A. et al. Transcriptomics in toxicogenomics, part III: data modelling for risk assessment. Nanomaterials (Basel) 10, (2020).
Valsami-Jones, E. & Lynch, I. NANOSAFETY. How safe are nanomaterials? Science 350, 388–389 (2015).
Shaw, S. Y. et al. Perturbational profiling of nanomaterial biologic activity. Proc Natl Acad Sci USA 105, 7387–7392 (2008).
Gallud, A. et al. Multiparametric profiling of engineered nanomaterials: unmasking the surface coating effect. Adv Sci (Weinh) 7, 2002221 (2020).
Lanone, S. et al. Comparative toxicity of 24 manufactured nanoparticles in human alveolar epithelial and macrophage cell lines. Part. Fibre Toxicol. 6, 14 (2009).
Bhattacharya, K., Kiliç, G., Costa, P. M. & Fadeel, B. Cytotoxicity screening and cytokine profiling of nineteen nanomaterials enables hazard ranking and grouping based on inflammogenic potential. Nanotoxicology 11, 809–826 (2017).
Manshian, B. B. et al. Single-walled carbon nanotubes: differential genotoxic potential associated with physico-chemical properties. Nanotoxicology 7, 144–156 (2013).
Catalán, J. et al. In vitro and in vivo genotoxic effects of straight versus tangled multi-walled carbon nanotubes. Nanotoxicology 10, 794–806 (2016).
Feliu, N. et al. Next-generation sequencing reveals low-dose effects of cationic dendrimers in primary human bronchial epithelial cells. ACS Nano 9, 146–163 (2015).
Mukherjee, S. P. et al. Next-Generation Sequencing Reveals Differential Responses to Acute versus Long-Term Exposures to Graphene Oxide in Human Lung Cells. Small 16, e1907686 (2020).
Scala, G. et al. Multi-omics analysis of ten carbon nanomaterials effects highlights cell type specific patterns of molecular regulation and adaptation. NanoImpact 11, 99–108 (2018).
Wang, B. et al. Similarity network fusion for aggregating data types on a genomic scale. Nat. Methods 11, 333–337 (2014).
Sayes, C. M., Reed, K. L. & Warheit, D. B. Assessing toxicity of fine and nanoparticles: comparing in vitro measurements to in vivo pulmonary toxicity profiles. Toxicol. Sci. 97, 163–180 (2007).
Duffin, R., Tran, L., Brown, D., Stone, V. & Donaldson, K. Proinflammogenic effects of low-toxicity and metal nanoparticles in vivo and in vitro: highlighting the role of particle surface area and surface reactivity. Inhal. Toxicol. 19, 849–856 (2007).
Shvedova, A. A. et al. Increased accumulation of neutrophils and decreased fibrosis in the lung of NADPH oxidase-deficient C57BL/6 mice exposed to carbon nanotubes. Toxicol. Appl. Pharmacol. 231, 235–240 (2008).
Rossi, E. M. et al. Airway exposure to silica-coated TiO2 nanoparticles induces pulmonary neutrophilia in mice. Toxicol. Sci. 113, 422–433 (2010).
Lenz, M., Müller, F.-J., Zenke, M. & Schuppert, A. Principal components analysis and the reported low intrinsic dimensionality of gene expression microarray data. Sci. Rep. 6, 25696 (2016).
Diaz-Uriarte, R. GeneSrF and varSelRF: a web-based tool and R package for gene selection and classification using random forest. BMC Bioinformatics 8, 328 (2007).
Fortino, V., Scala, G. & Greco, D. Feature set optimization in biomarker discovery from genome-scale data. Bioinformatics 36, 3393–3400 (2020).
Meyer, K. C. & Raghu, G. Bronchoalveolar lavage for the evaluation of interstitial lung disease: is it clinically useful? Eur. Respir. J. 38, 761–769 (2011).
Bornholdt, J. et al. Identification of gene transcription start sites and enhancers responding to pulmonary carbon nanotube exposure in vivo. ACS Nano 11, 3597–3613 (2017).
Kuempel, E. D. et al. Evaluating the mechanistic evidence and key data gaps in assessing the potential carcinogenicity of carbon nanotubes and nanofibers in humans. Crit. Rev. Toxicol. 47, 1–58 (2017).
Sayes, C. M. et al. Correlating nanoscale titania structure with toxicity: a cytotoxicity and inflammatory response study with human dermal fibroblasts and human lung epithelial cells. Toxicol. Sci. 92, 174–185 (2006).
Warheit, D. B., Webb, T. R., Reed, K. L., Frerichs, S. & Sayes, C. M. Pulmonary toxicity study in rats with three forms of ultrafine-TiO2 particles: differential responses related to surface properties. Toxicology 230, 90–104 (2007).
Tilton, S. C. et al. Three human cell types respond to multi-walled carbon nanotubes and titanium dioxide nanobelts with cell-specific transcriptomic and proteomic expression patterns. Nanotoxicology 8, 533–548 (2014).
Hartung, T. Perspectives on in vitro to in vivo extrapolations. Appl In Vitro Toxicol 4, 305–316 (2018).
Zhu, Z. et al. Acidic mammalian chitinase in asthmatic Th2 inflammation and IL-13 pathway activation. Science 304, 1678–1682 (2004).
Ober, C. et al. Effect of variation in CHI3L1 on serum YKL-40 level, risk of asthma, and lung function. N. Engl. J. Med. 358, 1682–1691 (2008).
Shurin, M. R. et al. Graphene oxide attenuates Th2-type immune responses, but augments airway remodeling and hyperresponsiveness in a murine model of asthma. ACS Nano 8, 5585–5599 (2014).
Zudaire, E. et al. The aryl hydrocarbon receptor repressor is a putative tumor suppressor gene in multiple human cancers. J. Clin. Invest. 118, 640–650 (2008).
Vogel, C. F. A. & Haarmann-Stemmann, T. The aryl hydrocarbon receptor repressor - More than a simple feedback inhibitor of AhR signaling: Clues for its role in inflammation and cancer. Curr. Opin. Toxicol. 2, 109–119 (2017).
Holt, B. D. et al. Carbon nanotubes reorganize actin structures in cells and ex vivo. ACS Nano 4, 4872–4878 (2010).
Dong, P. et al. Long Non-coding RNA NEAT1: A Novel Target for Diagnosis and Therapy in Human Tumors. Front. Genet. 9, 471 (2018).
Adriaens, C. et al. p53 induces formation of NEAT1 lncRNA-containing paraspeckles that modulate replication stress response and chemosensitivity. Nat. Med. 22, 861–868 (2016).
Kenneweg, F. et al. Long Noncoding RNA-Enriched Vesicles Secreted by Hypoxic Cardiomyocytes Drive Cardiac Fibrosis. Mol. Ther. Nucleic Acids 18, 363–374 (2019).
Yu, F., Jiang, Z., Chen, B., Dong, P. & Zheng, J. NEAT1 accelerates the progression of liver fibrosis via regulation of microRNA-122 and Kruppel-like factor 6. J. Mol. Med. 95, 1191–1202 (2017).
Liu, Y., Lu, F.-A., Wang, L., Wang, Y.-F. & Wu, C.-F. Long non–coding RNA NEAT1 promotes pulmonary fibrosis by regulating the microRNA–455–3p/SMAD3 axis. Mol. Med. Report. 23, (2021).
Zhang, Y., Yao, X. H., Wu, Y., Cao, G. K. & Han, D. LncRNA NEAT1 regulates pulmonary fibrosis through miR-9-5p and TGF-β signaling pathway. Eur. Rev. Med. Pharmacol. Sci. 24, 8483–8492 (2020).
Ryman-Rasmussen, J. P. et al. Inhaled carbon nanotubes reach the subpleural tissue in mice. Nat. Nanotechnol. 4, 747–751 (2009).
Nikota, J. et al. Meta-analysis of transcriptomic responses as a means to identify pulmonary disease outcomes for engineered nanomaterials. Part. Fibre Toxicol. 13, 25 (2016).
Wilkinson, M. D. et al. The FAIR Guiding Principles for scientific data management and stewardship. Sci. Data 3, 160018 (2016).
Franklin, J. The elements of statistical learning: data mining, inference and prediction. The Mathematical Intelligencer 27, 83–85 (2005).
Greco, D., Fortino, V., Serra, A. & Fratello, M. Data & code repository for the article “A predictive computational framework for the safety assessment of nanomaterials”. Zenodo (2020) doi:10.5281/zenodo.4247173.

There is NO Competing Interest.

NNANO20030636SupplementaryInformationFile24082021.docx
Supplementary methods and results

PDF

Version 1

posted 22 Apr, 2022

You are reading this latest preprint version

Biomarkers of nanomaterials hazard from multi-layer data

Status:

Version 1

Abstract

Figures

Introduction

Results

ENMs, apical toxicity and mechanistic data overview

Grouping of ENMs into hazard categories

Exploratory analysis based on PCA and univariate strategies

Selecting multivariate biomarkers for ENM safety assessment

Biomarkers to predict neutrophil infiltration

External validation of mRNA-based classifiers

External validation of selected molecular markers

Biological roles of the selected biomarkers

Conclusions

Methods

Synthesis and characterization of ENMs

Computational infrastructure

Grouping of the studied ENMs

Assessing the number of clusters for grouping

Training, validation and testing of RF-based classifiers

Exploratory analysis based on PCA and univariate strategies

Marker selection for the prediction of toxicity

Regression-based models for predicting BAL cell counts

External datasets focusing on mRNA-based features

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1