A transcriptomic overview of lung and liver changes one day after pulmonary exposure to graphene and graphene oxide

Hazard evaluation of graphene-based materials (GBM) is still in its early stage and it is slowed by their large diversity in the physicochemical properties. This study explores transcriptomic differences in the lung and liver after pulmonary exposure to two GBM with similar physical properties, but different surface chemistry. Female C57BL/6 mice were exposed by a single intratracheal instillation of 0, 18, 54 or 162 μg/mouse of graphene oxide (GO) or reduced graphene oxide (rGO). Pulmonary and hepatic changes in the transcriptome were profiled to identify commonly and uniquely perturbed functions and pathways by GO and rGO. These changes were then related to previously analyzed toxicity endpoints. GO exposure induced more differentially expressed genes, affected more functions, and perturbed more pathways compared to rGO, both in lung and liver tissues. The largest differences were observed for the pulmonary innate immune response and acute phase response, and for hepatic lipid homeostasis, which were strongly induced after GO exposure. These changes collective indicate a potential for atherosclerotic changes after GO, but not rGO exposure. As GO and rGO are physically similar, the higher level of hydroxyl groups on the surface of GO is likely the main reason for the observed differences. GO exposure also uniquely induced changes in the transcriptome related to fibrosis, whereas both GBM induced similar changes related to Reactive Oxygen Species production and genotoxicity. The differences in transcriptomic responses between the two GBM types can be used to understand how physicochemical properties influence biological responses and enable hazard evaluation of GBM and hazard ranking of GO and rGO, both in relation to each other and to other nanomaterials.


Introduction
Graphene-based materials (GBM) consist of a single or few layers of graphene sheets arranged in a 2-dimensional hexagonal lattice. Due to their nano-sized thicknesses and micron-sized lateral sizes, GBM fall under the high-aspect ratio nanomaterial category (Bianco et al., 2013;Wick et al., 2014). The simplest form, pristine graphene, is slow to produce (Smith et al., 2019) and can be challenging to work with it due heterogeneity in physicochemical properties across the GBM group makes hazard and risk assessment of GBM a cumbersome task (Park et al., 2017;Fadeel et al., 2018).
Graphene oxide (GO) and reduced graphene oxide (rGO) are two types of modified GBM that are industrially produced and used in a variety of applications (Dideikin and Vul, 2019). Created mainly through chemical oxidation of graphite flakes, the synthesis of GO is faster than that of pristine graphene. GO has a similar hexagonal carbon structure to graphene, but contains extensive oxidative modifications in its basal plane, such as hydroxyl, carbonyl, carboxylic acid and other oxygen-based functional groups (Smith et al., 2019). In addition to a higher oxygen content, resulting in C/O atomic ratios less than 3.0, the oxidation also generates defects by breaking up the sp 2 -hybridized structure of the stacked graphene sheets (Compton and Nguyen, 2010). The increased oxygen content reduces the electric conductivity of the material. This can be partly restored by chemical reduction, resulting in the formation of rGO. rGO differs from pristine graphene, as the procedures creating the material also introduce structural defects and residual oxygen.
Human exposure to GBM can either be intentional, e.g. injection for biomedical purposes, or unintentional, which can occur during manufacturing, handling, cleaning, packing, or disposal of GBM and GBM-containing products. For unintentional exposure, inhalation is the main route. A limited number of studies have investigated the toxic potential of GBM after pulmonary exposure (Duch et al., 2011;Schinwald et al., 2012;Li et al., 2013;Ma-Hock et al., 2013;Schinwald et al., 2014;Shurin et al., 2014;Han et al., 2015;Ma et al., 2015;Park et al., 2015;Shin et al., 2015;Wang et al., 2015;Kim et al., 2016;Roberts et al., 2016;Bengtson et al., 2017;Lee et al., 2017;Kim et al., 2018). These comprised of acute to sub-chronic inhalation and instillation/ aspiration studies, with large variety in the physicochemical properties of the GBM. The inhalation studies found minimal toxicity in male rats, despite exposing with different GBM, different doses, and study durations (highest dose: 10.1 mg/m 3 (6 h/day, 5 days)) (Ma-Hock et al., 2013;Han et al., 2015;Shin et al., 2015;Kim et al., 2016;Kim et al., 2018). In contrast, the instillation/aspiration studies, which primarily were conducted in mice, reported acute pulmonary inflammation in general and sub-chronic inflammation at the highest doses after exposure to both graphene and GO. However, differences in toxicity related to GBM size and types was also observed, such that GO appeared more toxic compared to graphene (Duch et al., 2011;Wang et al., 2015;Bengtson et al., 2017).
We have previously shown that commercially available micron-sized GO and rGO were not cytotoxic or genotoxic in lung epithelial cells in vitro but caused toxicity in vivo Bengtson et al., 2017). Mice exposed to GO and rGO by intratracheal instillation at occupationally relevant doses (18, 54 or 162 μg/mouse) displayed both short and long term toxicity, with respect to lung inflammation, acute phase response (APR) and genotoxicity (Bengtson et al., 2017). However, distinct differences between GO and rGO were observed, as GO generally was more inflammogenic than rGO, and strongly induced APR in lung and liver. In addition, increased APR protein levels were also observed in the blood of GO exposed mice. In contrast, rGO was found only to induce APR in lung tissue and the level was much lower than for GO. This indicates that the level of oxidation of these GBM affects the pulmonary toxicity following exposure. Both materials were present in the lungs of mice up to 90 days post exposure (Bengtson et al., 2017).
Detailed knowledge about mechanism of action in the lungs following pulmonary exposure to GBM has not yet been provided. We have previously investigated changes in the pulmonary and hepatic transcriptome after exposure to other carbonaceous nanomaterials such as carbon nanotubes (CNT) and nano-sized carbon black (CB) (Bourdon et al., 2012a;Poulsen et al., 2013;Poulsen et al., 2015a;Poulsen et al., 2015b;Halappanavar et al., 2019). As a continuation of these and our previous toxicological assessments of GO and rGO, the present study aims at mapping transcriptomic changes in the lung and liver after pulmonary exposure to GO and rGO in mice. This will allow us to directly compare mechanisms of toxicity that are distinctly associated with the two types of GBM. It will also enable comparison of transcriptomic changes found in the present study with those previously reported in the literature after exposure to other nanomaterials.

Materials
Tested materials in this study are two commercially available graphene materials, one graphene oxide (GO) and one reduced graphene oxide (rGO), manufactured and supplied by Graphenea (San Sebastian, Spain). GO were synthesized using liquid-phase exfoliation from graphite by a modified Hummers method. GO were then chemically reduced to form rGO using ascorbate. Details of the manufacturing has previously been published Bengtson et al., 2017).

