Systematic Identification of Genomic Markers for Guiding Iron Oxide Nanoparticles in Cervical Cancer Based on Translational Bioinformatics

Introduction

Conventional drug administration in cancer treatment is mainly nonselective, and many anti- cancer drugs are expelled out of cancer cells or cleared rapidly from the body, reducing chemotherapeutics efficiency and resulting in off-target side effects.^1,2 A promising strategy to deal with this challenge is the nano-drug delivery platform, which carries anti-cancer drugs and enhances their local concentration at the tumor site.^3–7 Magnetic iron oxide nanoparticles (MNPs) occupy a privileged position among other nano-drug delivery platforms,^8,9 given their properties in imaging contrast and magnetic drug targeting.^10–12 Although the MNP is designed as a carrier of the drugs and exhibits excellent biocompatibilities and low toxicities,^13–15 the effects of the MNP itself should not be ignored. However, currently, much effort is devoted to addressing several important issues such as MNP fabrication, functionalized modification, drug conjugation, and delivery selectivity,^16–19 and the focus on elucidating the molecular pathogenesis of MNP in cancer cells is lacking.

Magnetic iron oxide nanoparticles have been successfully applied in the in vivo delivery of doxorubicin to high-grade intramedullary spinal cord tumors in rats.²⁰ In in vitro models of hepatocellular carcinoma and prostate cancer, MNP delivery systems loaded with quercetin or docetaxel also exhibited potent anti-cancer efficacy.^21,22 However, it should be noted that a series of cascade reactions are supposed to be triggered by the MNP itself, which might impact the effects of anti-cancer drugs at the molecular level. However, despite the growing evidence for interaction among MNPs and bio-macromolecules,^23–25 information on the molecular mechanism of MNPs inside cancer cells and their interaction with anti-cancer drugs is very limited, which might be extremely important in clinic applications.

The systems pharmacology approach has been widely used in monitoring the effects of drugs and screening potential drug candidates.²⁶ This approach develops a global understanding of the drug’s pathophysiology and does not focus on the single interaction between a drug and its target or pathway. The Food and Drug Administration has recognized that gene expression profiling is a promising approach for advancing the development of medical products.²⁷ Owing to the unique property of MNPs, their molecular effects cannot be predicted, neither from the known chemical agents nor from iron metabolism. Therefore, this study rationalized that a non-hypothesis-driven systems pharmacology approach using whole genome expression profiling is likely to provide a broad insight into the biological effects of MNPs.

In the last few years, the progress in biological computational methods has made it possible to analyze genetic aberrations and molecular mechanisms and provide promising tools in medical oncology with great clinical applications.^28–30 However, these approaches have not been fully used in the field of nano-chemotherapy. In this study, bioinformatics data mining coupled with clinical data analysis were employed to achieve the systems pharmacology approach. The properties of MNPs were evaluated in cervical cancer cell lines, and four key genes relevant for MNP pharmacology, clinical features, somatic mutation, and immune infiltration in cervical cancer were identified. This study also showed that the four key genes were significantly correlated with the lysosome gene set. Potential targeting drugs of key genes were also presented.

This study validated the properties of MNPs in two cervical cancer cell lines, employed biological computational methods and identified protein tyrosine kinase 2 beta (PTK2B), liprin-alpha-4 (PPFIA4), mothers against decapentaplegic homolog 7 (SMAD7), and interleukin (IL) 1B as key genes of MNP onco-pharmacology in cervical cancer. Moreover, potential target drugs were also proposed for further investigation. This study can increase our understanding of the pharmacologic mechanism of MNP treatment in tumors, which is necessary when designing an effective MNP-based drug delivery system.

Materials and Methods

Chemicals

Chemicals, including iron (II) chloride tetrahydrate (FeCl₂.4H₂O), iron (III) chloride hexahydrate (FeCl₃.6H₂O), ethanol (EtOH), and ammonium hydroxide solution (NH₄OH), were obtained from Merck (Darmstadt, Germany).

The bare Fe₃O₄ nanoparticles (Fe₃O₄ NPs) were prepared via the co-precipitation method.³¹ Briefly, FeCl₃·6H₂O (0.02 mol/5.38 g) was mixed with FeCl₂·4H₂O (0.01 mol/1.98 g) in deionized water (100 mL) under a flow of nitrogen. After 10 min, NH₄OH (25%, 10 mL) was added dropwise into the above reaction with vigorous stirring at 80 °C for 30 min. The resulting black solution was filtered and washed with deionized water (100 mL) three times. The Fe₃O₄ NPs loaded with coumarin 6 (C6) used in this study were synthesized according to a previously published technique.³² For the cellular study, Fe₃O₄ NPs were sterilized with ethanol, and a bath sonicator was used to disperse the products into a colloid which was stored at 4 °C for further experiments.

Cell Culture and Treatment

