Transcriptome analysis and connectivity mapping of Cissampelos pareira L. provides molecular links of ESR1 modulation to viral inhibition

Bioactive fractions obtained from medicinal plants which have been used for the treatment of multiple diseases could exert their effects by targeting common pathways. Prior knowledge of their usage could allow us to identify novel molecular links. In this study, we explored the molecular basis of action of one such herbal formulation Cissampelos pareira L. (Cipa), used for the treatment of female hormone disorders and fever. Transcriptomic studies on MCF7 cell lines treated with Cipa extract carried out using Affymetrix arrays revealed a downregulation of signatures of estrogen response potentially modulated through estrogen receptor α (ERα). Molecular docking analysis identified 38 Cipa constituents that potentially bind (ΔG <  − 7.5) with ERα at the same site as estrogen. The expression signatures in the connectivity map (https://clue.io/;) revealed high positive scores with translation inhibitors such as emetine (score: 99.61) and knockdown signatures of genes linked to the antiviral response such as ribosomal protein RPL7 (score: 99.92), which is a reported ERα coactivator. Further, gene knockdown experiments revealed that Cipa exhibits antiviral activity in dengue infected MCF7 cells potentially modulated through estrogen receptor 1. This approach reveals a novel pathway involving the ESR1-RPL7 axis which could be a potential target in dengue viral infection.

www.nature.com/scientificreports/ shown to have antiviral and antimalarial activities, respectively 4,7 . Cipa methanolic extracts have also been shown to potentially modulate the estrus cycle in mice 8 . However, there is limited research on exploring whether there are shared pathways through which these extracts work in such diverse diseases. We used a genomics framework to explore the molecular signatures of the whole formulation of Cipa and to probe whether there are molecular players that connect the hormonal axis with antiviral inhibition. In this study, we carried out a transcriptomic analysis of Cipa in MCF7 cell lines and queried its effects on different pathways using gene ontology tools, Gene Set Enrichment Analysis, and probed its connectivity with drugs and other perturbagens including gene knockdowns using L1000 9,10 in CMAP (the connectivity map, the compendium of gene expression profiles 10 ). Pathways linked to lipid metabolism, viral transcription, translation, and estrogen axis seem to be modulated by Cipa. It also shows a high connectivity score with protein synthesis inhibitors that have demonstrated antiviral effect in in vitro studies. Besides the downregulation of ESR1, we also obtained a high connectivity score of Cipa with a knockdown signature of a selective estrogen receptor-alpha coactivator, RPL7 also involved in protein translation. Docking studies of Cipa constituents to estrogen receptor alpha show a high binding affinity of ~ 38 compounds to regions that bind known ERα modulators. Cipa inhibits dengue virus replication in MCF7 cell lines, and ESR1 knockdown in Cipa treated infected cells reduces this effect. Our study highlights a possible connection between ESR1 modulation and antiviral activity explored by the effect of the herb C. pareira.

Results
Cipa modulates estrogen receptor alpha in MCF7 cells. In order to probe the mechanism of action of C. pareira, we performed transcription profiling on Affymetrix microarray HTA 2.0. Genome wide expression assay of MCF7 cells provided 93 downregulated and 131 upregulated genes at 500 μg and 253 downregulated and 587 upregulated genes at 1000 μg, of Cipa treatment (Supplementary file 1). Gene ontology analysis revealed that the upregulated set had significant enrichment of transcription regulation by RNA polymerase II, tissue development, supramolecular fiber organization etc., whereas the downregulated genes were enriched in lipid metabolism processes (Fig. S3).
Gene Set Enrichment Analysis (GSEA) analysis using GSEA version 4.0.3, provided 59 positively and 12 negatively enriched gene sets (cut off set at p-value < 0.05 and FDR < 0.05). We observed a significant enrichment of estrogen response gene sets in both positive and negative enrichments (Fig. 1A). This indicates that Cipa may be modulating the expression of estrogen receptor alpha. Thus, we tested the mRNA expression of Estrogen receptor alpha (ESR1) by real-time quantitative analysis of mRNA levels in both above mentioned concentrations of Cipa. We observed that ESR1 levels were significantly downregulated by Cipa (Fig. 1B). Since estrogen receptor alpha (ERα), the protein encoded by ESR1 is an important transcription factor, we wanted ascertain whether the ERα binding sites, known as estrogen response elements (EREs) were differentially present in the promoter regions of the identified genes. Motif searches of estrogen response elements (ERE) of all the differentially expressed genes and 27 unchanged genes using PROMO tool reveal significant enrichments in ERE density in up and downregulated genes (Kruskal-Wallis, non-parametric test p value = 0.0071) in the 5 KB upstream region (Fig. 1C). There is no significant difference in the number of ESR1 elements between up regulated, downregulated and unchanged genes sets in the 5 KB upstream region. The regression analysis between fold change and the number of ERE sites revealed that the presence of ERE sites in the 5 KB upstream region have negative effects on expression. The differentially regulated genes also include well-characterized estrogen response genes with EREs in their promoters (Fig. S4).
Cipa shows high positive connectivity scores with protein synthesis inhibitors that are potential antivirals. A ranked gene signature of 69 down and 129 up genes was used to query the connectivity map (https:// clue. io/). At a cut-off score of 95 there were 69 positively connected small compounds. The topscoring compounds belonged to the class of protein translation inhibitors e.g., emetine (99.61), cycloheximide (99.86), verrucarin-a (99.47), cephaeline (99.51) ( Fig. 2A). Emetine, homoharringtonine, anisomycin, cephaeline, cycloheximide, and QL-XII-47 have reported antiviral activity (Table S2). In the genetic perturbations set, 77 genes had knockdown (KD) signatures score > 90, and all were involved in translation initiation, and viral transcription (Figs. 2B and S6). Since protein synthesis inhibition appears to be significant in our CMAP analysis, we decided to explore our expression data for genes involved in translation. The gene RPL7 (Large Ribosomal Protein 7) showed the highest knockdown connectivity score with Cipa. Although RPL7 was not identified among the differentially expressed genes, qPCR showed a significant downregulation of the gene mRNA levels (Fig. 2C). In vitro assays showed moderate inhibition of translation in response to Cipa (Fig. S7).

Virtual screening demonstrates Cipa constituents bind with estrogen receptor alpha.
To investigate the atomistic detail of possible Cipa constituent binding to the Estrogen receptor alpha, unbiased virtual screening of the constituent compounds against the receptor was employed. This revealed four preferential locations of binding for the 61 compounds. The largest cluster (pink) comprised 38 out of the 61 Cipa constituents, followed by the cluster of 8 (blue), 6 (orange), and 3 (green). As a positive control, we took previously reported SERMs, namely, tamoxifen and 17-β-estradiol. Our results show that both of these control molecules bind at the same location as the largest cluster in our findings with high concurrency to the reported binding energy (Fig. 3, Table 1). It is also noteworthy that more than 30 Cipa components had a minimum binding of < − 8.00 kJ/mol and at least 9 molecules have a minimum BE of < − 9 which is equal to or even higher than Tamoxifen binding energy − 9.2 (Table S1). Ligplot of the detailed binding sites for the controls and top binding molecules in Fig. S5 4A) with DENV-2 ( Fig. S9). Both cell lines were infected with DENV-2 at 10 MOI followed by treatment with Cipa. We observe that dengue inhibition was higher in MCF7 cells compared to MDA-MB-231 cells at both concentrations of Cipa (Fig. 4B). This differential inhibition could be a consequence of different cell types.
To further validate that the DENV inhibition by Cipa may be modulated by ESR1, we performed siRNA-mediated knockdown of ESR1 in MCF7 cells. Knockdown efficiency as confirmed by qRT-PCR 48 h post transfection show approximately 70-80% reduction in mRNA levels compared to that of non-targeting control (NTC) siRNA. After viral infection and subsequent Cipa treatment we observe that the inhibitory effect of Cipa was Graph is cumulative of three separate qPCR experiments each run with triplicates to ensure accuracy. Students t-test was performed for statistical significance ***p-value < 0.001, **p-value < 0.01. (C) Density of estrogen response elements 5 Kb upstream of differentially expressed genes identified using microarray at both concentrations of Cipa was calculated by dividing the number of EREs found with the length of the genomic region investigated (5000 bp each). The ERE elements density was found to be higher in down regulated genes than upregulated genes (Kruskal-Wallis test was applied to assess significance, p value < 0.05).

Discussion
Cissampelos Pareira or Patha, is a herb found in different parts of the world. Ayurveda describes its use in several formulations related to female hormonal disorders as well as in fevers. Recently the antiviral activity of CIPA also came to light. Despite this, the molecular mechanisms underlying such diverse activities of Cipa remain unexplored. We built our hypothesis around the Ayurvedic knowledge of its usage and the fact that there could    www.nature.com/scientificreports/ be possible crosstalk between the molecular mechanisms underlying stress response pathways implicated in these unrelated disorders. We carried out genome-wide expression profiling of MCF7 cell lines treated with Cipa to infer their plausible mechanisms of action as well as queried the signatures in connectivity map to infer similarities with the effect of other drugs, perturbagens, and gene-knockouts. Interestingly, the gene set enrichment analysis highlights an enrichment corresponding to the upregulation of genes upon inhibition of ESR1 and response to endocrine therapy in cancer (Fig. 1A). Upon qPCR analysis, we also find a downregulation of the ESR1 mRNA by Cipa. These results suggest that Cipa may play a role as a modulator of ESR1. When we searched for the estrogen response elements in the promoter region of the genes differentially expressed by Cipa, we found that the downregulated genes were highly enriched for such elements. Thus, highlighting that ESR1 downregulation is an important aspect of Cipa mechanism of action. Noteworthy, virtual screening of the constituent drug molecules against the estrogen receptor alpha revealed four distinct sites of binding. The possibility of multiple components binding at the same site shows inherent redundancy in the herbal extract. This could impart a synergistic activity in the whole extract of Cipa that might be lost in isolates. The analysis also revealed an abundance of the components in the formulation that have comparable binding energies to previously reported regulatory molecules. The inhibitory activity of Cipa and its constituents on estrogen and progesterone has been reported previously, however to our knowledge a molecular link has not been discovered 11 . For the first time we report that the antifertility and hormone regulatory activity of Cipa are perhaps modulated by regulating the levels of ESR1 and the binding activity of Cipa constituents with ERα.
CMAP analysis reveals that Cipa signature positively connects with protein synthesis inhibitors and the knockdown signature of genes are involved in translation and ribosomal protein processing that are needed for viral transcription and translation. The knockdown signature of the gene RPL7, a large ribosomal protein 7 and an ESR1 regulator 12,13 shows the highest positive connectivity with the Cipa gene signature. ESR1 is also known to regulate the immune response and estrogen regulators have been shown to possess antiviral activities [14][15][16][17] . Protein translation inhibition is a very well established first-line defense against viral infections 18,19 . Recently, many translation inhibitors with antiviral activity and antivirals as translation inhibitors have been reported 20,21 . All kinds of viruses have been seen to involve the endoplasmic reticulum and host translational machinery to form a membrane structure in which they reside and replicate [22][23][24] . Most interestingly, the biological pathways enriched for differential expression also include protein modification processes in upregulated and endoplasmic reticulum protein processing in downregulated genes. Cipa appears to inhibit translation without inducing cell death or inducing autophagy, hence could be potentially an effective antiviral agent.
Most striking, an HCV inhibitor Telaprevir has recently been observed to inhibit estrogen dependent proliferation of breast cancer cells by modulating levels and transcriptional activity of ERα 25 . Also, Emetine which is the highest scoring small compound in the CMAP and a potential viral inhibitor, has been shown to decrease the levels of ER alpha in MCF7 cells 26 .
Our results so far indicate an involvement of the estrogen axis in the antiviral activity of Cipa. Upon experimental validation we observed that Cipa induced reduction of viral titers was markedly lower in MDA-MB-231 cells that are naturally deficient in estrogen receptor alpha and also in MCF7 cells when ESR1 gene is knocked down. Thus, signifying the role of the estrogen axis in the antiviral action.
Mechanistically, ESR1 has recently been shown to modulate viral response against hepatitis D virus (HDV) via carbamoylphosphatesynthetase 2, aspartate transcarbamylase and dihydroorotase (CAD) enzyme modulation 27 . The pathway involves pyrimidine biosynthesis inhibition. Interestingly, among the pathways downregulated by Cipa, pyrimidine metabolism is significantly enriched (Fig. S3). Further, pyrimidine analogs have been observed to have inhibitory effects on various viruses 28,29 . It may be possible that Cipa inhibits dengue virus in MCF7 cells by modulating this pathway. This should however be validated with further mechanistic studies.
It is interesting to note that using a formulation with known and widespread usage, we were able to identify ESR1 as a factor involved in the antiviral activity of a therapeutic. This gives a hint of probable involvement of such connections in other viral infections as well, thus somewhat signifying a potential difference in the antiviral response between males and females.

Materials and methods
Cell culture and Cipa treatment. Cissampelos pareira L. whole plant extract was obtained commercially in lyophilized form, from a GMP certified manufacturer. Crude water extracts were reconstructed from powdered Cipa 30 . Fresh extracts were prepared right before each experiment. The chemical profiling of the extract was done using UPLC. The details are described in supplementary methods and the identified constituents are provided in supplementary Table S2. For subsequent experiments, guidelines set by CSIR, India were followed.
MCF7 (triple positive breast cancer) cell line was obtained from NCCS, Pune. It was maintained in DMEM high glucose media, supplemented with 10% FBS,1 M HEPES and 1X antibiotic antimycotics. The cells were kept at 37˚C and 5% CO 2 . MycoAlert (Cat no. LT07-318) was used to test the culture for mycoplasma.
MCF7 cells were treated with 1 µg, 10 µg, 100 µg, 500 µg and 1000 µg per ml of the aqueous extract from Cissampelos pareira L. The extract was prepared outside the cell culture hood and then filtered using 0.22 µm syringe filter before its use in cell culture. The MCF7 cells were seeded at a confluency of 60-70% 18-24 h. prior to treatment. The cells were treated with Cipa for 24 h for each experiment. A fresh extract was made for each set of experiments.
RNA isolation and whole transcriptome analysis. MCF-7 cells were seeded in 6 well plates at 60-65% confluency and then treated with vehicle, 1 µg, 10 µg, 100 μg, 500 µg and 1000 µg of Cipa for 24 h. The cells were then trypsinized and washed with PBS. Total RNA was isolated using TRIzol (Ambion, Cat no. 15596026) extraction method and its integrity was checked on 1% agarose gel which was followed by Nanodrop quantifica- www.nature.com/scientificreports/ tion (ND1000, Nanodrop Technologies, USA). Genome-wide transcriptome data was generated using 250 µg of the total RNA for each sample following the manufacturer's protocol. GeneChip (Affymetrix) Human Transcriptome Array 2.0 cartridges were used for this experiment. HTA 2.0 chip can capture 245,349 protein coding transcripts and 40,914 non protein coding transcripts in a 64-format. Briefly, the total RNA was prepared with poly A controls and first strand cDNA followed by second strand cDNA were synthesized. This was followed by in vitro transcription to form cRNA which was used as a template for the formation of single strand cDNA or ss-cDNA. Finally, the ss-cDNA was fragmented and labelled. At each step, purification and quantification of the samples were ensured. Approximately 200 µl of the sample, containing about 5.2 µg of the labelled ss-cDNA was loaded into the cartridges, which were then kept in the Affymetrix Hybridization Oven, set at 45 °C and RPM 60, overnight. The cartridges were registered on AGCC (Affymetrix GeneChip Command Console) and the fluidics of the experiment (washing and staining) was done using Affymetrix GeneChip Fluidics Station 450. The stained chips were then scanned using GeneChip Scanner 3000. The preliminary images were quality checked and .cell files were generated. The generated CEL files were background normalized using the RMA method. Batch effects were removed and differential gene expression analysis was done using the limma package in R. After background correction and RMA normalization, the probes that exceeded the p value < 0.05 were annotated and analyzed for differential expression. Since no probe could qualify the cut off at 1 µg, 10 µg and 100 µg, the differential gene expression was calculated using the dataset from 500 and 1000 µg treatments with the vehicle used as control. The genes with absolute log twofold change ≥ 1 and p value < 0.001 were taken as differentially expressed. The raw files and data have been submitted to GEO with the accession number GSE156445. cDNA synthesis and qPCR. cDNA was prepared from 1 μg DNase-treated RNA using Applied Biosys- Functional enrichment of the differentially expressed genes. For functional analysis we used g: profiler 31 (https:// biit. cs. ut. ee/ gprofi ler/ gost). After correcting for FDR < 5%, enrichments were analyzed for GO: Molecular Function, Cellular Compartment, Biological Processes, KEGG and Reactome.
In order to identify the specific sets of genes modulated by Cipa, we used Gene set enrichment analysis. We carried out a pre-ranked analysis for GSEA (Gene set enrichment analysis), in which the gene list was ordered from the highest positive gene expression to the lowest negative gene expression. The latest versions of gene sets databases were selected for query. Datasets with more than 500 and less than 15 genes were excluded from the query. The output was set to a minimum cut off of p value < 0.05 and corrected for FDR < 25%. An additional filter of FDR < 5% was applied for identification of significant enrichments.
Estrogen response elements (ERE) analysis. Sequence of 5 KB upstream region of genes were downloaded from UCSC genome browser for Human genome GRCh38 build. ERE sites were mapped to 5 KB upstream sequence using promo tool 32 with similarity cut off value of 8.7. High ERE density in 5 KB upstream leads to down regulation of genes as per regression analysis. ERE density's coefficient of slope of regressed line (β = − 0.0044) is negative and p-value is significant (0.018). Pair wise ERE density comparison of differentially expressed genes and unchanged genes was done using R software (version 3.2).
Virtual screening of constituent compounds for estrogen receptor binding affinity. Virtual screening of Cipa ligands was performed against the estrogen receptor crystal structure (PDBID: 3OS8) using Autodock Vina, a more accurate and faster version of Autodock 4. The Cipa ligands available in the PubChem database were downloaded from there, the 3D structures of the rest of the ligands were drawn using Marvin Sketch, a computational tool for drawing 3 and 2 dimensional chemical structures. The structures were randomized and minimized prior to docking. A blind docking study was performed for each ligand wherein the possible search-space was the complete receptor molecule. For each drug molecules, the docking parameters were as follows: center_x: 9.43, center_y: 22.811, center_z : 23.418, size_x: 60, size_y: 60, size_z: 60, exhaustiveness 5000, num_modes 50,000, energy_range: 20. Furthermore, for each ligand, 50 such runs were performed and the conformation with the minimum binding energy from among all the stable conformations of the 50 runs (a total of 1000 conformational possibilities) was selected for the cluster analysis.
The analysis was performed using ADT, UCSF Chimera and LigPlot + .
Connectivity map analysis for identifying similarities with known perturbations. The  www.nature.com/scientificreports/ the L1000 data as landmark or well inferred) was used to query the connectivity map using clue.io touchstone database.

Cipa treatment in DENV infection. For these experiments, institutional Biosafety clearance was obtained
at Translational Health and Science and Technology Institute, Faridabad, Haryana, India. MCF-7 were seeded at 100,000 cells per well in 24 well plate and were maintained for 24 h (37 degrees; 5% CO2). Virus challenge (at 10 MOI, see supplementary) was given for 1 h, here the media was supplemented with 2% FBS. After 1-h of viral adsorption, cells were washed with PBS and DMEM high glucose with 10% FBS with or without Cipa (50 μg and 100 μg) was added. DENV titers in the supernatant were determined by plaque assay at 24 h pi. DENV plaque assay was set up in BHK-21 cells (C-13), purchased from ATCC (Cat. no. CCL-10). 50,000 BHK-21 cells were plated per well in 24-well plate. Serial dilutions of virus were added in triplicates and allowed to adsorb for 1 h, followed by overlay with 0.5% carboxymethyl-cellulose (CMC; Sigma). After 72 h, cells were fixed in 3.7% formalin and plaques were visualized by staining with crystal violet.
p-value > 0.05 is considered non-significant (NS). Statistical significance of the differences was determined with a one/two-tailed Student's t test as per appropriate. Any additional steps are mentioned in the text.