Material dispersions
GO and rGO dispersion in the instillation media has previously been described in (Bengtson et al., 2017). Briefly, GO and rGO was added to 0.2 μm filtered, γ-irradiated Nano-pure Diamond UV water with 0.1% Tween80® (Sigma-Aldrich) to a final concentration of 3.24 mg/ml. The solution was sonication on ice for 16 min at 10% amplitude using a Branson Sonifier S-450D (Branson Ultrasonics Corp., Danbury, CT) equipped with disruptor horn. Dilutions (1.08 and 0.36 mg/ml) were further sonicated for 2 min. Vehicle control (VC) was prepared as above, but without GO and rGO. Intratracheal instillation in mice was conducted within 20 min following sonication.

Animal handling and exposure
Female C57BL/6J mice, 7 weeks upon arrival, were obtained from Taconic (Ry, Denmark) and were acclimatized for 1 week before the experiment. All mice were fed Altromin (no. 1324, Christian Petersen, Denmark) and had access to water ad libitum during the whole experiment. A total of 200 mice were used for the study.
All mice were randomly grouped according to graphene exposure, dose and day of euthanasia (n = 8 per dose and day for VC groups and n = 7 per dose and day for exposure groups), and they were housed in polypropylene cages with sawdust bedding and enrichment. Temperature and humidity was controlled at 21 ± 1 • C and 50 ± 10%, respectively, with a 12-h light and 12-h dark cycle. Daily observations of clinical signs of stress and discomfort were performed. The mean bodyweight at 8-weeks of age was 19.7 ± 1 g (Bengtson et al., 2017). The bodyweight was monitored during the study period and have been discussed previously (Bengtson et al., 2017). Mice were exposed to either VC, GO or rGO (18, 54 or 162 μg/mouse) by single intratracheal instillation (50 μl/mouse) under isoflurane sedation, as previously described (Jackson et al., 2010;Bengtson et al., 2017). VC groups included 8 mice per time point. GO and rGO exposure groups included 7 mice per dose and time point. Study design closely resembled that our other nanomaterial study designs in order to enable comparison between studies (Bourdon et al., 2012a;Husain et al., 2013;Poulsen et al., 2013;Poulsen et al., 2015b;Poulsen et al., 2016;Halappanavar et al., 2019;Knudsen et al., 2019).
Mice were euthanized at 1, 3, 28 or 90 days post-exposure by intraperitoneal administration of a 0.1 ml ZRF solution (Zoletil 250 mg, Rompun 20 mg/ml, Fentanyl 50 mg/ml in sterile isotone saline). Heart blood (800-1000 μl) was withdrawn via intracardiac puncture, stabilized with K 2 EDTA, and fractionated by centrifugation. The plasma was collected and stored at − 80 • C. Left lung was isolated and cut in 4 pieces (15-30 mg per piece). Collected samples were immediately snap-frozen in liquid nitrogen and then stored at − 80 • C until further analysis.
The study was in agreement with Directive 2010/63/EU of the

RNA isolation and quality
RNA from ~25 mg of lung tissue or 10-20 mg liver tissue was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and purified using RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) according to the manufactures protocol. RNA purity was measured using NanoDrop 2000c (ThermoScientific, USA). The mean RNA purity (A260/280 ratio) was determined to 2.07. Integrity of the isolated RNA was verified using Agilent Bioanalyzer 2100 with RNA Nano Chips (Agilent Technologies, USA) and conducted according to manufacturer's protocol. Total RNA was stored at − 80 • C until analysis.

Microarray hybridization
The following exposure groups were selected for the microarray analysis of lung samples: Day 1; all doses of both GO and rGO. Day 28; GO dose 18 μg and all doses of rGO. Day 90; GO dose 18 and 54 μg and all doses of rGO. The dose groups were selected as such, as the high dose of GO resulted in discomfort and weight loss in the mice. This was described in greater detail in a previous publication (Bengtson et al., 2017). VCs were included in each experimental group associated with time of exposure to account for the age associated changes. A total of 90 lung samples (five individual samples per experimental group) were analyzed using 12 individual microarray slides (8 samples per slide). For liver, day 1; 18 and 162 μg/mouse dose groups for both GO and rGO were analyzed (25 samples were analyzed using 4 microarray slides). Two hundred ng of total RNA from each tissue sample was used to perform microarray hybridization on Agilent Sureprint G3 Mouse GE 8x60K oligonucleotide arrays (Agilent Technologies, Mississauga, ON, CA) and scanned on an Agilent G2505B scanner at 5 μm resolution. Data were acquired using Agilent Feature Extraction software (v. 9.5.3.1). A detailed description of the DNA microarray analysis has been published previously (Poulsen et al., 2013).

