Gene expression analysis of MCF7 cell lines of breast cancer treated with herbal extract of Cissampelos pareira revealed association with viral diseases.

BACKGROUND: It is necessary to assess the cellular, molecular, and pathogenetic characteristics of COVID-19 and attention is required to understand highly effective gene targets and mechanisms. In this study, we suggest understandings into the fundamental pathogenesis of COVID-19 through gene expression analyses using the microarray data set GSE156445 publicly reachable at NIH/NCBI Gene Expression Omnibus database. The data set consists of MCF7 which is a human breast cancer cell line with estrogen, progesterone and glucocorticoid receptors. The cell lines treated with different quantities of Cissampelos pareira (Cipa). Cipa is a traditional medicinal plant which would possess an antiviral potency in preventing viral diseases such as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection.
METHODS: Utilizing Biobase, GEOquery, gplots packages in R studio, the differentially expressed genes (DEGs) were identified. The gene ontology (GO) of pathway enrichments employed by utilizing DAVID and KEGG enrichment analyses were studied. We further constructed a human protein-protein interaction (PPI) network and performed, based upon that, a subnetwork module analysis for significant signaling pathways.
RESULTS: The study identified 418 differentially expressed genes (DEGs) using bioinformatics tools. The gene ontology of pathway enrichments employed by GO and KEGG enrichment analyses of down-regulated and up-regulated DEGs were studied. Gene expression analysis utilizing gene ontology and KEGG results uncovered biological and signaling pathways such as "cell adhesion molecules", "plasma membrane adhesion molecules", "synapse assembly", and "Interleukin-3-mediated signaling" which are mostly linked to COVID-19. Our results provide in silico evidence for candidate genes which are vital for the inhibition, adhesion, and encoding cytokine protein including LYN, IGFBP5, IL-1R1, and IL-13RA1 that may have strong biomarker potential for infectious diseases such as COVID-19 related therapy targets.

Introduction

COVID-19 originated from China and is the cause of SARS-CoV-2 related to the human respiratory system. As of February 2021, a total of 104 million confirmed cases and 2.27 million deaths were reported worldwide, according to World Health Organization (WHO) estimation. While vaccines are developed and implemented at breakneck speed, the timeline for large-scale production and vaccination of the world's enormous population stays vague (McCallum et al., 2021). Moreover, although the primary endpoints are the vaccine's effectiveness and safety against laboratory-approved COVID-19, vaccine trials have not been strengthened to definitively evaluate efficacy by subgroup of any age with certain known or unknown medical conditions such as cancer, chronic kidney disease, COPD (chronic obstructive pulmonary disease), heart failure, coronary artery disease, or cardiomyopathies, sickle cell disease, smoking, type 2 diabetes (Garassino et al., 2020; Polack et al., 2020; Wu and McGoogan, 2020).

COVID-19 is the first global pandemic that humanity has encountered in this century, and the character of defensive immune responses is not fully established; it is still indefinite which vaccine strategy will be most effective (Jeyanathan et al., 2020). The success of any vaccine strategy on a population can only be observed if the vaccine is effective in high risk groups and is distributed extensively to the whole population (Hodgson et al., 2020). The above are all hazardous circumstances that can cause to severe COVID-19 infection and potentially lead to a worse prognosis and increased risk of death in risk groups. Older cancer patients had poor prognosis from COVID-19 (Chen et al., 2020; Wang and Zhang, 2020).

It is commonly considered that the fundamental relation between genetic and environmental factors is the reason for breast cancer (Brody and Rudel, 2003). Breast cancer is more common among women. In 2008, every minute 2.6 women are diagnosed with breast cancer and more than 52 women die every hour (Siegel et al., 2018). This study analyzed genome-wide expression profiling in conjunction with gene ontology (GO) enrichment analysis and protein-protein interaction (PPI) to provide us with more information MCF7 cell lines of breast cancer. Cipa transcription signatures overlap with related biological, molecular, and functional pathways and priority compound family predicted for COVID-19 from various sources (Amresh et al., 2007; Haider et al., 2021; Singh et al., 2013) Cipa is one such medicinal plant traditionally used to treat hormonal imbalance, fever, cough, skin ulcer, jaundice, and asthma. The most recent study suggested that herbal extracts of Cipa have antiviral and antimalarial properties. Biological and molecular pathways linked to lipid metabolism, viral transcription, translation, and estrogen appears to be regulated and altered by CIPA (Haider et al., 2021). The necessity for the identification of original targets and effective treatments is of paramount importance to defeat this pandemic and to improve survival rates among cancer patients is inevitable.

