Fractional Deletion of Compound Kushen Injection, A Natural Compound Mixture, Indicates Cytokine Signaling Pathways are Critical for its Perturbation of the Cell Cycle

We have used computational and experimental biology approaches to identify candidate mechanisms of action of a traditional Chinese medicine. Compound Kushen Injection (CKI), in a breast cancer cell line in which CKI causes apoptosis. Because CKI is a complex mixture of plant secondary metabolites, we used a high-performance liquid chromatography (HPLC) fractionation and reconstitution approach to define chemical fractions required for CKI to induce apoptosis in MDA-MB-231 cells. Our initial fractionation separated major from minor compounds, and showed that the major compounds accounted for most of the activity of CKI. By systematically perturbing the major compounds in CKI we found that removal of no single major compound could alter the effect of CKI on cell viability and apoptosis. However, simultaneous removal of two major compounds identified oxymatrine and oxysophocarpine as critical compounds with respect to CKI activity. We then used RNA sequencing and transcriptomics analysis to correlate compound removal with gene expression and phenotype data. We determined that many compounds in CKI are required for its effectiveness in triggering apoptosis but that significant modulation of its activity is conferred by a small number of compounds. In conclusion, CKI may be typical of many plant based extracts that contain many compounds in that no single compound is responsible for all of the bioactivity of the mixture and that many compounds interact in a complex fashion to influence a network containing many targets.


INTRODUCTION
Natural compounds are chemically diverse and have long served as resources for the identification of drugs (Harvey et al., 2015). However, the standard approach of fractionating natural product extracts to identify a single compound's biological activity can fail because the original activity of the mixture is not present in single compounds after fractionation. This failure to identify single compounds implies that some natural product mixtures derive their activity from the interaction of several bioactive compounds within the mixture. Characterising the mode of action of natural product mixtures has remained a difficult task as the combinatorial complexity of such mixtures makes it unfeasible to screen all combinations of the compounds in the mixture.
We introduce here a "subtractive fractionation approach" using high performance liquid chromatography (HPLC) that can pinpoint significant interacting compounds within a mixture when coupled with a suitable bioassay. We combined this approach with RNAseq (RNA sequencing) characterisation of our bioassay, correlating the removal of interacting compounds with concomitant alterations in gene expression. This allows us to identify specific combinations of compounds associated with specific pathways and regulatory interactions. In this report, we have applied this approach for the first time to a particular Traditional Chinese Medicine (TCM) formulation: Compound Kushen Injection (CKI), which is used to treat approximately 30,000 cancer patients/day in China in conjunction with Western chemotherapy.
Twenty-one chromatographic peaks have been identified from CKI with eight compounds being recognized as major components on the basis of their abundance (Ma et al., 2014).
The extract containing the most abundant compounds in CKI is derived from Kushen herb which has a long history in the treatment of cancer patients (Xu et al., 2005;Cheng et al., 2006).
The main component of CKI, macrozamin, is a derivative of baituling which has been a suggested therapeutic agent for the treatment of inflammatory disease (Jiang et al., 1997). Gao and colleagues showed that treatment with each of four of the main compounds of CKI (oxymatrine, matrine, sophoridine and N-methycytisine) at 4mg/mL significantly decreased cell viability (Gao et al., 2018). However, these concentrations are relatively high when compared to the contributing concentration of these four main compounds in CKI (Ma et al., 2014). The two main components of CKI, matrine and oxymatrine, may have significant anticancer activities in various types of solid tumors including cervical cancer, prostate cancer, synovial sarcoma, and hepatocellular carcinoma (Wu et al., 2015;Cai et al., 2016;Aung et al., 2017;Gao et al., 2018;Zhou et al., 2018). In contrast, toxicity of medicinal herbs containing matrine and oxymatrine as main components has been reported (Wang and Yang, 2003;Chan et al., 2014). Overall, understanding the effects of CKI based on the effects of single compounds present in CKI has been at best, partially successful.
Alternatively, by removing one, two or three compounds, we have been able to map the effects of these compounds and their interactions to effects on specific pathways based on altered gene expression profiles in a cell-based assay. This has illuminated the roles of several major compounds of CKI, which on their own have little or no activity in our bioassay. This approach can be used to dissect the roles and interactions of individual compounds from complex natural compound mixtures whose biological activity cannot be attributed to single purified compounds.

