Integrated analysis of transcriptomic and metabolomic data demonstrates the significant role of pyruvate carboxylase in the progression of ovarian cancer

Identifying candidate genes significantly associated with survival of EOC

To identify the candidate genes with prognostic value in ovarian cancer, we downloaded the gene expression profiles from the TCGA and GEO databases, which was followed by screening of DEGs and text mining.

Overall, we screened 1153 and 1022 genes related with prognosis from TCGA and GSE49997 dataset, respectively. After Venn diagram analysis, 66 overlapping genes were identified as candidate genes related to EOC progression (Figure 1).

After text mining, a total of 16 genes were found to be associated with EOC (Table 1). KEGG pathway analyses showed that the metabolic pathway (hsa01100) was the most significant pathway involved with GGPS1, ME1, DSE, NTPCR (also known as C1orf57, MGC13186), PPOX, and PC (Table 2). Among the metabolic pathway-related genes, ME1 and DSE have been reported to be associated with ovarian cancer. The KM survival curve for the remaining four genes is shown in Figure 2. High expression levels of GGPS1 and NTPCR correlated with a longer patient survival time, while high expression levels of PC and PPOX were closely related to a shorter patient-survival time.

The dysregulation of metabolic pathway-related genes in ovarian cancer cells

To evaluate the role of the metabolic pathway-related genes in ovarian cancer progression, we initially conducted real-time qPCR analyses to examine the expression levels of genes in different ovarian tumor cells. The results showed that PC was significantly upregulated in four cancer cell lines (SKOV3, CAOV-3, OV-1063, and OVCAR-3) compared to that in control IOSE80 cells (Figure 3A). Other genes were dysregulated in the different cell types. For example, GGPS1 was upregulated in the OV-1063 and OVCAR-3 cells and downregulated in CAOV-3 cells; PPOX was only upregulated in OV-1063 and OVCAR-3 cells. Considering these findings, the siRNA knockdown of PC was designed to explore its roles in ovarian cancer development. The knockdown of PC was confirmed using qRT-PCR analysis, which revealed that the expression of PC was significantly reduced in SKOV3 cells-transfected with shPC-1 and shPC-3, and shPC-2 presented off-target effects (Figure 3B). The degree of knockdown for shPC-1 and shPC-3 were 66.7% and 58%, respectively. Thus, shPC-1 transfected SKOV3 cells were used for subsequent experiments.

PC significantly inhibits proliferation, cell cycle progression, and cell migration and invasion in SKOV3 cells

The CCK-8 assay was conducted to evaluate the changes in proliferation of SKOV3 cells after PC knockdown. Depletion of PC significantly decreased proliferation in the shPC group compared to that of the shNC group (Figure 3C). The effect of PC on cell cycle progression was also examined using flow cytometry. The percentage of cells in the G1 phase was increased from 49.56% in the control group to 68.91% in shPC cells. In contrast, the percentage of cells in the S phase was decreased from 32.41% in the control group to 17.31% in the shPC group (Figure 3D). These results indicate that PC-knockdown significantly inhibited cell cycle transition from the G1 to S stage.

Invasion and metastasis are major events in the progression of cancers. Therefore, we further investigated the effect of PC on ovarian cell migration and invasion using the Transwell system. Results revealed that decreased expression of PC inhibited the invasive and metastatic abilities of ovarian tumor cells (Figure 3E).

Transcriptome sequencing reveals genes and pathways regulated by PC

After filtering out genes with low expression levels, a total of 22,838 gene expression matrices were obtained. The correlation between samples based on gene expression was measured using the Pearson’s correlation coefficient. A heat map of sample correlation is shown in Figure 4A. PCA analysis was performed for validation (Figure 4B) and the results showed that the samples in the different groups could be distinguished. Using thresholds of |log2FC| > 1.585 and FDR < 0.05, 1404 DEGs were identified in the shPC group compared to the shNC group, including 586 up-regulated genes and 818 down-regulated genes. Bidirectional clustering analysis illustrates that the expression profiles of DEGs are significantly different between the shPC and shNC groups (Figure 4C). We further analyzed the molecular functions of the genes and found that DEGs in shPC samples were mainly associated with steroid hormone biosynthesis and drug metabolism (Figure 4D). GSEA showed that the down-regulation of genes in the shPC group was more associated with pathogenesis of ovarian cancer compared to that in the shNC group (Figure 4E).

A PPI network was constructed, which consisted of 1495 edges connecting 431 nodes (Figure 5A). In the PPI network, STAT1 (degree = 54) and TP53 (degree = 46) had higher degree values than other genes, which indicated the significant role of STAT1 and TP53 in ovarian cancer progression. In the miRNA-TF-gene network, several miRNAs and TFs were identified as hub factors associated with ovarian cancer. Hsa-miR-124-3p, EZH2, and EP300 were the significant nodes in the miRNA-TF-gene network (Figure 5B).