Human cervical cancer HeLa and SiHa cells were obtained from the China Academy of Chinese Medical Sciences (Beijing, China) and were cultured under ideal growth conditions (37 °C, 5% CO₂, and 95% air atmosphere) using RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc., USA) containing 10% fetal bovine serum (FBS, Gibco; Thermo Fisher Scientific, Inc., USA), 100 U/mL penicillin, and 100 µg/mL streptomycin. HeLa and SiHa cells were incubated in the presence of 0, 50, 100, 200, and 400 μg/mL MNP. Primarily, events under 100 μg/mL MNP exposure over 24 h were chosen for the following experiments.

MTT Assays

Cytotoxicity assessment was performed using the 3-[4,5-dimethylthiazol-2-yl]-2,5 diphenyl tetrazolium bromide (MTT) assay. HeLa and SiHa cells were plated in 96-well plates at a density of 1×10⁴ cells/well and cultured with RPMI-1640 medium. After 24 h of incubation, the cells were supplemented with MNP (10 μg) at a final concentration of 100 μg/mL and incubated for a further 24 h. Control cells were incubated in RPMI-1640 medium with 10% FBS, and H₂O₂ (0.5 mM) was employed as the positive control. After incubation, 10 μL of 5 mg/mL MTT in 90 μL of RPMI-1640 medium was added to each well. After 4 h, the formazan crystals were dissolved by adding 150 μL of dimethylsulfoxide. The absorption of the formazan solution was measured using a Fluostar Omega microplate reader (BMG Labtech, Germany) at a wavelength of 570 nm.

Cell Fluorescence Imaging

Magnetic iron oxide nanoparticles were labeled with C6 for cell imaging, and LysoTracker Red (LTR) was utilized to indicate lysosomes. Bright yellow fluorescence indicated the colocalization of the MNP and lysosome in HeLa cells. Quantitative colocalization analysis was performed using ImageJ software 1.52a (http://rsb.info.nih.gov/ij/). The LTR (Beyotime, China) probes are fluorescent acidotropic probes used for labeling and tracking acidic organelles in live cells; the variance in fluorescence intensity may also imply lysosomal pH changes. Acridine orange (AO; Merck, Germany), a lysosomotropic weak base that accumulates at high concentrations in intact lysosomes because of proton trapping was also used as the shift of granular red fluorescence to diffused green fluorescence may monitor either lysosomal leakage or change in lysosomal pH. Briefly, HeLa and SiHa cells seeded (5 × 10⁴ cells/well) on coverslips placed in 24-well plates were treated with or without 100 μg/mL MNP (24 h), followed by loading with 5 μg/mL AO for 30 min and Hoechst 33,342 for 5 min, or 100 nM LTR for 30 min under specific growth conditions (37 °C, 5% CO₂). Then, the slides were rinsed twice with warm (37 °C) PBS (1 mL) and observed under a Zeiss Axio Imager Z2 (Carl Zeiss, Germany).³³

Cathepsin-B (CTSB) Activity Assay

The lysosomal function was evaluated by detecting CTSB-specific activity using Z-FR-AMC (R&D Systems, USA), a synthetic fluorogenic peptide substrate for CTSB. The cathepsin-B-specific activities of the HeLa and SiHa cells treated with MNPs were measured in an assay buffer containing 200 mM NaOAc, 0.1% Brij 35, 0.1% Triton X-100, and 1 mM EDTA at pH 6.7. Briefly, 9 μL of cellular extraction and 1 μL of the CTSB substrate Z-FR-AMC (100 mg/mL) were added to 90 μL of the assay buffer. Then, the samples were incubated in each well of a black 96-well plate at 37 °C for 1 h. The relative fluorescence units, corresponding to the CTSB cleavage of the substrates, were measured at an excitation wavelength of 355 nm and emission wavelength of 460 nm in a Fluostar Omega microplate reader (BMG Labtech, Germany).

Apoptosis Analysis by Flow Cytometry

Apoptosis was determined by Annexin-V and propidium iodide (PI) staining, according to the manufacturer’s protocol (556,547; BD Biosciences). Briefly, HeLa and SiHa cells were seeded (2 × 10⁵ cells/well) into six-well plates. After treatments, the cells were collected and washed and then stained with Annexin-V and PI. Cells were quantified using a guava easyCyte flow cytometer (Millipore).

Quantitative PCR

Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA). Then, 2 μg of total RNA was used to synthesize cDNA using the Omniscript Reverse Transcription kit (TaqMan, Applied Biosystems, USA). Quantitative reverse transcription PCR (qRT-PCR) was performed to estimate the mRNA levels of PPFIA4, PTK2B, IL1B, and SMAD7 using SYBR SuperMix (TransGen Biotech, CA, USA). All reactions were performed in duplicate on a Bio-Rad CFX96 Touch system (Bio-Rad, Hercules, CA, USA). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal control, and the relative expression level for each gene was calculated using the 2^−ΔΔCt relative quantification method and normalized to that of endogenous housekeeping gene. All reactions were performed using the following cycling parameters: initial denaturation at 94 °C for 30s; followed by 45 cycles of 94 °C for 5s, 60 °C for 15s, and 72 °C for 10s. The primers used were as follows: PTK2B gene, 5′-GAAGTCCGATGAGATCCACTGG-3′ (sense) and 5′-CTTCTGGCAAGTAGCGGATTT-3′ (anti-sense); PPFIA4 gene, 5′-CTCTGCGGATGTTGTCTCCC-3′ (sense) and 5’-ATGCTGCCACTGGTTACACG-3′ (anti-sense); IL1B gene, 5′-ATGATGGCTTATTACAGTGGCAA-3′ (sense) and 5’-GTCGGAGATTCGTAGCTGGA-3′ (anti-sense); SMAD7 gene, 5′-GGACGCTGTTGGTACACAAG-3′ (sense) and 5′-GCTGCATAAACTCGTGGTCATTG-3′ (anti-sense); GAPDH gene, 5′-ACAACTTTGGTATCGTGGAAGG-3’ (sense) and 5′-GCCATCACGCCACAGTTTC-3′ (anti-sense); and PIK3R2 gene, 5′-AATGCAGCAAGGAATACCTGG-3′ (sense) and 5′-GCTCTCATGGATCTCGGCAA-3′ (anti-sense).