Subtractive fractionation overview
Well resolved chromatographic separation of CKI was used to collect all of the major components of CKI as individual fractions ( Figure 1A). We then reconstituted all of the separated fractions except for those we wished to subtract. We tested the reconstituted combination of compounds/peaks to see if removal of a single (N-1) or multiple compounds, (N-2 or N-3, where N represents the number of compounds in CKI), or removal of all major peaks (minor, MN) or depletion of all minor peaks (major, MJ) significantly altered the effect of CKI in our cell based assays. Our cell based assays (Qu et al., 2016) measured MDA-MB-231(human breast adenocarcinoma) cell viability, cell cycle phase and cell apoptosis. A summary of the subtractive fractions used in the cell based assays is shown in Table 1. We then carried out RNA isolation of cells treated with CKI, individual compounds or CKI deletions for RNAseq. Differentially expressed (DE) genes in these samples allowed the association of specific compounds with cell phenotype and underlying alterations in gene regulation. By comparing DE genes across treatment combinations we identified specific candidate pathways that were altered by removal of single or multiple compounds, as detailed below.

HPLC fractions and content identification using LC-MS/MS
HPLC fractionation and reconstitution was used to generate a number of N-1, N-2, N-3, MJ and MN mixtures, ( Figure 1A, B, C and Supplementary Figure 1) with specific combinations and their components shown in Table 1. The concentrations of known compounds in CKI and reconstituted subtractive fractions were determined from standard curves (Supplementary Data1) for nine available reference compounds, using cytisine as an internal standard ( Table 2). The combined concentration of 9 reference compounds from CKI was approximately 10. 461 mg/mL, whereas subtractive fractions N-OmtOspc and N-MacOmtOspc had concentrations of reference compounds of 3.045 mg/mL and 2.335 mg/mL which were equivalent to the concentrations of these compounds in unfractionated CKI. The depleted OmtOspc and MacOmtOspc were not observed in the N-OmtOspc and N-MacOmtOspc respectively. These collectively suggested any effects observed after the treatments of N-OmtOspc and N-MacOmtOspc were not influenced by the concentrations. A total of 9 (N-1), 4 (N-3) and 9 (N-2) combinations, along with MJ and MN deletions were tested in our cell based assays (Table 1). The MJ subtractive fraction contained a total of 9 compounds including eight previously identified MJ peaks (Ma et al., 2014) and adenine (unpublished data from Ma Yue) ( Figure   1A) and the MN fraction contained the remaining peaks ( Figure 1C). MJ had no effect on cell viability, while MN reduced cell viability to the same extent as CKI (Figure2A).

Phenotypic changes associated with compound deletion
The 9 major compounds were individually depleted from CKI and tested as 9 (N-1) subtractive fractions, with no significant alterations in cell viability compared to CKI Sets of 3 compounds from the 9 major/standard compounds of CKI were depleted to generate 3 (N-3) subtractive fractions. The nine reference compounds were allocated into three groups, one of which contained structurally similar compounds (Omt, Ospc, Spc) and two other groups ([Mac, Ade, Tri] and [Nme, Mt, Spr]) that contained structurally different compounds. Of these three fractions, N-OmtOspcSpc decreased the number of viable cells significantly more than CKI while none of the sets of three compounds on their own had any effect on cell viability (Supplementary Figure 2). We then generated 9 (N-2) subtractive fractions based on the N-3 subtractive fractions (Table 1). Out of 9 (N-2) subtractive fractions (Supplementary Figure 2), only N-OmtOspc significantly decreased proliferation compared to CKI (P<0.05) ( Figure 2C).
We then depleted macrozamin, the only major compound derived from Baituling, from N-OmtOspc as N-3 (N-MacOmtOspc) in order to determine if there was an additional effect when compared to CKI. N-OmtOspc and N-MacOmtOspc both decreased cell proliferation to the same extent ( Figure 2C, D).

