Exposure to graphene oxide sheets alters the expression of reference genes used for real-time RT-qPCR normalization

Studies unraveling the interactions between graphene oxide (GO) and the biological milieu, including cells and tissues, are multiplying quickly as the biomedical applications of this and other 2D materials continue to be explored. Many of such studies rely on real-time RT-qPCR as a powerful yet simple technique to assess gene expression. However, a systematic investigation of potential GO-induced changes in the expression of reference genes, crucial for appropriate qPCR data normalization, is still lacking. We aimed to cover this gap investigating the stability of the expression of ten candidate reference genes upon exposure to increasing, but subtoxic, GO concentrations, with two established algorithms (Bestkeeper and NormFinder). The study was performed in a human cancer cell line (MCF7) and in mouse, non-cancerous, primary cells (mouse embryonic fibroblasts, MEFs), to assess different behaviors between cell types. Both algorithms evidenced significant deviations in the expression of various reference genes. Ribosomal proteins scored among the most significantly dysregulated in both cell types. ACTB and GAPDH, the most frequent calibrators in real-time RT-qPCR, were also affected, although differences existed between cell lines. This study illustrates the need to validate reference genes for appropriate real-time RT-qPCR normalization, according to specific experimental conditions, when GO-cell interactions occur.

economical demands of deep sequencing 9 . However, correct normalization of relative gene expression data is critical for the reliability of the technique and has been the topic of numerous studies covering both mathematical approximations and the adequacy of calibrators [10][11][12] . To account for possible differences in the amount of starting material, RNA recovery and integrity, efficiency of cDNA synthesis and overall transcriptional activity of the cell or tissue of interest, all of which can easily compromise the accuracy of the analysis if not properly corrected, so-called "reference" or "housekeeping genes" are used as calibrators to normalize real-time RT-qPCR data. Traditionally, those are selected among genes involved in very fundamental cellular functions (e.g. transcription, translation, cytoskeletal structure and basic metabolic pathways) under the assumption that they are constitutively expressed across different cell types, tissues and conditions 13 . However, several reports have unveiled the lack of stable expression of many of such traditionally considered housekeeping genes under a number of circumstances, including different cell types and tissues 14,15 , physiological and diseased states [16][17][18] and experimental conditions 19,20 . The MIQE guidelines, released in 2009 to promote transparency and good practice in RT-qPCR studies, strongly advise against the normalization of gene expression data with a single, non-validated, reference gene 21 . Several algorithms and software tools, including Bestkeeper 22 , NormFinder 23 and GeNorm 24 , have been developed to facilitate the validation of intended calibrators.
It is of special concern that the impact of GO exposure in the expression of commonly used reference genes remains, to our knowledge, unexplored. Even more alarming is the absence of such validation in a number of studies, reviewed elsewhere 25 , where GO was used as a component of siRNA and mRNA delivery vectors and that relied on RT-qPCR to evaluate the efficacy of the silencing or forced gene expression achieved.
In this study, we aimed to assess the stability of the expression of ten candidate reference genes (Tables 1, 2) upon in vitro exposure to sub-cytotoxic concentrations of highly characterized, endotoxin-free, GO flakes. We performed the study in the human cancer cell line MCF7 and in murine, non-cancerous, primary cells (mouse embryonic fibroblasts, MEFs), to investigate if the impact of GO in gene expression depended on cell type. Real-time RT-qPCR analysis in this work was performed in strict compliance with MIQE guidelines to ensure the reliability of the results, and Bestkeeper and NormFinder algorithms were used to rank the performance of the candidates as stable reference genes. A more in-depth analysis of the impact of GO exposure was performed for the most dysregulated genes in each cell type.

Results
Selection of candidate reference genes. We based our selection of candidate reference genes on a literature search to identify those most commonly used in real-time RT-qPCR normalization and on a meta-analysis conducted by de Jonge et al. 26 that analyzed 13,629 human and 2,543 mouse gene arrays from previous publications. De Jonge et al. identified various novel housekeeping genes-including RPS13, RPL27, RPL30 and OAZ1-whose expression proved significantly more stable than that of more commonly used calibrators, such as GAPDH and ACTB. We paid special attention to select candidates with different cellular functions, to minimize