Statistical analysis of microarray data
A reference randomized block design was created (Kerr, 2003;Kerr and Churchill, 2007), with the samples labeled with Cy5 and the Universal Reference labeled with Cy3, to analyze gene expression microarray data. LOcally WEighted Scatterplot Smoothing (LOWESS) (Cleveland, 1979) regression modeling was used to normalize data and statistical significance of the differentially expressed genes was determined using MicroArray ANalysis Of VAriance (MAANOVA) (Wu et al., 2003) in R statistical software (http://www.r-project.org). A shrinkage estimator for the gene-specific variance components was used to test the treatment effects (Fs statistic (Cui et al., 2005)). A permutation method (30,000 permutations with residual shuffling) was used to estimate the p-values for all the statistical tests. These p-values were then adjusted for multiple comparisons by using the false discovery rate multiple testing correction (Benjamini and Hochberg, 1995). Fold change calculations were based on the least-square means. Genes with a fold-change in expression of >1.5 in either direction (False Discovery Rate, FDR ≤ 0.05) when compared to VC were considered as differentially expressed genes (DEG) and were used in the downstream analysis. All microarray data have been deposited in the NCBI Gene Expression Omnibus database at NCBI (http://www.ncbi.nlm.nih.gov/geo/) accession number: GSE159707.

Analysis and interpretation of microarray data
Functional gene ontology analysis of the top 10 regulated genes from each exposure group were analyzed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 (Huang et al., 2009a;Huang et al., 2009b). Benjamini-Hochberg corrected gene ontology biological processes with a Fisher's exact p ≤ 0.05 were considered to be significantly enriched. In order to understand and interpret data from the microarray experiments, we used Ingenuity Pathway Analysis, (IPA) version 01-07 (Ingenuity systems, Redwood City, CA, USA) as previously described (Poulsen et al., 2015b). The reference in IPA was set to Agilent Sureprint G3 Mouse GE 8x60K Microarray and only direct relationships were considered. Unless specifically mentioned, only canonical pathways, diseases and disorders, upstream regulators and networks with Benjamini-Hockberg multiple testing corrected p-value (B -H p-value) <0.05 were considered for discussion. To infer the activation state of (predicted), canonical pathways and upstream regulators we incorporated the activation z-score from IPA in the analysis (z-score ≥ 2 or ≤ − 2 was considered as a predicted activation or inhibition, respectively). For the functional analyses in IPA, the individual enriched functions ("Molecular and Cellular Functions" and "Diseases and Disorders") were filtered by: 1) removing redundant functions with overlapping genes and annotations. This was done by importing the gene lists behind the function from IPA to excel and then compare genes. If all genes from one list could be found in a gene list from another function, then it was removed. 2) removing functions that were not directly relevant to the present study (e.g. dermal diseases and ophthalmic diseases).

Validation of DNA microarray results using qRT-PCR
From the microarray results, 12 differentially expressed genes were selected for validation by qRT-PCR. The genes (Il6,Timp1,Tgfβ1,Noxo1,Myd88,Cxcr5,Mt2,Apoa2,Ldlr,Tlr2,Saa1,Il1β) showed high magnitude of differential expression at least in one dose or time point group in the microarray experiment, and were included in the analysis because of their involvement in pathogen recognition, immune response, lipid metabolism, oxidative stress or fibrosis. Lung samples from 1 day exposure to 18 and 162 μg/mouse GO and rGO were considered in the analysis.
cDNA synthesis was prepared from isolated total RNA using TaqMan Reverse Transcription Reagent Kit (ThermoFischer Scientific, Denmark) according to manufacturer's protocol. The validation was conducted using custom-made TaqMan Array 96-well plates according to the manufacturer's instructions (Life Technologies, Denmark). Hprt, Actb, and Rn18s were used as reference genes for normalization. Following probe numbers were used: Rn18s (Hs99999901_s1), and Actb (Mm00607939_s1). Calculated fold changes of the 12 genes from the qRT-PCR analyses were compared to fold changes from the microarray analysis by linear regression and by a Pearson Correlation analysis performed in SAS version 9.4 (SAS Institute Inc., Cary, NC). Fold changes were log-transformed prior to the statistical analyses.

Hierarchical clustering
Hierarchical clustering is a statistical technique that groups experimental conditions based on their similarity or dissimilarity and measured by distances. The analysis was conducted on three data sets. First on the gene expression data from the current study. Then twice on historical gene expression data from previous publically available microarray studies on exposure to different nanomaterials such as multiwalled carbon nanotubes (MWCNT), CB, and nano-sized TiO 2 , combined with the present data set. One hierarchical clustering analysis used all time points, whereas the other only focused on post-exposure day 1, generating two clusters. These analyses were relevant to perform as procedures, doses and type of animals used were close to identical across the different studies. The following studies were included: GSE29042 (Guo et al., 2012), GSE35193 (Bourdon et al., 2012a), GSE41041 , GSE47000 (Poulsen et al., 2013), GSE60801 (Halappanavar et al., 2015), GSE61366 (Poulsen et al., 2015b). All of these studies are conducted on the same Agilent microarray platform. The software and gene annotation for this platform are continuously updated and regularly incorporated in the analysis. Data from all studies were merged using the Agilent probe ID, which is a unique identifier. Data were then collapsed to the gene symbol using the gene annotation from GPL10787-9758, downloaded from GEO. Each condition is represented by the log2 fold change. The fold change is a normalized value such that any batch related effects were subtracted out and thus controlled for. A 1-Pearson correlation dissimilarity metric was used for the hierarchical clustering analyses. The average linkage function, which determines how distances between sets of observations are calculated, was used and the analysis was conducted using the hclust() function in the R software package.

Plasma HDL, LDL/VLDL and total cholesterol
Plasma levels of high-density lipoprotein (HDL), low-density lipoprotein/very low-density lipoprotein (LDL/VLDL) and total cholesterol were determined with the EnzyChrom™ AF HDL and LDL/VLDL assay kit (EHDL-100, BioAssay Systems) according to the manufacturer's instructions. All time points and doses were evaluated. Within each dose and time point, animals were pooled two and two, such that n = 3 for each group. This choice was made to reach the sample volume needed for the kit. Remaining samples were not included in this analysis. After fractionation, plasma cholesterol, HDL, and LDL/VLDL isolations and a standard cholesterol reference supplied by the manufacturer were placed in 50 μl aliquots as duplicates in a 96-well plate. After addition of NAD-enzyme buffer mix and incubation, OD values were read at 340 nm on an Epoch microplate spectrophotometer (BioTek Instruments, USA). Concentrations were determined by comparison to the standard sample.

Results and discussion
In order to ensure the occupational safety of employees working with GBM, a deeper toxicological understanding of these materials is needed. Here, we investigated pulmonary and hepatic transcriptomic changes after intratracheal exposure to a type of GO and a type of rGO, with the aim of identifying key genes and pathways involved in the toxicity of these materials.

Materials
GO and rGO were chosen as they are produced at an industrial scale and as their dimensions are similar those of other GOs investigated in the literature (Park et al., 2015;Lee et al., 2017;Kim et al., 2018). A detailed physicochemical characterization of GO and rGO, both as raw materials and in relevant suspensions, has been published previously and is summarized in Table 1 Bengtson et al., 2017).
Briefly, GO and rGO were highly similar in their physical properties with comparable lateral sizes and number of layers (Table 1). This is probably a direct consequence of rGO being a reduced form of GO. Inorganic impurities were found to be present at relatively low levels (< 1.5%) for both GBM, but highest in rGO. Further, the level of endotoxin was found to be below the level required to induce pulmonary inflammation (Dong et al., 2009). The major difference in the physicochemical properties of GO and rGO lies with their chemical surface composition, with C/O and C/H ratios of 1.4 and 1.7 for GO and 8.5 and 13.2 for rGO (Table 1). This result in a much greater proportion of -COOH and -OH side groups on GO compared to rGO. As the two GBM resemble each other in all other physicochemical properties, differences in toxic potential could likely be attributed to their surface chemistry.

Microarray analysis
The present study assessed global transcriptomic changes in lung and liver using DNA microarray. All microarray data have been deposited in the NCBI Gene Expression Omnibus (GSE159707).