Metabolomics sequencing reveals the dysregulated metabolic profile induced by PC knockdown

After normalization, the data were subjected to PCA and PLS-DA discriminant analyses under ESI+ and ESI- modes, separately. PCA showed that the samples in the different groups could be clustered together, and there were small offsets for the QC samples, which suggested data stability (Figure 6A). PLS-DA score plots showed obvious group separation between samples, indicating that the data were stable and reliable (Figure 6B). The cumulative R2Y and Q2Y values were both close to 1.0 (Figure 6C), which indicated that the PLS-DA models were stable, credible, and supported the differences between groups.

As for the shPC vs. shNC groups, a total of 44 metabolites met the cutoff value of FDR < 0.01, |log2FC| > 1, and VIP > 1, including 22 up-regulated and 22 down-regulated metabolites. Bidirectional clustering heat maps and volcano plots were used to visualize the results (Figure 7A and 7B). The differential metabolites in the shPC samples were mainly enriched in the metabolic pathways, ABC transporters, and nicotinate and niacinamide metabolism (Figure 7C).

Dysregulated pathways at both transcriptomics and metabolomics levels

Using an integrated analysis of transcriptomic and metabolomics data, eight pathways were found to be dysregulated by DEGs and differential metabolites in SKOV3 cells induced by PC knockdown, such as neuroactive ligands receptor interaction; glycine, serine, and threonine metabolism; and ABC transporter pathways (Figure 7D).

Data sources

The mRNA expression profiles associated with ovarian cancer were downloaded from The Cancer Genome Atlas (TCGA) database, which were derived from 299 tissue samples. The microarray dataset GSE49997 was downloaded from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/) based on the platform of GPL2986 ABI Human Genome Survey Microarray Version 2. The GSE49997 dataset includes 204 EOC tissue samples, of which 194 samples contain corresponding clinical information.

Identification of survival-related genes

To identify crucial genes in ovarian cancer development, the gene expression profiles of samples with clinical information from TCGA and GSE49997 datasets were further analyzed. The raw count data were downloaded and normalized using the count per million (CPM) method. Genes with low expression levels were removed, and the remaining 15,683 genes with high expression levels were subjected to survival analyses. Prognostic risk assessment was based on survival information and normalized gene expression profiles. Patients were divided into high-expression and low-expression groups based on the median value of gene expression. Kaplan-Meier (KM) curves were generated using the Survival package in R and statistically tested using the log-rank test. Genes with a p-value of less than 0.05 were considered as potential genes associated with prognosis in EOC patients. Similarly, after normalization, the gene expression datasets from EOC tissues in the GEO dataset were subjected to survival analyses. The genes closely related with prognosis based on the two databases were analyzed using Venn diagrams. The overlapping genes were then selected for further analysis.

Pathway enrichment analysis and literature search

Text mining was performed for the overlapping genes using Perl code. The published genes that are closely related with ovarian cancer were searched in the PubMed database. In addition, the overlapping genes were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool. Pathways with p < 0.05 and counts ≥ 2 were considered significant.

Cell culture

Four human ovarian cancer cell lines (SKOV3, Caov-3, OV-1063, OVCAR-3), normal epithelial ovarian cells (IOSE80), and human embryonic kidney (293T) cells were purchased from the cell bank of China Academic of Science. The SKOV3 cells were cultured in McCoy’s 5A Media (modified with tricine) supplemented with 10% fetal bovine serum (FBS). The Caov-3 and 293T cell lines were maintained in 90% DMEM with 10% FBS. The OV-1063 and IOSE80 cell lines were maintained in 90% RPMI 1640 with 10% FBS, and the OVCAR-3 cells were cultured in 80% RPMI 1640 with 20% FBS, sodium pyruvate, and 0.01 mg/ml bovine insulin. All the cell lines were cultured in an atmosphere of 5% CO₂ and 95% air at 37 °C.

Real-time qPCR analysis

Total RNA was extracted from cells using the RNAiso Plus (Trizol) reagent (TaKaRa, Japan) and cDNA was synthesized using the PrimeScript™RT Master Mix (Perfect Real Time) kit (TAKARA, Japan) according to the manufacturer’s instructions. Real-time PCR was performed to evaluate the expression levels of GGPS1, NTPCR, PPOX, PC, PRICKLE2, TCF7L1, and PPP3CA in tumor cells. A total of 8 μl of cDNA was used as template in a final 20 μl PCR volume containing 1 μl forward primer, 1 μl reverse primer, and 10 μl SYBR Premix EX Taq (2x). PCRs were run as follows: 50.0 °C for 3 min, 95.0 °C for 3 min, followed by 40 cycles of 95.0 °C for 10 s and 60.0 °C for 30 s. Following PCR, a melting curve was obtained at temperatures from 60 °C to 95 °C, at increments of 0.5 °C for 10 s. Primer sequences are listed in Table 3.