This study addresses cancer patients in particular, breast cancer patients in the high-risk group who need an urgent vaccination strategy. One of the aims of this study is deciphering hub genes which perform critical functions in cytokine production in breast cancer. There are various cytokines that controls the inflammatory tumor microenvironment. Tumor progression is advanced by Interleukin (IL)-1, IL-6, IL-11, and transforming growth factor-β (TGF-β) that induces tumorigenesis, and penetration and cytokine receptor trigger. TGF-β is one of the most investigated cytokines stimulating cancer cell expansion, progression, and invasion (Nicolini et al., 2006). One of the recent results showed IL-1R1 signaling represses breast cancer cell growth early in tumor formation, and reduced breast cancer growth and pulmonary invasion (Dagenais et al., 2017). On the other hand, IL-1R blockade is suggested in cases with severe COVID-19 and hematological malignancies since the IL-1/IL-6 pathway is exceedingly regulated in patients with serious infection (van de Veerdonk and Netea, 2020). Therefore, researching the unique properties of SARS-CoV-2 in breast cancer patients is essential in developing treatments to improve patient outcomes.

Consequently, such therapies can be beneficial in the case of unvaccinated persons or poor responders to vaccination, as well as when immunity is reduced or antigenically different strains appear (McCallum et al., 2021). Presently, with genomic data such as The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) available to the public, many bioinformaticians studied gene expression data with clinical data to try to predict prognosis and find biomarkers for treatment (Chin et al., 2011). Microarrays have become a significant tool in investigating cancer genes and targeting therapeutic drugs. Extracting these gene expression patterns and demonstrating their molecular mechanisms and signaling pathway degrees will assist clinicians and researchers to subgroup and handle COVID-19 more precisely for breast cancer patients.

The analysis procedure proceeds by detecting DEGs from which the functions and pathways of the genes involved are studied. A variety of research in breast cancer have been performed previously (AbuHammad and Zihlif, 2013; Deng et al., 2019). Even though most of this research only emphasized the detection of hub genes in breast cancer MCF7 cell lines, examining Cipa transcription signatures with related molecular pathways and priority compound family prediction for COVID-19 would be the first among other studies.

The present study seeks to detect the DEGs and to conduct Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis carried out to explore the biological functions and pathways. Consequently, this study aim to investigate whether there is a detectable link between MCF7 cell lines treated with Cipa and COVID-19 at the gene expression level. A human protein-protein interaction (PPI) network was constructed and a module analysis was studied to find the significant genes in the data set. Our study offers original insight of the mechanism of infectious diseases such as COVID-19 formation in breast cancer MCF7 related hub genes. The related pathways further may function as promising targets for therapy of breast cancer cases with COVID-19.

Materials & methods

Preprocessing of the data set

Our overall study design is demonstrated in Fig. 1 . Affymetrix microarray data of the gene expression profile of GSE156445 was obtained from the NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/) which is based on the Affymetrix Human Transcriptome Array 2.0. Here, 3 replicates of experiments of vehicle treated cells replicates of samples (GSM4731601, GSM4731602, GSM4731603), and 15 cells treated with Cissampelos pareira replicates of 3 experiments as the following quantities; 1 μg (GSM4731604, GSM4731605, GSM4731606), 10 μg (GSM4731607, GSM4731608, GSM4731609), 100 μg (GSM4731610, GSM4731611, GSM4731612), 500 μg (GSM4731613, GSM4731614, GSM4731615), and 1000 μg (GSM4731616, GSM4731617, GSM4731618) samples. Genomic information ranging from cell lines of gene expression data was obtained. Fig. 2 demonstrates uniform manifold approximate and projection (UMAP) GSE154245 data set for quality control (QC) (Diaz-Papkovich et al., 2021).

Fig. 1

Design and workflow of the computation steps. DEGs, differentially expressed genes; p, p value; FC, fold change.

Fig. 2

Uniform manifold approximate and projection (UMAP) of the raw GSE154246 microarray data set to explore structure in high-dimensional data sets to perform quality control (QC).

The GSE156445 data set was analyzed by using the GEOquery package in Bioconductor following standard procedures in R studio (Durinck et al., 2009, Durinck et al., 2005; Tarca et al., 2006; Warnes et al., 2009). The other packages are used in R studio are as the following; Biobase, biomaRT, umap and gplots packages (Davis and Meltzer, 2007; Durinck et al., 2005; Konopka and Konopka, 2018; Warnes et al., 2009). To estimate the adjusted p value and avoid Type I errors, the Benjamini–Hochberg Procedure is used to correct multiple testing. In order to adapt the statistical tests locally, a hypergeometric model was used for both of the down-regulated and up-regulated DEGs in the functional GO and pathway enrichment analysis, and computing a false discovery rate (FDR) were computed (Benjamini and Hochberg, 1995; Dudoit et al., 2003; Hochberg and Tamhane, 1987).

Experimental data and analysis codes

Analyses were conducted using the R statistical environment. Sample codes and analysis of the data can be found in a GitHub repository. Prior to the analyses, the low quality reads and genes with very low number of reads were removed and gene expression values converted to base-2 logarithmic scale using R language. Samples are separated into six groups: vehicle treated cells and cells treated with 1 μg, 10 μg, 100 μg, 500 μg, 1000 μg Cipa. We then employed a one way ANOVA test to normalize the gene expression data set and supply predictions of alterations in expression values that are adjusted for probability of extraneous determinant effects between each category (Kerr et al., 2000). The outcomes were stated as the means±standard deviation (S.D.). Once the significant differences of variances were found using the anova() function, statistical differences between categories of the data samples were determined utilizing the TukeyHSD() function to perform Tukey's post hoc test which are available in base R programming. Differences in mean levels were taken significant at the level of 5%.

