Amorphous silica nanoparticles induce tumorigenesis via regulating ATP5H/SOD1-related oxidative stress, oxidative phosphorylation and EIF4G2/PABPC1-associated translational initiation

Background Recent studies indicate amorphous silica nanoparticles (SiNPs), one of the widely applied nanomaterials, have potential toxicity in humans and induces cell malignant transformation. However, its carcinogenic mechanisms remain poorly understood. This study’s purpose was to investigate the underlying toxic mechanisms of amorphous SiNPs on human lung epithelial cells model by using microarray data. Methods Microarray dataset GSE82062 was collected from Gene Expression Omnibus database, including three repeats of Beas-2B exposed to amorphous SiNPs for 40 passages and three repeats of passage-matched control Beas-2B cells. Differentially expressed genes (DEGs) were identified using linear models for microarray data method. Protein–protein interaction (PPI) network was constructed using data from the STRING database followed by module analysis. The miRwalk2 database was used to predict the underlying target genes of differentially miRNAs. Function enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool. Results A total of 323 genes were identified as DEGs, including 280 downregulated (containing 12 pre-miRNAs) and 43 upregulated genes (containing 29 pre-miRNAs). Function enrichment indicated these genes were involved in translational initiation (i.e., eukaryotic translation initiation factor 4 gamma 2 (EIF4G2), poly (A) binding protein cytoplasmic 1 (PABPC1)), response to reactive oxygen species (i.e., superoxide dismutase 1 (SOD1)) and oxidative phosphorylation (i.e., ATP5H). PABPC1 (degree = 15), ATP5H (degree = 11) and SOD1 (degree = 8)] were proved to be hub genes after PPI-module analyses. ATP5H/SOD1 and EIF4G2/PABPC1 were overlapped with the target genes of differentially expressed pre-miR-3648/572/661 and pre-miR-4521. Conclusions Amorphous SiNPs may induce tumorigenesis via influencing ATP5H/SOD1-related oxidative stress, oxidative phosphorylation and EIF4G2/PABPC1-associated translational initiation which may be regulated by miR-3648/572/661 and miR-4521, respectively.


INTRODUCTION
Nanomaterials refer to those with sizes ranging from 1 to 100 nm in at least one dimension. Nanosized particles possess a number of superior physicochemical properties compared with the same material fabricated in a conventional manner, including high thermal and chemical stability, hydrophobicity, heat and electrical insulation, resistance to oxidation, good biocompatibility and minimal immunogenicity, which make them as attractive and promising candidates for a wide range of advanced applications (Liang et al., 2008;Peng et al., 2014). Silica is the most frequently used material to create nanoparticles due to its higher abundance on earth and thus a relatively low-cost for preparation. It is estimated that approximately one million tons of synthetic amorphous silica nanoparticles (SiNPs) may be produced per year worldwide to act as additives to cosmetics, drug tablets, paints, varnishes, food or deliveries for gene, protein and drugs (Fruijtier-Pölloth, 2012). The widespread application of amorphous SiNPs leads to frequent human exposure and therefore its safety is of high concern. Although amorphous SiNPs are supposedly non-carcinogenic for humans according to the classification of the International Agency for Research of Cancer (Group 3), recently published studies proposed long-term exposure of amorphous SiNPs may induce cell malignant transformation and tumorigenesis (Fontana et al., 2017;Guo et al., 2017). Therefore, how to early diagnose and prevent the development of amorphous SiNPs-induced cancer may be an underlying challenge that needs to be solved. This resulted in the requirements for understanding the molecular mechanisms of the tumor-promoting effects.
In a recent study, Guo et al. (2017) used a microarray analysis to investigate the genes significantly changed by amorphous SiNPs in human lung epithelial cells and found amorphous SiNPs may trigger tumorigenesis by influencing 821 significant differentially expressed genes (DEGs) (five upregulated and 816 downregulated) to regulate oxidative stress, oxidative phosphorylation, DNA damage, p53 and MAPK signaling pathways. However, the related mechanism leading to the development of cancer by amorphous SiNPs still remains unclear. In present study, we aim to re-analyze the microarray data established by Guo et al. (2017) via different bioinformatic approaches: the DEGs were identified by the linear models for microarray data (LIMMA) method, but not random variance model; Compared with study of Guo et al. (2017), a strict threshold was used for screening DEGs (|logFC(fold change)| > 2 & P < 0.05 vs FC > 2 & P < 0.05), which may be beneficial to obtain more crucial and verifiable genes associated with amorphous SiNPs; the whole protein-protein interaction (PPI) network for all DEGs were established, but not signal-net analysis network PPI; In addition, the miRNA-target genes interaction network was also analyzed to explore the regulatory mechanisms of DEGs and then screen key upstream targets for amorphous SiNPs, which had not previously been performed.

