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.
Materials and Methods
The microarray data were extracted from the gene expression omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE82062 (Guo et al., 2017), in which three repeats of human lung epithelial cells, Beas-2B continuously exposed to five μg/mL amorphous SiNPs for 40 passages (BeasSiNPs-P40 group; Supplemental Information 1.1, 1.2, 1.3) and three repeats of passage-matched control Beas-2B cells (Beas-P40 group) (Supplemental Information 1.4, 1.5, 1.6) were available for the analysis.
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.
miRNA-target genes interaction prediction
Because there were some pre-miRNAs included in the DEGs, the pre-miRNA–mRNA interaction was also investigated to explore the potential regulatory mechanisms of mRNAs. mRNA targets of differentially expressed miRNAs (DE-miRNAs) were predicted using the miRWalk database (v2.0; http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) (Dweep & Gretz, 2015), which were then overlapped with the DEGs to obtain DE-miRNA–DEGs interaction relationships. Then, the DE-miRNA–DEG interaction network was constructed and visualized using Cytoscape software (Kohl, Wiese & Warscheid, 2011).
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).
|Down||GO:0006334∼nucleosome assembly||2.32E-16||HIST1H2BB, HIST1H3J, HIST1H4L, HIST1H1E, HIST1H1C, HIST1H2BF, HIST1H2BO, HIST1H2BK, HIST1H4B, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H4C, HIST1H3D, HIST1H4D, HIST1H3F|
|GO:0051290∼protein heterotetramerization||5.54E-15||HIST1H3J, HIST1H4L, HIST1H4B, XRCC6, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, ANXA2|
|GO:0032200∼telomere organization||1.51E-13||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0045815∼positive regulation of gene expression, epigenetic||5.82E-13||HIST1H3J, HIST1H4L, POLR2L, POLR2K, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0006335∼DNA replication-dependent nucleosome assembly||8.76E-13||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0000183∼chromatin silencing at rDNA||3.77E-12||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0045814∼negative regulation of gene expression, epigenetic||7.03E-11||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0031047∼gene silencing by RNA||4.04E-10||HIST1H3J, HIST1H4L, POLR2L, POLR2K, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO:0044267∼cellular protein metabolic process||7.84E-10||HIST1H3J, HIST1H4L, HIST1H4B, TGFBI, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, B2M|
|GO:0006614∼SRP-dependent cotranslational protein targeting to membrane||2.29E-08||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO:0098609∼cell-cell adhesion||7.12E-08||HIST1H3J, CHMP5, EEF2, RPL24, PDLIM1, SFN, ANXA2, KRT18, HIST1H3A, HIST1H3B, HIST1H3D, CNN2, HIST1H3F, ENO1|
|GO:0019083∼viral transcription||1.06E-07||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO:0000184∼nuclear-transcribed mRNA catabolic process, nonsense-mediated decay||1.80E-07||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO:0006413∼translational initiation||5.98E-07||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO:0060968∼regulation of gene silencing||7.87E-07||HIST1H3J, HIST1H3A, HIST1H3B, HIST1H3D, HIST1H3F|
|GO:0006364∼rRNA processing||3.14E-06||RPS28, RPS29, RPL41, RRP36, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO:0006303∼double-strand break repair via nonhomologous end joining||5.72E-06||HIST1H4L, HIST1H4B, XRCC6, HIST1H4E, BABAM1, HIST1H4C, HIST1H4D|
|GO:0045653∼negative regulation of megakaryocyte differentiation||7.02E-06||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO:0034080∼CENP-A containing nucleosome assembly||1.33E-05||HIST1H4L, HIST1H4B, HIST1H4E, CENPW, HIST1H4C, HIST1H4D|
|GO:0006412∼translation||1.37E-05||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, MRPS7, RPS5|
|GO:0006342∼chromatin silencing||1.67E-05||HIST1H2AB, HIST2H2AB, HIST2H2AA4, HIST2H2AC, HIST1H2AE, HIST1H2AM|
|GO:0016233∼telomere capping||1.98E-05||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO:0006336∼DNA replication-independent nucleosome assembly||3.28E-05||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO:0006352∼DNA-templated transcription, initiation||1.22E-04||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO:1904837∼beta-catenin-TCF complex assembly||2.47E-04||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO:0000028∼ribosomal small subunit assembly||3.17E-04||RPS28, RPS15, MRPS7, RPS5|
|GO:0006123∼mitochondrial electron transport, cytochrome c to oxygen||3.71E-04||COX7B, COX8A, COX6B1, COX6A1|
|GO:1902600∼hydrogen ion transmembrane transport||9.45E-04||UQCR10, COX7B, COX8A, COX6B1, COX6A1|
|GO:0043154∼negative regulation of cysteine-type endopeptidase activity involved in apoptotic process||1.50E-03||PRDX5, SFN, THBS1, IFI6, DHCR24|
|GO:0019731∼antibacterial humoral response||3.81E-03||HIST1H2BK, HIST1H2BF, HIST1H2BI, B2M|
|GO:0042274∼ribosomal small subunit biogenesis||5.69E-03||RPS28, RRP36, RPS15|
|GO:0002576∼platelet degranulation||6.37E-03||PSAP, THBS1, SOD1, SRGN, FN1|
|GO:0002227∼innate immune response in mucosa||1.36E-02||HIST1H2BK, HIST1H2BF, HIST1H2BI|
|GO:0043086∼negative regulation of catalytic activity||1.66E-02||OAZ1, NQO1, ANXA2, ANXA2P2|
|GO:0050434∼positive regulation of viral transcription||1.81E-02||POLR2L, POLR2K, SUPT4H1|
|GO:0050830∼defense response to Gram-positive bacterium||2.30E-02||HIST1H2BK, HIST1H2BF, HIST1H2BI, B2M|
|GO:0006368∼transcription elongation from RNA polymerase II promoter||2.37E-02||POLR2L, POLR2K, ELOF1, SUPT4H1|
|GO:0033490∼cholesterol biosynthetic process via lathosterol||2.83E-03||DHCR7, DHCR24|
|GO:0033489∼cholesterol biosynthetic process via desmosterol||2.83E-02||DHCR7, DHCR24|
|GO:0007568∼aging||3.07E-02||CCL2, EEF2, SOD1, NQO1, APEX1|
|GO:0000302∼response to reactive oxygen species||3.16E-02||TXN, PRDX5, SOD1|
|GO:0001895∼retina homeostasis||3.31E-02||ACTG1, SOD1, B2M|
|GO:0045471∼response to ethanol||3.94E-02||CCL2, EEF2, SOD1, NQO1|
|GO:0006356∼regulation of transcription from RNA polymerase I promoter||4.21E-02||POLR2L, POLR2K|
|GO:0098532∼histone H3-K27 trimethylation||4.21E-02||HIST1H1E, HIST1H1C|
|GO:0007596∼blood coagulation||4.30E-02||HIST1H3J, HIST1H3A, HIST1H3B, HIST1H3D, HIST1H3F|
|GO:0009615∼response to virus||4.43E-02||IFITM2, IFITM3, RPS15A, ENO1|
|GO:0043066∼negative regulation of apoptotic process||4.48E-02||KRT18, PLK2, PLK1, TMBIM4, PRDX5, THBS1, NQO1, DHCR24|
|GO:1900121∼negative regulation of receptor binding||4.90E-02||ANXA2, B2M|
|Up||GO:0045727∼positive regulation of translation||1.57E-02||EIF4G2, PABPC1|
|GO:0006413∼translational initiation||4.01E-02||EIF4G2, PABPC1|
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).
|Down||hsa05322: Systemic lupus erythematosus||2.06E-16||HIST1H2AB, HIST1H2BB, HIST1H3J, HIST1H4L, HIST2H2AA4, HIST1H2BF, HIST1H2AE, HIST1H2BO, HIST2H2AB, HIST1H2BK, HIST1H4B, HIST2H2AC, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, HIST1H2AM|
|hsa05034: Alcoholism||5.44E-14||HIST1H2AB, HIST1H2BB, HIST1H3J, HIST1H4L, HIST2H2AA4, HIST1H2BF, HIST1H2AE, HIST1H2BO, HIST2H2AB, HIST1H2BK, HIST1H4B, HIST2H2AC, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, HIST1H2AM|
|hsa03010: Ribosome||9.39E-06||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, MRPS7, RPS5|
|hsa05016: Huntington’s disease||1.80E-04||UQCR10, NDUFS5, POLR2L, POLR2K, COX7B, COX8A, COX6B1, COX6A1, SOD1, NDUFA1, ATP5H|
|hsa00190: Oxidative phosphorylation||3.14E-04||UQCR10, NDUFS5, ATP6V0E1, COX7B, COX8A, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa05203: Viral carcinogenesis||1.32E-03||HIST1H2BO, HIST1H2BB, HIST1H4L, HIST1H2BK, HIST1H2BF, HIST1H4B, HIST1H2BI, HIST1H4E, HIST1H4C, HIST1H4D|
|hsa05010: Alzheimer’s disease||1.47E-03||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1, PSENEN, NDUFA1, ATP5H|
|hsa05012: Parkinson’s disease||2.41E-03||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa04932: Non-alcoholic fatty liver disease (NAFLD)||1.36E-02||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1, NDUFA1|
|hsa04260: Cardiac muscle contraction||1.62E-02||UQCR10, COX7B, COX8A, COX6B1, COX6A1|
|hsa03008: Ribosome biogenesis in eukaryotes||2.64E-02||SNORD3A, SNORD3C, SNORD3B-1, SNORD3B-2, NOP10|
h_eif4Pathway: Regulation of eIF4e and p70 S6 Kinase
KEGG, Kyoto encyclopedia of genes and genomes.
PPI network construction
After mapping the DEGs into the PPI data, a DEG-related PPI network was constructed (Fig. 2A), including 92 nodes (86 downregulated and six upregulated) and 291 edges (interaction relationships, such as EIF4G2-PABPC1, SOD1-ATP5H) (Supplemental Information 4). After calculating the degree, PABPC1 (degree = 15), ATP5H (degree = 11) and SOD1 (degree = 8) were suggested to be hub genes (Fig. 2B).
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).
|1||hsa03010: Ribosome||8.25E-17||RPS28, RPL41, RPS29, RPS15, RPL26, RPS4Y1, RPS13, RPS15A, RPL24, RPS5|
|2||hsa05322: Systemic lupus erythematosus||9.23E-09||HIST1H2BO, HIST1H2BB, HIST1H2BK, HIST2H2AC, HIST1H3B, HIST1H2AM, HIST3H3|
|3||hsa05012: Parkinson’s disease||5.44E-12||NDUFS5, UQCR10, COX8A, COX7B, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa00190: Oxidative phosphorylation||6.08E-12||NDUFS5, UQCR10, COX8A, COX7B, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa05010: Alzheimer’s disease||3.06E-11||NDUFS5, UQCR10, COX8A, COX7B, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa05016: Huntington’s disease||6.21E-11||NDUFS5, UQCR10, COX8A, COX7B, COX6B1, COX6A1, NDUFA1, ATP5H|
|hsa04260: Cardiac muscle contraction||1.73E-06||UQCR10, COX8A, COX7B, COX6B1, COX6A1|
DEGs, differentially expressed genes; KEGG, Kyoto encyclopedia of genes and genomes.
miRNA–mRNA interaction network
A total of 18,896 target genes were predicted for the 41 DE-miRNAs by using the miRWalk database, 130 of which (including 124 downregulated and six upregulated) were found to be shared with the DEGs. Then, a miRNA–mRNA interaction network was constructed by using the 130 target genes and their corresponding 24 miRNAs, resulting in 1,615 and 17 interrelationship pairs for the 18 upregulated miRNAs and six downregulated miRNAs, respectively (Fig. 4).
Function enrichment analysis also demonstrated these 130 target genes were involved in GO terms of response to reactive oxygen species (hsa-miR-3155b/3648/4522/4523/4533/4634/4734/572/661/933/1470-SOD1), and translational initiation (miR-5095/548ai/622/4521-EIF4G2, miR-4521-PABPC1) (Table 4) and KEGG pathways of Oxidative phosphorylation (hsa-miR-1976/3155/3648/4461/4522/4523/4533/4634/4734/572/661/933-ATP5H) (Table 5).
|GO: 0006334∼nucleosome assembly||1.02E-17||HIST1H2BB, HIST1H3J, HIST1H4L, HIST1H1E, HIST1H1C, HIST1H2BF, HIST1H2BO, HIST1H2BK, HIST1H4B, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, HIST3H3|
|GO: 0051290∼protein heterotetramerization||1.02E-16||HIST1H3J, HIST1H4L, XRCC6, ANXA2, HIST1H4B, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H4C, HIST1H3D, HIST1H4D, HIST3H3, HIST1H3F|
|GO: 0032200∼telomere organization||1.63E-13||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO: 0045815∼positive regulation of gene expression, epigenetic||6.38E-13||HIST1H3J, HIST1H4L, POLR2L, POLR2K, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO: 0006335∼DNA replication-dependent nucleosome assembly||9.45E-13||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO: 0000183∼chromatin silencing at rDNA||4.07E-12||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO: 0031047∼gene silencing by RNA||2.43E-11||HIST1H3J, HIST1H4L, POLR2L, POLR2K, HIST1H4B, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H4C, HIST1H3D, HIST1H4D, PABPC1, HIST1H3F|
|GO: 0045814∼negative regulation of gene expression, epigenetic||7.57E-11||HIST1H3J, HIST1H4L, HIST1H4B, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F|
|GO: 0044267∼cellular protein metabolic process||8.57E-10||HIST1H3J, HIST1H4L, HIST1H4B, TGFBI, HIST1H3A, HIST1H3B, HIST1H4E, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, B2M|
|GO: 0006413∼translational initiation||4.23E-09||EIF4G2, RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, PABPC1, RPS5|
|GO: 0098609∼cell-cell adhesion||9.43E-09||HIST1H3J, CHMP5, PDLIM1, EEF2, RPL24, SFN, ANXA2, EIF4G2, KRT18, HIST1H3A, HIST1H3B, HIST1H3D, CNN2, HIST1H3F, ENO1|
|GO: 0000184∼nuclear-transcribed mRNA catabolic process, nonsense-mediated decay||1.42E-08||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, PABPC1, RPS5|
|GO: 0006614∼SRP-dependent cotranslational protein targeting to membrane||2.47E-08||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO: 0019083∼viral transcription||1.14E-07||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO: 0006303∼double-strand break repair via nonhomologous end joining||3.36E-07||HIST1H4L, HIST1H4B, XRCC6, HIST1H4E, BABAM1, HIST1H4C, HIST1H4D, HIST3H3|
|GO: 0016233∼telomere capping||5.42E-07||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D, HIST3H3|
|GO: 0060968∼regulation of gene silencing||8.14E-07||HIST1H3J, HIST1H3A, HIST1H3B, HIST1H3D, HIST1H3F|
|GO: 0006364∼rRNA processing||3.38E-06||RPS28, RPS29, RPL41, RRP36, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, RPS5|
|GO: 0045653∼negative regulation of megakaryocyte differentiation||7.26E-06||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO: 0034080∼CENP-A containing nucleosome assembly||1.38E-05||HIST1H4L, HIST1H4B, HIST1H4E, CENPW, HIST1H4C, HIST1H4D|
|GO: 0006412∼translation||1.48E-05||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, MRPS7, RPS5|
|GO: 0006342∼chromatin silencing||1.73E-05||HIST1H2AB, HIST2H2AB, HIST2H2AA4, HIST2H2AC, HIST1H2AE, HIST1H2AM|
|GO: 0006336∼DNA replication-independent nucleosome assembly||3.39E-05||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO: 0006352∼DNA-templated transcription, initiation||1.26E-04||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO: 1904837∼beta-catenin-TCF complex assembly||2.55E-04||HIST1H4L, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D|
|GO: 0000028∼ribosomal small subunit assembly||3.25E-04||RPS28, RPS15, MRPS7, RPS5|
|GO: 0006123∼mitochondrial electron transport, cytochrome c to oxygen||3.80E-04||COX7B, COX8A, COX6B1, COX6A1|
|GO: 0045727∼positive regulation of translation||5.72E-04||EIF4G2, EEF2, QARS, PABPC1, THBS1|
|GO: 1902600∼hydrogen ion transmembrane transport||9.75E-04||UQCR10, COX7B, COX8A, COX6B1, COX6A1|
|GO: 0043154∼negative regulation of cysteine-type endopeptidase activity involved in apoptotic process||1.55E-03||PRDX5, SFN, THBS1, IFI6, DHCR24|
|GO: 0019731∼antibacterial humoral response||3.90E-03||HIST1H2BK, HIST1H2BF, HIST1H2BI, B2M|
|GO: 0042274∼ribosomal small subunit biogenesis||5.79E-03||RPS28, RRP36, RPS15|
|GO: 0002227∼innate immune response in mucosa||1.39E-02||HIST1H2BK, HIST1H2BF, HIST1H2BI|
|GO: 0050434∼positive regulation of viral transcription||1.84E-02||POLR2L, POLR2K, SUPT4H1|
|GO: 0050830∼defense response to Gram-positive bacterium||2.35E-02||HIST1H2BK, HIST1H2BF, HIST1H2BI, B2M|
|GO: 0006368∼transcription elongation from RNA polymerase II promoter||2.42E-02||POLR2L, POLR2K, ELOF1, SUPT4H1|
|GO: 0033489∼cholesterol biosynthetic process via desmosterol||2.85E-02||DHCR7, DHCR24|
|GO: 0033490∼cholesterol biosynthetic process via lathosterol||2.85E-02||DHCR7, DHCR24|
|GO: 0007568∼aging||3.15E-02||CCL2, EEF2, SOD1, NQO1, APEX1|
|GO: 0000302∼response to reactive oxygen species||3.21E-02||TXN, PRDX5, SOD1|
|GO: 0001895∼retina homeostasis||3.36E-02||ACTG1, SOD1, B2M|
|GO: 0002576∼platelet degranulation||3.84E-02||PSAP, THBS1, SOD1, SRGN|
|GO: 0043488∼regulation of mRNA stability||3.84E-02||PSMB3, HSPA1B, PABPC1, APEX1|
|GO: 0045471∼response to ethanol||4.03E-02||CCL2, EEF2, SOD1, NQO1|
|GO: 0006356∼regulation of transcription from RNA polymerase I promoter||4.25E-02||POLR2L, POLR2K|
|GO: 0098532∼histone H3-K27 trimethylation||4.25E-02||HIST1H1E, HIST1H1C|
|GO: 0007596∼blood coagulation||4.41E-02||HIST1H3J, HIST1H3A, HIST1H3B, HIST1H3D, HIST1H3F|
|GO: 0009615∼response to virus||4.52E-02||IFITM2, IFITM3, RPS15A, ENO1|
|GO: 0043066∼negative regulation of apoptotic process||4.65E-02||KRT18, PLK2, PLK1, TMBIM4, PRDX5, THBS1, NQO1, DHCR24|
|GO: 1900121∼negative regulation of receptor binding||4.94E-02||ANXA2, B2M|
|hsa05322: Systemic lupus erythematosus||5.48E-18||HIST1H2AB, HIST1H2BB, HIST1H3J, HIST1H4L, HIST2H2AA4, HIST1H2BF, HIST1H2AE, HIST1H2BO, HIST2H2AB, HIST1H2BK, HIST1H4B, HIST2H2AC, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, HIST1H2AM, HIST3H3|
|hsa05034: Alcoholism||2.00E-15||HIST1H2AB, HIST1H2BB, HIST1H3J, HIST1H4L, HIST2H2AA4, HIST1H2BF, HIST1H2AE, HIST1H2BO, HIST2H2AB, HIST1H2BK, HIST1H4B, HIST2H2AC, HIST1H2BI, HIST1H3A, HIST1H4E, HIST1H3B, HIST1H3D, HIST1H4C, HIST1H4D, HIST1H3F, HIST1H2AM, HIST3H3|
|hsa03010: Ribosome||6.90E-06||RPS28, RPS29, RPL41, RPS15, RPL26, RPS13, RPS15A, RPS4Y1, RPL24, MRPS7, RPS5|
|hsa05016: Huntington’s disease||6.46E-04||UQCR10, NDUFS5, POLR2L, POLR2K, COX7B, COX8A, COX6B1, COX6A1, SOD1, ATP5H|
|hsa05203: Viral carcinogenesis||1.03E-03||HIST1H2BO, HIST1H2BB, HIST1H4L, HIST1H2BK, HIST1H2BF, HIST1H4B, HIST1H2BI, HIST1H4E, HIST1H4C, HIST1H4D|
|hsa00190: Oxidative phosphorylation||1.36E-03||UQCR10, NDUFS5, ATP6V0E1, COX7B, COX8A, COX6B1, COX6A1, ATP5H|
|hsa05010: Alzheimer’s disease||5.06E-03||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1, PSENEN, ATP5H|
|hsa05012: Parkinson’s disease||8.75E-03||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1, ATP5H|
|hsa04260: Cardiac muscle contraction||1.45E-02||UQCR10, COX7B, COX8A, COX6B1, COX6A1|
|hsa04932: Non-alcoholic fatty liver disease (NAFLD)||4.12E-02||UQCR10, NDUFS5, COX7B, COX8A, COX6B1, COX6A1|
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 SiNPs-mediated 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, scaffolding protein eIF4G and ATP-dependent RNA helicase eIF4A (Gingras, Raught & Sonenberg, 1999). Hereby, the expression of translation protein related genes may be upregulated in malignant transformation, which have been confirmed in several cancers because downregulation of eIF4GII was reported to decrease cell proliferation, but induces cellular senescence (Emmrich et al., 2016; Xie et al., 2017). During translation process in eukaryotic cells, nonsense mutations (premature stop codon) may be present, which disrupts production of full-length, functional proteins (such as tumor suppressor gene p53) and thus may induce the development of various diseases, including cancers (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 (Liu et al., 2017). 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.
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.