Number of differentially expressed genes
In general, pulmonary gene expression changes in mice exposed to GO and rGO varied across dose and post-exposure day. A descriptive analysis of the number of differentially expressed genes (DEG) showed a robust response at day 1, at which a total of 1363, 3302 and 2343 DEG were observed at 18, 54 and 162 μg of GO, respectively, whereas rGO exposure resulted in 805 and 860 DEG at doses 54 and 162 μg, respectively ( Fig. 1A). No DEG were observed after exposure to 18 μg of rGO on post-exposure day 1. In contrast, at day 28, only 1 DEG was observed in the 18 μg GO group, whereas 15, 2 and 0 DEG were observed for rGO dose 18, 54 and 162 μg, respectively. On post-exposure day 90, rGO solely induced 217 DEG at dose 54 μg. The time-dependent variations in DEG reflect previously observed variations in studied endpoints from the same mice, such as neutrophil influx, which were greatly increased on post-exposure day 1 and 3, but had returned to baseline level at day 28, and slightly increased on post-exposure day 90 (Bengtson et al., 2017).
Most DEG were upregulated across all exposure groups (Fig. 1A). In GO exposed samples, only ~8% of DEG were common across all the dose groups; however, ~35% of DEG were common to both middle and high dose groups (Fig. 1B). For rGO, the number of common DEG across the dose groups was ~58% (Fig. 1C). Surprisingly, a larger overlap in DEG (~24%) between the GO 18 μg dose group and the rGO 162 μg dose group was identified. In contrast, the high dose groups of GO and rGO shared only 5% DEG) (Fig. 1D-E). This indicates distinct differences in pulmonary responses after high dose GO and rGO exposure, as well as dose-dependency in GO-induced toxicity, since low dose GO exposure shared approximately the same number of DEG with GO dose 162 μg as it did with rGO dose 162 μg.
A similar analysis was conducted for DEG observed in the liver tissue 1 day post exposure. The GO 162 μg exposure group showed a larger number of DEG compared to the other exposure groups (Fig. 1A), with larger fold changes. Interestingly, in contrast to lung tissue, more DEG Lateral size, specific surface area, C/O and C/H ratios were determined using transmission electron microscopy (TEM), Brunauer-Emmet-Teller (BET) and combustion elemental analysis, respectively. Data were adapted with permission from Bengtson et al., 2017). a Determined at particle concentration 3.24 mg/ml in 0.1% TW80. were found to be downregulated in the liver (Fig. 1A). The most differentially expressed genes in the liver across all analyzed groups are presented in Supplementary Table 1.

Hierarchical cluster analysis
Hierarchical cluster analysis was conducted on lung samples from all time points to explore the grouping of the samples ( Supplementary  Fig. 1). Two separate clusters were observed: Cluster 1, which comprised of GO exposure at dose 54 and 162 μg, 1 day post exposure, and Cluster 2, which comprised the remaining samples ( Supplementary Fig. 1). This indicates distinct differences between GO exposure at dose 54 and 162 μg and the remaining exposure groups. Lowest dose of GO at day 1 was located in Cluster 2, albeit separated from all other groups (Supplementary Fig. 1). This is in agreement with the previous observation that the gene expression pattern after low dose GO exposure at day 1 resembled that of high dose rGO exposure ( Fig. 1B and D). The main part of Cluster 2 could be separated into two large clusters: An Effect cluster and a No-effect cluster. In the Effect cluster, exposure conditions leading to DEG were located, such as the rGO exposures at day 1 and the 54 μg rGO exposure at day 90. However, samples from exposure groups with few or no DEG were also present in this group, indicating low variance among these samples. The No-effect cluster contained the majority of VC samples and samples from exposure condition leading to few or no DEG. The last cluster was the 1 day VC cluster. The VC samples from day 1 clustered together, because the exposure method, intratracheal instillation, induces low-level inflammation (Jackson et al., 2010). Such lowlevel inflammation in the VC group was also detected in our previous publication using the same mice (Bengtson et al., 2017). Inflammation had returned to baseline levels 3 days post exposure.
This hierarchal cluster analysis supported the conclusions made based on number of DEG from the different exposure groups. Although many samples did not cluster distinctly with samples from their exposure groups, we consider this a natural consequence in a dataset including many samples with few or no DEG.

qRT-PCR
In order to validate results obtained through the microarray analysis, we conducted RT-qPCR analysis on 12 selected genes involved in pathogen recognition (Tlr2, Myd88), immune response (Il6, Saa1, Cxcr5, Tgfβ1), lipid metabolism (Ldlr, Apoa2) and oxidative stress (Timp1, Mt2). mRNA fold changes were determined 1 day after exposure to 18 and 162 μg GO and rGO, and were compared to the fold changes obtained in the microarray analysis ( Supplementary Fig. 2). In general, high consensus in fold changes was observed between qRT-PCR and microarray data. This trend was verified in both a strong significant correlation and linear regression (p < 0.0001), validating the microarray results.

BAL fluid composition and cytokine expression
In addition to the qRT-PCR validation, BAL fluid cell composition of the GO and rGO exposed mice (neutrophil levels at post-exposure day 1, lymphocyte levels at post-exposure day 3 and eosinophil levels at postexposure day 3) were compared by regression analysis to gene expression levels of several of their known chemoattractants 1 day after exposure to GO and rGO at all doses (Table 2). In addition, Il1β, Tnfα and IL6 were included for each analyzed cell type. BAL fluid cell composition was reported previously (Bengtson et al., 2017). For neutrophil cells in the BAL fluid, increasing Tnfα, Cxcl5, Saa3, and Ccl7 expression levels were identified as associated with increased neutrophil influx, which add additional validation to the microarray assay. In agreement with this, Saa3 and Ccl7 were both identified as correlated with neutrophil influx in a recent large comparison across several nanomaterials of very different physicochemical properties (Hadrup et al., 2020). Il1β and Ccl2 expression levels were identified as positively associated with lymphocyte cells in the BAL fluid, whereas Il16 was identified as negatively associated with lymphocyte cell numbers. For eosinophil levels, only Ccl24 was identified as associated with increased levels ( Table 2).

Lung functional analyses
Due to the high number of DEG 1 day post exposure to GO and rGO, compared to the other post-exposure days, further analyses were only conducted for this time point. The implications of the results obtained in these sections will be further discussed in the following sections.

