Bioinformatics analysis and verification of molecular targets in ovarian cancer stem-like cells.

BACKGROUND: Epithelial ovarian cancer (EOC) is a lethal and aggressive gynecological malignancy. Despite recent advances, existing therapies are challenged by a high relapse rate, eventually resulting in disease recurrence and chemoresistance. Emerging evidence indicates that a subpopulation of cells known as cancer stem-like cells (CSLCs) exists with non-tumorigenic cancer cells (non-CSCs) within a bulk tumor and is thought to be responsible for tumor recurrence and drug-resistance. Therefore, identifying the molecular drivers for cancer stem cells (CSCs) is critical for the development of novel therapeutic strategies for the treatment of EOC.
METHODS: Two gene datasets were downloaded from the Gene Expression Omnibus (GEO) database based on our search criteria. Differentially expressed genes (DEGs) in both datasets were obtained by the GEO2R web tool. Based on log2 (fold change) >2, the top thirteen up-regulated genes and log2 (fold change) < -1.5 top thirteen down-regulated genes were selected, and the association between their expressions and overall survival was analyzed by OncoLnc web tool. Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways analysis, and protein-protein interaction (PPI) networks were performed for all the common DEGs found in both datasets. SK-OV-3 cells were cultured in an adherent culture medium and spheroids were generated in suspension culture with CSCs specific medium. RNA from both cell population was extracted to validate the selected DEGs expression by q-PCR. Growth inhibition assay was performed in SK-OV-3 cells after carboplatin treatment.
RESULTS: A total of 200 DEGs, 117 up-regulated and 83 down-regulated genes were commonly identified in both datasets. Analysis of pathways and enrichment tests indicated that the extracellular matrix part, cell proliferation, tissue development, and molecular function regulation were enriched in CSCs. Biological pathways such as interferon-alpha/beta signaling, molecules associated with elastic fibers, and synthesis of bile acids and bile salts were significantly enriched in CSCs. Among the top 13 up-regulated and down-regulated genes, MMP1 and PPFIBP1 expression were associated with overall survival. Higher expression of ADM, CXCR4, LGR5, and PTGS2 in carboplatin treated SK-OV-3 cells indicate a potential role in drug resistance.
CONCLUSIONS: The molecular signature and signaling pathways enriched in ovarian CSCs were identified by bioinformatics analysis. This analysis could provide further research ideas to find the new mechanism and novel potential therapeutic targets for ovarian CSCs.

Introduction

Epithelial ovarian cancer (EOC) is a lethal and aggressive gynecologic malignancy with poor long-term prognosis. According to the American Cancer Society, 21,750 new cases and 13,940 deaths are predicted from ovarian cancer (OC) in 2020 [1]. Despite recent advances, standard therapies have provided only a marginal survival rate in the past decade. According to the data published by the SEER 18 (2010–2016), the current 5-year relative survival rate in the US is approximately 48.6% [2]. A crucial factor limiting the efficacy of current treatments is that most patients experience delayed onset of signs and symptoms until the disease has locally advanced [2]. The primary treatment consists of surgical removal of the tumor, followed by chemotherapeutics and radiation therapy [3, 4]. Although patients are highly responsive to the standard first-line chemotherapy, the vast majority of patients (approximately 85%) present with recurrent disease and drug-resistance [5, 6].

EOC tumors consist of a heterogeneous population of cells, including a small population known as CSCs. CSCs are pro-tumorigenic and have been identified in several human solid malignancies, including OC [7, 8]. Therefore, eradication of CSCs is a potential therapeutic avenue to increase the survival rate of EOC patients. Recently, studies have shown that targeting key molecular players and genetic aberrations in CSCs using small molecule inhibitors, including rapamacyin and triciribine can eradicate CSCs in brain tumors and improve overall survival in glioblastoma and neuroblastoma patients [4, 8, 9, 10, 11, 12]. Identifying the molecular signature of ovarian CSCs will provide a novel therapeutic strategy to selectively target this subpopulation of cells and improve overall survival. The primary challenge is to identify and isolate CSCs [13]. Due to the heterogeneous nature of the tumor [6, 14], it is critical to differentiate CSCs from cancer cells (CCs). Although, various cell surface and non-surface markers are currently in use for the identification and isolation of ovarian CSCs [6, 15, 16], the phenotypic and functional heterogeneity in CCs within the tumor require further dissection to understand intra-tumor cell population complexity to eradicate CSCs more effectively [11, 17].

It is widely believed that the origin of CSCs is attributed to cellular reprogramming during oncogenic transformation [18]. CCs and CSCs exist in a dynamic equilibrium that is tightly regulated by cellular signaling pathways associated with tumor microenvironments [19]. Any disturbance in this dynamic equilibrium facilitates CCs to acquire stem cell-like properties by undergoing a differentiation process [18]. The regulation of such equilibrium involves a complex crosstalk between CSCs and CCs in CSCs niche. Therefore, differential gene expression profiling of CCs and CSCs is critical for the identification of novel molecular targets that can be exploited to improve OC patient outcomes.

This present study assesses the DEGs in CCs (monolayer cancer cells) and CSLCs (spheroids) using two microarray datasets. All DEGs were further analyzed for gene ontology to explore their biological processes and molecular functions. Moreover, the signaling pathway enrichment analysis and PPI network of DEGs were carried out. Based on the fold change value of genes, the top thirteen up-regulated and thirteen down-regulated DEGs were selected commonly from both datasets and further validated by q-PCR. Survival analyses of all thirteen up-regulated and thirteen down-regulated genes were performed to investigate the association between their expressions and clinical prognosis. To determine the potential role of DEGs in carboplatin sensitivity/resistance, we investigated the IC50 value of carboplatin in SK-OV-3 cells and explored the changes in gene expression of all twenty-six DEGs induced by carboplatin treatment in SK-OV-3 cells following carboplatin treatment. Taken together, we performed a preliminary analysis of the DEG's role in OC prognosis and drug sensitivity.

Materials and methods

Gene expression datasets collection from GEO database