We further used the transformed data set and separated into two groups provided that vehicle treated cells and cells treated with Cipa. The data set was normalized by computing the means of the samples of each group in R programming language. The analysis on separated samples which is grouped by categories was performed as computing fold-change (biological significance) difference between the means of the vehicle treated and treated groups. A broadly performed statistical model is the t-distribution and its variants. A t-test compares the discrepancy of the average gene expression levels between the two samples or subgroups, given the noisiness of the data, i.e., the discrepancy in means between samples divided by the standard deviation. The genes are filtered in accordance with both fold change and p value criterion. Despite the fact that methods correcting for multiple comparisons have been applicable for a long time such as the Bonferroni correction, most of these methods are not appropriate to analyze gene expression data sets (Tarca et al., 2006). We highlight statistical significance performing t-test by taking p value cut-off <0.05 and |log2(FC)| > 1 to identify down and up-regulated DEGs between each group understudy.

Differentially expressed genes and clustering analysis

Using the GEOquery package in Bioconductor, gene expression values were pulled out for each sample and converted to base-2 logarithmic scale using R language. We used the gplots package of R to create heatmaps of DEGs using the heatmap.2() function. Clustering analysis of DEGs was performed to compare the expression pattern of DEGs between vehicle treated cells and cells treated with Cipa.

Gene ontology (GO) terms and analysis of the pathway enrichments

Expression measurement annotations for up-regulated and down regulated DEGs for each group probe were mapped to gene names using the Ensemble Biomart package in R. All of the DEGs were characterized by their biological processes, molecular functions, and cellular components of GO enrichments were studied using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang et al., 2007; Sherman & Lempicki, 2009). All classified genes were examined and further parts like the Universal Protein resource, physical properties like GO, and annotation types were taken using DAVID and KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa et al., 2016).

Human protein-protein interaction (PPI) network of DEGs

NetworkAnalyst (Zhou et al., 2019), publicly accessible on the web, provides analysis of PPI networks for single gene lists using STRING Interactome (Szklarczyk et al., 2016). To comprehensively decipher the regulatory mechanisms in the gene expression data set of MCF7 cell lines of breast cancer patients, DEGs from vehicle treated cells and cells treated with Cipa were analyzed to form a PPI network with previously reported GO classifications and enrichments.

Results

Experimental data analysis

In Fig. 3A is a boxplot of the non-normalized base-2 logarithmic scale of raw GSE156445 gene expression data with samples of vehicle treated and treated with Cipa. We detected DEGs in a total of 418 genes from vehicle treated cells and cells treated with Cipa groups which was demonstrated in the volcano plot (Fig. 3B). The volcano plot indicated that compared to vehicle treated, there were less upregulated genes than downregulated genes in treated with Cipa MCF7 cancer cell lines. We identified 202 down-regulated and 216 up-regulated DEGs, Fig. 4 is displaying top 30 down-regulated and up-regulated DEGs of each.

Fig. 3

(A) Boxplot showing gene expression of each sample of the base-2 logarithmic scale GSE156445 gene expression data. (B) Volcano plot of comparing all of the significant DEGs from vehicle treated cells and cells treated with Cipa samples. The red and blue dots show DEGs that were significant at p value < 0.05 and |log2(FC)| > 1.

Fig. 4

Heat map of the top 60 differentially expressed genes (DEGs) of vehicle treated and treated with Cipa samples. The heat map presents 30 upregulated genes (red) and 30 downregulated genes (green). Base-2 logarithmic values of the gene expression data are calculated. Each column represents samples of vehicle treated and treated with Cipa, and each row presents genes. The progressive color changing from red to green represents the ranging from up to down-regulated DEGs.

GO enrichment analysis and KEGG pathway analysis

Gene expression studies using GO enrichment analysis on MCF7 cell lines treated with Cipa extract reveal first set of descriptors in base-2 logarithmic scale of GSE156445 is shown in Fig. 5 . Table 2A demonstrates the significant enrichments of down-regulated DEGs using a biological processes (BP) enzyme linked receptor protein signaling pathway, cellular protein modification process, protein modification process, and gene expression. The significant enrichment of down-regulated DEGs in cellular component (CC) contains RNA binding, poly(A) RNA binding, and enzyme inhibitor activity. Finally, the significant enrichments GO terms in molecular function (MF) is revealed mitochondrial membrane, nucleoplasm, and mitochondrial envelope. KEGG pathway study outcomes illustrated that the down-regulated DEGs were considerably enriched in hsa03013:RNA transport.

Fig. 5

The significant enrichments of DEGs using (A) KEGG pathways, gene ontology (GO) enrichments of (B) biological processes (BP), (C) cellular component (CC), and (D) molecular function (MF).