Data normalization and DEGs identification
The raw data (CEL files) of GSE82062 were downloaded from the Affymetrix Human Transcriptome Array 2.0 platform GPL17586. The raw data were preprocessed, background-corrected and summarized using robust multichip average algorithm (Irizarry et al., 2003) in the "affy" package of Bioconductor R (v3.4.1; http://www.bioconductor.org/ packages/release/bioc/html/affy.html). The DEGs between BeasSiNPs-P40 and Beas-P40 groups were identified using the LIMMA method (Ritchie et al., 2015) in the Bioconductor R package (v3.4.1; http://www.bioconductor.org/packages/release/bioc/ html/limma.html). The t-test was used to identify the P-value and FC was calculated. Genes were considered to be differentially expressed with the threshold value setting to |logFC| > 2 and P < 0.05.
To determine the specificity of DEGs between BeasSiNPs-P40 and Beas-P40 groups, bidirectional hierarchical clustering analysis with Euclidean distance (Szekely & Rizzo, 2005) was performed for the top 50 DEGs using the pheatmap package in R (version 1.0.8; Kolde, 2015).

PPI network construction and module analysis
The DEGs were imported into STRING database (v10.0; Search Tool for the Retrieval of Interacting Genes; https://string-db.org/) (Szklarczyk et al., 2015) to obtain the PPI data. The PPI network was constructed and visualized using Cytoscape software (v2.8; www.cytoscape.org/) (Kohl, Wiese & Warscheid, 2011). The hub genes with more interactions with other partners (degree) were selected and plotted with ggplot2 in R package (v3.4.1; R Core Team, 2017).
Furthermore, the Molecular Complex Detection (MCODE) plugin of Cytoscape software was also employed to identify functionally related and highly interconnected clusters from the PPI network with degree cutoff of 6, node score cutoff of 0.2, k-core of 5 and maximum depth of 100 (ftp://ftp.mshri.on.ca/pub/BIND/Tools/MCODE) (Bader & Hogue, 2003). Sub-modules were defined to be significant with MCODE score 4 and nodes 6.

Function enrichment analysis
Gene ontology (GO) biological process terms, Kyoto Encyclopedia of Genes and Genomes (KEGG), and BioCarta pathways were also enriched for all DEGs, DEGs in PPI and miRNA-DEG interaction networks by the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (v6.8; http://david.abcc.ncifcrf.gov) (Dennis et al., 2003). P-value < 0.05 was chosen as the threshold value for functional enrichment analyses.

Identification of DEGs
After preprocessing and data normalization (gene expression in all samples are shown in Supplemental Information 2), DEGs between BeasSiNPs-P40 and Beas-P40 groups were identified by the LIMMA method. As a result, 323 genes were considered to be differentially expressed due to exceeding the difference threshold (|logFC| > 2 and P < 0.05), including 280 downregulated (containing 12 pre-miRNAs) and 43 upregulated genes (containing 29 pre-miRNAs) (Fig. 1), which was fewer than the study results of Guo et al. (2017). All the DEGs are listed in in Supplemental Information 3.

Function enrichment analysis of DEGs
The above differential genes were subjected to the online tool DAVID for function enrichment analysis. As a result, 49 GO biological process terms were enriched for downregulated DEGs, including response to reactive oxygen species (ROS) (i.e., superoxide dismutase 1 (SOD1)). Only two GO biological process terms were enriched for upregulated DEGs, mainly involving in protein translation (i.e., eukaryotic translation initiation factor 4 gamma 2 (EIF4G2); poly (A) binding protein cytoplasmic 1 (PABPC1)) ( Table 1).
Kyoto Encyclopedia of Genes and Genomes pathway analyses were also performed by DAVID software. In line with the GO terms, the downregulated DEGs were predicted to participate in oxidative stress related metabolism pathways, such as oxidative phosphorylation (i.e., ATP5H, ATP synthase, H+ transporting, mitochondrial F0 complex, subunit D), while upregulated DEGs were enriched for Regulation of eIF4e and p70 S6 Kinase (i.e., EIF4G2, PABPC1) ( Table 2).
Following cluster analysis using the MCODE algorithm, three significant modules were obtained (Fig. 3). Function enrichment analysis showed that ATP5H was included in module 3 and participated in Oxidative phosphorylation pathway (Table 3).

DISCUSSION
After PPI hub genes and module analysis as well as miRNA target gene prediction, our present study preliminarily suggests downregulated ATP5H, SOD1 and upregulated EIF4G2, PABPC1 may be especially important genes involved in amorphous SiNPsmediated tumor initiation. EIF4G2 and PABPC1 may be involved in the progression of cancer by affecting translational initiation, while SOD1 and ATP5H may participate in carcinogenesis via influencing oxidative stress and oxidative phosphorylation. Although genes were not similar, these enriched pathways seemed to be in line with the study of Guo et al. (2017), further indicating these biological processes may be crucial. Recently, accumulating evidence has indicated exposure to amorphous SiNPs induces oxidative stress in cells or organs (Guo et al., 2015;Wu et al., 2016;Nemmar et al., 2016). For example, Jiang et al. (2016) incubated the erythrocytes with SiNPs and found the oxidative damage biomarker malondialdehyde was significantly increased, while the activity of antioxidant superoxide dismutase (SOD) was decreased. Di Cristo et al. (2016) showed exposure of SiNPs to macrophages elicited a greater oxidative stress than was assessed from heme oxygenase-1 induction and ROS production. Nemmar et al. (2016) observed intraperitoneal administration of SiNPs induced significantly increased lipid peroxidation in the lung, liver, kidney and brain of mice, displaying reduced SOD and catalase activities. Consistent with these studies, our study also found anti-oxidative gene SOD1 was significantly downregulated after SiNPs treatment. The increased oxidative stress was reported to induce the switch of glucose metabolism from oxidative phosphorylation to aerobic glycolysis (the Warburg Effect) to promote excessive proliferation and growth of cells, leading to the development and progression of cancer (Pani, Galeotti & Chiarugi, 2010;Molavian, Kohandel & Sivaloganathan, 2016).  ATP synthase is an enzyme to be responsible for the synthesis of abundant ATP in oxidative phosphorylation process. The expression of ATP synthase may be reduced due to the decreased oxidative phosphorylation. As expected, Feichtinger et al. (2018) found significantly reduced levels of ATP Synthase Subunit ATP5F1A was correlated with earlier-onset prostate cancer. Shin et al. (2005) used the two-dimensional gel electrophoresis data to demonstrate the expression and activity of ATP synthase were lower in 5-FU-resistant cells compared with parent cancer cells. Inhibition of ATP synthase by oligomycin A or siRNA transfection strongly antagonized 5-FU-induced suppression of cell proliferation and increased cell viability. Similarly, the study of Song et al. (2018) indicated loss of ATP5H conferred a stem-like, invasive phenotype to tumor cells as well as multimodal resistance to immunotherapy, chemotherapy and radiotherapy. Also, reduced levels of ATP synthase were associated with poor prognosis in cancer patients (Song et al., 2018). In agreement with these studies, we also found ATP5H could interact with SOD1 and was downregulated in malignant transformation after SiNPs treatment.
In accordance with our hypothesis that the increased oxidative stress and aerobic glycolysis may accelerate the cell proliferation, we also found translational initiation was abnormally activated in human lung epithelial cells after SiNPs exposure, which may, on one hand, promote cell division and on the other hand, maintains cell survival (Pyronnet & Sonenberg, 2001). Translation process requires the protein complex known as eukaryotic initiation factor 4F (eIF4F) which consists of cap-binding protein eIF4E,   (Kashofer et al., 2017). Nonsense-mediated mRNA decay (NMD) represents a surveillance mechanism that eliminates transcripts with nonsense mutations and prevents cancer development. In contrast, inhibition of NMD may result in the initiation of cancer (Cao et al., 2017). Recent studies showed PABPC1 can interact with eIF4G to inhibit NMD (Fatscher et al., 2014;Peixeiro et al., 2012). Thus, the upregulation of PABPC1 may be one underlying reason for tumor formation. This hypothesis has been demonstrated in the gastric carcinoma (Zhu et al., 2015) and hepatocellular carcinoma (Zhang et al., 2015b) samples. As expected, our present study also demonstrated SiNPs may induce the tumorigenesis of Beas-2B cells by upregulating EIF4G2 and PABPC1. These two genes could interact with each other. More interestingly, our study showed the downregulation of ATP5H and SOD1 may be resulted from the upregulation of their common upstream miR-3648/572/661, while the downregulated miR-4521 may lead to the upregulation of the transcription of EIF4G2 and PABPC1. Although all these were the potential mechanisms firstly obtained for the carcinogenicity of amorphous SiNPs, recent studies on the miR-3648/572/661 and miR-4521 may indirectly demonstrate their importance in cancer development. For example, Rashid et al. (2017) demonstrated overexpression of miR-3648 promoted the growth of HeLa cells, while opposite results were obtained when miR-3648 was inhibited by antagomir. The mechanism of miR-3648 for cancer progression was to suppress the expression of a tumor suppressor gene adenomatous polyposis coli 2. miR-572 was also found to be highly expressed in human ovarian cancer tissues and cell lines. Ectopic overexpression of miR-572 promoted ovarian cancer cell proliferation and cell cycle progression in vitro and tumorigenicity in vivo by inhibiting its direct target suppressor of cytokine signaling 1 and cyclin-dependent kinase inhibitor 1A (p21KIP). Kaplan-Meir analysis indicated that high level expression of miR-572 was associated with poorer overall survival (Zhang et al., 2015a). miR-572 also can induce proliferation, invasion and inhibit apoptosis of nasopharyngeal carcinoma cells by targeting protein phosphatase 2 regulatory subunit Bgamma (Yan et al., 2017). miR-661 was observed to be upregulated in non-small cell lung cancer (NSCLC) tissues as compared to paired adjacent tissues. Furthermore, miR-661 promoted proliferation, migration and metastasis of NSCLC by regulating RB1 and mediating epithelial-mesenchymal transition process in NSCLC . Yamaguchi et al. (2017) observed miR-4521 was downregulated in chemotherapy-resistant renal cell carcinoma. However, the studies on miR-3648/572/661 and miR-4521 remains rare and their regulatory relationship with our predicted target genes have not been investigated. However, there are some limitations in this study. First, this is a preliminary study to identify the potential carcinogenic mechanisms of amorphous SiNPs. Additional wet experiments are necessary to confirm the expressions of identified genes and miRNAs (i.e., PCR), their interaction relationships (i.e., dual-luciferase, overexpression or knockout in vitro and in vivo) as well as the influence on the cell proliferation, apoptosis, migration and invasion. Second, our obtained miRNAs from the transcriptome array may be pre-miRNAs. Thus, further miRNA microarray (such as Affymetrix GeneChip miRNA v4) or sequencing analysis is essential to screen more crucial mature miRNAs.

CONCLUSION
Our findings reveal amorphous SiNPs may exert a carcinogenic effect by targeted regulating of miR-3648/572/661 and miR-4521 followed by influencing their downstream target genes (ATP5H/SOD1 and EIF4G2/PABPC1, respectively). These target genes may be involved in cancer development by promoting oxidative stress, translational initiation, while also inhibiting oxidative phosphorylation and NMD. Accordingly, the four miRNAs and their target genes may be underlying biomarkers for prediction of carcinogenesis when exposed to SiNPs and potentially other targets for cancer treatment.