Gene expression datasets were manually searched and gathered from the NCBI GEO database (www.ncbi.nlm.nih.gov/geo), using “ovarian spheroids” as keywords and selected two datasets, i.e., GSE28799 and GSE80373. Dataset GSE28799 consists of the gene expression profile of six samples, expression (in triplicate) of ovarian carcinoma OVCAR-3 cell cultured in adherent culture medium (2D), and spheroid derived cells (3D) cultured in suspension culture. The platform for the experiments in dataset GSE28799 was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. Gene expression profile of eight samples was available in dataset GSE80373, which contains a quadruplicate of each HEY ovarian cancer 2D monolayer sample and 3D spheroid sample, where the experimental platform is used as GPL13667 ([HG-u219] Affymetrix Human Genome U219 Array).

Identification of DEGs using GEO2R

We defined the sample groups as “ovarian cancer cells (CCs)” for cells cultured in adherent culture medium (2D) and “ovarian cancer stem-like cells (CSLCs)” for cells cultured in suspension culture (3D spheroids). The value distribution was checked for all samples in both datasets. The value data was median-centered across the samples, which indicates that the datasets were suitable for cross-comparison among their defined grouped samples. The raw p-values of both the datasets were adjusted to the FDR using the Benjamini & Hochberg method [20]. Fold change of the genes present in CSLCs and CCs were analysed using the limma package in R (Version: R 3.2.3 https://doi.org/10.1093/nar/gkv007). We determined the fold change of individual genes using the GEO2R web tool from GEO database. To determine the genes that are significantly expressed in CSLCs compared to CCs, the FDR adjusted p-values were kept to less than 0.05. Based on the fold change, genes were categorized into two classes, up-regulated genes (log2FC > 1) and down-regulated genes (log2FC < -1). Unknown genes in the experimental data were removed, and genes with multiple expression data were adjusted to single expression data based on their higher fold change value. Genes that were commonly up-regulated or down-regulated in both datasets were used for further analysis.

Volcano plot and heatmap analysis

The GraphPad Prism v8 was used for the generation of volcano plot and heatmap analysis. Volcano plots were generated using whole gene expression data obtained after GEO2R analysis for both datasets separately. Fold change values of the top thirteen up-regulated and thirteen down-regulated genes present commonly in both datasets were standardized in excel and estimated as z score values.

GO and pathways enrichment analysis of DEGs

All the common DEGs (117 up-regulated genes and 83 down-regulated genes) in both datasets were analyzed using the DAVID v6.8 (https://david.ncifcrf.gov/) online database. Using a functional annotation tool in DAVID, GO analysis and signaling pathway enrichment analysis were achieved, where the p-value was kept less than 0.05. KEGG and Reactome pathways were used for signaling pathway enrichment analysis with a p-value of less than 0.05. For GO and pathway enrichment analysis, the annotations and the background species were kept as Homo sapiens in DAVID.

Protein class identification and PPI network

All the common DEGs (117 up-regulated genes and 83 down-regulated genes) were used to identify protein classes and the PPI network. PANTHER Classification System online database was used to identify DEG's protein classes. The statistical analysis tool used for the analysis was the PANTHER Overrepresentation Test (Released 20190711) with Annotation Version and Release Date PANTHER version 14.1 Released 2019-03-12. The reference list used for the DEGs (analyzed list) comprises all Homo sapiens genes present in the PANTHER database. Annotation dataset was used as PANTHER protein classes, and test type was kept as Fisher's exact test with FDR correction.

STRING v11.0 (https://string-db.org/) online database was used to find the interaction between translated proteins of the DEGs. A confidence score higher than 0.4 was kept as the cut-off value for the interaction. Cytoscape software was used to visualize and construct the PPI network, where nodes in the network represent translated protein from the identified DEGs and edge in the network represents the interaction.

Cell and spheroid culture

SK-OV-3 (ATCC® HTB-77™) cells were purchased from ATCC and cultured in Mc Coy's 5A medium (Thermo Fisher Scientific) supplemented with 10% of FBS (Sigma) and 1% penicillin-streptomycin (Thermo Fisher Scientific). Cells were cultured as an attached monolayer at 37 °C in 5% CO2. Cells in passages 3 to 8 were used for all studies. Mycoplasma contamination of cells was screened periodically by using Mycoplasma Detection Kit (Lonza). Spheroids were generated from SK-OV-3 cells in ultra-low attachment petri-dishes (Corning) and cultured in knockout DMEM (Thermo Fisher Scientific)/F12 medium supplemented with 20% knockout serum replacement (Life Technologies), 20 ng/mL epidermal growth factor (EGF), 10 ng/mL basic fibroblast growth factor (bFGF), 1% L-glutamine, and 1% penicillin-streptomycin.

Growth inhibition assay

Growth inhibition assay was performed in accordance with the method published by Smith et al. [21] with slight modifications. SK-OV-3 monolayer cells were plated at a density of 5000 cells/well in 96 well plates. After 24 h, cells were either treated with diluents (control) or treated with varying concentrations of carboplatin (Sigma). After 72 h, 100 μl of MTT (Millipore) solution was added (final concentration of 0.5 mg/ml) and kept for 4 h. The supernatant was replaced with 100 μl of DMSO and absorbance was read at 570 nm. The IC50 of carboplatin was calculated by using GraphPad Prism version 8.0. All experiments were done in triplicates.

Experimental verification of selected DEGs by q-PCR analysis

To examine whether the spheroids were enriched in CSCs, gene expression of the stem cell markers such as CD-133, NANOG, NESTIN, OCT-4, and SOX-2 were performed by q-PCR. Briefly, RNA was extracted from both SK-OV-3 monolayer and spheroid cells with Trizol reagent (Thermofisher Scientific). 2 μg of RNA was converted into cDNA by cDNA Reverse Transcription Kit (Bio-rad). Gene expression was quantified through the use of SYBR green q-PCR (Bio-rad) assay and gene-specific primers (Table S1). Relative gene expression was quantified by normalizing the gene-specific amplification to that of housekeeping gene 18S in the same sample. The relative expression of each gene was calculated using (ΔΔCT) methods [22]. Experiments were performed in triplicates.

For experimental verification of selected top thirteen up-regulated and top thirteen down-regulated DEGs, total RNA was isolated from both SK-OV-3 monolayer and spheroid cells, and converted into cDNA. q-PCR analysis was performed using the methods mentioned above. The Log2(Fold Change) of both cell populations was calculated using GraphPad Prism version 8.0. Data were shown Log2(Fold Change) relative transcriptional changes in Y-axis vs. genes in X-axis. Experiments were performed in triplicate.

Carboplatin treatment and q-PCR analysis

A total of 0.6 × 106 SK-OV-3 monolayer cells were plated in a 60 mm dish. After 24 h, cells were treated with either diluents (control) or carboplatin (IC50 concentration) for 12, 24, and 48 h. Total RNA was extracted from both control and carboplatin-treated cells with Trizol reagent (Thermofisher Scientific). q-PCR analysis was performed as previously described using gene-specific primers (Table S1). Experiments were performed in triplicates.

Survival analysis using OncoLnc platform

OncoLnc (www.oncolnc.org) is an online interactive tool that explores survival correlations and downloads clinical data coupled with expression data for mRNAs, miRNAs, and long non-coding RNAs [23]. OncoLnc provides a broad set of survival data from 8,647 patients across 21 cancer studies from the TCGA and MiTranscriptome beta analysis. Using the OncoLnc platform, Cox regression data were obtained for MMP1 and PPFIBP1 in ovarian serous cystadenocarcinoma (OV). Patients were sorted by the expression of gene of interest. We kept 25 as lower and upper percentile as recommended by OncoLnc. Upon submission, Kaplan-Meier plots with a logrank p-value were plotted for further analysis.

Statistical analyses

The data are represented as mean ± SD. An unpaired student's t-test was performed to assess the significant difference between groups. The level of significance was based on a p-value < 0.05. Survival curves were plotted using OncoLnc. All results are represented with p-values obtained from a logrank test.

Results

Screening of DEGs between ovarian CCs and CSCs

Two microarray datasets (GSE28799 and GES80373) that compare ovarian CCs with CSLCs were used. GSE28799 dataset consisted of 2,388 DEGs, including 1,255 up-regulated genes and 1,133 down-regulated genes. GSE80373 dataset consisted of 2,100 DEGs, including 1,104 up-regulated genes and 996 down-regulated genes. Volcano plots were generated for each microarray data set to show DEGs (Figure 1A, B). By comparing both datasets, we found a total of 117 genes were commonly up-regulated, and 83 genes were down-regulated based on their gene expression-fold change (|log2 (fold change)| ≥ 1) in both datasets (Figure 1C). Based on log2 (fold change)>2, the top thirteen up-regulated genes and log2 (fold change) < -1.5 for the top thirteen down-regulated genes were selected. An expression heat map of all twenty six commonly selected up-regulated and down-regulated genes among all DEGs from both datasets was generated (Figure 1D, E).

Figure 1

Identification of DEGs between ovarian CCs and CSLCs. (A & B) The DEGs between ovarian CCs and CSLCs in GSE28799 (A) and GSE80373 (B) datasets were presented in the volcano plots. The black nodes represent genes that are not differentially expressed, and the green and red dots represent down-regulated and up-regulated genes in ovarian CSLCs respectively. Any points on or above the grey horizontal line are significant (p < 0.05). (C) Venn diagrams show up-regulated and down-regulated genes from both datasets. Intersected areas represent the common DEGs in both datasets. (D & E) Heat maps show the selected top 13 up-regulated and down-regulated DEGs profiles between CCs and CSLCs commonly found in both datasets. The heatmaps indicate up-regulation (red), down-regulation (green), and mean gene expression (black). CCs, Cancer cells, CSLCs, Cancer stem- like cells, DEGs, Differentially expressed genes.

Experimental validation of DEGs

To validate the expression of top up-regulated and down-regulated DEGs, SK-OV-3 cells were either cultured adherently or in suspension to form a monolayer or spheroids, respectively (Figure 2A, B). To examine whether spheroids were enriched in CSCs, we examined the gene expression of stem cell markers, including CD-133, NANOG, NESTIN, OCT-4, and SOX-2 by q-PCR. The relative transcript levels of all stem cell markers were significantly elevated in SK-OV-3 spheroids than SK-OV-3 monolayer cells (Figure 2C). To validate the expression of selected DEGs, total RNA from SK-OV-3 monolayer and spheroids were extracted and q-PCR analysis was performed. Out of the 13 up-regulated genes, 11 genes were over-expressed in spheroids compared to monolayer cells. Interestingly, two genes, AKR1C1 and RRAD, were down-regulated (Figure 2D). Moreover, out of the 13 down-regulated genes, 7 were significantly down-regulated and few were over-expressed in SK-OV-3 spheroids (Figure 2E).

Figure 2

(A) SK-OV-3 cells cultured in an adherent culture medium. (B) SK-OV-3 spheroids are cultivated in a suspension culture medium. (C) SK-OV-3 spheroids are enriched in stem cell markers compared to monolayer cells. Relative gene expression of CSC markers in SK-OV-3 spheroids compared to monolayer cells performed by q-PCR. Results are represented as mean ± SD from triplicates, ∗∗p < 0.01, ∗∗∗p < 0.001 compared with monolayer cells. (D) Relative gene expression of thirteen up-regulated and, (E) thirteen down-regulated genes in SK-OV-3 spheroids compared to monolayer cells performed by q-PCR. Fold changes of genes were calculated using (ΔΔCT) methods, and 18S gene expression was used to normalize their expressions. The Log2(fold change) of SK-OV-3 monolayer and spheroid cells were calculated using GraphPad Prism version 8.0. Data were shown Log2 relative transcriptional changes (fold changes) in Y-axis vs. genes in X-axis. Results were represented as mean ± SD from triplicates, ∗∗p < 0.01, ∗∗∗p < 0.001 compared with monolayer cells.

Functional annotation for DEGs

Based on the three GO term classes, (molecular functions, biological processes, and cell components), biological processes was the most enriched in DEGs groups. Some of the up-regulated DEGs in CSLCs were enriched in the extracellular matrix (ECM) of the cells (GO:0031012). Most were enriched in biological processes including cell proliferation (GO: 0042127), tissue development (GO:0009888), and response to lipid (GO:0033993), chemical stimulus (GO:0070887), and organic substances (GO:0010033). Down-regulated DEGs were enriched in the biological processes, including molecular function regulation (GO:0065009) (Table 1).

Table 1

Gene ontology enrichment analysis of up-regulated and down-regulated genes in ovarian cancer stem like cells found in both datasets.

Up-regulated genes
Category	Term	Count	p-value	Genes
Biological Process	Tissue development (GO:0009888)	30	2.88E-06	S100A4, XDH, MTSS1, PTGS2, PGF, SOX4, HSBP1, LGR5, CALB1, TGF-β1, AKR1C3, AKR1C2, SAP30, ISG15, SERPINE2, HEY1, SORBS2, COL6A1, AKR1C1, ETV4, BMP4, ARID5B, PDE4D, EPHA4, ADM, PECAM1, TGFBR3, PRDM1, BAMBI, LUC7L
Biological Process	Response to lipid (GO:0033993)	20	3.80E-06	BMP4, CXCL1, PTGS2, NR4A2, PDE4D, MBD2, TGF-β1, IFNAR1, AKR1C3, ZFP36L2, PRMT2, INSIG2, ADM, HEY1, SLPI, FBXO32, PRDM1, LOX, F2R, NR2F1
Biological Process	Cellular response to chemical stimulus (GO:0070887)	38	4.47E-06	CXCL1, XDH, LTBP2, IFITM1, PTGS2, PGF, SOX4, CALB1, TGF-β1, AKR1C3, ZFP36L2, AKR1C2, PRMT2, ISG15, HEY1, CXCR4, PDE4A, COL6A1, AKR1C1, NR2F1, BMP4, MT1M, VAV3, NR4A2, PDE4D, MBD2, KCNK3, IFNAR1, IFIT3, IFIT2, RAB31, IFIT1, CD58, BNIP3L, TGFBR3, FBXO32, PRDM1, BAMBI
Biological Process	Response to organic substance (GO:0010033)	39	7.62E-06	XDH, CXCL1, LTBP2, IFITM1, PTGS2, PGF, SOX4, CALB1, TGF-β1, AKR1C3, ZFP36L2, AKR1C2, PRMT2, INSIG2, ISG15, HEY1, CXCR4, GPX3, COL6A1, LOX, AKR1C1, NR2F1, BMP4, NR4A2, PDE4D, MBD2, IFNAR1, IFIT3, IFIT2, RAB31, IFIT1, ADM, CD58, TGFBR3, SLPI, FBXO32, PRDM1, BAMBI, F2R
Biological Process	Regulation of cell proliferation (GO:0042127)	26	2.54E-05	CXCL1, XDH, MTSS1, IFITM1, PTGS2, PGF, SOX4, LGR5, TGF-β1, AKR1C3, AKR1C2, SERPINE2, DPP4, ETV4, BMP4, VAV3, MBD2, TNFSF9, IFIT3, FBLN1, ADM, TGFBR3, PRDM1, BAMBI, EMP3, F2R
Cellular Components	Extracellular matrix (GO:0031012)	14	3.00E-05	BMP4, LTBP2, SPOCK1, NID2, TGF-β1, MMP1, FBLN1, SERPINE2, PI3, SLPI, TGFBR3, COL6A1, SERPINA1, LOX
Down-regulated genes
Biological Process	Regulation of molecular function (GO:0065009)	29	1.16E-05	IL18, GPR87, TPM1, SHISA9, WNK4, GIT2, ARHGAP11A, CCNA2, CRIM1, KCNG1, CYR61, PRKCA, EGFR, UBXN1, SOCS2, CCNL1, TRIO, ANXA3, NCAM1, STOM, DUSP1, CHML, ID1, EIF4A2, MCPH1, SPTBN1, ZNF462, ID3, ADD2

Using information from KEGG and Reactome databases, significantly enriched pathways of both up-regulated and down-regulated DEGs were found from the DAVID database (Table 2). Biological pathways including interferon-alpha/beta signaling (IFIT3, IFIT2, IFIT1, ISG15, IFITM1, IFNAR1), elastic fiber signaling (BMP4, FBLN1, LTBP2, TGF-β1), and bile acid and bile salt synthesis (AKR1C3, AKR1C2, AKR1C1) were significantly enriched in up-regulated DEGs. Biological pathways, including the African trypanosomiasis pathway (PRKCA, IL18, IL12A), and EGFR transactivation by gastrin pathway (PRKCA, EGFR) were significantly decreased in down-regulated DEGs in CSLCs.

Table 2

Signaling pathways enrichment analysis of up-regulated and down-regulated genes in ovarian cancer-stem like cells found in both datasets.

Up-regulated genes
Pathway	Term	Count	p-value	Genes
Reactome	R-HSA-909733: Interferon alpha/beta signaling	6	1.92E-04	IFIT3, IFIT2, IFIT1, ISG15, IFITM1, IFNAR1
Reactome	R-HSA-2129379: Molecules associated with elastic fibres	4	3.44E-04	BMP4, FBLN1, LTBP2, TGF-β1
Reactome	R-HSA-193775: Synthesis of bile acids and bile salts via 24-hydroxycholesterol	3	5.45E-04	AKR1C3, AKR1C2, AKR1C1
Reactome	R-HSA-193807: Synthesis of bile acids and bile salts via 27-hydroxycholesterol	3	6.26E-04	AKR1C3, AKR1C2, AKR1C1
Reactome	R-HSA-193368: Synthesis of bile acids and bile salts via 7alpha-hydroxycholesterol	3	1.57E-03	AKR1C3, AKR1C2, AKR1C1
Down-regulated genes
KEGG	hsa05143: African trypanosomiasis	3	0.008357	PRKCA, IL18, IL12A
Reactome	R-HSA-2179392: EGFR transactivation by Gastrin	2	0.039951	PRKCA, EGFR

Protein classes identification

Various classes of translated proteins from DEGs were identified by using the PANTHER classification system online database (Table 3). The top five protein classes that were enriched in up-regulated DEGs in CSLCs were signaling molecules, enzyme modulators, transcription factors, oxidoreductases, and hydrolases. In contrast, the top five protein classes that were enriched in down-regulated DEGs were enzyme modulators, cytoskeleton proteins, transferases, transcription factors, and nucleic acid-binding proteins (Figure 3). Protein classes such as ECM proteins were found only in up-regulated DEGs. However, down-regulated DEGs consisted of protein classes such as adapter proteins, defense/immunity proteins, etc.

Table 3

Protein classes of up-regulated and down-regulated genes identified using PANTHER database.

Protein Classes	Up-regulated Genes	Down-regulated Genes
Signaling molecule	ADM, ANGPTL2, CXCL1, C3, S100A4, ETV1, ETV4, PGF, TGF-β1, LTBP2, SPRY4, BMP4	GDF15, CYR61, TRIO, OSMR
Enzyme modulator	CCNG2, SLPI, C3, PI3, SERPINA1, SPOCK1, DPP4, CCNL2, SERPINE2, DBI	CCNA2, CCNL1, RND3, RIMS2, STOM, SOCS2, CHML
Transcription factor	SOX4, ARID5B, ETV4, ETV1, NR2F1, NR4A2, BHLHE41, CCNL2, SPRY4, HEY1	ZNF462, CCNL1, ID1, ID3
Oxidoreductase	GPX3, IDH3A, ECHDC1, XDH, CD163L1, PTGS2, AKR1C1, AKR1C2, AKR1C3	EHHADH
Hydrolase	PRRG4, GLS, BACE1, PFKFB4, ECE1, CD163L1, MMP1, DPP4	CCNL1, ASNS, USP25
Nucleic acid binding	SAP30, RPS23, ETV4, ETV1, ZFP36L2, NR2F1, ISG15, CCNL2	CCNL1, RIMKLB, MCPH1
Receptor	IFNAR1, MMD, NR2F1, CD163L1, NR4A2, COL6A1, LGR5	GPR87
Calcium-binding protein	S100A4, LTBP2, NID2, CALB1	PRKCA
Membrane traffic protein	SNAP29, SYNGR1, TMED5, SEC24A	Not available
Transfer/carrier protein	PLEKHA3, TMED5, DBI	PRKCA
Extracellular matrix protein	LTBP2, LGR5	Not available
Cytoskeletal protein	MTSS1, EMP3	TPM1, PKP2, KIF23, STOM, ADD2
Cell adhesion molecule	PCAM1, LTBP2	Not available
Defense/immunity protein	IFNAR1, C3	HEATR6, CRISPLD1, OSMR
Transporter	SLC6A6	SLC16A7, SCN8A
Transferase	ECHDC1	WNK4, UGCG, PRKCA, CHML
Lyase	ECHDC1	EHHADH, PCK2
Ligase	ECHDC1	ASNS
Isomerase	ECHDC1	EHHADH
Adaptor protein	Not available	GRB10

Figure 3

Protein classes of DEGs enriched in up-regulated (N = 117) and down-regulated genes (N = 83), respectively. The protein of DEGs is identified using the PANTHER classification system online database and classified according to its function. Annotation dataset was used as PANTHER protein classes, and test type was kept as Fisher's exact test with FDR correction. N, Number of genes.

PPI network analysis

Up-regulated and down-regulated DEGs in CSLCs were used to find protein-protein interactions in the STRING database separately. Out of 117 up-regulated DEGs in CSLCs, only 71 proteins were found in the PPI network (Figure 4A). A total of 71 nodes and 120 edges were found in the network with an average local clustering coefficient of 0.382 and a PPI enrichment p-value of 7.64 × 10−9. Out of 83 down-regulated DEGs in CSLCs, 55 proteins were present in the PPI network (Figure 4B). A total of 55 nodes and 57 edges were found in the network, with an average local clustering coefficient of 0.463 and PPI enrichment p-value of 2.48 × 10−5. Nodes with different colors represent the respective protein classes (Figure 4A, B). Proteins that regulate the interferon-alpha/beta signaling pathway interact with one another and form one cluster in the PPI network of up-regulated DEGs. Oxidoreductases AKR1C1, AKR1C2, and AKR1C3 interact with one another, and CXCR4 was found to interact with many other up-regulated genes. Of the proteins that were down-regulated in CSCs, EGFR was found to interact with a variety of different proteins.

Figure 4

(A) Protein-protein interaction (PPI) network of up-regulated genes (B) Protein-protein interaction network of down-regulated genes. A confidence score higher than 0.4 was kept as the cut-off value for the interaction. Different colors represent different protein classes.

Effect of carboplatin treatment in DEGs in SK-OV-3 cells

To find a potential effect of carboplatin on DEGs expression, we first performed a growth inhibition assay to detect the IC50 of carboplatin in SK-OV-3 cells. Experiments were performed using SK-OV-3 monolayer cells. The IC50 was analyzed using a dose-response curve, and the IC50 level of carboplatin was determined to be 78.55 μM (Figure 5A). Further SK-OV-3 monolayer cells were treated with IC50 concentration of carboplatin at different time intervals and cell morphology was observed at each time point. Intriguingly, we noticed a significant change in cell morphology with increasing treatment time (supplementary figure S1). In addition, we checked the expression of DEGs at different time points. Among the top 13 up-regulated genes, ADM, CXCR4, KCNK3, LGR5, and PTGS2 were significantly increased in time-dependent manner. KCNK3, MT1M, PRDM1, RRAD, and SERPINA1 exhibited an increase in gene expression after 24 h of treatment. Following 24 h, their expression was decreased. AKR1C1 and IFIT2 expressions were decreased after 48 h treatment (Figure 5B). Unlike up-regulated DEGs, gene expressions of down-regulated DEGs were reduced following 48 h of treatment. CCDC80, CENPF, DTNA, EGFR, IL12A, PKP2, PPFIBP1, RIMS2, TPM1, CYR61, ID3, and ZMAT1 gene expression were increased in the first few hours of carboplatin treatment, and reduced after that (Figure 5C).

Figure 5

(A) Dose-response curve of SK-OV-3 monolayer cells to carboplatin (IC50 = 78.55 μM) determined by MTT assay. The IC50 of carboplatin was calculated using GraphPad Prism version 8.0. (B & C) Relative gene expression of up-regulated (B) and down-regulated (C) genes in carboplatin treated SK-OV-3 monolayer cells for the different time intervals by q-PCR. Fold changes of genes were calculated using (ΔΔCT) methods, and 18S gene expression was used to normalize their expression. Results were represented as mean ± SD from triplicates ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 and non-significant (ns) p > 0.05 compared with control cells.

Overall survival (OS) analysis

The association of patient survival rate and all thirteen selected up-regulated and thirteen down-regulated DEG's expression was investigated. Kaplan plot (OncoLnc (www.oncolnc.org)) demonstrated that MMP1 expression was significantly correlated with overall patient survival (Figure 6A; Logrank p-value = 0.0473) (http://www.oncolnc.org/search_results/?q=MMP1). The prognosis of OC patients was closely associated with MMP1 expression levels. High expression of PPIFIB1 was correlated with poor overall survival rate (Figure 6B; Logrank p-value = 0.0756) (http://www.oncolnc.org/search_results/?q=PPFIBP1). The prognosis of OC patient overall survival with high expression of PPIFIBP1 was worse compared to low expression. Few other selected DEGs expression was associated with poor prognosis in other types of cancer. High expression of ADM was associated with a low survival rate in esophageal carcinoma (Logrank p-value = 0.05) and glioblastoma (Logrank p-value = 0.03). Similarly, high expression of CXCR4 was associated with a low survival rate in kidney renal clear cell carcinoma (Logrank p-value = 0.00027).

Figure 6

Association between DEG expression and overall survival. (A) Kaplan plot analysis for overall survival of ovarian cancer patients with high and low expression of MMP1. (Logrank p-value = 0.047) (B) Kaplan plot analysis for overall survival of ovarian cancer patients with high and low expression of PPFIBP1 (Logrank p-value = 0.07).

Discussion

EOC is a principal cause of death among women globally. Despite recent development, current treatments are restricted due to disease recurrence and chemoresistance. Emerging evidence indicates that this is largely due to CSCs, located along with CCs within the bulk tumor. The molecular signatures that characterize CSCs can be exploited to uncover novel therapeutic targets for the treatment of EOC. Analysis of differential gene expression is a commonly used method for identifying abnormally expressed genes in two different conditions. In this present study, we integrated two microarray datasets, and bioinformatics tools to screen DEGs between CCs and CSLCs. A total of 200 DEGs were screened, consisting of 117 up-regulated genes and 83 down-regulated genes. These genes were then evaluated for GO, and KEGG and reactome pathway enrichment analysis. GO analysis disclosed that most of the commonly up-regulated DEGs were enriched in the ECM of cells as well as in biological processes, including cell proliferation, tissue development, and response to lipids, chemical, and organic substances. Down-regulated DEGs were enriched in biological process involved in molecular function regulation. Commonly up-regulated DEGs in CSLCs play a pivotal role in various cell signaling pathways. Genes associated with interferon-alpha/beta signaling (IFIT3, IFIT2, IFIT1, ISG15, IFITM1, IFNAR1), fibers elasticity (BMP4, FBLN1, LTBP2, TGF-β1), and bile acids and bile salts synthesis (AKR1C3, AKR1C2, AKR1C1) were significantly enriched in up-regulated DEGs. Very few DEGs that were down-regulated in CSLCs were enriched in EGFR transactivation by gastrin pathway (PRKCA, EGFR). Understanding of the molecular signatures and signaling pathways that characterize CSLCs will advance the development of targeted therapy.

On further examination of DEGs, we selected the top thirteen commonly up-regulated (ADM, AKR1C1, CXCR4, IFIT2, KCNK3, LGR5, MMP1, MT1M, PRDM1, PTGS2, RRAD, S100A4, SERPINA1) and thirteen down-regulated DEGs (CCDC80, CENPF, CHRDL1, CYR61, DTNA, EGFR, ID3, IL12A, PKP2, PPFIBP1, RIMS2 TPM1, ZMAT1) in CSCs based on fold change values in both datasets. The identification of such significant findings provides an opportunity to examine for common operational genes affected by two subpopulations, CCs and CSCs. These DEGs may prove to be specific therapeutics targets. Moreover, they could also be used as a diagnostic marker and for the ultimate management of OC.

ADM is a multi-functional peptide involved in cellular proliferation, tumor angiogenesis, and inhibition of programmed cell death. Previous studies have reported that ADM interacts with VEGF and increases its expression to promote angiogenesis in EOC [24]. VEGF-mediated signaling has known to regulate the self-renewal and survival of CSCs [25]. The VEGFR/ADM axis might be a potential marker and therapeutic target for ovarian CSCs. IFIT2 (interferon-induced protein with tetratricopeptide repeats) is an interferon's stimulated gene. A recent study has identified its involvement in drug and radiation resistance through transcriptome analysis. In breast cancer, lower IFIT2 expression is associated with poor prognosis [26]. Baicalein, a traditional medicinal herb, suppresses stem-like characteristics in radio and chemo-resistant breast cancer cells by up-regulating IFIT2 expression [26]. Given its pro-tumorigenic role in breast cancer, IFIT2 likely has a similar role in other types of cancer including EOC. Therefore, additional studies investigating its specific role in EOC might establish the potential therapeutic effects of IFIT2 in ovarian CSCs. AKR1C1 belongs to Aldo-keto reductase enzyme family, known to catalyze NADPH dependent reductions and have an essential role in biosynthesis, metabolism, and detoxification [27]. In head and neck squamous cell carcinoma (HNSCC), AKR1C1 expression is a poor prognostic and recurrent biomarker, and a crucial regulator of cisplatin-resistance [28]. AKR1C1 might be a useful predictive factor for EOC patients treated with platinum-based therapy. CXCR4 is a chemokine receptor, belongs to the serpentine family of the G-protein coupled receptor, found up-regulated in high-grade serous EOC, and promotes metastasis [29]. CXCR4 and its ligand SDF-1 (CXCL12) have been demonstrated to be functional in CSCs subpopulations of various types of cancer [30]. Our findings also indicate its potential use as a biomarker of ovarian CSCs. CXCR4-CXCL12 axis may be a specific therapeutic target of ovarian CSCs. KCNK3 gene encodes K2P3.1, which belongs to the potassium channel superfamily proteins. Reduced expression of KCNK3 has been observed in many cancer types and linked to alteration in membrane potentials. Altered membrane potentials are known to induce cellular proliferation, migration, and apoptosis [31]. LGR5 is a member of glycoprotein hormone receptor sub-family, which has been demonstrated to regulate tumorigenesis and metastasis in EOC through Notch1 signaling [32]. The notch pathway is known to regulate the self-renewal and survival of CSCs. It is an attractive target for treatment as they eliminated both CCs and CSCs [33]. LGR5-Notch1 could be targeted as potential therapeutic strategies for EOC management. LGR5 promotes CSCs traits and chemoresistance through Wnt/β-catenin signaling pathway in cervical cancer [34]. LGR5 up-regulation is correlated with OC stemness properties and conferring chemoresistance. In osteosarcoma, MMP1 up-regulation is associated with tumor progression, metastasis, and cancer stem-like properties [35]. MMPs are involved in the breakdown of ECM in a normal physiological process as well as in diseased state. In our study, MMP1 expression was significantly correlated with OC survival and predicted prognosis. The mechanism underlying its up-regulation in ovarian CSCs and possible links with chemoresistance warrant further investigation. MT1M expression is involved in thyroid cancer progression and acts as tumor suppressor [36]. Our study indicates the potential role of MT1M in biological processes. PRDM1 has been associated with cancer progression, metastasis, and altered expression correlated with poor prognosis in lung cancer [37]. In ovarian CSCs, PRDM1 might be used as prognostic tool and serve as a new therapeutic target. PTGS2 is one of the rate-limiting enzymes of synthesis of prostaglandins from arachidonic acid. In radiation-resistant glioma cells, PTGS2 expression is up-regulated and increased radio-tolerance by activating downstream signaling effectors NF-κB signaling [38]. NF-κB signaling is known to maintain the characteristics of CSCs associated with cancer progression. After cisplatin treatment, the dormant cells displayed an increased expression of cells with CSCs markers in OC, indicating the resistance of this cell fraction [39]. Targeting the PTGS2-NF-κB axis in CSCs will provide new insight into OC therapy. In lung cancer, RRAD binds with the p65 subunit of NF- κB and inhibits its nuclear translocation and consequently, inhibits GLUT1 transport to the plasma membrane [40]. In bladder cancer stem cells, S100A4 is up-regulated and associated with cellular proliferation by binding to IKK and activating the NF- κB signaling pathway [41]. Thus, it can also be implicated in the characteristics of ovarian CSCs associated with tumor progression. SERPINA1 is a protease inhibitor that inhibits serum trypsin and is generally found to be higher in serum concentration of cancer patients than healthy persons. SERPINA1 is associated with distant metastasis of various cancers and correlated with lymph nodes in colorectal cancer [42, 43]. Our study identified that SERPINA1 is up-regulated in OC, suggesting that it could serve as a biomarker.

Abnormal CCDC80 expression inhibits migration and growth proliferation of melanoma cells [44]. CENPF is a cell cycle-associated nuclear protein that is over-expressed in breast and lung cancer and is correlated with poor prognosis [45]. CHRDL1 is a secretory protein and functions as an antagonist of bone morphogenetic protein (BMP). In gastric cancer, hypermethylation of CHRDL1 promoter induces low expression of CHRDL1 and promotes proliferation and metastasis [46]. Cyr61 is a matricellular protein present in ECM. In pancreatic cancer cells, Cyr61 is expressed in CSCs and activates PI3k signaling to increase chemotherapy resistance [47]. Our study found that Cyr61 expression was reduced in CSLCs compared to CCs. DTNA expression has been found to increase in hepatitis B-virus (HBV)-induce hepatocellular carcinoma cells and promote cellular proliferation and inhibit apoptosis. DTNA binds with STAT3 and induces TGF-β expression, and reduces p53 expression [48]. This gene might be associated with inhibition of apoptosis in OC. In HNSCC, EGFR overexpression is correlated with poor prognosis. It contributes to the development of CSCs by stabilizing Sox2. Sox2 is an established stem cell marker that serves as a substrate for EGFR. Inhibition of EGFR signaling has been shown to reduce Sox2 expression [49]. Targeting EGFR with an existing therapy might be a potential strategy for EOC treatment through the eradication of CSCs. In intrahepatic cholangiocarcinoma, ID3 promotes CSCs by activating β-catenin and is correlated with poor prognosis [50]. This gene expression may be associated with OC prognosis. IL12 has been known to inhibit the survival of CSCs [51]. PKP2 is a regulator of EGFR, and increased PKP2 activates EGFR signaling and increased cellular proliferation and migration of CCs [52]. PKP2-EGFR might be a potential target for OC treatment. The reduced expression of miR-194-5p promotes IGF1R and PPFIBP expression, promoting OC progression. NF-κB binds to the miR-194-5p promoter and negatively regulates miR-194-5p expression [53]. TPM1 overexpression has been demonstrated to promote cell apoptosis, inhibit migration, and its expression is correlated to OSCC prognosis [54]. TPM1 expression might be linked to OC prognosis.

Biological pathways such as interferon-alpha/beta signaling were enriched in ovarian CSLCs. Activation of interferon-alpha signaling is known to contribute to cellular senescence, cell death, and attributed to increasing migratory cell rate and drug resistance depending on the genes transcribed after interferon-stimulation [55]. Type I Interferon (IFN-1) has been demonstrated to control tumor initiation, progression, and immune surveillance [56]. Previous studies have shown that the impairment of IFN-1 signaling alters the expression of CSCs marker ALDH and affects breast CSCs population in tumors [57]. In oral squamous cell carcinoma, IFN-α has shown to activate the transcriptional expression of CSCs marker CD44 and ALDH1A1, and IFN-α-primed enhances cytotoxic inhibition of multiple therapeutic drugs [58]. Combinational treatment of IFN-α and IFN-γ inhibit growth and induces OC cells [59]. IFN-α signaling is mediated through a receptor complex of interferon-alpha/beta receptor-1 (IFNAR1) and IFNAR2, and these two subunits are significantly associated with an immunosuppressive environment in cancer [60]. IFITs are interferon's stimulated genes and identified as IFIT1, IFIT2, IFIT3, and IFIT5 in humans [61, 62]. These genes may be associated with OC prognosis. The Interferon-stimulated gene (ISG15) is an ubiquitin-like protein and is a known prognostic marker in OC. While ISG15 plays an equivocal role in cancer progression and treatment response to solid tumors, its underlying mechanism is unexplored [63]. Higher expression of IFIT2 in ovarian spheroids in both datasets indicates its potential roles in OCSCs. This gene might be used as a diagnostic marker of ovarian CSCs.

Molecules associated with elastic fibers were enriched in ovarian CSCs. After secretion from cells, TGF-β1 is deposited in the ECM with its LAP. TGF-β1 is critical for the proper mobilization of the latent cytokine and its activation. TGF-β1 binds and forms a complex with LTBP isoform -1, -3, and -4. LTBP2 is missing the binding site for LAP is secreted without latent TGF-β. LTBP2 has been demonstrated to manifest tumor-suppressing or tumor-promoting functions in multiple cancers, including OC [64, 65]. Bone morphogenetic protein 4 (BMP4) is a family of signaling molecules that belong to TGF-β super-family. Higher expression of BMP4 in NSCLC patients is associated with therapeutic resistance [66]. In gastric cancer, BMP4 is identified as a modulator of cisplatin sensitivity [67]. BMP4 is known to regulate ovarian tumor microenvironment and cancer progression. A stem cell factor, Lin28, interacts with BMP4, and promotes BMP4 expression in EOC and forms a complex with Oct4 to regulate cell proliferation [68]. Fibulin-1 (FBLN-1) is a member of the extracellular glycoprotein family, binds to ECM proteins, and promotes tumor progression. In bladder cancer, Fibulin-1 deregulation is associated with angiogenesis, migration, and disease recurrence. This gene expression may be correlated with EOC progression and disease recurrence.

Bile acid and salts such as AKR1C3, AKR1C2, and AKR1C1 are known to catalyze NADPH dependent reductions and have essential roles in biosynthesis, metabolism, and detoxification [27]. In colon cancer, inhibition of AKR1C3 and AKR1C1 is associated with increased cisplatin sensitivity [69], and overexpression of AKR1C1 has been shown to induce cisplatin resistance in OC cells [70]. This gene may be a potential therapeutic target for OC.

Comparative expression analysis of CSC markers between monolayer and spheroid demonstrated that spheroids are enriched in CSCs. Validation of the top thirteen up-regulated and down-regulated genes was examined by q-PCR and demonstrated robust correlations between both expression data. CCs are frequently altered at the cellular and molecular levels and acquire resistance to chemotherapeutics. Carboplatin is among the leading treatments for various cancer types, including OC. Carboplatin treatment induces molecular alterations in DEGs. Gene expression study of DEGs in SK-OV-3 cells treated with carboplatin revealed that up-regulated DEGs such as ADM, CXCR4, LGR5, and PTGS2 might be potential candidates for drug-resistant in OC. The silencing of the ADM gene has shown to inhibit proliferation and increases the chemo-sensitivity of the OC cells [71]. CXCR4 is a critical molecule in cisplatin treatment for EOC patients and its inhibition can be a potential strategy to address chemoresistance [72]. LGR5 promotes ovarian cell proliferation, progression, and EMT [32] and is associated with chemoresistance in cervical cancer [34]. In bladder cancer, PAX5/PTGS2 cascade is known to induce cisplatin resistance [73]. MMP-1 and PPFIBP1 correlate with OC patient survival. Overexpression of MMP-1 in osteosarcoma cells contribute to cancer progression, metastasis, and stem-like properties [35]. In addition, a recent study has demonstrated that PPFIBP1 contributes to the tumorigenesis of OC [53]. MMP1 and PPIFIB1 expression have decreased with increasing drug time exposure that indicates MMP1 and PPIFIB1 expression might be associated with carboplatin sensitivity. In carcinoma-associated fibroblasts (CAFs), drug sensitivity was increased when MMP1 expression was inhibited, chemotherapy-induced collagen IV over-expression, and protects cancer cells. Cui et al. [74] demonstrated that GM6001 (a specific MMP inhibitor) treatment induces reduced expression of collagen IV as well as MMP-1. Thus, Collagen IV reduces the therapeutic effect of Taxotere on breast cancer via TGF-β signaling [74]. Gemcitabine treatment induces higher expression of collagen IV and promotes OC cells resistance to this drug [75]. Collagen is an ECM component that interacts with the integrin receptor of cancer cells and may induce cell-adhesion-mediated drug resistance and reduce the chemotherapy effects [74]. Carboplatin may suppress tumor progression via down-regulating MMP1 expression. However, MMP1 is a known tumor promoter as well as a tumor suppressor. Further studies are needed to understand their role in carboplatin sensitivity.

In conclusion, comprehensive data mining and bioinformatics analysis of DEGs that may affect both CCs and CSLCs population in OC was conducted. A total of 200 DEGs were screened. ADM, AKR1C1, CXCR4, IFIT2, KCNK3, LGR5, MMP1, MT1M, PRDM1, PTGS2, RRAD, S100A4, and SERPINA1 were top thirteen up-regulated. CCDC80, CENPF, CHRDL1, CYR61, DTNA, EGFR, ID3, IL12A, PKP2, PPFIBP1, RIMS2 TPM1, ZMAT1 were top thirteen down-regulated genes commonly identified in both dataset. Furthermore, the expression of all DEGs was validated by q-PCR. The most significant enrichments are interferon-alpha/beta signaling, molecules associated with elastic fibers, and bile acids and bile salts synthesis. Also, gene expression analysis of DEGs was performed following carboplatin treatment in SK-OV-3 cells, and higher expressions of ADM, CXCR4, LGR5, and PTGS2 were noted. The patients with high expression levels of PPFIBP1 and low levels of MMP1 expression exhibited worse prognosis. Overall this work highlights the functional difference caused by DEGs that may affect CCs and CSLCs. However, a further understanding could be useful for designing effective therapeutic strategies to increase overall patient survival.

Study limitations

In this study, rigorous bioinformatics analysis was performed, however, there were shortcomings. The sample size in the microarray data set was small, and further increment of sample size is required to attain more accuracy. This study lacks clinical samples and relevant animal experiments to conclusively verify the role of DEGs in ovarian CSCs maintenance. Mechanistic understanding of these up-regulated and down-regulated genes in OC recurrence and drug resistance is needed. In addition, validation of the DEGs was performed only in SK-OV-3 cells. Follow-up experiments will be performed after acquiring more ovarian cell lines for both monolayer and spheroids. We will also expand our investigations into sorted CSCs instead of using spheroids only.

Conclusion

Bioinformatics analysis identified the molecular signatures and signaling pathways that are enriched in ovarian CSLCs. DEGs were found between ovarian CCs and CSLCs, which could help in the development of therapeutic strategies to target both cell populations. These DEGs may be used as targets for specific treatment and could provide further research idea to find the new mechanism and potential therapeutic targets for ovarian CSCs.

Availability of data and materials

The data and materials used to support the findings of this study are available from the corresponding author upon request.

Declarations

Author contribution statement

S. Kumar: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

A. Behera: Performed the experiments; Wrote the paper.

R. Ashraf: Performed the experiments.

A. Srivastava: Analyzed and interpreted the data.

Funding statement

This work was supported by Department of Biotechnology, Government of India (IN) (BT/RLF/Re-entry/13/2016), and (CRG/2019/002104).

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.