Table 2

GO functional enrichment analyses of differentially expressed genes of MCF7 cell line of breast cancer associated with COVID-19.

Category	Term/gene function	Gene count	%	p-Value
A, Down-regulated
GOTERM_BP_FAT	GO:0007167 ~ enzyme linked receptor protein signaling pathway	24	0.09330534	6.49E−06
GOTERM_BP_FAT	GO:0006464 ~ cellular protein modification process	57	0.22160019	9.28E−06
GOTERM_BP_FAT	GO:0036211 ~ protein modification process	57	0.22160019	9.28E−06
GOTERM_BP_FAT	GO:0010467 ~ gene expression	71	0.2760283	1.47E−05
GOTERM_BP_FAT	GO:0043624 ~ cellular protein complex disassembly	9	0.0349895	2.83E−04
GOTERM_MF_FAT	GO:0003723 ~ RNA binding	37	0.14384574	5.80E−07
GOTERM_MF_FAT	GO:0044822 ~ poly(A) RNA binding	28	0.10885623	1.04E−05
GOTERM_MF_FAT	GO:0004857 ~ enzyme inhibitor activity	13	0.05054039	2.47E−04
GOTERM_MF_FAT	GO:0003676 ~ nucleic acid binding	54	0.20993702	0.00284797
GOTERM_MF_FAT	GO:1901363 ~ heterocyclic compound binding	72	0.27991603	0.00354973
GOTERM_CC_FAT	GO:0031966 ~ mitochondrial membrane	16	0.06220356	0.00102392
GOTERM_CC_FAT	GO:0005654 ~ nucleoplasm	42	0.16328435	0.00133637
GOTERM_CC_FAT	GO:0005740 ~ mitochondrial envelope	16	0.06220356	0.00181153
GOTERM_CC_FAT	GO:0030864 ~ cortical actin cytoskeleton	5	0.01943861	0.00279154
GOTERM_CC_FAT	GO:0031967 ~ organelle envelope	20	0.07775445	0.00514539
KEGG_PATHWAY	hsa03013:RNA transport	5	0.01943861	0.05072706
B, Up-regulated
GOTERM_BP_FAT	GO:0016339 ~ calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules	8	0.04620005	9.84E−10
GOTERM_BP_FAT	GO:0007156 ~ homophilic cell adhesion via plasma membrane adhesion molecules	10	0.05775006	1.35E−06
GOTERM_BP_FAT	GO:0007416 ~ synapse assembly	9	0.05197505	2.86E−06
GOTERM_BP_FAT	GO:0098742 ~ cell-cell adhesion via plasma-membrane adhesion molecules	10	0.05775006	2.60E−05
GOTERM_BP_FAT	GO:0050808 ~ synapse organization	9	0.05197505	1.78E−04
GOTERM_MF_FAT	GO:0005509 ~ calcium ion binding	16	0.09240009	1.19E−04
GOTERM_MF_FAT	GO:0015459 ~ potassium channel regulator activity	4	0.02310002	0.00417623
GOTERM_MF_FAT	GO:0004602 ~ glutathione peroxidase activity	3	0.01732502	0.00912591
GOTERM_MF_FAT	GO:0016247 ~ channel regulator activity	5	0.02887503	0.01431307
GOTERM_MF_FAT	GO:0017080 ~ sodium channel regulator activity	3	0.01732502	0.02428745
GOTERM_CC_FAT	GO:0031988 ~ membrane-bounded vesicle	33	0.19057519	0.03817299
GOTERM_CC_FAT	GO:0070062 ~ extracellular exosome	27	0.15592516	0.04033309
GOTERM_CC_FAT	GO:1903561 ~ extracellular vesicle	27	0.15592516	0.04255514
GOTERM_CC_FAT	GO:0043230 ~ extracellular organelle	27	0.15592516	0.04271734
GOTERM_CC_FAT	GO:0031226 ~ intrinsic component of plasma membrane	18	0.1039501	0.05286412
KEGG_PATHWAY	hsa00480:Glutathione metabolism	4	0.02310002	0.00622493
KEGG_PATHWAY	hsa04146:Peroxisome	4	0.02310002	0.02328098
KEGG_PATHWAY	hsa04142:Lysosome	4	0.02310002	0.06008537

BP, biological processes; FAT, functional annotation tool; MF, molecular function; CC, cellular component.

Similarly, Table 2B demonstrates the significant enrichments of up-regulated DEGs using biological processes (BP) calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules, homophilic cell adhesion via plasma membrane adhesion molecules, and synapse assembly. The significant enrichment of up-regulated DEGs in molecular function (MF) contains calcium ion binding, potassium channel regulator activity, and glutathione peroxidase activity. Finally, the significant enrichments GO terms in cellular component (CC) is revealed membrane-bounded vesicle, extracellular exosome, and extracellular vesicle. KEGG signaling pathway study outcomes demonstrated that the up-regulated DEGs were significantly enriched in hsa00480:glutathione metabolism, hsa04146:peroxisome, and hsa04142:lysosome.