Cell cycle
While no change in cell viability was found across all N-1 treatments, cell cycle analysis was performed to identify more subtle differences. There was no statistically significant difference in phases of the cell cycle of MDA-MB-231 cells for any of the N-1 treatments compared to CKI ( Figure 3A). However, N-OmtOspc treatment significantly altered the cell cycle for MDA-MB-231 cells and induced significant higher apoptosis from 0.25mg/mL through 2mg/mL treatments as compared to CKI at both time points ( Figure 3B  These results indicated that the N-OmtOspc and N-MacOmtOspc subtractive fractions induced apoptosis not only in cancerous cells but also in non-cancerous cell lines. In contrast to this, no significant apoptosis was triggered by CKI on HFF cells. A small but significant apoptotic induction was observed for HEK-293.

Cytotoxicity
Because of the significant decreased viability accompanied by increased apoptosis triggered by subtractive fractions, cytotoxicity tests were carried out for HEK-293 and HFF cells using N-OmtOspc and N-MacOmtOspc subtractive fractions at concentrations equivalent to CKI 2mg/mL . N-OmtOspc and N-MacOmtOspc at equivalent concentration to CKI 2mg/mL were significantly cytotoxic to both non-cancerous cell lines ( Figure 3G).
Overall, these results indicate that removal of combinations of specific compounds from CKI had unpredictable effects on the ability of CKI to kill cells. While removal of all major compounds from CKI caused no loss of activity and removal of all minor compounds caused total loss of activity, removal of selected major compounds (N-OmtOspc) paradoxically caused major, significant increases in the ability of CKI to reduce viability and kill cells.

Differential gene expression
In order to understand the interactions of the components in CKI as a result of depletion, we The number of differentially expressed (DE) genes associated with each treatment was calculated using pair-wise comparative analysis. CKI treatment was used as a baseline to compare all other treatments in order to emphasize the effect of depleted compounds and CKI treatment was compared to UT.
There were thousands of upregulated and downregulated genes at 24 and 48 hours in most pairwise comparisons ( Figure 4B). However DE genes between OmtOspc and MacOmtOspc treatments were not observed and there were almost no DE genes between N-Mac, N-Nme and N-Tri treatments ( Figure 4B) indicating that these three subtractive fractions had very similar effects on gene expression.
When we compared the DE genes found between treatments, there were a large number of DE genes (~71.3%) shared between all four (N-1) treatments (Supplementary Figure 8,9 and Supplementary Table3). A similar number of shared DE genes (~24.6%) between four (N-1), OmtOspc and MacOmtOspc and between four (N-1), N-OmtOspc and N-MacOmtOspc as compared to CKI at 48 hours indicated that gene expression patterns from N-1 treatments were mostly different from N-OmtOspc, N-MacOmtOspc, OmtOspc and MacOmtOspc treated cells. 55% of the DE genes between UT, OmtOspc and MacOmtOspc were shared. When CKI treatment was compared to the four (N-1) treatments 42.8% of DE genes were shared, and when CKI was compared to N-OmtOspc and N-MacOmtOspc treatments, 50.1% DE genes were shared, indicating that N-OmtOspc and N-MacOmtOspc treatments appeared to be more similar to CKI than N-1 treatments.
The overall levels of similarity in DE genes were as follows: 1)All N-1 treatments had approximately 70 % similar gene expression patterns, 2) OmtOspc and MacOmtOspc treatments were approximately 50% similar to UT and 33% similar to N-1 treatments, 3) downstream gene expression patterns between N-1, N-OmtOspc and N-MacOmtOspc were approximately 37% similar.