Microarray Data of HeLa Cells Treated with MNPs

The microarray datasets of GSE15248 were downloaded from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). These RNA profiles were based on the GPL6102 platforms (Illumina human-6 v2.0 expression beadchip) containing a total of three MNP-treated HeLa cell samples and three control samples.

Identification of Differentially Expressed Genes (DEGs) and Hub Genes

The downloaded matrix files were converted and processed in R (version 3.5.3) software with the “Affy” package to read the expression data for the CEL files. Gene differential expression analysis was performed using the “limma” package, after which the genes were sorted according to the log fold-change values. Finally, the genes with P < 0.05 and |logFC| > 0.5 were considered as DEGs. The “pheatmap” and “ggplot2” packages were used to draw the heatmap and volcano plot of the top ranked genes and DEGs, respectively.

The protein–protein interaction (PPI) network was also constructed by uploading all DEGs to the STRING database (https://string-db.org/). Then, the PPI network was further analyzed in Cytoscape software (www.cytoscape.org/). CytoHubba, a Cytoscape plug-in, sorts the genes by analyzing 12 parameters, including MCC, DMNC, MNC, DEGREE, EPC, BOTTLENECK, EcCentricity, CLOSENESS, RADIALITY, BETWEENNESS, STRESS, and ClusteringCoefficient. The top 50 genes ranked by each parameter were recorded. The genes sorted by eight or more parameters were selected as the hub genes being highly essential in the functional network.

Survival Analysis and Efficacy Evaluation

Expression data and survival information of 310 patients with cervical cancer from The Cancer Genome Atlas (TCGA) were obtained using the online software SangerBox (http://sangerbox.com). A Kaplan–Meier survival plot was constructed using the “survival”, “survminer”, and “ggplot2” packages in R software. Differential expression of the key genes among the clinical features was also plotted as a dot-box plot graph with the Kruskal–Wallis test.

Somatic Mutation of Key Genes

Genetic alterations of the key genes across the patients with cervical cancer were explored in the CBio Cancer Genomics Portal (http://www.cbioportal.org/). Oncoplots for the mutation spectrum visualization and Lollipop plots for the amino acid changes were processed with the “maftools” package in R software. Mutually exclusive or co-occurring sets of genes were also detected using somatic interaction function in the “maftools” package, which performs the pair-wise Fisher’s exact test to detect significant pairs of genes.

Correlation of Key Genes with Immune Infiltration

Systematic analysis of the immune infiltrates was conducted through Tumor Immune Estimation Resource (TIMER), which is a comprehensive resource with diverse cancer types (https://cistrome.shinyapps.io/timer/). The TIMER applies a previously published statistical algorithm, that is, deconvolution, to estimate the abundance of immune infiltrating levels from gene expression profiles. Based on 310 patients with cervical cancer from TCGA, the correlation of the expression of four key genes with the abundance of immune infiltrates was analyzed, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, dendritic cells, and tumor purity via gene modules. Finally, Spearman correlation and the estimated statistical significance were displayed as scatter-plots. The gene expression level was displayed with log2 RSEM.

Correlation of Key Genes with Lysosomes

Patients with cervical cancer were divided into two groups according to the median expression values of each key gene. First, Gene Set Enrichment Analysis (GSEA) was applied using the Java GSEA implementation to determine whether the defined set of lysosomal genes shows statistically significant concordant differences. The number of permutations was set as 1000, and c2.cp.kegg.v6.2.symbols.gmt was downloaded as the reference gene set. A false discovery rate q value < 0.05 and nominal P < 0.25 were set as the cut-off threshold. Second, correlations among key genes and individual lysosomal genes were calculated as Pearson’s correlation coefficients and processed using the “ggstatsplot” package in R software.

Identification of Potential Anti-Cancer Drugs Suitable for MNP Delivery

CBioPortal was used to explore the correlation among the four key genes, their most frequently altered neighbor genes, and anti-cancer drugs to screen new potential anti-cancer drugs suitable for MNP delivery.

Chemotherapy Resistance Gene Informatics Analysis

The Cancer Genome Atlas gene expression data were used for chemotherapy resistance analysis. As shown in Table 1, most of the patients received cisplatin, whereas the sample size of patients that received other drug regimens was too small to perform gene expression analysis. Finally, 106 patients that received cisplatin were selected for further analysis. Neoplasm cancer status was selected as the indicator of drug resistance. Expression data were transformed into transcripts per kilobase million before analysis in SangerBox software (Mugu Tech, Hangzhou, China). Finally, the genes with P < 0.05 and |logFC| > 0.5 were considered as DEGs. Hub genes were also explored with the PPI network and CytoHubba.

Table 1 Correlation Between the Clinicopathologic Variables and the Four Key Genes mRNA Expression in Cervical Cancer

Statistical Analyses

GraphPad Prism 8 statistical software (GraphPAD Software Inc., San Diego, CA, USA) was used for experimental data statistical analyses and image processing. Differences between groups were determined by the unpaired Student’s t-test.

Results

Cytotoxic Effect of MNP

The transmission electron microscopy images of MNPs in Figure 1A–C show a spherical morphology with a size of 20−100 nm. Magnetic iron oxide nanoparticles had a good magnetic response in the external magnetic field (Figure 1D). No cytotoxicity was observed even at the highest particle loading level (400 μg/mL) in the culture medium (Figure 1E, F, H, I). The effects of MNPs on apoptosis of cervical cancer cells were further investigated. HeLa (Figure 1G) and SiHa (Figure 1J) cells treated with MNP for 24 h and 72 h were stained with fluorescein isothiocyanate-labeled Annexin-V/PI and subsequently subjected to flow cytometry analysis. The results revealed that MNP (100 μg/mL) did not significantly increase the number of apoptotic cervical cancer cells.

Figure 1 Properties of MNP in cervical cancer cells. (A and B) TEM images of MNPs in 25,000× (A), 80,000× (B), 600,000× (C). (D) MNP exhibited room temperature superparamagnetism with good magnetic field responsivity towards the external magnet. Effect of the MNP on the viability of the HeLa (E and F) and SiHa (H and I) cells via the MTT assay. Cells were treated with MNP for 0–72h at a concentration range 0–400μg/mL, H2O2 (0.5mM) were employed as positive control. HeLa (G) and SiHa (J) cells treated with MNP (100μg/mL) for 24 and 72 h were stained with FITC-labeled annexinV/PI and subsequently subjected to flow cytometry analysis.

MNP Induced Inhibition of Lysosomal Function

To investigate the successful cell internalization of MNPs, colocalization experiments were conducted using fluorescence microscopy. The colocalization correlation coefficient was determined to be 0.95 (Figure 2A), indicating the internalization of MNP in the lysosome. A decrease in relative fluorescence units of Z-FR-AMC was observed, as a sign of decline in CTSB enzymatic activity (Figure 1B and C), indicating decreased lysosomal function. Lysosomotropic probes were applied to evaluate the effect of MNP on lysosomal pH in the HeLa and SiHa cells. LysoTracker Red exhibits red fluorescence in the lysosome in a pH-dependent manner, and a decrease in red fluorescence indicates an increase in lysosomal pH. The 24-h treatment with 100 μg/mL MNP neutralized the lysosomal pH, which was displayed as decreased red fluorescence of LTR compared with that of the control cells (Figure 1B and C). In addition, AO staining presents red fluorescence when it is accumulated within lysosomes owing to it being highly protonated but shows green fluorescence after diffusion into the cytosol. Thus, an increase in diffuse green fluorescence with a decrease in punctate red fluorescence suggests an elevated lysosomal pH. The 24-h treatment of HeLa and SiHa cells with 100 μg/mL MNP led to a decrease in the red fluorescence and increase in the green fluorescence with AO staining (Figure 1F–G).

Figure 2 Biological effects of MNP in lysosome of cervical cancer cells. (A) Colocalization of MNP-C6 with lysosomes in HeLa cells. Gray-value calculation and quantitative colocalization analysis was carried out using ImageJ software. CTSB activity in HeLa (B) and SiHa (C) cells was measured as RFU with CTSB substrate Z-FR-AMC. After MNP treatment, CTSB activity decreased significantly (p<0.001). Fluorescence value of lysotracker decreased obviously after MNP treatment in HeLa (D) and SiHa (E) cells. HeLa (F) and SiHa (G) cells treated with MNP exhibit a shift of granular red fluorescence to diffused green fluorescence of AO. Images at 1000× magnification were also shown (row 2, 4). Scale bar = 20 μm. ***p < 0.001. Abbreviations: MNP, magnetic iron oxides nanoparticle; C6, coumarin 6; CTSB, cathepsin B; RFU, relative fluorescence units; AO, acridine orange.

Identification of DEGs After MNP Treatment

Microarray datasets GSE15248 obtained from the Gene Expression Omnibus database were based on the GPL6102 platforms, which included three MNP-treated HeLa cell samples and three control samples. The results of the GSE15248 provided insights into the potential mechanisms underlying the effects of MNP on cervical cancer cell function. After normalization, 179 DEGs were identified, of which 67 were significantly upregulated genes and 112 were downregulated genes. The results are displayed in the heatmap and volcano plots for DEG expression (Supplementary Figure S1).

Identification of Key Genes Associated with the Survival Time and Clinical Features in Patients with Cervical Cancer

To evaluate the involvement of the above-identified hub genes in cervical cancer, the survival data of 310 patients from the TCGA were used for Kaplan–Meier survival analysis using the Log rank test. Four out of 44 hub genes demonstrated significant associations with overall survival (Figure 3). Patients with high expression of receptor-like protein tyrosine phosphatases-interacting protein (liprin) alpha 4 (PPFIA4) (P = 0.0031; Figure 3A), SMAD7 (P = 0.044; Figure 3B), and IL1B (P = 0.0085; Figure 3C) had a poor overall survival, whereas patients with high expression of PTK2B exhibited a better performance in survival time (P = 0.017; Figure 3D). The four genes were identified as key genes for further analysis. The relationship between the clinical features and the expression of the four key genes was also analyzed (Table 1).

Figure 3 Four key genes were associated with overall survival time and clinical features in cervical cancer patients. Survival analysis of PPFIA4 (A), SMAD7 (B), IL1B (C), and PTK2B (D) expression levels associating with overall survival time in 310 TCGA cervical cancer patients using Kaplan–Meier curves. The patients were divided into 2 groups, based on median gene expression value as a cut-off point. (E–G) The expression of PPFIA4 was highly associated with clinical stage with a p value of 2e-7 (E), and tumor (T) classification with a p value of 7.7e-06 (F). The expression of IL1B was significantly associated with tumor (T) classification with a p value of 0.0072 (G). The association between genes expression level and clinical features was estimated by Kruskal–Wallis test. (see also Figure S2).

Further subgroup analysis was conducted among patients with different clinical stage, histologic grade, and tumor (T) and lymph node (N) classification. High PPFIA4 expression correlated with a poor overall survival in patients with clinical stages I/II (P = 0.029), T1/T2 (P = 0.021), and N0 (P = 0.002). High SMAD7 expression correlated with a poor overall survival in patients with clinical stages I/II (P = 0.048), histologic grade G3/G4 (P = 0.034), and T1/T2 (P = 0.031). High IL1B expression correlated with a poor overall survival in patients with clinical stage I/II (P = 0.036), histologic grade G1/G2 (P = 0.027), T1/T2 (P = 0.028), and N0 (P = 0.033). High PTK2B expression correlated with a better overall survival in patients with histologic grade G3/G4 (P = 0.016) and T1/T2 (P = 0.013) (Supplementary Figure S2).

Somatic Mutations of Key Genes with Various Mutation Patterns in Patients with Cervical Cancer

Genetic alteration of the oncoplot in CBioPortal was used to show the alteration status of the four key genes in patients with cervical cancer. As shown in Figure 4, the four key genes were altered in a total of 71 (22.9%) out of 310 patients (Figure 4A), among which PPFIA4 was altered the most, with an alteration frequency of 11%. Mutational hotspots of PPFIA4 and PTK2B showed that PPFIA4 has two mutation types (S743F/R1007M/G1134D-Missense, X1087-Splice) (Figure 4B) and PTK2B has three mutation types (R159Q/E494K/R586H/R936L-Missense, X369-Splice, I597Hfs*2-Frameshift) (Figure 4C). The mutation spots of PPFIA4 were highly clustered, which is the feature of a proto-oncogene. The mutation spots were discrete in PTK2B, as expected in a tumor suppressor gene. Mutation interaction of the key genes with the top mutated genes was also explored (Figure 4D) in patients with cervical cancer. PPFIA4 and PTK2B had a significant co-occurrence of mutation, with seven of the top mutated genes (Figure 4E).

Figure 4 Somatic mutations of four key genes that were detected in the TCGA cervical cancer patients. (A) Oncoplot showed somatic mutation patterns of the four key genes. Each column in the figure represented an individual cancer sample. (B and C) Plots showed mutation distribution and the protein domains for the mutated proteins of PPFIA4 (B) and PTK2B (C). (D) Top 30 most recurrent mutated genes in cervical cancer patients were represented. The bar plot on the top panel represented the total number of mutations. (E) PPFIA4 and PTK2B were in significant co-occurrence of mutation with 7 of the top mutated genes (□). p < 0.05; *p < 0.01.

Key Genes Play a Role in Immune Infiltration in Cervical Cancer

Immune infiltration was investigated using the TIMER deconvolution algorithm, which is more reliable than other methods. Analysis showed that PTK2B expression level had significant positive correlations with infiltrating levels of neutrophils (r = 0.38, P = 6.04e-11), dendritic cells (r = 0.317, P = 7.2e-08), CD8⁺ T cells (r = 0.217, P = 3.01e-04), and CD4⁺ T cells (r = 0.18, P = 2.7e-03) in patients with cervical cancer (Figure 5A), whereas there were significant negative correlations of PPFIA4 expression with infiltrating levels of CD8⁺ T cells (r = –0.201, P = 7.96E-04) and neutrophils (r = –0.172, P = 4.09e-03) (Figure 5B). In addition, expression levels of PTK2B and PPFIA4 had significant negative (r = –0.283, P = 1.59E-06) and positive correlations (r = –0.20, P = 8.01E-04), respectively, with tumor purity in patients with cervical cancer (Figure 5A and B). Significant correlations between immune infiltration and the expression levels of SMAD7 and IL1B were also observed (Figure 5C and D). These findings demonstrated that the four key genes, especially PTK2B and PPFIA4, play a specific role in immune infiltration in cervical cancers.

Figure 5 Correlations of key genes expression with immune infiltration level in TCGA cervical cancer patients. (A) PTK2B expression was significantly negatively related to tumor purity, and had significant positive correlations with infiltrating levels of CD8+ T cells, neutrophils, and dendritic cells in cervical cancer patients. (B) PPFIA4 expression had significant negative correlations with infiltrating levels of CD8+ T cells. (C) SMAD7 expression had weak positive correlations with infiltrating levels of neutrophils and dendritic cells, and a weak negative correlation with tumor purity. (D) IL1B expression had significant positive correlations with infiltrating levels of neutrophils and dendritic cells. Each dot represented a single tumor sample.

Key Genes Significantly Correlated with the Lysosome

To determine whether lysosomes are associated with key genes in patients with cervical cancer, GSEA associated with the expression of key genes was conducted. The GSEA revealed significant differences (nominal P-value < 0.05, false discovery rate q value < 0.25) in the enrichment of the lysosome gene set for the expression of PTK2B and PPFIA4 (c2.cp.kegg.v6.2.symbols.gmt) (Figure 6A and B). The correlations among PTK2B, PPFIA4, and lysosomal genes were then evaluated by two methods. First, mutation interactions of these genes in patients with cervical cancer were explored, which showed significant co-occurrence of mutations in 26 lysosomal genes and in PPFIA4 and PTK2B (Figure 6C). Second, correlation analyses of gene expression levels were also performed, and strong correlations between key genes and certain lysosomal genes were observed. The results are displayed as Pearson’s correlation coefficient bubble-plots (Figure 6D), heatmap (Figure 6E), and scatter-plots (Figure 6F).

Figure 6 Correlations of key genes with lysosome. GSEA results showed KEGG term “lysosome” is differentially enriched for the expression of PTK2B (A) and PPFIA4 (B) in cervical cancer. (C) PPFIA4 and PTK2B were in significant co-occurrence of mutation with 26 lysosomal genes in KEGG term “lysosome” gene set. (D–F) Correlation analyses of gene expression levels were also performed, and strong correlations between key genes and certain lysosomal genes were observed. The results were displayed as Pearson’s correlation coefficient bubble-plots (D), heatmap (E) and lysosomal genes with top highest Pearson’s r with key genes were shown in scatter-plots (F). Abbreviations: GSEA, Gene Set Enrichment Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Network of Key Genes and Anti-Cancer Drugs

Figure 7 shows the relationship between the key genes and the other 46 most frequently altered neighbor genes. The correlations between anti-cancer drugs and the key genes were also illustrated. PTK2B and IL1B were the direct targets of anti-cancer drugs (genistein, leflunomide, and gallium nitrate). Although there was no anti-cancer drug directly targeting PPFIA4 and SMAD7, many of their most frequently altered neighbor genes have been identified as the targets of anti-cancer drug amounts.

Figure 7 Network of key genes and anti-cancer drugs. The network contained 50 genes (4 key genes and their 46 most frequently altered neighbor genes). Relationships between these genes and anti-cancer drugs (hollow hexagon) were also illustrated. Yellow hexagons indicate the FDA approved anti-cancer drugs. PTK2B and IL1B were the targets of anti-cancer drugs including genistein, leflunomide and gallium nitrate.

Relationship Between Clinical Features and Key Gene Expression

The relationships between the expression of the four key genes and clinical features were also explored. The expression of PPFIA4 was highly associated with clinical stage (P = 2e-7; Figure 3E) and tumor (T) classification (P = 7.7e-06; Figure 3F). The expression of IL1B was significantly associated with tumor (T) classification (P = 0.0072; Figure 3G). In addition, the expression of PPFIA4 was associated with age and radiation treatment and was thus unsuitable for further analysis (Table 1).

Chemotherapy Resistance Analysis

A total of 1288 DEGs were identified, of which 319 were significantly upregulated genes, and 969 were downregulated genes. The top genes sorted by six or more parameters in the CytoHubba plug-in were identified as hub genes. PIK3R2 was identified as a common hub gene in both MNP treatment and chemotherapy resistance (Figure 8A).

Figure 8 PIK3R2 as a common hub genes in MNP treatment and cisplatin resistance (A). mRNA expression levels of key genes in HeLa (B) and SiHa (C) cells by qRT-PCR. Synergistic killing effects of MNP with cisplatin in HeLa (D) and SiHa (E) cells using MTT assay. (F) mRNA expression levels of key genes after combination of cisplatin and MNP in HeLa cells by qRT-PCR. Data were presented as mean ± SEM. *p < 0.05; **p < 0.01; ***p < 0.001.

The validation of the key genes was performed at the mRNA level in both the HeLa and SiHa cells by qRT-PCR. The MNP treatment significantly downregulated the mRNA expression levels of PTK2B, IL1B, SMAD7, and PIK3R2 (Figure 8B and C).

The MTT assay was performed in HeLa and SiHa cells to explore whether MNP enhanced the killing effects of cisplatin in cervical cancer. The results show that MNP increased the sensitivity of HeLa and SiHa cells to cisplatin after 24 h of treatment (Figure 8D and E). The qRT-PCR demonstrated that a combination of MNP with cisplatin significantly reduced the expression levels of SMAD7 and PIK3R2 in HeLa cells (Figure 8F).

Discussion

In the present study, we showed that MNP has no cytotoxicity in cervical cancer cells; however, changes in the ratio of 650 to 540 nm fluorescence of AO were observed. The LTR red fluorescence decreased, indicating a decline in the lysosomal proton gradient. The altered fluorescence could also reflect lysosomal membrane permeabilization, besides an elevated lysosomal pH indicator.³³ However, the loss of the ability to retain lysosomotropic dyes can be attributed to either lysosomal membrane permeabilization or proton gradient dissipation. Nonetheless, it indicates lysosomal destabilization.³⁴ Overall, these results demonstrated that MNP interferes with lysosomes but has no toxicity in cancer cells, which are the two most important properties of an MNP as a carrier in drug delivery systems. A previous study demonstrated that MNP was transported into lysosomes, where it interfered with the lysosomal hydrolases and caused lysosome instability,^35,36 which is consistent with our findings. Lysosomes are thought to play a crucial role in regulating the biological activity of MNPs. Accumulating data shows that malignant transformation is associated with dramatic alterations in lysosomal structure and functions to meet increased metabolic demands. Furthermore, lysosomal cathepsins are involved in tumor progression and invasion. Overall, these changes greatly sensitize cancer cells to lysosomal destabilization.^37,38 Thus, the intrinsic properties of the MNP involved in lysosomes make it a promising material as a chemotherapeutic carrier. The kinetics of nanoparticle dissolution were found to be strongly dependent on the lysosomal pH and lysosomal membrane permeabilization, only in acidic conditions the MNP significantly released iron ions in solution.^39–41 The effects of lysosomal destabilization on the dissolution of iron oxide nanoparticles should be emphasized.⁴²

A systems approach to address the issues of onco-pharmacologic mechanism of drugs is recognized to have advantages over the conventional hypothesis-driven methodology, which dominantly identifies a single target or its associated pathway.^{23,24,43–45} Here, we followed a data-driven approach, which provides an unbiased view of alterations in cancer cells in response to MNP exposure by statistically analyzing whole-genome expression profiling without a previously defined hypothesis. Jameel et al studied transcriptional alterations following MNP exposure in cervical cancer cells and showed that MNP interferes with transforming growth factor-β (TGF-β) signaling by inhibiting the expression of ID and SMAD genes.⁴⁶ The identified genes were partially different from our results, and notably, the authors selected the genes for further investigations by ranking the false discovery rate, which would lead to the loss of functional genes. In the present study, we employed CytoHubba to explore important nodes or hubs by ranking 12 parameters and excluded genes without clinical relevance in cervical cancer. Navin et al proposed a systems pharmacology approach to identify genes possibly relevant for cisplatin resistance,⁴⁷ which provided a broad insight into the molecular mechanisms. We made improvements by analyzing multi-omics data and associating key genes with clinical phenotypes.

This work identified PPFIA4, PTK2B, SMAD7, and IL1B as key players in MNP onco-pharmacology. PPFIA4, also known as liprin-α4, is a member of the liprin gene family.⁴⁸ PPFIA4 was found to contribute to the proliferation of SBC-5 small-cell lung cancer cells both under normoxia and hypoxia.^49,50 The protein tyrosine kinase PTK2B, also known as PYK2, is a non-receptor tyrosine kinase acting as a signaling mediator involved in tumorigenesis, invasion, and metastasis.^51–53 A completely opposite role of PTK2B in cancer was also reported in ovarian carcinoma cells and prostate cancer cells as a tumor suppressor,^54,55 which is in agreement with our observation that high expression levels of PTK2B are correlated with improved survival in patients with cervical cancer. Interleukin-1β genes encode an inflammatory cytokine IL-1β, which favors the maintenance of inflammatory tumor microenvironment that contributes to the acquisition of the aggressive phenotype in cancer cells.⁵⁶ SMAD7 is an intracellular protein that allows transmembrane delivery of TGF-β. Similar to PTK2B, SMAD7 also possesses dual roles as a tumor suppressor and activator depending on tumor types.^57,58

The findings of drug interaction analysis showed that PTK2B was the target of genistein and leflunomide, and IL1B was the target of gallium nitrate. Considering the upregulation of these two genes upon MNP treatment in cervical cancer cells, genistein, leflunomide, and gallium nitrate were recommended as potential partners of MNP in drug delivery systems for further investigation. No drug targeting PPFIA4 and SMAD7, which might be the promising targets of new anti-cancer drugs, was identified. Genistein as an active soy isoflavone has been found to mediate anti-cancer effects. Hua-Yan et al showed that compared to free genistein, the MNP-genistein delivery system exhibits a significantly enhanced inhibition effect on the gastric cancer cell line SGC-7901.^59,60 Leflunomide, as an inhibitor of dihydroorotase dehydrogenase, was initially approved for the treatment of rheumatoid arthritis. Recently, leflunomide was also identified as a potent anti-cancer drug.⁶¹ Pund et al explored the delivery of leflunomide by nanoemulgel in melanoma and observed a significantly enhanced therapeutic response.⁶²

The four key genes identified have been reported to affect the sensitivity of cancer cells to several anti-cancer drugs. Mendoza et al reported that the IL-1β-induced upregulation of genes potentiates the acquisition of cisplatin resistance in breast cancer cells.⁵⁶ PTK2B has been demonstrated to play a pivotal role in integrin signaling in resistance to metformin, cisplatin, and proteasome inhibitors in multiple cancers, including breast cancer, ovarian cancer, myeloma, and non-small-cell lung cancer.^47,53,63,64 SMAD7 has also been found to be associated with chemoresistance to HDAC4/5 inhibitor, cisplatin, and sorafenib in leukemia, lung cancer, and hepatocellular carcinoma.^65–67 In a previous study by Hideya et al, PPFIA4 inhibition sensitized SBC-5 cells to cisplatin under hypoxia.⁵⁰ However, evidence is lacking for the involvement of the four key genes in cervical cancer chemoresistance, which might provide critical information on drug selection in the MNP delivery system in cervical cancer treatment.

To achieve personalized medicine, the systematic investigation of the interface between the MNP and biological environments of different cancer types is necessary owing to the highly heterogeneous tumor biology and microenvironments.^68–70 It should be noted that key genes of MNP in cancer cells and their roles may vary among different cancer types. A full understanding of the genomic alterations upon MNP treatment in different tumor types might have implications for specific cancer treatments. Furthermore, it would be interesting to investigate the genomic alterations in non-cancer cells. Consequently, delivery drug selection decisions will be based on not only the tumor molecular profiles in response to drugs but also the delivery carriers.

Conclusion

The results demonstrated that treatment of bare MNPs in cervical cancer cells significantly reduced the mRNA expression levels of PTK2B, PPFIA4, SMAD7, and IL1B, which play crucial roles in cervical cancer development and progression. The significance of the four key genes was validated by performing multi-omics data analysis, clinical relevance estimation, and qRT-PCR. Our study suggested that interactions of the key genes with anti-cancer drugs must be considered in the rational design of an MNP drug delivery system. Thus, this study could pave a novel way toward a systems pharmacology approach, characterizing the most potential molecular reaction to MNP in cervical cancer, which may help in further research on drug selection in MNP delivery systems to exert synergetic or additive anti-cancer effects.