Biological functions
The functional significance of the gene expression changes was determined using Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, www.ingenuity.com). Comparison analyses showed that GO and rGO perturbed or activated several "Diseases and disorders" (Fig. 2A) and "Molecular and cellular functions" (Fig. 2B). In general, GO and rGO exposure affected both common and unique functions. For commonly affected functions, the total number of DEG was consistently greater for GO-exposed groups.
The top regulated function in the category "Diseases and Disorders" was Organismal Injury and Abnormalities ( Fig. 2A). Annotations of this function after exposure to GO was related to inflammation and cancerous changes, which is a broad term related to genes often involved in the cell cycle processing and maintenance. For rGO exposure, this function was only annotated related to inflammation, not cancerous changes. In general, inflammation was greatly represented among top regulated "Diseases and Disorders", including functions: Immunological Disease, Inflammatory Response, Infectious Diseases and Hypersensitivity Response. Both GO and rGO exposure resulted in predicted activation of hypersensitive reaction in the Immunological Disease function, and activation of annotations related to mast cells and eosinophils in the  Although not immediately apparent, many of the top regulated "Molecular and cellular functions" were closely related to inflammation (Fig. 2B). As an example, the function Cellular Movement primarily consisted of annotations related to movement of inflammatory cells (both for GO and rGO exposure). Similarly, annotations related to the function Cellular Growth and Proliferation mainly involved proliferation or stimulation of inflammatory cells. Although GO and rGO exposures affected similar functions related to inflammation, some distinctions in annotation were observed. Whereas GO exposure primarily induced activation of cell viability and survival in the function Cell Death and Survival, rGO exposure also activated cell death and cytotoxicity. The only functions not directly involved in inflammation were Lipid Metabolism, Molecular Transport, and Free Radical Scavenging. Molecular Transport was affected by all exposure types, however, differences in annotations were observed. Whereas annotations related to GO exposure at dose 54 and 162 μg involved lipid metabolism within the function, low dose GO and rGO exposure primarily activated annotations related to mobilization and flux of Ca 2+ . The function Free Radical Scavenging was activated across all exposure types with the same annotations related to production, synthesis and metabolism of reactive oxygen species and production of superoxide. The function Lipid Metabolism was almost exclusively regulated by GO exposure, with activation of synthesis of lipids for all GO doses.

Pathway analysis
To further broaden our understanding of the potential biological effects of the identified DEG, an analysis of perturbed pathways was also conducted in IPA. Interestingly, although the total number of DEG for GO-exposed pathways were consistently higher, their significance level were in some cases found to be lower than rGO-affected pathways, indicating a greater clustering of genes in the perturbed pathways of the rGO-affected functions. An overview of significantly regulated pathways and their molecular functions are depicted in Fig. 3. Large differences were observed between the different exposure groups, especially between GO exposure at dose 162 μg and the rGO exposures. As observed in previous analyses, GO exposure at dose 18 μg perturbed pathways in a similar manner to both GO high dose exposure and to rGO exposure (Fig. 3). All perturbed pathways are listed in Supplementary Table 2. The pathways with the highest significance levels one day after exposure to GO and rGO are shown in Table 3    in the lung across all groups related to pathogen recognition, immune response, cellular signaling and cellular growth. The heatmap was generated to visualize the perturbed pathways, organized by classification categories. Identified pathways and categories were derived from the Ingenuity Pathway Analysis. cholesterol homeostasis (Acute Phase Response Signaling, LXR/RXR Activation, and Superpathway of Cholesterol Biosynthesis) (Table 3). Although IL-10 is considered an anti-inflammatory cytokine, analysis of the IL-10 Signaling pathway revealed that the DEG resulted in a strong predicted inhibition of IL-10. Both Granulocyte and Agranulocyte Adhesion and Diapedesis were perturbed by GO exposure, however as the pathways did not distinguish between different immune cells within the two categories, the analyses could not differentiate between the innate and the adaptive immune response.
Top regulated pathways after rGO exposure were primarily involved in inflammation (Table 3) and overlapped quite heavily across the two doses. The only non-inflammatory top-regulated pathway was Tec Kinase Signaling, which is involved gene expression, apoptosis and Ca 2+ mobilization. Although it appears as though rGO exposure resulted in uniquely regulated pathways, a more detailed analysis revealed that for the majority of these pathways the same genes were differentially expressed across both GO and rGO datasets and the predicted activation based on this was almost identical.