Protein-protein interaction (PPI) network

The human PPI network was built to identify the most significant biological components and proteins and that may be key functions in the progress of in breast cancer, in accordance with the predicted interactions of detected DEGs. 1219 nodes and 10,553 edges were filtered from the PPI network in total (Fig. 6 ). 10 core nodes with the highest degrees in MCF7 cell lines gene expression data set DEGs are selected in Table 3. These were RPL7L1, NIFK, LYN, MRPS5, MRPL1, IRS1, MET, CASP8, HCK, and TNKS. Genes are given in a sequence showing the number of neighboring nodes and changes in proteins/genes in the PPI network. The PPI network hub genes had the following cytokine-related genes and receptors which are also key regulators of COVID-19; IL-13RA1, IL-4R, IL-13, IL4, IL-6ST, CXCR4, TNFSF10, TNFRSF10A, TNFRSF1A, TNFRSF10B, CSF2RA, and CSF2RB.

Fig. 6

Human PPI network of DEGs of gene expression in MCF7 cell lines of breast cancer identified by NetworkAnalyst. The network displayed is “continent” (Subnetwork 1). The colors represent the expressions of nodes. Specifically, “green” and “purple” indicate the nodes are up- and down-regulated, respectively. The progressive color changing represents the expression levels. The areas of the nodes indicate the degrees that the nodes connect to others. The node sizes represent ranking of significant genes in terms of degree centrality i.e., the greater quantity of neighbors a node has.

Table 3

Top 15 core genes of PIP network of DEGs of with MCF7 cell line of breast cancer.

Gene ID	Genes	Node degree	Betweenness centrality	Expression
285855	RPL7L1	98	133,959.29	2.91
84365	NIFK	93	87,919.46	3.09
4067	LYN	91	149,289.72	2.97
64969	MRPS5	77	28,828.99	3.16
65008	MRPL1	70	19,893.02	2.75
3667	IRS1	59	110,871.06	3.10
4233	MET	47	72,644.15	2.88
841	CASP8	47	47,694.94	3.04
3055	HCK	45	14,205.28	2.52
8658	TNKS	41	44,879.61	2.70
1399	CRKL	39	39,138.52	3.04
331	XIAP	38	104,019.94	2.76
84135	UTP15	37	14,354.31	2.79
56902	PNO1	32	11,100.54	2.84
7316	UBC	30	430,261.22	2.72

Module analysis of the Interleukin-3-mediated signaling pathway

The core of this project due to it is close association with breast cancer within all of the significantly (p value < 0.05) enriched pathways of DEGs is the Interleukin-3-mediated signaling pathway (GO:0038156 in QuickGO) in which IL-3 was noted as predictive marker related to severity of SARS-CoV-2 (Bénard et al., 2020). There were 29 DEGs particularly engaged in this pathway, containing IL-1R1, IL-13RA1, FOXC1, TRIO, PPM1L, IRS1, XIAP, PTPRK, RASSF2, FNTA, FLRT3, HES1, SPTBN1, KLF10, LYN, EPHA4, CDKN2B, XBP1, RBFOX2, IGFBP5, TIPARP, CRIM1, EIF2AK3, PTPN12, SULF1, ACVR2A, HCK, FGF18, ANOS1, BMPR1B and MET (Fig. 6 and Table 1 ). This study used the primary DEGs associated with Interleukin-3-mediated signaling pathway analysis. This task is performed to confirm that the IL3 pathway is enriched by the DEGs found in Cipa treated group. We implemented Interleukin-3-mediated signaling pathway genes to construct Fig. 7A showing subnetwork 1 and Fig. 7B showing subnetwork 2 protein-protein interaction networks.

Table 1

GO functional enrichment and KEGG pathways analyses of differentially expressed genes of MCF7 cell line of breast cancer.