RnA samples and Rt-qpcR reactions.
To study the impact of the material on the expression of candidate reference genes, MCF7 cells and MEFs were exposed to increasing concentrations of endotoxin-free GO (0, 5, 10 and 50 μg/ml). This range was selected based on previous studies from our laboratory that confirmed lack of toxicity in a number of cell lines, including MCF7 27 and MEFs 28 used in this work. Full characterization of the material, produced in house by a modified Hummer's method as previously described 29,30 , has been reported in a previous publication 31 and the most relevant parameters are summarized in Table S1. In brief, lateral dimensions did not surpass the 2 μm threshold and thickness corresponded to 1-2 single GO layers. Functionalization degree was estimated as 41%. Exposure took place in the absence of FBS for the first 4 h, to mimic conditions commonly used when testing nanomaterials-cell interactions in vitro. Gene expression was assessed 24 h after the initial exposure by real-time RT-qPCR. All procedures were performed in strict compliance with MIQE guidelines 21 , to ensure reliability of the results. To provide the transparency required by these recommendations, a MIQE guidelines checklist is provided in Table S2. To avoid error introduced by poor RNA quality, A 260/A280 and A 260/230 ratios of all RNA samples included in the study ranged between 1.70 and 2.1, their RIN values were >8.9, and their 28 S/18 S ratios ranged between 2.0-3.9 (Tables S3, S4). Primers, designed in house for the study, amplified all transcription variants of each target with equal product length. Their details are given in Tables S5, S6, including the efficiencies of qPCR reactions (E), determined by serial dilution of the cDNA template.
Expression stability of candidate reference genes in MCF7 cells exposed to GO. We first used Bestkeeper software 22 to obtain preliminary information regarding the expression of each candidate reference gene in MCF7 cells treated with GO. Bestkeeper provides the standard deviation (SD) and the coefficient of variance (CV) of the quantification cycles (Cq, expressed as crossing point (CP) in the original software), for each gene. CV is calculated as the percentage of the Cq SD to the Cq mean. Genes with SD >1 are considered to have an unacceptable range of variation 22 . As reported in Table 3, all gene candidates showed Cq SD <1, thus none had to be excluded from further analysis on such grounds. However, β-Actin (ACTB) clearly scored as the less stably expressed gene among all candidates, with SD = 0.58 and CV = 4.15 (Table 3) Table 3). Bestkeeper assesses stability of expression based on the two parameters above (SD and CV) but does not offer a strict ranking of the most stable housekeeping genes. It provides a normalization index (Bestkeeper index) to normalize each sample, calculated as the geometric mean of Cqs from the best performing candidates 22 . We used the model-based NormFinder approach 23 to obtain a precise ranking. In this case, the algorithm accounts for all different experimental groups included in the study and considers intra-and inter-group variation to provide a direct measure of stability, defined as stability value. The lower this value, the higher the stability of expression of the candidate 23 . In agreement with Bestkeeper data, ACTB and the three ribosomal proteins (RP27, RPL30 and RPS13) scored the highest (less stable) values (0.367, 0.264, 0.230 and 0.136, respectively) ( Fig. 1). GAPDH was ranked as the most stably expressed gene (stability value = 0.029). NormFinder also defines the stability value for the best combination of two reference genes. In this case, the combination of GAPDH and HMBS did not improve, but obtained the same stability value as GAPDH alone, 0.029 (Fig. 1).
Overall, these data highlighted differences in the stability of expression among various common reference genes utilized for RT-qPCR data normalization, when MCF7 cells were treated with different concentrations of GO. Expression stability of candidate reference genes in MEFs exposed to GO. To confirm whether the observations above where exclusive to MCF7, a human cancer cell line, we repeated the study with  Table 3. Descriptive statistics from Cq values extracted from Bestkeeper algorithm, MCF7 study. Geometric mean (geoMean), standard deviation (SD) and coefficient of variance (CV) of Cq data from ten candidate reference genes in MCF7 cells exposed to increasing concentrations of GO (n = 12). CV is calculated as the percentage of the Cq SD to the mean Cq.
www.nature.com/scientificreports www.nature.com/scientificreports/ non-cancerous, primary mouse cells. MEFs were isolated from E12.5 embryos from the CD1 background, cultured for no more than three passages, and treated under the same conditions used in the MCF7 study. The homologue genes were evaluated as candidate references, since they are conserved in the mouse and human genomes ( Table 2), but specific primer pairs were designed (Table S6). Descriptive statistics from Cq values retrieved from Bestkeeper software are summarized in Table 4. The mean Cq values of several candidates differed from those registered in MCF7 cells (Table 3), which underlines the need to validate the calibrator of choice when comparing gene expression in different cell types. As in the MCF7 cell line, Actb Cqs showed the highest variation among all candidates in MEF samples (SD = 0.21, CV = 1.47). However, SD and CV values were generally lower that those observed in MCF7 cells, which evidenced an overall higher stability of expression of the candidates upon GO treatment. The same observation was confirmed by NormFinder analysis, whereby the stability values of all candidates ranged between 0.172 and 0.071 in MEFs (Fig. 2), while they spanned from 0.367 to 0.029 in the MCF7 study (Fig. 1). This algorithm identified Rpl27 as the most unstable gene in MEFs exposed to GO, with a stability value of 0.172, and closely followed by Gapdh (0.154). On the opposite side of the ranking, Hmbs was the most stable candidate (stability value = 0.071), but the combination of Rpl30 and Tbp provided an even lower stability value (0.041) (Fig. 2).
Go treatment induces dose-dependent RPL27 downregulation in both cell types. Moved by the poor stability scores of RPL27 in both MCF7 cells and MEFs, we decided to investigate more closely the relationship between the GO treatment and the mRNA levels of this gene. In MCF7 cells, RPL27 expression was normalized to the geometric mean of GAPDH and HMBS Cqs, the best combination of two reference genes inferred from NormFinder analysis (Fig. 3a). Following the same algorithm, Rpl30 and Tbp were used as calibrators to normalize Rpl27 data in MEFs (Fig. 3b). In both cases, we found a dose-dependent and statistically significant Figure 1. Stability values of ten reference genes in MCF7 cells exposed to GO, calculated with NormFinder algorithm. Candidate reference genes are ranked according to decreasing stability values (i.e. increasing stability of expression). GAPDH was ranked as the most stable gene. The most stable combination of two genes (GADPH + HMBS) is also represented. See Table 1 for gene names.  Table 4. Descriptive statistics from Cq values extracted from Bestkeeper algorithm, MEF study. Geometric mean (geoMean), standard deviation (SD) and coefficient of variance (CV) of Cq data from ten candidate reference genes in MEFs exposed to increasing concentrations of GO (n = 12). CV is calculated as the percentage of the Cq SD to the mean Cq.
downregulation of the target in the presence of GO, which confirmed the inadequacy of RPL27 as reference gene under such experimental conditions. The dysregulation was more pronounced in MCF7 cells, where exposure to the lowest concentration of GO tested (5 μg/ml) already induced an 18% downregulation of RPL27 (p = 0.03). At 50 μg/ml, the downregulation reached 56% (p = 0.00005). In MEFs, however, the latter concentration was the only one to produce a statistically significant change in Rpl27 expression (25% downregulation, p = 0.02). Data normalization using the Bestkeeper index confirmed these results ( Figure S2).
Archetypal housekeeping genes are dysregulated in the presence of GO. ACTB and GAPDH are without doubts the two genes most frequently utilized to normalize real-time RT-qPCR data 13 . However, various studies have called to question their status as housekeeping genes, based on significant fluctuation of their expression under specific treatment or experimental conditions and even among different tissues of the same organism 16,32,33 . Indeed, in de Jonge et al. 's ranking of stable human reference genes, ACTB and GAPDH occupied positions 57 and 139, respectively 26 . In our study, ACTB was ranked the most unstable candidate reference gene in MCF7 cells treated with GO, by both Bestkeeper and NormFinder algorithms (Fig. 1, Table 3). In this cell line, ACTB mRNA levels were significantly upregulated in the presence of 50 μg/ml GO (p = 0.04, Fig. 4a). In MEFs, Actb showed the highest Cq CV (Table 4), albeit this variation did not translate into significant changes in the expression of the transcript (Fig. 4c). GAPDH returned even more discrepant results between the two cell types. It proved the most stable candidate in MCF7 cells, according to NormFinder results (Fig. 1), and consequently its mRNA relative levels remained stable across all experimental groups (Fig. 4b). However, in MEFs, the same gene scored the second highest stability value (0.154) upon NormFinder analysis (Fig. 2) and showed the second highest Cq CV in Bestkeeper (Table 4). Indeed, 50 μg/ml GO induced a 1.44-fold Gapdh upregulation (p = 0.001, Figure 2. Stability values of ten reference genes in MEFs exposed to GO, calculated with NormFinder algorithm. Candidate reference genes are ranked according to decreasing stability values (i.e. increasing stability of expression). Hmbs was ranked the most stably expressed single gene, but the combination of (Rpl30 + Tbp) provided an even lower stability value. See Table 2 for gene names. www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 4d). All such results were confirmed when Bestkeeper index was used as normalization factor ( Figure S3). Overall, these findings corroborate the inadequacy of ACTB and GAPDH as reference genes under particular experimental conditions in which GO is involved.