Gene ontology and pathway annotation of DE genes
DE genes were analysed for over-representation in our data sets with respect to biological function using Gene Ontology (GO) annotation. We looked for shared DE genes between treatments and identify over-represented genes in these shared genes. The only common function enriched across all comparisons was for "cell cycle checkpoint". This confirmed earlier results (Qu et al., 2016) and was consistent with the phenotype data for CKI.

Subtracted fractionation altered pathways
We also used pathway based analysis to look for pathway level perturbation by comparing DE genes within Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between treatments. We used Signaling Pathway Impact Analysis (SPIA) to identify pathways with statistically significant perturbation values expected to alter pathway flux. We identified 86 pathways (SuppFig 11) with statistically significant (P< 0.05) perturbations of gene expression and of these, 15 pathways were most obviously linked to our phenotypes of cell viability, cell cycle and apoptosis ( Figure 5). By comparing the pathway gene expression perturbation scores (pG) between treatments we can see that three specific observations can be made: 1) N-1 fractional deletions vs CKI had significant effects on flux in some pathways without Overall our results support the concept of multi-compound/multi-target interactions for plant extract based drugs that contain many plant secondary metabolites. Biological effects of complex plant extracts may result from interactions of multiple compounds, with negligible effects from single compounds alone. This has implications for how we assess the functional evidence for such extracts.

DISCUSSION
Previous studies have demonstrated that CKI can alter the cell cycle and induce apoptosis in cancer cell lines (Xu et al., 2011;XU et al., 2013;Qu et al., 2016). CKI arrest of the cell cycle at checkpoints is correlated with the induction of double strand breaks by CKI treatment (Cui et al., 2018). In contrast to our experiments reported above, oxymatrine was previously shown to arrest the cell cycle and induce apoptosis in primary human myeloma cells through caspasedependent pathway (Wu et al., 2014) and inhibit the proliferation of laryngeal squamous cell carcinoma Hep-2 cells (Ying et al., 2015). As shown in this report, oxymatrine or oxysophocarpine or combined OmtOspc treatment caused no significant change in cell viability, the cell cycle or apoptosis, in agreement with prior work showing that oxymatrine and oxysophocarpine showed no significant effect on apoptosis, cell cycle or cell proliferation in HCT116 human colon cancer cells . However, at the level of gene expression, gene ontology analysis indicated that "cell cycle checkpoint" was significantly enriched in genes expressed in cells treated with all fractionated mixtures or mixtures of Omt and Ospc. Our results demonstrate that these compounds had little or no phenotypic effect on their own, but that when both were deleted, the remaining compounds unexpectedly had significantly greater effects on phenotype and gene expression.
When examined in the context of specific pathways, treatment with OmtOspc or N-OmtOspc which had strikingly different effects on phenotype, had similar effects on the perturbation of the "cytokine-cytokine receptor" pathway, the most commonly perturbed pathway seen in our analysis, that interestingly did not show up when comparing CKI to UT. This is consistent with previous information showing that CKI induced cytokines IL4 and IL10 (Tu et al., 2016). The observation that many genes in the "Cytokine-Cytokine Receptor Interaction" pathway were not affected by OmtOspc and MacOmtOspc compared to deletion fractions confirmed that removal of compounds rather than treatment with single or a few compounds can be more informative of the role and significance of individual compounds as part of mixtures/extracts.
The paradoxical result that removal of OmtOspc caused a striking increase in apoptosis is most simply explained by a model that assumes an integreation of effects of multiple compounds on many targets. Our approach allowed the identification of both synergistic and antagonistic interactions within the drug mixture. Viewed as a network where the compounds and the targets are nodes and the interactions between compounds and targets, and between targets are edges, it is clear that the edges (interactions) determine the overall effect of the compound mixture.
By removing one or two compounds from a mixture, we can potentially perturb the target network(s) to either reduce the effect of the mixture for some outcome or potentiate it for another. We believe this approach may be of general use for the study of herbal medicines/extracts, avoiding failures that stem from exclusive reliance on the identification of a single compound that accounts for most of the biological activity in mixtures. HPLC separation was achieved using a reparative C18 column (5 µm, 250 x 10 mm, Phenomenex) with the following mobile phases: 0.01 M ammonium acetate (adjusted to pH 8.0, solvent A) and acetonitrile + 0.09% trifluroacetic acid (solvent B).. The flow rate was 2 mL/min and linear gradient was adopted as follows; 0 min, 100% A; 60 min, 65% A, 70 min, 100%A. The chromatogram was recorded from 200 nm to 280 nm, with monitoring at 215 nm.