Category	Term	Count	p-Value	Genes
GOTERM_BP_FAT	GO:0016339 ~ calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules	8	2.89E−07	PCDHB2, PCDHB6, PCDHB5, PCDHB4, PCDHB3, PCDHB11, PCDHB10, PCDHB9
GOTERM_BP_FAT	GO:0007416 ~ synapse assembly	10	1.88E−04	PCDHB2, FLRT3, TPBG, PCDHB6, PCDHB5, PCDHB4, PCDHB3, PCDHB11, PCDHB10, PCDHB9
GOTERM_BP_FAT	GO:0007156 ~ homophilic cell adhesion via plasma membrane adhesion molecules	10	8.24E−04	PCDHGB5, PCDHB2, PCDHB6, PCDHB5, PCDHB4, PCDHB3, PCDHB11, PCDHB10, PCDHB9, PCDHB8
GOTERM_BP_FAT	GO:0038156 ~ interleukin-3-mediated signaling pathway	29	8.33E−04	IL-1R1, IL-13RA1, FOXC1, TRIO, PPM1L, IRS1, XIAP, PTPRK, RASSF2, FNTA, FLRT3, HES1, SPTBN1, KLF10, LYN, EPHA4, CDKN2B, XBP1, RBFOX2, IGFBP5, TIPARP, CRIM1, EIF2AK3, PTPN12, SULF1, ACVR2A, HCK, FGF18, ANOS1, BMPR1B, MET
GOTERM_CC_FAT	GO:0005740 ~ mitochondrial envelope	23	0.002144	LYN, PTCD3, GSTK1, EPHA4, ANXA1, SIRT5, CHCHD10, CPOX, MRPL47, DUSP18, MAIP1, GDAP1, MRPS5, MRPL32, HK2, PARL, COQ4, CASP8, MRPL1, PMPCB, EFHD1, SLC25A6, ABCG2
GOTERM_CC_FAT	GO:0031966 ~ mitochondrial membrane	22	0.002365	LYN, PTCD3, GSTK1, EPHA4, ANXA1, SIRT5, CPOX, MRPL47, DUSP18, MAIP1, GDAP1, MRPS5, MRPL32, HK2, PARL, COQ4, CASP8, MRPL1, PMPCB, EFHD1, SLC25A6, ABCG2
GOTERM_CC_FAT	GO:0031967 ~ organelle envelope	31	0.003082	NUP205, TNKS, NUDT1, CPOX, MGST2, DTX2, RGPD8, DUSP18, MAIP1, MRPL32, HK2, PARL, GLE1, CASP8, MRPL1, PMPCB, EFHD1, LYN, PTCD3, GSTK1, EPHA4, ANXA1, CHCHD10, SIRT5, MRPL47, RANBP6, GDAP1, MRPS5, COQ4, ABCG2, SLC25A6
GOTERM_CC_FAT	GO:0031975 ~ envelope	31	0.003282	NUP205, TNKS, NUDT1, CPOX, MGST2, DTX2, RGPD8, DUSP18, MAIP1, MRPL32, HK2, PARL, GLE1, CASP8, MRPL1, PMPCB, EFHD1, LYN, PTCD3, GSTK1, EPHA4, ANXA1, CHCHD10, SIRT5, MRPL47, RANBP6, GDAP1, MRPS5, COQ4, ABCG2, SLC25A6
GOTERM_MF_FAT	GO:0005509 ~ calcium ion binding	24	0.00137609	PLA2G12A, ANXA1, PCDHGB5, ANXA3, SULF1, PCDHB11, PCDHB10, AIF1L, FKBP9P1, PCDHB2, PCDHB6, PLS3, EFHD1, PCDHB5, PCDHB4, ESYT2, MAN1A1, PCDHB3, PAM, SPTAN1, PCDHB9, PCDHB8, PLS1, MATN2
GOTERM_MF_FAT	GO:0004725 ~ protein tyrosine phosphatase activity	7	0.006801	PTPN18, MTMR1, EYA2, DUSP18, PTPRK, PTPN12, PTPN3
GOTERM_MF_FAT	GO:0004672 ~ protein kinase activity	20	0.009606	LYN, EPHA4, TRIO, CRIM1, EIF2AK3, STK39, LMTK2, PRKX, PHKA2, ACVR2A, CLK1, HCK, RASSF2, C8ORF44-SGK3, STK17B, FGF18, SGK3, TRIB1, BMPR1B, MET
GOTERM_MF_FAT	GO:0004857 ~ enzyme inhibitor activity	14	0.01064	CAST, CDKN2B, ANXA1, ANXA3, CRIM1, LMTK2, XIAP, TFPI, FLRT3, ANOS1, PDE6D, PHACTR2, FNIP1, TRIB1
KEGG_Pathway	hsa03013:RNA transport	8	0.01820	NUP205, SMN2, RGPD5, RGPD8, PABPC1L2B, PABPC1L, PABPC1L2A, SMN1
KEGG_Pathway	hsa05100:Bacterial invasion of epithelial cells	5	0.03378	ARHGAP10, CLTCL1, FN1, MET, CRKL
KEGG_Pathway	hsa00480:Glutathione metabolism	4	0.04551	GSTK1, GPX3, IDH1, MGST2
KEGG_Pathway	hsa03015:mRNA surveillance pathway	5	0.05434	FIP1L1, PABPC1L2B, PABPC1L, PABPC1L2A, WDR33

BP, biological processes; FAT, functional annotation tool; MF, molecular function; CC, cellular component.

Fig. 7

Hub genes associated with Interleukin-3-mediated signaling pathway. (A) Subnetwork 1 and (B) Subnetwork 2 protein-protein interaction network. Specifically, “green” and “purple” indicate the nodes are up- and down-regulated, respectively. The progressive color changing represents the expression levels. The areas of the nodes indicate the degrees that the nodes connect to others. The node sizes represent ranking of significant genes in terms of degree centrality i.e., the greater quantity of neighbors a node has.

