An ancestral molecular response to nanomaterial particulates

Knowledge assortment and pre-processing of nanomaterial datasets

The meta-analysis of ENMs toxicogenomics research in vitro and in vivo can determine frequent molecular MOA impartial from the organic system below analysis. To this finish, we carried out a meta-analysis of 66 transcriptomic datasets derived from the general public knowledge assortment curated by ref. 8 (, supplemented with knowledge beforehand revealed in ref. 10 (GSE157266) and ref. 9 (GSE148705) (Supplementary Desk 8). From the unique assortment, we excluded rats’ datasets and those based mostly on outdated microarray platforms, as they shared little probes with the newer variations (Supplementary Desk 8).

The datasets GSE148705 and GSE157266 had been pre-processed utilizing the eUTOPIA software program (model commit December 2021), as beforehand described8,65. Briefly, we filtered probes with a price larger than the 0.8 quantile in opposition to the damaging management in at the very least 75% of the samples. Knowledge had been normalized between arrays utilizing quantile normalization66. No batch correction was wanted for the dataset GSE148705, whereas GSE157266 was corrected for technical variation related to variables ‘dye’ and ‘slide’ utilizing the ComBat technique67. Lastly, we used the limma package deal (model 3.52.4) to compute the gene expression distinction between every publicity within the dataset and the corresponding controls, correcting the P worth utilizing the Benjamini–Hochberg process. The aggregated, normalized and corrected expression matrix was then exported with no extra filtering.

On this research, we chosen all pairwise comparisons between time–dose exposures and their respective controls. The ultimate dataset comprised 584 particular publicity situations (therapy, publicity time and dose, and organic system) and three,676 genes, ranging throughout varied human and mouse tissues and cell varieties.

The pool of three,676 genes represents the intersection of genes current in all of the experiments, limiting the mouse–human conversion to 1:1 orthology relationships (that’s, the place each genes within the pair have just one ortholog within the different species). The ortholog genes had been transformed utilizing the getLDS perform of the biomaRt R package deal (model 2.52.0)68. All knowledge can be found within the on-line Zenodo repository (

Assortment and pre-processing of transcriptomics knowledge for drug publicity and chemical compounds

To evaluate the specificity of the ENM signature, uncooked microarray knowledge for 158 medicine was downloaded from the Open Toxicogenomics Undertaking-Genomics Assisted Toxicity Analysis System database69. We utilized the identical pipeline to the in vivo exposures of rats to 3 dose ranges of every drug.

Uncooked knowledge had been imported into R utilizing the justRMA perform from the R library Affy (model 1.60.0)70. Probe annotation to Ensembl genes was carried out by the customized annotation recordsdata rat2302rnensgcdf (v. 22.0.0), downloaded from the mind array web site, resulting in 12,153 genes. The expression values had been quantile normalized by the use of the normalizeQuantile perform from the R limma library (model 3.52.4)66. Differential expression evaluation was carried out for every drug, for every mixture of dose degree and time level (1,839 comparisons). The analyses had been carried out evaluating the handled samples to the matched management samples of the identical time level. Because of this, log2 fold modifications, P values and adjusted P values (by the use of the false discovery fee (FDR) correction) had been retrieved for all genes for every comparability.

To make a comparable analysis, we solely chosen genes that had been included within the meta-analysis rank. As for chemical compounds, we downloaded, from the CTD, 142 gene signatures of ionic and covalent compounds71. We discarded gene units smaller than 15 and greater than 1,000, as it might have altered the statistics of the check, in addition to hydrocarbons, alcohols and ethers, as they’re labeled as natural. For every of the remaining gene units, we carried out a gene set enrichment evaluation (GSEA) of the ENM-associated rank and regarded a significance threshold of P worth adjusted to 0.05.

Assortment of ecotoxicological transcriptomics knowledge

To check the translatability of our mannequin to non-mammal species of eco-toxicological curiosity, we verified whether or not genes altered in response to ENM in different non-mammal species are regulated by C2H2–ZNF transcription components. To this purpose, we downloaded from the Gene Expression Omnibus (GEO) seven datasets masking 17 exposures to well-known eco-toxicological mannequin organisms (D. rerio, C. elegans, E. albidus and A. thaliana) and report the lists of differentially expressed genes (GSE80461, GSE32521, GSE70509, GSE73427, GSE77148, GSE41333, GSE47662)11,51 (Supplementary Desk 11). The genes have been used to carry out promoter evaluation as beforehand described for the invention assortment (evaluate with ‘Promoter evaluation’).

Characterization of the 584 experimental situations

From the preliminary assortment, we recognized 584 experimental situations which might be distinctive due to the organic system, ENM or experimental situation used. We manually annotated every experimental situation by curating the data within the unique publications (Supplementary Desk 3). First, we grouped samples in accordance with the organic system and publicity setting. We additionally included the publicity period by grouping samples into quick, intermediate and lengthy exposures. Totally different thresholds had been outlined for in vivo and in vitro. Intimately, for in vitro experiments, exposures have been thought-about quick at 24 h, intermediate between 24 and 72 h and lengthy after 72 h, respectively, as we anticipated the (sub-)acute toxicity to be noticed throughout the first days of the publicity. Certainly, ref. 72 just lately reported that in vitro techniques examined between 6 and 72 h reproduce a situation accounting for acute toxicity evaluation of chemical compounds.

For in vivo exposures, thresholds had been 3 days (quick), between 3 days and 1 month (intermediate) and greater than 1 month (long-term). Within the well being chemical analysis, the 28 day publicity Organisation for Financial Cooperation and Growth protocols are thought-about because the preliminary checks to evaluate long-term toxicity, whereas persistent toxicity research ought to have a size of 12 months73,74,75,76 (Supplementary Fig. 10).

As for doses, solely nominal doses for every experiment had been obtainable. It’s noteworthy that nominal doses usually are not comparable between the exposures, additional sophisticated by the heterogeneity of the measuring items reported (Supplementary Fig. 2b). We investigated that almost all of the research included right here examined sub-toxic dose exposures through a semi-automatic pipeline to scan the unique manuscripts. We retrieved all of the PubMed Central identifiers of the unique articles. When the experimental particulars weren’t obtainable, we thought-about the cited protocol in its references.

We used the BioPython Entrez (model 1.81)77 api to retrieve the paperwork by way of the PubMed Central IDs. Lastly, we parsed the ensuing XML paperwork and searched the summary and/or the entire article for key phrases (Supplementary Desk 8) referring to sub-toxic doses. Every optimistic consequence was returned with its context and manually validated. For a lowered variety of datasets, the corresponding articles couldn’t be retrieved routinely, so that they had been manually checked. We had been capable of finding indications of sub-toxic doses in 47 out of 58 publications.

The nanomaterials used within the experiment had been labeled in accordance with the chemical traits and the presence or absence of functionalized teams.

A panel of data was extracted from the unique publications (when doable), masking crystal part, purity, absence of endotoxins, coating, stabilizer and provider info, in addition to protocol info.

Lastly, when characterised, knowledge had been reported relating to the nominal diameter, size and particular floor space; transmission electron microscopy diameter, width and size; Brunauer-Emmett-Teller floor space; variety of partitions; dynamic mild scattering imply diameter and polydispersity index in water and medium; and zeta potential in water and medium.

Computation of molecular descriptors

A set of 159 ENM descriptors masking each molecular and digital construction properties was computed. Liquid drop mannequin molecular attributes78 are calculated assuming that ENM will be represented as a spherical drop, the place elementary molecules are tightly packed, whereas the density of clusters is the same as the particle mass density78,79. We computed the Wigner–Seitz radius (rw), the variety of ENMs within the analysed agglomerate (n), the variety of floor parts (S), the floor–quantity ratio (SV) and the aggregation parameter (AP)78. The Wigner–Seitz radius characterizes the minimal radius of interactions between particular person molecules and is represented by equation (1):

$$r_{mathrm{w}} = left( {frac{{3M}}{{4pi rho N_{mathrm{A}}}}} proper)^{frac{1}{3}}$$


the place M is the molecular weight, ρ is mass density, and NA is the Avogadro fixed.

The variety of ENMs within the agglomerate (n) is represented utilizing equation (2):

$$n = left( {frac{{r_0}}{r_{mathrm{w}}}} proper)^3$$


the place r0 is the radius of every ENM.

The variety of floor parts (S) is represented by equation (3):

$$S = 4n^{ – frac{1}{3}}$$


the place S describes the ratio of floor molecules to molecules within the quantity (or floor ENMs in agglomerates).

The floor–quantity ratio (SV) is represented utilizing equation (4):

$$SV = frac{S}{{1 – S}}$$


the place SV is the characteristic that describes the ratio of floor molecules to molecules in quantity (or floor ENMs in agglomerates).

Measurement-dependent interfacial thickness (h) was calculated with equation (5)

$$h = 0.01 occasions (T – 273) occasions r^{0.35}$$


the place r is the nominal dimension of the ENM and T is temperature80.

The ENM digital construction descriptors had been computed by density useful concept and semi-empirical quantum chemical strategies, whereas the Hamaker constants had been evaluated from atomistic drive fields and a continuum technique81,82. ENMs work together through long-range van der Waals interplay, which is a serious contribution to calculating the adsorption energies of biomolecules in water. Due to this fact, Hamaker constants are evaluated to explain bio–nano interactions in water by way of an atomistic drive subject method and through Lifschitz concept82. Within the Lifschitz concept82 two supplies are interacting by way of a medium; the Hamaker fixed for the ENM and a biomolecule in water is calculated from optical parameters which might be experimentally decided (Supplementary Desk 13), whereas within the drive subject method long-range dispersion interplay is calculated utilizing the Lorentz–Berthelot guidelines for sigma (atom dimension) and epsilon (atom–atom interplay amplitude)83,84. For steel ENMs, we used CHARMM drive subject parameters85. For steel oxides and carbon ENMs in addition to amino acids, lipids and sugars, the drive fields are reported in ref. 86. Contemplating all atom–atom interactions between two molecular entities, the Hamaker fixed is derived by an approximation of the mixed sigma and epsilon dispersion parameters87. On this work, we additionally thought-about the interplay between two ENM items in water. The geometric constructions of the majority ENMs had been optimized with density useful concept and the Perdew-Burke-Ernzerhof useful88 utilizing the SIESTA code89. The band gaps had been additionally calculated by PBE88, whereas the warmth of formation, electronegativity, absolute hardness, dispersion vitality per atom, dipole second and static polarizability descriptors81 had been obtained on the self-consistent subject degree by way of the semi-empirical code MOPAC (http://OpenMOPAC.web) utilizing the PM6-D390 parametrization. Lastly, ionization potentials, electron affinities, and the worldwide electrophilicity index had been computed by way of self-consistent cost calculations (ΔSCC calculation) for the digital states of the impartial and ion ENMs through the GFN1-xTB parameterization of the GFN-xTB code91,92,93,94,95.

We additional included a set of atomistic descriptors which might be based mostly on the chemical composition, potential vitality, lattice vitality, topology, dimension and drive vectors96,97. Constitutional descriptors are the counts of atoms of various identification and/or location. Potential vitality descriptors are derived from the force-field calculations, comparable to the arithmetic technique of the potential energies for particular atom varieties and/or places within the ENM. Lattice energies are based mostly on the identical potential energies however offered as per steel oxide nominal items (MxOy) and describe the vitality wanted to tear away stated unit from the ENM floor. The coordination variety of atoms is outlined because the depend of the neighbouring atoms that lie contained in the radius,

$$R = 1.2 occasions (R_{mathrm{M}};{mathrm{and}};R_{mathrm{O}})$$


the place RM and RO are the ionic radii of steel and oxygen ions, respectively. A low coordination quantity signifies that some atoms have lacking neighbours and thus makes the ENM extra unstable. The dimensions was derived from the precise calculated ENM diameter. The drive vector lengths have been derived from the construction optimization. For instance, to derive the common size (V) of the floor regular element of the drive vector for a shell area atom, its coordinates (x, y, z), drive vector elements (fx, fy, fz) and distance from the centre of mass (d) are used:

$$V = frac{{(xf_x + yf_y + zf_z)}}{d}$$


Pattern values for TiO2 (10 nm) are reported (Supplementary Desk 14). For amorphous ENMs a multi-step process was used requiring the simulation of bulk steel or steel oxide supplies above their melting temperature, the extraction of ENMs with the specified dimension and form, and their subsequent cooling on the temperature of curiosity with a prescribed fee. Such a process has been utilized to construct spherical amorphous ENMs with assistance from the automated Enalos Demokritos KNIME nodes. All the info are hosted on the NanoPharos database ( and had been transformed right into a ready-for-modelling format.

Meta-analysis implementation

We carried out a consensus of three algorithms for meta-analysis to prioritize the shared 3,676 genes.

As beforehand proposed98,99, our pipeline is predicated on the effect-size, P value-based and rank-product strategies. Often, meta-analysis frameworks are based mostly on effect-size strategies, assessing within- and between-study variations throughout a number of research. These strategies outperform others when there’s massive between-study variation and small pattern sizes. To implement it, the ‘effect_sizes’ perform from the esc R package deal (model 0.5.1) was used with the P-value argument and the ‘chi_esc’ perform100. The Fisher’s sum of logs technique combines particular person P values. Fisher’s sum of logs technique was carried out by utilizing the ‘sumlog’ perform of the R package deal metap (model 1.8), giving as enter the P values of every gene101. Lastly, the rank product is a non-parametric statistical technique to mix differential gene expression evaluation outcomes from particular person research based mostly on the within-study gene ranks. To this finish, the genes in every experiment had been ranked based mostly on the relevance of their related P values, and the ‘RP.advance’ perform of the RankProd R package deal (model 3.24.0) was used to merge them through one-class evaluation of the rank-product technique102,103. This perform permits combining knowledge coming from completely different research, resembling within the case of datasets generated by completely different laboratories. For every technique, a rank was generated. Lastly, all of the ranks had been mixed by way of the Borda perform of the TopKlists R package deal (model 1.0.8)104. The ultimate imply rank is reported in Supplementary Desk 1.

GSEA and have choice step

To pick out probably the most biologically related portion of the rank, we carried out a GSEA of the meta-analysis rank on 5 databases (Wikipathways105, Gene Ontology106, Reactome107, Kyoto Encyclopedia of Genes and Genomes108 and MsigDB109). In every case, the ‘fgsea’ perform from the fgsea R package deal (model 1.22.0) was used110. For every check, we recognized the place of the rank having the best peak worth of cumulative enrichment statistics. We created a lowered illustration of the meta-analysis gene rank by setting as a threshold the highest tenth percentile of such values (1,873 genes).

Computation of the frequency rating and hierarchical clustering

To seek out genes related to in vitro and in vivo long-term exposures for every gene of the lowered meta-analysis rank, we calculated a frequency rating as the proportion of samples during which the gene was statistically important.

The genes had been clustered in accordance with the Euclidean distance of their frequency scores. The hierarchical clustering algorithm, with Ward linkage technique, carried out into the ‘hclust’ perform of the R package deal ‘stats’ (model 4.2.0) was used.

For every kind of publicity system, we chosen the cluster with probably the most incessantly deregulated genes. To functionally annotate them, we carried out a pathway enrichment evaluation by way of the EnrichR on-line software (accessed in 2021), utilizing the MsigDB and Reactome databases107,109,111,112,113.

Computation of the molecular descriptors–gene expression correlation

To determine associations between ENM chemical properties and molecular alterations induced in cells and organisms by their publicity, the Pearson correlation coefficient was computed for every pair of gene and molecular descriptor, after a pre-processing step. Particularly, a Winsorize perform of the DescTools R package deal (model 0.99.43)114 was used to switch excessive values of log2 fold modifications with much less excessive ones. Furthermore, a dice root transformation was utilized to the molecular descriptor values.

Because the molecular descriptor knowledge layer accommodates lacking knowledge, the Pearson correlation was computed for the subset of samples the place values had been obtainable.

For every descriptor the highest 10% of probably the most correlated genes had been chosen. First, the gene units had been enriched in opposition to the Kyoto Encyclopedia of Genes and Genomes pathways by the use of the FunMappOne software (model commit December 2021)115. Solely pathways with FDR-corrected P < 0.05 had been thought-about considerably enriched. The molecular descriptors had been additional clustered in 9 teams based mostly on the Jaccard Index similarity of the shared enriched pathways. Lastly, the fgsea R package deal (model 1.22.0)110 was used to carry out a GSEA evaluation and determine the molecular descriptors whose set of related genes is enriched on the highest of the ranked listing of genes recognized with the meta-analysis method. Solely molecular descriptors with an adjusted P < 0.01 had been chosen.

Promoter evaluation

For every gene within the subset of curiosity, the sequence of the promoter area [−500 bp, +100 bp] across the TSS was downloaded utilizing the biomart package deal and the getSequence perform in ‘coding_gene_flank’ mode68. On this modality the perform returns the flanking area of the gene together with the untranslated areas.

Motif discovery was performed with the MEME software program suite (model 5.5.1)116. The motif website distribution was set as any variety of repetitions; the search was restricted to motifs ranging between 6 and 15 bases and the P-value threshold to 0.05.

For every consequence, the Factorbook database was interrogated to discover the TFBS ( Factorbook is a transcription factor-centric web-based repository related to ENCODE ChIP-seq knowledge, in addition to a number of databases of TFBSs. We chosen the TFBS that greatest matches the question in accordance with the software.

For every organism (Homo sapiens and Mus musculus within the discovery set, and D. rerio, C. elegans, E. albidus and A. thaliana within the eco-toxicological comparability) we annotated whether or not the TFBS could be acknowledged by a C2H2–ZNF member. To judge the statistical significance of C2H2–ZNF overrepresentation, we carried out a Fisher check with the fisher.check perform within the stats R package deal. We used as a background the set of non-redundant transcription issue binding profiles offered within the JASPAR database117. The contingency matrix was constructed by utilizing the set of TFBS of the C2H2–ZNF members of the family and all of the others, respectively.

Twin luciferase assay

BEAS-2B cells (ATCC, CRL-9606) had been grown in BEGM (Lonza, CC-3170). Cells had been cultivated in 75 cm2 tradition flasks at 37 °C with a humidified ambiance of 5% CO2. For all experiments, 500 µl of cells was seeded at a density of three.75 × 103 cells per ml in 48-well plates. Cells had been then left to relaxation in a single day earlier than transfections. One hour earlier than transfection, the media was changed with 225 µl of contemporary BEGM per effectively. Two vectors had been used for transfection: human cytomegalovirus (CMV) (optimistic vector management) and the ZNF362, created with VectorBuilder (Supplementary Fig. 7b). Per effectively, 0.25 µg of DNA vector in Opti MEM lowered serum medium (Gibco, 31985062) was added 1:1 with Lipofectamine 3000 reagent 6% V/V (ThermoFisher Scientific, 15338030) and P3000 enhancer reagent 4% V/V in Opti MEM lowered serum medium (Gibco, 31985062). About 25 µl of this transfection resolution was added per effectively and combined by light agitation of the plate. Twenty-four hours post-transfection, cells had been uncovered to one of many following: NM401 (JRC MWCNTs-NM401-JRCNM04001) or carbon black (CB, Orion Engineered Carbons, Printex 90) nanomaterials at both low (20 µg ml−1) or excessive (100 µg ml−1) focus; or nefazodone hydrochloride (Sigma-Aldrich, N5536) at low (25 µM) or excessive (50 µM) focus. NM401 and CB nanomaterials had been ready in accordance with the nanogenotox protocol ( Briefly, in a glass vial, 0.5% of ultimate inventory quantity of ethanol was added to the preliminary weighed nanomaterial powder; 0.05percentW/V BSA-BEGM was added for a ultimate inventory focus of 0.2 mg ml−1. The vial was then sonicated 2 occasions for 15 min in a water tub; this inventory resolution was then diluted in 0.05percentW/V BSA-BEGM to create ultimate options, and ultimate options had been sonicated for 15 min earlier than addition to the wells. Car management (VC) for nanomaterials was 0.05percentW/V BSA-BEGM. Nefazodone hydrochloride was ready in 0.05percentW/V BSA-BEGM and DMSO (ultimate DMSO focus in effectively of 0.5%). VC for nefazodone hydrochloride was 0.5% DMSO in 0.05percentW/V BSA-BEGM. Exposures had been carried out for twenty-four and 48 h.

The Twin-Glo luciferase Assay System (Promega, E2920) was used as per the producer’s pointers to measure firefly and renilla luciferase exercise, on a Spark multiplate reader (Tecan).

There have been three samples measured for every vector and automobile management, for every publicity. The imply sign of three background wells (cells solely) was used to subtract background from luciferase measurements. The firefly luciferase exercise was normalized to renilla luciferase and energy reworked. When current, outliers had been eliminated with the boxplot perform in R. The t-test was used to analyze the variations between the experimental and management samples.

Dose-dependent gene evaluation

To confirm if the C2H2–ZNF mannequin can clarify the dose-dependent portion of the ENM response, a dose–response evaluation of 62 research derived from 33 datasets (all initially included within the evaluation however GSE146708) was carried out following the technique carried out within the BMDx software (model commit February 2022)118. Briefly, a number of fashions had been fitted, and the optimum mannequin was chosen because the one with the bottom Akaike info criterion. The efficient doses (BMD, BMDL and BMDU) had been estimated below the idea of fixed variance. The benchmark response was recognized by the use of the usual deviation method with a benchmark response issue (BMRF) of 1.349, comparable to a minimal of 10% distinction with respect to the controls. Solely genes with lack-of-fit P > 0.01 and with estimated benchmark dose (BMD), benchmark dose decrease confidence restrict (BMDL) and benchmark dose higher confidence restrict (BMDU) values had been deemed related. Genes with BMD or BMDU values larger than the best publicity dose had been eliminated. Moreover, genes whose ratio between the expected doses is larger than the urged values (BMD/BMDL > 20, BMDU/BMD > 20 and BMDU/BMDL > 40) had been faraway from the evaluation.

AOP enrichment evaluation

To counterpoint key occasions and AOPs, we exploited the just lately curated annotation from ref. 48. For every gene set, we carried out an enrichment of the C2H2–ZNF targets as derived from ref. 119. For particular person occasions, the importance threshold was set at P worth adjusted to 0.05.

As for the entire AOPs, we first enriched as described the genes annotated to every pathway and discarded P values larger than 0.05. In a subsequent step, we evaluated if at the very least one-third of the occasions contained in it handed the identical threshold of significance (0.05).

Reporting abstract

Additional info on analysis design is out there within the Nature Portfolio Reporting Abstract linked to this text.

Related Articles


Please enter your comment!
Please enter your name here

Latest Articles