MDA
Samples were frozen and lyophilised using a Christ Alpha 1-2 LD lyophilizer (Martin Christ Gefriertrocknungsanlagen GmbH, Germany). Several cycles of lyophilisation and resuspension were used to remove all remaining HPLC solvents and final reconstitution carried out using MilliQ water and buffered with 10mM Hepes pH 6.8-7.0 (Gibco, Life technologies, USA). Lyophilised samples were resuspended to create an equivalent dilution for compounds in the sample compared to CKI.

Identification of reconstituted mixtures by LC-MS/MS
Agilent 6230 TOF mass spectrometer was used to determine the concentration of the known compounds from the CKI and reconstituted N-OmtOspc and N-MacOmtOspc mixtures. 10uL sample was injected with the flow rate of 0.8 mL/min, a gradient program of 0 min, 100% A; 40% B; 25 min, 60% B, 35 min, and solvents H2O + 0.1% formic acid (solvent A) and acetonitrile + 0.1% formic acid (solvent B). The column used was C18 (5 , 150 x 4.6 mm, Diamosnsil, Dkimatech). The recovered contents of the samples was measured by spike-in compound cytosine. Gas phase ions were generated with an electrospray source, with with key instrument parameters: gas temperature, 325; sheath gas temperature, 350; vCap, 3500; fragmentor, 175; acquisition range (m/z) 60-17000. Calibration curves for 9 standard compounds containing various concentrations were shown in Supplementary Data1.

Cell viability Assay
XTT [2,3-bis-(2-methoxy-4-nitro-5-sulfophenyl)-2H-tetrazolium-5-carboxanilide] and PMS (N-methyl dibenzopyrazine methyl sulfate) (50:1, Sigma-Aldrich) assay was used to assess cell viability following. Briefly, 8,000 cells in 50ul of medium were plated in 96 wells trays overnight prior to drug treatments in triplicate. Cells were subsequently treated with 50ul of drug mixtures to provide final concentrations of 0.25, 0.5, 2 and 2mg/mL in 100uL. Cell viability was then measured at 24 and 48 hours after drug treatment by the addition of 50 ul of XTT:PMS (50:1 ratio). An equal volume of medium and treating agents plus XTT:PMS was used to subtract the background optical density (OD). The absorbance of each well was recorded using a microplate reader at 492 nm.

Annexin V/PI apoptosis assay
Apoptosis, or programmed cell death, resulting from treatment was determined using an Annexin V-FITC apoptosis detection kit (eBioscience™ Annexin V-FITC Apoptosis Detection Kit, Thermofisher Scientific) according to the manufacturer's protocol. Briefly, 4x10 5 cells were seeded in 6-well plates in triplicate overnight prior to treatment. On the following day, cells were treated with the agents as described for 24 and 48 hours. Data was acquired with a BD LSR-FORTESSA (NJ, USA) flow cytometer, and FlowJo software (TreeStar Inc., OR, USA) was used to analyse the acquired data and produce percent apoptosis values.