Associated genes with the DEGs of the data set enriched with the Interleukin-3-mediated signaling pathway would be associated with severity of SARS-CoV-2 and other infectious diseases. In subnetwork 1 of Interleukin-3-mediated signaling pathway (Fig. 7A), LYN gene is in the center and the most significant gene in terms of betweenness centrality (Table 3). In a previously reported study (Weisberg et al., 2020), we focused our attention on the member of the Src protein family, such as the LYN gene, including kinase inhibitors involved in and targeting these proteins in SARS-CoV, MERS-CoV and associated virus infections.

As previous studies confirmed that IGFBP5 plays role in key positions such as inducing cell adhesion, increasing cell survival and inhibiting cell migration in MCF-7 human breast cancer cells (Sureshbabu et al., 2012) and reported as one of the significant genes of COVID-19 (Vastrad et al., 2020). These outcomes confirm the key task of the Interleukin-3-mediated signaling pathway engaged in breast cancer, COVID-19, and treatments of similar infectious diseases offering new molecular therapeutic targets to improve fundamental drug agents.

Effects of Cissampelos pareira (Cipa) treatment on hub genes

We searched to see if any statistically significant differences existed in DEGs considering the means of 6 categories of vehicle treated and treated groups. We employed a post-hoc test, i.e., Tukey-HSD (Honestly Significant Difference) after one-way ANOVA was performed. There were vehicle_1, 1μg_1, 10μg_1, 100μg_1, 500μg_1, and 1000μg_1 replicates of 3 were used to study the related genes. So there are 6 treatments - which gives 15 possible pairwise comparisons between them: tr1-vhc_tr, tr10-vhc_tr, tr100-vhc_tr, tr500-vhc_tr, tr1000-vhc_tr, tr10-tr1, tr100-tr1, tr500-tr1, tr1000-tr1, tr100-tr10, tr500-tr10, tr1000-tr10, tr500-tr100, tr1000-tr100, and, tr1000-tr500.

Using Tukey's range test the LYN, IGFBP5, IL-1R1, and IL-13RA1 (Fig. 8 ) has all these comparisons displayed at different heights with the label on the left y axis side. The extended lines show the 95% confidence intervals. Based on this, the LYN, IGFBP5, IL-1R1, and IL-13RA1 genes treated with 500 μg and 1000 μg amounts of Cipa was statistically significant relative to the rest of the pairwise comparisons of treatments. Furthermore, most of the DEGs associated with the Interleukin-3-mediated signaling pathway treated with 500 μg and 1000 μg amounts of Cipa were statistically significant from the rest of the pairwise comparisons of treatment. XIAP and PTN treated with 1000 μg amounts of Cipa was also statistically significant from the rest of the pairwise comparisons of treatments.

Fig. 8

LYN, IGFBP5, IL-1R1, and IL-13A1 genes Tukey-HSD test plots. Abbreviations: vhc_tr, vehicle treated; tr1, treatment 1 μg; tr10, treatment 10 μg;tr100, treatment 100 μg, tr500, treatment 500 μg; tr1000, treatment 1000 μg.

Discussion

Knowledge of the pathogenesis of COVID-19 is vital for developing therapies particularly important for risk groups such as cancer patients for unvaccinated persons or poor responders to vaccination. Here, we analyzed gene expression profile of MCF7 cell lines of breast cancer treated with Cipa to determine COVID-19 associated genes, proteins, biological, and molecular pathways. The results of down-regulated DEGs showed three metabolic KEGG pathways “hsa00480:Glutathione metabolism”, “hsa04146:Peroxisome” and “hsa04142:Lysosome”. A recent study, which has confirmed our findings, has shown the target human genes were down-regulated which account for a significant percentage of COVID-19 (Verma et al., 2020).

Cancer patients with SARS-CoV-2 have higher mortality rates among normal cancer cases (Curigliano, 2020; Wang and Zhang, 2020). Standard independent causes related with increased risk of mortality are those involving older age, history of smoking, active performance status. In a recent study in the U.K., 13% of the 800 patients with COVID-19 had breast cancer. Among the 800 patients, the breast cancer mortality rate for those with COVID-19 is 8%. Breast cancer ranked 3rd among other cancer types with those ill from COVID-19 in terms of mortality rates. In another recent study from the US, it was reported that 12% of 1482 cases with COVID-19 had a family history or history of comorbidity, of which 48% were obese (Garg, 2020). The roles of peroxisome metabolism in breast cancer patients were studied in the obesity framework previously (Carter and Church, 2009). In this framework, the link between breast cancer and COVID-19 cannot be underestimated. This current analysis demonstrated COVID-19 primarily has a negative impact on calcium-dependent cell-cell adhesion via plasma membrane cell adhesion molecules, synapses assembly, homophilic cell adhesion via plasma membrane adhesion molecules, and interleukin-3-mediated signaling. Similarly, in a recent study in China with 32 control and 39 COVID-19 patients, high expression of cell adhesion molecules is directly correlated with the burden of COVID-19 disease and may affect coagulation dysfunction (Tong et al., 2020).