Additional analyses
In addition to the functional and pathway analyses, we also conducted analyses of the 10 most regulated (up and down) genes in each exposure and upstream analyses. The results from these analyses were in line with the functional and pathway analyses and can be found in the supplementary material ( Supplementary Figs. 3-4

Pulmonary inflammatory response
As described in the overall analyses of perturbed functions, perturbed pathways, and upstream regulators, changes related to inflammation and inflammatory processes were the most common observations in the lung. This was true for both GO and rGO exposure, with many overlapping functions and pathways as described above. The inflammogenic response was stronger for GO compared to rGO, both in terms of the number of DEG and fold change. For example, a total of 26 cytokines and chemokines were differentially regulated after exposure to GO, compared to 11 cytokines after rGO exposure (Supplementary  Table 5). Cytokines Ccl2 and Il6, which were differentially expressed solely after GO exposure, have recurrently been reported upregulated in the literature after both inhalation and instillation exposures to GBM (Duch et al., 2011;Ma-Hock et al., 2013;Park et al., 2015;Roberts et al., 2016). The main differences between the inflammatory responses elicited by GO and rGO exposure were related to the innate immune response and activation of lymphocytes. GO exposure affected the innate immune system to a much larger degree than rGO exposure, with activation of key innate pathways such as IL-6 Signaling, IL-8 Signaling, HMGB1 signaling pathways, and LXR/RXR Activation (Supplementary  Table 2). This links GO exposure with neutrophil influx as IL-6 Signaling and LXR/RXR Activation are involved in neutrophil recruitment (Zelcer and Tontonoz, 2006;Hunter and Jones, 2015). In addition, the GOexposed groups uniquely increased expression levels of a collection of cytokines (Ccl2, Ccl4, Cxcl2, Cxcl6) linked to phagocyte interaction during onset or resolution of inflammation (Supplementary Table 5) (Soehnlein and Lindbom, 2010).
In contrast to the strong association between GO exposure and the innate immune response, rGO exposure uniquely upregulated expression of Ccl5. This cytokine is linked to decreased neutrophil infiltration (Ariel et al., 2006). In addition, rGO exposure uniquely activated several pathways involved in activation and maturation of the adaptive immune system, such as CD28 Signaling in T Helper Cells and FcγRIIB Signaling in B Lymphocytes (Supplementary Table 2). This indicates that the inflammatory response was further in the progress towards resolution 1 day after exposure to rGO compare to GO exposure. It should be noted, that the observed increased maturation and activation of lymphocytes after rGO exposure was not reflected in the BAL fluid cellular composition (Bengtson et al., 2017). However, this analysis did not distinguish between different types of lymphocytes, which may hide the effect seen at mRNA level. Similar difference in resolution progress between GO and rGO was observed in our previous study, with neutrophilic levels in BAL fluid already at or close to baseline levels 3 days after exposure to rGO, whereas GO exposure resulted in greatly increased neutrophilic influx at the same time point (Bengtson et al., 2017). As in vivo comparisons between GO and rGO (or graphene) are very scares in the literature (Duch et al., 2011;Wang et al., 2015), this difference in resolution progress has not been observed previously. However, other observations challenge this hypothesis. Although rGO uniquely regulated pathways involved in the adaptive immune system, a detailed look at the genes involved in these pathways revealed that they generally were similarly regulated after rGO and GO exposure. However, this was not detected by the bioinformatics tools used, due to the much larger innate response quenching the smaller adaptive response. This would suggest that GO and rGO induced similar baseline responses and that GO induces a much larger innate inflammatory response on top of this. As GO 18 μg exposure show similarities with both high dose GO and high dose rGO responses across all analyses, it suggests that the lower dose of GO induces an intermediate response that contains a strong inflammatory response as the higher GO doses, but which also appears closer to resolution of inflammation similar to the response seen after rGO exposure. This highlights how transcriptomic analyses can provide important insight regarding dose dynamics, which could be used in future hazard assessment for GBM. The large observed differences in inflammogenic potential could probably be attributed to the surface chemical composition of the GBM, namely a much greater proportion of -COOH and -OH side groups on GO compared to rGO, as they have quite similar dimensions (Table 1). A previous study in mice exposed to nanocellulose showed a hydroxyl group-dependent induction of inflammation and the APR (Hadrup et al., 2019). These effects were significantly reduced when -OH groups were blocked by carboxylation. It is therefore possible that the large increase in inflammation seen in the present study after GO exposure can be attributed to the presence of hydroxyl groups. In addition, we observed an upregulation of genes related to toll-like receptor 2 (Tlr2) and toll-like receptor 4 (Tlr4) pathway signaling (Lbp,Cd14,Tirap,Tram,Myd88,Mapk and Tbk1) (Tsukamoto et al., 2018) after exposure to GO, whereas only two of these genes were upregulated for rGO (Tlr1 and Tlr2). This suggests that the abundance of -OH groups on the GO surface could result in GO being recognized as bacteria or LPS, due to structural similarities between GO and the lipid A part of LPS and with lipoteichoic acid. Due to its different surface composition, rGO exposure does not seem to trigger the same immunological mechanisms as GO, thus promoting a less aggressive immune response.

Cholesterol metabolism and the acute phase response in the lung
One of the main differences between GO and rGO exposure was their activation of functions and pathways related to lipid metabolism, cholesterol homeostasis and the APR. In "Diseases and Disorders", the function Cardiovascular disease was perturbed solely after GO exposure ( Fig. 2A), with Atherosclerosis predicted as activated for all doses. In similar fashion, Lipid metabolism was almost exclusively perturbed after exposure to GO (93, 361, 302 DEG for GO compared to 18 DEG for rGO) (Fig. 2B). Synthesis of lipids was activated for all GO doses, as well as Concentration of lipid, Release of lipid, and Release of fatty acid for GO exposure at 18 μg. Several pathways related to cholesterol homeostasis were also perturbed after GO exposure only, namely Acute Phase Response Signaling, LXR/RXR Activation, and Superpathway of Cholesterol Biosynthesis (Table 3). Acute Phase Response Signaling was among the top regulated pathways for all GO exposures, and although some DEG related to Acute Phase Response Signaling were also observed for rGO exposure ( Supplementary Fig. 5), the numbers (9-14) were much lower than the numbers observed for the GO-exposed groups (31-51). The large difference in responses after GO and rGO exposure in regards to cholesterol metabolism and the APR may have severe implications, as these pathways have previously been linked to increased risk of developing cardiovascular disease (Dong et al., 2011;Saber et al., 2013;Saber et al., 2014;Poulsen et al., 2015a;Thompson et al., 2018;Vogel and Cassee, 2018;Hadrup et al., 2020).
The atherogenic potential of GO is strongly linked to the induction of the APR and specifically to the APR genes of the serum amyloid A family (Saa). In the present pulmonary dataset, Saa3 was the overall most differentially regulated gene, with isoforms Saa1 and Saa2 among the top regulated genes (Supplementary Table 3). However, they were almost exclusively differentially expressed after GO exposure. This is in concordance with the qRT-PCR results from our previous study, which revealed that increased Saa3 expression was also present at postexposure day 3, but not at day 28 and 90 (Bengtson et al., 2017). In addition, increased levels of SAA3 protein were detected in the plasma of GO exposed mice 3 days after exposure as compared to the vehicle exposed. Such increase was not seen following rGO exposure. Studies in APOE − /− mice have shown that overexpression of SAA leads to increased plaque progression and inhibition of SAA synthesis leads to lowered plaque progression (Dong et al., 2011;Thompson et al., 2018).
Increased Saa expression may lead to atherogenic changes, as SAA is incorporated into HDL lipoproteins during an APR (Feingold and Grunfeld, 2016). This results in alteration to the cholesterol flow, leading to accumulation of cholesterol in peripheral tissue and macrophages, thereby interrupting reverse cholesterol transport (Ono, 2012). In support of this, a significant downregulation of Abca1 and Abcg1 was observed in the LXR/RXR Activation pathway after GO exposure. These genes are involved in cholesterol efflux out of cells. In addition, analysis of this pathway also revealed that the LDL receptor (Ldlr) was upregulated after exposure to GO at doses 54 and 162 μg. This receptor mediates the endocytosis of LDL particles from the bloodstream into the cells. Abca1, Abcg1 and Ldlr are all important players in the reverse cholesterol transport, and their expression suggest that cholesterol is maintained in the lung. In addition, SAA and HDL-SAA have the ability to turn macrophages foam cells (Lee et al., 2013). These combined effects may lead to plaque progression and atherosclerosis. The transcriptomic analysis of the present paper therefore indicates that pulmonary exposure to GO, but not rGO, could promote SAA-mediated development of atherosclerosis.
Interestingly, the lower exposure doses induced greater Saa3 levels compared to the higher doses, similar to that seen in the qRT-PCR analysis presented in (Bengtson et al., 2017). Although Saa3 has previously been identified as the most differentially regulated gene in the lungs after pulmonary exposure other carbon-based nanomaterials, the same trend of the greatest expression levels seen at the lowest doses has not previously been observed (Bourdon et al., 2012a;Poulsen et al., 2015b).

Free radical scavenging in the lung
GO and rGO exposure (all doses) induced activation of the function Free Radical Scavenging (Fig. 2B). Annotations of this function were to a larger extent predicted as activated after rGO exposure compared to GO exposure, although the annotations contained the most DEG in the GOexposed groups. Activation of the annotation "production of superoxide" was found for both GO and rGO, with the highest number of DEG observed after GO exposure, and superoxide was also identified as an activated upstream regulator for GO. Indeed, we have previously shown that GO and rGO generate reactive oxygen species (ROS) in an acellular assay, similar to other carbon based nanomaterials, such as CB, CNT and diesel exhaust particles (Jacobsen et al., 2008a;Jacobsen et al., 2008b;Poulsen et al., 2015b;Bengtson et al., 2016). In concordance with this, we also observed increased genotoxicity, in terms of increased levels of DNA strand breaks by the Comet assay, in the lungs of mice following pulmonary exposure to GO and rGO (Bengtson et al., 2017), as well as CB and diesel exhaust particles (Bourdon et al., 2012b;Kyjovska et al., 2015b;Kyjovska et al., 2015a). However, despite GO being more potent in inducing acellular ROS, no difference was observed in their ability to induce in vivo genotoxicity. Similar observations were made at the mRNA level in the present study. The lack of difference between GO and rGO regarding the genotoxic potential and transcriptomic changes related to ROS-mediated toxicity, suggests that the observed pulmonary genotoxicity is not related to the surface composition of the GBM and the inflammatory response. Similar observations was seen for MWCNT (Poulsen et al., 2016). Here carbon nanotube diameter predicted DNA damage, rather than the level of surface functionalization.

Pulmonary graphene-induced fibrosis markers
The fibrotic potential of GBM has been investigated in several animal studies. While some reported pulmonary fibrotic reactions in rodents following inhalation, intratracheal instillation or pharyngeal aspiration of graphene or GO (Duch et al., 2011;Wang et al., 2015;Lee et al., 2017), other studies reported no fibrotic findings (Ma-Hock et al., 2013;Kim et al., 2016;Mao et al., 2016). In this study several genes, which have been proposed as key mediators in the onset of fibrosis, were identified in the upstream analyses, mainly as regulators in the GO groups. This included the pro-fibrotic growth factors PDGF-complex and TGFb1, which were identified as activated for GO exposure only. Both GO and rGO activated upstream regulation of Th2 cytokines (IL4 and IL13) and pro-fibrotic cytokines (TNFα, IL6, and IL1β)(Supplementary Table 4). Several metalloproteinases were significantly expressed after GO exposure only (Mmp3, Mmp8, Mmp9, Mmp14 and Mmp19). The proteins encoded by the Mmp genes are generally associated with matrix degradation, but have also been linked with both pulmonary and hepatic fibrosis (Giannandrea and Parks, 2014). Similarly, inhibitor of metalloproteinase (Timp1) was highly upregulated for GO (14-24 fold) compared to rGO (2-fold), and the pathway "Inhibition of Matrix Metalloproteases" was perturbed for GO high dose exposure (Supplementary Table 2). These changes indicate that exposure to GBM, in particular GO, leads to changes at the transcriptional level that could potentially lead to fibrosis. Induction of epithelial-mesenchymal transition has been proposed as possible mechanism behind the fibrotic potential of graphene materials (Liao et al., 2018), however, in the present study we found no up-regulation of key epithelial-mesenchymal transition genes. This is most likely due to the analyzed time point, as epithelial-mesenchymal transition is a late event in the fibrotic process.
A fibrosis biomarker gene set of 17 genes (PFS17) were recently proposed by Rahman and colleagues (Rahman et al., 2020) for assessing the fibrogenic potential of nanomaterials using an ex vivo lung culture and exposure system. A comparison between this list and the transcriptomic changes after GO and rGO exposure revealed that following GO exposure 14, 10 and 4 genes of the 17 genes in PFS17 were found upregulated at dose 18, 54 and, 162 μg, respectively. One gene of the 17 was downregulated at the 54 and 162 μg doses (Supplementary Table 6). For rGO exposure, 6 out of 17 genes were upregulated at both dose 54 μg and 162 μg. These results suggest that GO and rGO has the potential to induce lung fibrosis, with GO being more potent than rGO. This potential may be linked to the surface chemisty of the GBM, via a strong induced inflammatory response, as several of the identified key regulated genes related to fibrogenesis were cytokines (Il4,Il13,Tnfα,Il6,and Il1β). As reported previously, the inflammatory response was almost back to baseline level 28 days after exposure to GO and rGO, and no fibrotic lesions were detected from the histological analyses 90 days after exposure to the GBM, despite the presence of particular agglomerates in the lung tissue (Bengtson et al., 2017). This indicates that sustained inflammation is crucial for GBM-induced fibrosis. However, when looking at our previous MWCNT studies, we found no fibrosis in the lungs of mice 90 days after exposure to a thin and entangled type of MWCNT at a dose of 54 μg, despite sustained inflammation throughout the entire post-exposure period (Poulsen et al., 2016). In contrast, two other MWCNT (one thin/entangled (NRCWE-026) and one thick/rigid (NM-401)) both promoted pulmonary fibrosis 28 days after exposure at dose 162 μg, with sustained inflammation throughout the post-exposure period (Poulsen et al., 2015b). This indicates that more than inflammation is needed to induce fibrosis, at least at the investigated subacute/sub-chronic time points. The dispersion state has previously been highlighted as important for the fibrotic potential of GBM, as welldispersed graphene and GO did not induce fibrosis, whereas aggregated graphene did (Duch et al., 2011).

Liver functional analyses
Pulmonary exposure to rGO had very limited effect on transcriptomic changes in the liver, whereas GO exposure affected hepatic gene expression in similar fashion to MWCNT (Poulsen et al., 2015a). Due to this highly variable number of DEG found in the liver 1 day post exposure (Fig. 1A), a contextual analysis was only conducted on the GO 162 μg exposure group. As previously described, this was conducted in IPA.

Biological functions
Although several functions in each category ("Molecular and Cellular Functions" and "Diseases and Disorders") were perturbed, predicted activation of annotations was scarcer and mostly decreased (Fig. 4). This is in line with a higher number of downregulated genes in the hepatic tissue compared to the pulmonary results (Fig. 1A). Similar to the pulmonary functional analysis, the top regulated function in "Diseases and disorders" was Organismal Injury and Abnormalities (Fig. 4A). However, in contrast to the pulmonary results, this was mainly driven by a predicted decreased in the annotation: Inflammation of liver. Same annotation was also present in the regulated functions: Hepatic System Disease and Inflammatory Response. In Inflammatory Response, this decrease was flanked by a predicted decrease in the annotation: Immune response of cells.
Among the top-regulated "Molecular and Cellular Functions" (Fig. 4B), Small Molecule Biochemistry comprised of several annotations related to synthesis and accumulation of lipids, which were predicted as decreased. In addition, a tendency for activation of lipid and fatty acid transport was also observed, albeit no direction of the lipid transport could be concluded from the analysis. The gene list for the Small Molecule Biochemistry function revealed that the expression of several genes encoding ATP-binding cassette transporters proteins involved in lipid transport were perturbed. However, no consistency in regulation was observed. Similar to "Diseases and Disorders" several annotations related to inflammation were predicted as decreased across "Molecular and Cellular Functions".

Pathway analysis
Several pathways were perturbed after exposure to high dose GO. However, Unfolded protein response was the only pathway with a significance level lower than 0.01 (Supplementary Table 7). A large proportion of the perturbed pathways were involved in lipid synthesis and homeostasis, and they were predicted to be activated. An important pathway for cholesterol homeostasis, LXR/RXR Activation pathway, was significantly perturbed after exposure to high dose GO. A closer look at the pathwayrevealed that cholesterol metabolism and biosynthesis were predicted as increased, whereas lipogenesis and cholesterol transport was predicted as inhibited ( Supplementary Fig. 6). In addition, the entire Superpathway of Cholesterol Biosynthesis was upregulated. This indicates that the liver decreases its cholesterol uptake, increases its cholesterol synthesis and increases its efflux of cholesterol. Similar to that seen for the pulmonary microarray analysis, this indicates enhanced hepatic regulation of cholesterol metabolism after GO, while these pathways were unaffected by rGO exposure. These changes could result in a net increase in plasma cholesterol levels, similar to that seen 3 days after MWCNT exposure (Poulsen et al., 2015a). In general, a large overlap in perturbed hepatic pathways and functions between GO and MWCNT exposure was discovered. These were mainly involved in lipid homeostasis.
Among the regulated genes in the LXR/RXR pathway was 2 cytochrome P450 (Cyp) genes. In total, 19 Cyp genes were differentially regulated in the GO 162 μg exposure group (15 upregulated and 4 downregulated). This large superfamily of hemeproteins is involved in drug, foreign material and lipid metabolism (Nebert and Russell, 2002), and hepatic regulation of several members was also observed after MWCNT exposure (Poulsen et al., 2015a). Interestingly, Cyp7a1, which was differentially upregulated after GO exposure (4.24) was found to be strongly downregulated after exposure to thin entangled type of MWCNT 3 days after exposure (− 20.16). This indicates time-and physicochemical dependent variations in Cyp expression.

Plasma lipid composition
As both lung and liver transcriptomic analyses showed changes to cholesterol metabolism, plasma concentrations of HDL, LDL/VLDL and total cholesterol were investigated at all available time points. Although trends were observed, blood lipid composition were not significantly different from those of the VCs (Supplementary Fig. 7). The only significant difference compared to the vehicle exposed mice was lowered total cholesterol levels 90 days after exposure to 162 μg rGO. This is in contrast to previous observed changes is lipid composition after exposure to MWCNT and CB (Bourdon et al., 2012a;Poulsen et al., 2015a). However, the fact that the cholesterol experiment was conducted with low group numbers (n = 3) may have resulted in too low statistical power to detect a possible effect.

Comparison with other nanomaterials
The experimental setup in the current study is very similar to that of our previous transcriptomic studies using MWCNT, CB and nano-sized TiO 2 , in regards to animals, doses, exposure route, and time points (Bourdon et al., 2012a;Husain et al., 2013;Poulsen et al., 2013;Halappanavar et al., 2015;Poulsen et al., 2015b;Rahman et al., 2016). It was therefore possible to compare the pulmonary expression profiles from the different nanomaterial exposures through hierarchal clustering. First, all time points were included ( Supplementary Fig. 8). The data set visibly divided into two different clusters: the TiO 2 cluster and the MWCNT cluster. rGO (all doses) and the lowest dose of GO (18 μg) clustered with TiO 2 , whereas middle and high doses of GO clustered with CB and MWCNT in the MWCNT cluster ( Supplementary Fig. 8), albeit closer to MWCNT than CB. In the second cluster analysis involving only post-exposure day 1 samples, four clusters were identified: Two TiO 2 clusters, a MWCNT cluster and a rGO cluster. GO (all doses) clustered with middle and high doses of MWCNT and CB, whereas rGO clustered distinctly on a separate branch with low dose CB (Supplementary Fig. 9). This supports the contextual analyses presented in this paper, which concluded that GO and rGO exposure induced very different responses at many levels, and that GO exposure to a higher degree induces gene expression changes similar to MWCNT, compared to rGO.
Transcriptomic data also enables other ways of comparing and ranking nanomaterial responses. Here, we compared the number of upand downregulated genes 1 day post-exposure to 162 μg GO, rGO, nano-CB, and different types of MWCNT and nano-TiO 2 (Supplementary Table 8). At this time point and dose, GO and NRCWE-026 display the largest potency with most DEG, followed by Mitsui7, NM401 and rGO. CB and some of the TiO 2 types grouped together in potency and then lastly the remaining TiO 2 types displayed very few DEG. This is only one of multiple ways the data sets can be utilized for comparison and ranking. Bench Mark Dose response analyses on transcriptomic data has previously been applied to rank nanomaterials (Halappanavar et al., 2019). They highlighted the long, stiff MWCNT as especially potent compared to the other nanomaterials and of high priority for further toxicity testing. Transcriptomic data pave the way for very detailed comparisons across nanomaterial exposures, and future comparisons and ranking would be both interesting and could prove very valuable for identification of high priority nanomaterials.

Conclusion
Here, we show that pulmonary exposure to GO and rGO results in different transcriptomic responses in the lung and liver 1 day postexposure. Overall, GO exposure was more potent in inducing DEG, and affecting functions and canonical pathways compared to rGO exposure, both in lung and in liver. Although both GBM induced transcriptomic changes associated with pulmonary inflammation, the response after GO exposure was stronger and specifically involved the innate immune system. In contrast, the response to rGO was weaker, with a larger proportion of DEG related to the adaptive immune system. In addition, only GO exposure activated the APR in the lung and affected lipid homeostasis in the liver, which are changes that have been linked to increased risk of atherosclerosis development. GO exposure also resulted in more transcriptomic changes related to fibrosis compared to rGO, although this was not reflected in the previously conducted histopathological analysis. GO and rGO share similar physical properties, but differ greatly in their surface chemistry. Thus, the results suggest that a larger number of hydroxyl groups on the surface of GO may be causative to the robust transcriptomic responses observed. In contrast, GO and rGO induced similar changes related to ROS production and genotoxicity, indicating that these changes are unrelated to surface chemistry.
GO and rGO were compared to other nanomaterials and the hierarchal cluster analyses revealed that GO clustered with MWCNT, whereas rGO tended to cluster separately or with TiO 2 . This highlights the difference in potency of the two GBM. Conducting such comparisons is an important first step in the ranking of nanomaterials, and could prove a valuable tool for future risk assessment of nanomaterials.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.