Cell cycle assay
Cell culture and drug treatments were performed as described above for cell cycle analysis. A Propidium Iodide (PI) staining protocol (Riccardi and Nicoletti, 2006) was used to detect the changes in cell cycle as a result of treatment after 24-and 48-hours. The characteristics of stained cells were measured using a BD LSR-FORTESSA (NJ, USA) flow cytometer, and acquired data were analysed using FlowJo software.

Cytotoxicity assay
MDA-MB-231, HEK-239 and HFF cells were seeded in 96-well plates at a density of 2.5x10 3 cells per well in triplicate. CKI and fractionated mixtures to produce a final concentration of 1mg/mL and 2mg/mL were added to each well and after 24 hours of incubation, viable cells were measured using the AlamarBlue assay (Molecular Probes, Eugene, OR). Mercuric chloride (Sigma-Aldrich , St. Louis, MO) (5µM) was used as a positive control and wells without cells were set as a negative control in the same plate.

Sample preparation and RNA sequencing
Cells were plated in 6 well plates with a density of 2x10 5 cells/mL overnight prior to drug treatments. On the following day, CKI (to give a final concentration of 2mg/mL) and fractionated mixtures (equivalent dilutions of CKI) were added. Total RNA was isolated by using an RNA extraction kit (Thermo Fisher Scientific) according to the manufacturer's instructions and RNA samples were quantified and quality determined using a Bioanalyzer at the Cancer Genome Facility of the Australian Cancer Research Foundatin (SA, Australia).
RNA samples with RNA integrity number (RINs) > 7.0 were sent to be sequenced at Novogene (China). Briefly, after QC procedures were performed, mRNA was isolated using oligo(dT) beads and randomly fragmented by adding fragmentation buffer, followed by cDNA synthesis primed with random hexamers. Next, a custom second-strand synthesis buffer (Illumina) , dNTPs, RNase H and DNA polymerase I were added for second-strand synthesis After end repair, barcode ligation and sequencing adaptor ligation, the double-stranded cDNA library was size selected and PCR amplified. Sequencing was carried out on an Illumina HiSeq X platform with paired-end 150 bp reads.
Removal of unwanted variance (RUVs) package in R was applied to two different batches of transcriptome datasets to eliminate batch variance (Risso et al., 2014). APE was used to cluster the treatments (Paradis et al., 2004) followed by RUVs. Gene Ontology (GO) overrepresentation analyses were performed using clusterProfiler with the parameter ont = "BP"(Biological Process), pAdjustMethod = "BH", pvalueCutoff = 0.01, and qvalueCutoff = 0.05 (Yu et al., 2012). Signalling Pathway Impact Analysis (SPIA) was carried out to identify the commonly perturbed pathways within the treatments using the SPIA R package (Tarca et al., 2008). Significantly perturbed pathways were visualized using Pathview package in R (Luo and Brouwer, 2013).

Statistical analysis
Statistical analyses were carried out using GraphPad Prism 8.0 (GraphPad Software Inc., CA, USA). Student's t-test or ANOVA (one-way or two-way) was used when there were two or three groups to compare respectively. Post hoc "Bonferroni's multiple comparisons test" was performed when ANOVA results were significant. Statistically significant results were represented as p<0.05 (*) or p<0.01 (**) p<0.001 (***), or p<0.0001 (****); ns (not significant). All data were shown as mean ± standard deviation (SD).  Effect of N-MacOmtOspc subtractive fraction and MacOmtOspc, compared to CKI.
Statistically significant results relative to CKI treatment shown as p<0.05 (*) or ns (not significant), all data were shown as mean ± standard deviation (SD).

Conflict of interest
The authors declare no conflicts of interest.

Declaration of transparency and scientific rigour
This declaration acknowledges that this paper adheres to the principles for transparent reporting and scientific rigour of preclinical research recommended by funding agencies, publishers and other organisations engaged with supporting research.