Drugs targeting cell-cell adhesion molecules, RNA transport, bacterial invasion of epithelial cells, and mRNA surveillance pathways may be beneficial in case of decreasing immunity or emerge as antigenically distinct strains in the near future. Our results also confirmed that the importance of regulating metabolisms such as protein kinase activity is significant mediators of several viral infections, specifically SARS-CoV and MERS-CoV (Weisberg et al., 2020).

The consequent structure of the PPI network detected 10 DEGs; RPL7L1, NIFK, LYN, MRPS5, MRPL1, IRS1, MET, CASP8, HCK, and TNKS as promising biomarkers involved in COVID-19 particularly MCF7 cell lines of breast cancer. Ribosomal Protein L7 Like 1 (RPL7L1) plays a key role in “RNA-binding”. Nucleolar protein interacting with the FHA domain of MKI67 (NIFK) is involved in mitosis and cell cycle progression as well as “RNA-binding”. It was previously suggested that Src family tyrosine kinase (LYN) gene kinases might be essential for effective MERS-CoV replication (Shin et al., 2018). Thus, LYN can be a potential target gene to produce a significant decline in COVID-19 titers. MRPS5 and MRPL1 are ribosomal proteins that bond to proactively translated mRNA that if detected may serve as potential inhibitor parts or mechanisms to significantly reduced SARS-CoV-2 replication (Bojkova et al., 2020; Wang et al., 2020). Insulin receptor substrates (IRS1) which plays a role directing signal transduction whereas significant inflammatory pathways suppress signal propagation by serine kinase activity (Olefsky and Glass, 2010). For instance, metformin which is the medication for the treatment of type-2 diabetes indirectly reduces Akt activity through phosphorylation of IRS, making it an option of the usage as anti-COVID19 (Sharma et al., 2020). MET is defined as proto-oncogene (RTK) and its alterations are drivers of human cancer which can also serve as a potential biomarker for the therapy of COVID-19 (Rezaei et al., 2020). PPI network intersects with the following cytokine-related genes and receptors; IL-13RA1, IL-4R, IL13, IL-4, IL-6ST, CXCR4, TNFSF10, TNFRSF10A, TNFRSF1A, TNFRSF10B, CSF2RA, and CSF2RB. IL-3 is a glycoprotein and a growth factor for B lymphocytes which participates in the process of hematopoiesis, spread, and differentiation of hematopoietic stem cells. IL-3, moreover, is known to activate several tyrosine kinases such as LYN, FYN, SRC, SYK, TEC1 and HCK.

Furthermore, module analysis of the human PPI network of the Interleukin-3-mediated signaling pathway is involved in the process of hematopoiesis. In subnetwork 2, IGFBP5 centered PPI network of the Interleukin (IL)-3-mediated signaling detected genes are IGF1, IGF2, VTN, THBS1 PAPPA2, IGFALS, SERPINE1, and FHL2. Suppression of insulin-like growth factor (IGF1) is recommended to treat SARS-CoV-2 related adult respiratory distress syndrome (Winn, 2020). In the case of low expression levels of Vitronectin (VTN), a gene from the BTK (Bruton tyrosine kinase) family, it is a potential biomarker in the progression of COVID-19 disease whereas we would not be able to detect VTN's status of upregulation or down-regulation status (Treon et al., 2020). In a recent study, Insulin-like growth factor binding protein, acid labile subunit (IGFALS) was detected to be a possible novel biomarker of the COVID-19 prognosis (Kimura et al., 2020). In another study conducted on colorectal cancer (CRC) patients with COVID-19 identified Plasminogen Activator Inhibitor 1 (SERPINE1) was identified one of the pharmacological targets of niacin (Li et al., 2020). The limitations of the study are Cipa treated MCF7 cell lines of breast cancer in vitro tested against COVID-19 that impacted only some types of breast cancer. Further, the study is designed experimentally in Ayurgenomics framework to explore the molecular signatures of the whole formulation of Cipa and to probe whether there are molecular biomarkers that connect MCF7 cell of breast cancer with antiviral inhibition. This study is a pure analysis of gene expression profiles of from single data set of the molecular mechanism of Cipa as an antiviral potential. Validation should be performed by using data sets from different sources and by using more samples of MCF7 cell lines in-vivo or cells as well.

Conclusions

Here, we focused on publicly available GSE156445 data set of Gene Expression Omnibus (GEO) database, i.e., gene expression profile of MCF7 breast cancer cell lines exposed to Cipa treatment with control group to reveal association with COVID-19. Taken together, this current study sheds light into potential host metabolic and immune pathways vulnerable to breast cancer patients with viral infections such as SARS-CoV-2, SARS-CoV, and MERS-CoV. Further, our results suggest a list of potential key target genes such as LYN, IGFBP5, IL-1R1 and IL-13RA1 which may facilitate unvaccinated groups or poor responders to vaccination.

CRediT authorship contribution statement

Emine Güven: Conceptualization, Methodology, Software, Data curation, Writing – original draft, Visualization, Investigation, Supervision, Validation, Writing – review & editing.

Declaration of competing interest

The authors state that there's no financial/personal conflict of interest.