Discussion
We have shown here that GO induces significant changes in the expression of common reference genes, at the mRNA level. This is, to our knowledge, the first study to systematically explore the impact of GO exposure in the stability of traditionally considered housekeeping genes, a gap that needed to be filled considering the ever-growing number of studies exploring the GO-cell biology interface. Strikingly, the validation of suitable reference genes was so far absent even in studies in which real-time RT-qPCR data constitutes a central readout of the GO-cell interaction. Most of such works relied on GAPDH, exclusively, as calibrator, without data that supported the stability of its expression under their specific experimental conditions [34][35][36][37][38] . Another study did not even specify the reference used for normalization, nor whether it had been validated 39 .
The above is a matter of concern because (a) MIQE guidelines strongly discourage the use of a single, non-validated gene for standardization 21 and (b) the inadequacy of GADPH as calibrator has been evidenced by a number of studies that detected large variations of its mRNA levels across tissue types and under different experimental conditions 16,32 . Our results indeed support these observations and the current thought that a specific gene may behave as stable reference in some, but not all, scenarios. Gapdh mRNA levels were significantly upregulated by high (but subtoxic) concentrations of GO in MEFs (Fig. 4d, Figure S3d), rendering it of no use as calibrator. However, this was not the case in MCF7 cells, in which GAPDH scored as the most stable candidate reference gene (Fig. 1). A similar scenario was found for ACTB, also a very popular and widespread calibrator that has received increased attention due to the inconsistency of its expression 16 . ACTB showed the highest stability value, according to NormFinder analysis, in MCF7 cells (Fig. 2) and its expression significantly increased with GO exposure (Fig. 4a, Figure S3a).
One can also infer from our results that ribosomal proteins are, overall, not appropriate calibrators to normalize real-time RT-qPCR data in the presence of GO. All such candidates included in our study, RPL27, RPL30 and RPS13, scored the highest instability in MCF7 cells, in that order and only surpassed by ACTB (Fig. 1, Table 3). RPL27, in particular, experienced significant downregulation in the presence of GO, which was dose-dependent and confirmed in both cell types, MCF7 and MEFs (Fig. 3). These results oppose those reported in de Jonge et al. 's meta-analysis, where RPS13, RPL27 and RPL30 occupied top positions (first, second and fourth, respectively) in the stability ranking of candidate housekeeping genes 26 . This discrepancy highlights even more the necessity to www.nature.com/scientificreports www.nature.com/scientificreports/ validate the reference genes of choice in the presence of GO, as interaction with the nanomaterial seems to dysregulate genes that are otherwise rather stably expressed across an ample variety of conditions. Although RPL27 mRNA levels followed the same trend in MCF7 cells and MEFs exposed to GO, we also found important differences in the impact of the material in the gene expression profiles of both cell types. First, instability introduced in the expression of candidate reference genes was more pronounced in MCF7 cells compared to MEFs, in which Cq variation was narrower (see SD and CV values in Tables 3, 4, and stability values in Figs 1, 2). Second, we found differences in the ranking of most stable housekeeping genes-even when the same algorithm was used-with stable references in one cell type scoring very poorly in the other (see for instance GAPDH data in Figs 1, 2, as discussed above). The discrepancies between MCF7 and MEFs may be due to a combination of cell-intrinsic differences in basal mRNA expression and differences in the way the GO material interacts with each cell type, but overall stresses the absolute requirement to validate housekeeping genes under specific experimental conditions. Indeed, inconsistencies in the expression of traditionally considered reference genes between different cell types have already been described, even when such cells were not subjected to any exogenous treatment 14,15 . In addition, GO, as well as other nanomaterials, is known to interact differently with different cell types, including in what concerns to gene expression 8 . Given the many parameters that can alter gene expression, our study does not intend to provide a list of reference genes to be used in in vitro GO studies, but to highlight the necessity of their systematic validation under specific experimental conditions.
We also found differences within the rankings provided by Bestkeeper and NormFinder (Figs 1, 2, Tables 3, 4). Such discrepancies were expected, based on the different algorithms that support each software, and have indeed been reported by other studies prior to ours 14,15,18,40 . However, we show here that choosing NormFinder or Bestkeeper index for normalization did not change the results obtained for the relative expression of various targets (Figs 3, 4, S2 and S3), which grants additional significance to the results shown here.
In spite of the GO-induced changes in gene expression shown here, it is nevertheless to be said that the instability inferred by the material was less pronounced than that triggered by other physiological or experimental conditions. In our study, the highest stability value recorded was 0.367 for ACTB in MCF7 cells and 0.172 for Rpl27 in MEFs (Figs 1, 2 18 . This comparison, however, does not eliminate the need for appropriate validation of reference genes, since the GO-induced changes in their expression that we have reported here would introduce statistically significant errors in the normalization of relative gene expression data.
The work described here tackles a very specific issue, the common lack of and absolute need for validation of reference genes utilized in gene expression studies. Therefore, the investigation of changes at the protein level is out of the scope of this study and it is indeed not required by the MIQE guidelines to ensure reproducibility and high quality of gene expression data. Besides, the sensitivity of common tools utilized to assess the expression of specific proteins, such as the Western blot technique, is known to be significantly lower than that of real-time RT-qPCR. Therefore, it is anticipated that very small changes in mRNA levels that have been shown here to significantly affect real-time RT-qPCR normalization, will not be detected at the protein level.

conclusion
We have shown here that in vitro exposure to GO sheets alters the expression of various candidate reference genes used to normalize real-time RT-qPCR data, including very frequently used calibrators such as GAPDH and ACTB. We have also demonstrated that the magnitude and nature of the changes induced vary between different cell types and therefore reference gene validation cannot be extrapolated but must be specifically determined according to defined experimental conditions. Those may include, but may not be limited to, physicochemical characteristics of the material, dose, exposure conditions and cell type of interest. Using stable reference genes is imperative to obtain reliable gene expression data.

Methods
Graphene oxide. GO was produced in house following a modified Hummer's method as previously described 29,30 . Full characterization of this material was reported in a previous publication 31 , where it was termed small GO (s-GO) in order to differentiate it from larger GO flakes that have not been used in the present study. A summary of this information is provided in Table S1. In brief, the lateral dimensions of GO flakes do not surpass 2 µm and their thickness corresponds to 1-2 layers of GO (1-2 nm). The degree of functionalization was estimated as 41% by thermogravimetric analysis (TGA). Presence of oxygenated functionalities in form of hydroxyls, carboxyls and epoxides was confirmed via x-ray photoelectron spectroscopy (XPS). Surface charge was strongly negative (ζ = −55.9 ± 1.4 mV).
Primary cell extraction, cell lines and culture. The MCF7 human breast cancer cell line was obtained from the American Tissue Culture Collection (ATCC, HTB-22 ™ ) and cultured in Eagle's Minimum Essential Medium (MEM, M4655, Sigma) supplemented with 10% fetal bovine serum (FBS, 10500, Gibco, Lot 08G3057K) and 1% antibiotics (PenStrep, P4333, Sigma). Cells were maintained in a 5% CO 2 atmosphere at 37 °C. Authenticity was verified at the DNA sequencing facility of The University of Manchester by Short Tandem Repeat (STR) analysis. Mouse embryonic fibroblasts (MEFs) were extracted from E12.5 CD1 embryos following a standard protocol 41 and maintained in Dubelco's Modified Eagles Medium (DMEM, D6429, Sigma) supplemented with 10% FBS and 1% antibiotics, in a 5% CO 2 atmosphere at 37 °C. MEFs were used for a maximum of three passages. primer design. Details of the primers used in this study are reported in Tables S4 and S5. Primers were designed with Primer Basic Local Alignment Search Tool (Primer BLAST) adhering to the following requisites: amplicon size 75-200 bp, GC content 50-65%, ≤3 G or C repetitions, ≤4 base repetitions, melting temperature (Tm) 55-65 °C. When gene targets had several splicing variants (including predicted variants), primer pairs were designed to amplify all of them at the same product length. Each primer pair was verified with Blast tool (NCBI) to confirm the its specificity for the desired target. Primers were synthetized by Sigma (UK) and purified by desalting.
Real-time quantitative Polymerase Chain Reaction (qPCR). 2 ul of cDNA were used for each 20 μl qPCR reaction, which was performed with SYBR green chemistry (PowerUp TM SYBR TM Green Master Mix, A25742, Applied Biosystems) that contains Dual-lock ™ Taq DNA polymerase. Reactions were loaded manually (in technical duplicates) on white-wall, clear-well, hard-shell 96-well plates (HSP9601, Biorad) and run on a CFX96 Touch real-time PCR detection system (Biorad) according to the following protocol: 2 min at 50 °C, 2 min at 95 °C, (15 sec at 95 °C and 1 min at 60 °C) x 40 cycles. The reaction was followed by a melt curve analysis, as specified by the manufacturer, to confirm amplification of a single product. No amplification was detected in non-template controls (NTC) and non-reverse transcription (NRT) controls. Cq <40 for all detections, which evidenced establishment of the limit of detection. Each gene was assessed in the same run for the totality of the samples, to avoid inter-run variability. Cq values were determined with the Single Threshold mode in the CFX Manager software (Biorad). PCR efficiencies (E) were calculated from the slope of a linear regression of the Cq values obtained from a dilution series of the starting cDNA, following the equation = − ( ) E 10 slope 1 Descriptive statistics of Cq values (Bestkeeper algorithm). Cq values from technical duplicates were averaged and the Bestkeeper Excel tool was used to retrieve descriptive statistics (n = 12). Variation was expressed as standard deviation (SD) and coefficient of variance (CV), the latter calculated as the percentage of the Cq SD to the Cq mean. Bestkeeper index for data normalization was calculated as the geometric mean of the most stable candidate genes, identified by repeated pair-wise correlation analysis, as described by Pfaffl et al. 22 .
Assessment and ranking of reference gene stability (NormFinder algorithm). The model-based NormFinder algorithm, as described by Andersen et al. 23 , was used to assess the stability of the expression of ten candidate reference genes across all experimental groups. Cq values (n = 12) were transformed to relative quantities, according to the formula: E (lowest Cq − Cq) , that takes into account the efficiency of the PCR reaction (E) and uses the lowest Cq as a calibrator. Stability values were calculated for each candidate housekeeping gene, taking into account intra and intergroup variation. The stability value of the best combination of two genes was also calculated.
www.nature.com/scientificreports www.nature.com/scientificreports/ Relative gene expression analysis. Relative gene expression was calculated following the Livak method 10 .
As calibrator, either the geometric mean of Cqs of two reference genes indicated by NormFinder algorithm, or the Bestkeeper index, was used. All data was normalized to the control (untreated) group. Error was propagated according to the formula: Error a Error b ( ) ( ) ( ) 2 2 Statistical analysis of relative gene expression data. Differences in gene expression among experimental groups were assessed with OriginPro v9.1, using ∆Cq data. First, normal distribution and homogeneity of variances were confirmed (Levene's test). Differences between groups were assessed by one-way ANOVA and Tukey's post-hoc test. P values ˂0.05 were considered significant.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.