Cell transfections

Based on designing sigma software (https://www.sigmaaldrich.com/life-science/functional-genomics-and-rnai/sirna/mission-predesigned-sirna.html), the shRNA sequence of PC was obtained. Three shRNAs sequences in CDS region for PC, GCCCAGTTTATGGTGCAGAAT (shPC-1), GCCAAGGAGAACAACGTAGAT (shPC-2), and ATGGGCATCCGCCTGGATAAT (shPC-3) were selected to knock down endogenous PC, and messy sequence (NC- GTTCTCCGAACGTGTCACGTC) was used as control sequence. Then, the primers were designed and synthesized by Suzhou Jinweizhi Biological Technology Co. LTD. (Table 3). Inducible shRNA construction was performed by ligating annealed oligonucleotides into the pLKO.1 Puro vector digested with EcoRI and AgeI. The ligated products were transformed into chemically competent Stbl3 E. coli using Lipofectamine 3000 reagent (Invitrogen) according to the manufacturer’s instructions and grown in LB supplemented with ampicillin. After validation of shRNA insertion using gene sequencing, the plasmids were packaged into lentiviral particles using the third-generation lentiviral production system. SKOV3 cells were divided into three groups: black (without any treatment), shNC (transfected with negative control plasmid), shPC (transfected with shRNA-PC). Stable SKOV3 cells were selected with 0.5 μg/ml of puromycin over 3 days. Lastly, the efficiency of PC knockdown in SKOV3 cells was confirmed using real-time qPCR analysis, and the degree of knockdown [(1- the relative mean expression of PC in shPC group/shNC group) × 100%] for three PC shRNAs were calculated. The sable shRNA-PC transfected SKOV3 cells with highest knockdown degree were used for following experiments.

Cell viability and cell cycle analyses

Cell viability assays were performed using the cell counting kit-8 (CCK-8; Beyotime Biotechnology, China). Briefly, the SKOV3 cells in three groups (3 × 10⁴ cells/ml) were seeded in 96-well plates at a density of 3,000 cells/well, respectively. Cells were incubated at 37 °C in 5% CO₂ atmosphere conditions for 24, 48, 72, and 96 hours, and then 10 μl of CCK-8 solution was added into each well for 1.5 hours. The OD value of each well sample was measured at 450 nm using the Multiscan Spectrum System (BD Biosciences, USA). All assays were repeated in triplicate.

Flow cytometry was used for cell cycle analyses. shPC and control SKOV3 cells were digested with trypsin and collected using centrifugation at 300 × g, 4 °C for 5 min. Then, cells were fixed in 70% ethanol overnight at 4 °C. Subsequently, the samples were treated with RNase A for 30 min and stained with propidium iodide (BD Biosciences, USA) for 15 min. Flow cytometry was immediately conducted using a FACSCAN flow cytometer (BD Biosciences, USA) following the manufacturer’s instructions. All assays were repeated in triplicate.

Transwell migration and invasion assays

Invasion and migration assays were performed using BD Matrigel culture inserts. First, 8 μm pore 24-well Transwell inserts (Corning, USA) were coated with 20 μl of Matrigel (BD Biosciences) and incubated for 30 min at 37 °C in a 5% CO₂ incubator to allow gel formation. Control shNC and shPC SKOV3 cells were suspended at equal cell densities in serum-free medium. A total of 100 μl of cells (2 × 10⁵ cells/ml) were seeded in the upper chamber of untreated Transwells for cell migration analyses and in the upper chamber of Transwells treated with Matrigel for cell invasion analyses. The lower chambers were supplemented with 500 μl of 10% FBS-containing medium. Then, cells were incubated at 37 °C in 5% CO₂ atmosphere conditions for 12 and 24 h, respectively. The cells in Transwells were treated with 4% paraformaldehyde for 30 min and stained with crystal violet for 10 min. Subsequently, the upper chambers were washed three times in PBS and non-invading cells on the inner surface were carefully removed using cotton swabs. Finally, the invading cells on the reverse side of the upper chamber were counted using a microscope at three different fields to assess cell migration and invasion abilities. All assays were repeated in triplicate.

Metabolomics data acquisition and identification

Metabolomics data acquisition was performed using an ultra-performance liquid chromatography-quadrupole-time of flight liquid chromatography/mass spectrometry (UPLC-Q-TOF LC/MS) system [32]. The samples included 10 shPC-infected SKOV3 cell samples (1 × 10⁷ cells for each sample), 10 shNC SKOV3 cell samples, and 5 quality control (QC) samples. Cells in each sample were mixed with 500 μl of methyl alcohol/acetonitrile/distilled water (2:2:1, v/v/v) solution and disrupted using ultrasonication for 30 min. After incubation at -20 °C for 1 h, the supernatant was collected using centrifugation at 13,000 rpm at 4 °C for 15 min followed by freeze-drying. For metabolomic analyses, samples were re-dissolved in 100 μl of acetonitrile solution (1:1 ratio of acetonitrile and water, by volume) and centrifuged at 14,000 × g at 4 °C for 15 min. The supernatant was collected for LC/MS analysis.

The metabolic products acquisition and identification were both performed in positive-ion (ESI+) and negative-ion (ESI-) modes and analyzed using a Triple-TOF 5600 mass spectrometer. Data normalizations were performed on LC/MS data based on internal standards using ESI+ and ESI- modes.

Principal component analysis (PCA), extensively used in the statistical learning field, and the partial least squares discriminant analysis (PLS-DA) are commonly applied to evaluate differences between groups [33, 34]. In the present study, we used two-dimensional PCA to ensure detection stability and three-dimensional PLS-DA for cross validation. In addition, univariate nonparametric analyses and multivariate PLS-DA were performed for metabolic profiling using the ropls software, version 1.6.2 (http://bioconductor.org/packages/release/bioc/html/ropls.html). The thresholds for potential biomarker selection were set as follows: a false discovery rate (FDR) of <0.01, |log2 fold change (log2FC)| > 1 in univariate analysis, and Variable Importance in the Projection (VIP) > 1 in multivariate analysis. Bidirectional hierarchical clustering analyses were performed to assess data classification ability and concentration levels of the screened metabolites. For differential metabolites obtained from tumor cells, we transformed these data into the KEGG ID format using the MetaboAnalyst online tool [35] and performed KEGG pathway analyses using MBROLE 2.0 (http://csbg.cnb.csic.es/mbrole2/analysis.php) with a p-value < 0.05.

Transcriptome sequencing

Total mRNA extractions and cDNA library preparations were performed for shPC and control shNC cell samples. Then, two-paired end sequencing was performed using the Illumina platform (Illumina, San Diego, CA, USA).

Clean reads data were obtained using the Trimmomatic tool (version 3.6) [36]. Then, the reads were mapped to the human reference genome (GRCH38, Gencode) [37] with the Hisat 2 software, version 2.05 [38]. Gene expression levels were evaluated by counting reads mapped to protein-coding regions using FeatureCounts tools (v1.6.0) [39], and expression values were normalized using Fragments Per Kilobase of Exon Per Million Fragments Mapped (FPKM) method. Genes with an FPKM value ≥ 0.1 in at least three samples were further analyzed. According to filtered gene abundance expression profiles, the correlation of gene expression levels between samples was analyzed using PCA with the ggord package (https://zenodo.org/badge/latestdoi/35334615) [40].

DEGs between the shPC and shNC groups were identified using the quasi-likelihood F-tests method of the edgeR software [41]. Genes with |log2FC| > 1.585 and FDR < 0.05 were considered to be significant, and two-dimensional clustering heatmaps were used to visualize the gene expression profiles of DEGs. Subsequently, gene ontology (GO) function and KEGG pathway enrichments were performed for genes of interest using the Clusterprofiler package [42]. A p-value < 0.05 was considered as a significant difference.

Ovarian carcinoma-related genes were retrieved from the DisGeNET database [43] (http://www.disgenet.org/web/DisGeNET/menu/home). Differentially expressed genes were subjected to gene set enrichment analysis (GSEA) using the DOSE package [44].

The PPI pairs with required confidence (combined score > 0.7) were obtained with the STRING online tool [45] (https://string-db.org/). Cytoscape software [46] was utilized to construct a PPI network and the network topological properties were analyzed based on degree [47], betweenness [48], and closeness [49] using the cytoscape CytoNCA plugin [50].

Moreover, the significant modules in the PPI network were mined using the MCODE tool [51] with a screening score > 10. For the intriguing DEGs, we performed microRNA and transcription factor (TF) predictions by using Enrichr tools [52]. The corresponding miRNA-gene pairs and TF-gene pairs were also identified from the miRTarBase and ENCODE databases [53]. Finally, we integrated these miRNAs, TFs, and candidate genes to construct a systemic regulatory network.

Integrated pathway analysis of transcriptomics and metabolomics data

Integrated analyses are conducted in order to understand the biological function of post-genomic data at a higher level than individual biomolecules. The “IMPaLA” web tool has been used to integrate more than one type of omics data for pathway analysis [54]. In this study, the joint pathway analysis was performed for DEGs and differential metabolites based on the KEGG database. Lastly, pathways with number_of_overlapping_metabolites/gene > 0 and metabolite with p < 0.05 were considered to be significant.