High-throughput transcriptional profiling of perturbations by Panax ginseng saponins and Panax notoginseng saponins using TCM-seq

Junyun Cheng , Jie Chen , Jie Lio , Tinho Wng , Xin Sho , Jinbo Long ,Penghui Yng , Anyo Li , Zheng Wng , Xioyn Lu ,d,e,*, Xiohui Fn ,d,e,**

a Pharmaceutical Informatics Institute, College of Pharmaceutical Sciences, Zhejiang University, Hangzhou, 310058, China

b Future Health Laboratory, Innovation Center of Yangtze River Delta, Zhejiang University, Jiaxing, Zhejiang, 314100, China

c State Key Laboratory of Medical Molecular Biology, Department of Biochemistry and Molecular Biology, Institute of Basic Medical Sciences, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing,100005, China

d Innovation Center in Zhejiang University, State Key Laboratory of Component-Based Chinese Medicine, Hangzhou, 310058, China

eJinhua Institute of Zhejiang University, Jinhua, Zhejiang, 321016, China

Keywords:

Panax ginseng

Panax notoginseng

Sample multiplexing

RNA sequencing

ABSTRACT Panax ginseng (PG) and Panax notoginseng (PN) are highly valuable Chinese medicines (CM).Although both CMs have similar active constituents, their clinical applications are clearly different.Over the past decade,RNA sequencing(RNA-seq)analysis has been employed to investigate the molecular mechanisms of extracts or monomers.However, owing to the limited number of samples in standard RNA-seq, few studies have systematically compared the effects of PG and PN spanning multiple conditions at the transcriptomic level.Here, we developed an approach that simultaneously profiles transcriptome changes for multiplexed samples using RNA-seq (TCM-seq), a high-throughput, low-cost workflow to molecularly evaluate CM perturbations.A species-mixing experiment was conducted to illustrate the accuracy of sample multiplexing in TCM-seq.Transcriptomes from repeated samples were used to verify the robustness of TCM-seq.We then focused on the primary active components, Panax notoginseng saponins (PNS) and Panax ginseng saponins (PGS) extracted from PN and PG, respectively.We also characterized the transcriptome changes of 10 cell lines, treated with four different doses of PNS and PGS,using TCM-seq to compare the differences in their perturbing effects on genes,functional pathways,gene modules, and molecular networks.The results of transcriptional data analysis showed that the transcriptional patterns of various cell lines were significantly distinct.PGS exhibited a stronger regulatory effect on genes involved in cardiovascular disease, whereas PNS resulted in a greater coagulation effect on vascular endothelial cells.This study proposes a paradigm to comprehensively explore the differences in mechanisms of action between CMs based on transcriptome readouts.

1.Introduction

Panax ginseng (PG) and Panax notoginseng (PN) are valuable Chinese medicines(CM)that have been widely used for centuries in clinical therapies.The global PG market size is estimated at USD 8 billion in 2022 and is expected to reach USD 11.7 billion by 2026, and the Chinese market for PN reached USD 1.9 billion in 2019 and is still increasing annually.Pharmacological studies have demonstrated that both PG and PN have pharmacological properties of anti-Alzheimer's [1,2], cardiovascular protection[3,4], booster immunization [5,6], and hypoglycemic agents [7,8].As systematically summarized in our previous review [9],phytochemical studies have identified that PG and PN have similar chemical compositions, including saponins, polysaccharides,amino acids,and volatile oils,whereas their primary clinical indications vary widely [10-12].Therefore, systematically investigating the underlying mechanisms of action of the active components in PG and PN is essential to harness their therapeutic potential.

Previous methods for studying the active ingredients in PG and PN, such as chemical activity identification [13], in vitro screening[14], and in vivo animal modeling [15], have been limited to singular phenotype readout and have impeded the uncovering of PG and PN mechanisms of action.With the rapid development of sequencing technologies, RNA sequencing (RNA-seq) has been predominantly employed as a routine analysis for investigating the pharmacological effects of CM containing complex components[16-18].RNA-seq analysis has demonstrated that 20(S)-ginsenoside Rh2,a well-known protopanaxadiol-type ginsenoside from PG,has inhibited HepG2 cell proliferation via the p53 signaling pathway and DNA replication [19].Transcriptome analysis has revealed that ginsenoside-Rg1, one of the principal ingredients of PG, has regulated several important hub genes in rats with nonalcoholic fatty liver disease [20].Using RNA-seq, ginsenoside Rf has also been shown to reduce human tau proteotoxicity in a Caenorhabditis elegans tauopathy model[21].According to RNA-seq data,notoginsenoside R1 also enhances wound healing in diabetic rats by enhancing extracellular matrix remodeling and inflammation[22].However, these studies have focused on the mechanisms of action of some monomer saponins while disregarding the overall regulation by other main components of PG and PN.Furthermore,the workflow of standard RNA-seq is labor-intensive and costprohibitive, which results in low-throughput samples and impedes systematic investigation and comparison of PG and PN under multiple conditions.In the last decade, many sample multiplexing strategies have been developed to increase the number of RNA-seq samples for cost-efficient library generation [23].The core of multiplexed RNA-seq is “early barcoding” with sample-specific DNA tags in the primers used during complementary DNA (cDNA)generation, which is promising in high-throughput studies[24-27];however,no approaches have been developed specifically for simultaneous transcriptome profiling of a large number of CM samples.

To overcome these challenges,in this study,we have developed a multiplexed CMs profiling technology, namely profiling transcriptome changes for multiplexed samples using RNA-seq(TCM-seq),based on a DNA barcoding strategy.This can specifically label cDNAs of up to 468 samples, treated with CMs, for pooled library construction and sequencing, hence, improving the throughput of parallel sequenced samples and reducing the consumption of reagents.We focused on investigating the differences and similarities in the efficacy of the main active ingredients in PN and PG, namely Panax notoginseng saponins (PNS) and Panax ginseng saponins (PGS).With the throughput advantage of TCMseq, we compared the transcriptional responses of 10 cell lines treated with four different doses of PGS and PNS,and compared the differences in the efficacy of both saponins in terms of differentially expressed genes, biological functions, gene co-expression, and molecular network regulation to shed light on the efficacy and potential clinical applications of PG and PN.

2.Materials and methods

2.1.Materials and chemicals

PG was collected from Jilin Province in China and purchased from Zhejiang Yingte Pharmaceutical Co., Ltd.(Hangzhou, China).PN was collected from Yunnan Province in China and purchased from KPC Pharmaceuticals,Inc.(Kunming,China).Both root plants were identified by Xingchu Gong, an associate professor at the College of Pharmaceutical Sciences at Zhejiang University, according to the Chinese Pharmacopoeia.Reference ginsenosides Rg1,Re,Rf, Rb1, Rc, Ro, Rb2, Rb3, Rd, and notoginsenoside R1 were purchased from Shanghai Yuanye Bio-Technology Co., Ltd.(Shanghai,China).Analytical grade ethanol was purchased from Shanghai Lingfeng Chemical Reagent Co., Ltd.(Shanghai, China).Chromatographic-grade methanol and acetonitrile were purchased from Merck(Darmstadt,Germany)and that of phosphoric acid was purchased from the Tedia Company (Fairfield, IA, USA).Deionized water was obtained using a Milli-Q academic water purification system(Millipore, Milford, MA, USA).

2.2.Cell culture and treatment

Rat cardiomyocyte(H9c2),human malignant melanoma(A375),human lung adenocarcinoma (A549), human hepatocellular carcinoma (HepG2), human breast adenocarcinoma (MCF7), human glioma (U251), and human cardiomyocyte (AC16) cell lines were purchased from the Cell Bank of Chinese Academy of Sciences(Shanghai, China), and mouse hippocampal neuronal (HT22), human umbilical vein endothelial cell (HUVEC), human prostate cancer (PC3), and human colon cancer (HT29) cell lines were purchased from American Type Culture Collection (Rockville, MD,USA).Human breast adenocarcinoma (MCF7), HT29, and PC3 cell lines were cultured in minimum essential medium (12561072;Gibco,Grand Island,NY,USA),McCoy's 5A(16600082;Gibco),and F12K (21127022; Gibco) medium, respectively, containing 100 U/mL penicillin and 100 μg/mL streptomycin(15140122;Gibco) with 10% heat-inactivated fetal bovine serum (FBS, 10099141; Gibco).The other eight cell lines were cultured in dulbecco's modified eagle medium (11995073; Gibco) containing 100 U/mL penicillin-100 μg/mL streptomycin and 10% heat-inactivated FBS, and incubated at 37°C with 5% CO2.

2.3.Multiplexed RNA-seq by TCM-seq

Different cell lines were seeded into 384-well plates at a cell density of 5 × 103/well.After incubation for 24 h, cell lysis in situ and mRNA purification were performed using TurboCapture 384 mRNA capture plate according to the manufacturer's protocol for TurboCapture 384 mRNA Kit (72271; QIAGEN, Dusseldorf, Germany).Then,15.5 μL of TCE buffer of the capture plate in each well was transferred from well-to-well to pre-assigned polymerase chain reaction (PCR) plate wells containing 2.5 μL of reverse transcription (RT) mix,1 μL of barcoded primer (Table S1), and 1 μL of 1:1000 ERCC RNA Spike-In Mix1 (00718793; Invitrogen, Carlsbad,CA, USA) per well for RT reaction.RT mix contains 200 U Proto-Script®II Reverse Transcriptase (2690A; Takara, Kyoto, Japan),1× RT buffer (00788967; Thermo Fisher Scientific Inc., Waltham,MA, USA), 7.5 U RNaseOUT (10777019; Thermo Fisher Scientific Inc.), 1 M betaine (B0300; Sigma-Aldrich, Saint Louis, MO, USA),8 mM MgCl (AM9530G; Invitrogen, Dusseldorf, Germany), 0.8 μM template switch oligo(TSO) (Table S1),and 0.08 μM dNTP(R0192;Thermo Fisher Scientific Inc.).After incubation at 42°C for 90 min,the RT reaction products from each well were pooled and purified using a DNA Clean&Concentrator-100 Kit(D4029;Zymo Research,Orange, CA, USA) to obtain 150 μL of clean eluate.The eluate(150 μL)was concentrated,using a DNA Clean&Concentrator-5 Kit,into 32 μL of eluate for the nucleic acid endonucleation reaction to remove excess primers by adding 4 μL of exonuclease I(Exo I)and 4 μL of Exo I buffer (EN0581; Thermo Fisher Scientific Inc.) and incubating at 37°C for 30 min,followed by incubation at 85°C for 15 min to inactivate the nuclease enzyme, and kept at 4°C.

Next, 50 μL of 2× Kapa HIFI PCR ReadyMix (KK2602; Kapa Biosystems, Cape Town, South Africa) and 10 μL of 10 μM preamp_PCR primer(Table S1)were added to the above 40 μL of eluate mix to perform the first PCR thermocycling as follows: 98°C for 3 min,then six cycles at 98°C for 20 s,65°C for 45 s,72°C for 6 min,followed by a single 10 min at 72°C,and held at 4°C.The first PCR product was purified using VAHTS DNA clean beads (N411-01;Vazyme, Nanjing, China) with a 1:0.8 (PCR product:beads, V/V)bead ratio according to the manufacturer's protocol and then eluted in 40 μL of nuclease-free water.Then, 50 μL of 2× Kapa HIFI PCR ReadyMix and 10 μL of 10 μM preamp_PCR primer were added to the above 40 μL of eluate mix to perform the second PCR thermocycling as follows:98°C for 3 min,then ten cycles at 98°C for 20 s,67°C for 20 s,72°C for 6 min,followed by a single cycle at 72°C for 10 min and holding at 4°C.The second PCR product was purified using VAHTS DNA clean beads at a 1:0.7 (PCR product:beads, V/V)bead ratio, and then eluted in 15 μL of nuclease-free water.cDNA products (2 μL) were run on a Tapestation 4200 (Agilent Technologies, Santa Clara, CA, USA) to determine cDNA quality.

Ten microliters of cDNA product were mixed with 10 μL of nuclease-free water, 25 μL of Tagment DNA buffer, and 5 μL of Tagment DNA Enzyme (20034210; Illumina, San Diego, CA, USA),followed by incubation at 55°C for 5 min, holding at 10°C, purification using a MinElute PCR Purification Kit (28004; QIAGEN),and eluted in 25 μL of nuclease-free water.Then, 25 μL of fragmented cDNA was mixed with 15 μL of NPM(15027873;Illumina,San Diego,CA,USA),5 μL of P5_PCR(5 μM)(Table S1),and 5 μL of P7_PCR mix (5 μM) (Table S1) to run PCR thermocycling as follows:72°C for 3 min,95°C for 30 s,followed by 15 cycles at 95°C for 10 s,55°C for 30 s,72°C for 30 s,followed by a single cycle at 72°C for 5 min,and held at 4°C.After PCR amplification,the PCR products were purified sequentially using 0.6×and 0.2×ratios of VAHTS DNA clean beads to generate an Illumina-compatible cDNA library.After quality control by qPCR and Agilent 2100 Bioanalyzer, sequencing was performed using custom Read1 sequencing primers (Table S1) and universal Nextera Read2 sequencing primers on the Illumina®Xten (Illumina, San Diego,CA, USA) platform.

2.4.TCM-seq data analysis

TCM-seq data demultiplexing was performed using the barcode-splitter (v0.18.6, https://pypi.org/project/barcode-splitter/) Python package by providing the barcode list and running the default program.The Read2 sequence from each sample was then aligned to the corresponding genome reference (GRch38 or Rn6 reference) based on the species of samples with hisat2.After sorting the indexed SAM file into a BAM file with samtools[28]and counting the reads of uniquely mapped alignment for each sample via the featureCounts function in the subread package[29],a gene expression matrix of all samples was generated and used in downstream analysis.Principal component analysis was conducted to remove apparent outliers from each cell line.Then,the updated gene expression matrix was employed to analyze differentially expressed genes (DEGs) using the DESeq2 [30] R package for each treatment group, where triplicate samples represent a unique combination of CM and doses with medium treatment as a control.A comparative analysis of the similarity of expression profiles for replicate samples or samples from different species was performed by computing the Pearson correlation coefficient with R.

Unsupervised clustering analysis was performed using t-distribution stochastic neighbor embedding in the R package Seurat(v4.0.1) [31].Gene Ontology (GO) functional annotation and enrichment analysis were performed using clusterProfiler(v3.18.1)[32] via the compareCluster function.

Weighted correlation network analysis (WGCNA) was performed to characterize distinct gene modules across various cell lines using the WGCNA [33] R package with a matrix containing samples only treated with PGS/PNS, and the top 5000 DEGs were sorted by standard deviation.

2.5.Species-mixing experiment

H9c2 and AC16 cell lines were trypsinized at the logarithmic growth stage, seeded into 84 wells of a 384-well plate at a cell density of 5×103/well followed by incubation at 37°C for 24 h.The 168 wells/samples then were subjected to multiplexed RNA-seq,where the processes of cDNA barcoding and pooled cDNA library construction were the same as those described in Section 2.3.(sample well positions were detailed in Table S2).The sequencing data of individual wells obtained through demultiplexing according to predefined well-specific barcode sequences were separately aligned to both rat and human genome references,and the number of aligned reads was determined for each well/sample.

2.6.Reproductivity of TCM-seq

Ten replicates in each of the 11 cell lines (H9c2, HT22, A375,A549, MCF7, HepG2, HUVEC, PC3, U251, AC16, and HT29 cells)across three species seeded in 384-well plates for 24 h with 5 × 103cells per well were used to measure the reproducibility of TCM-seq.The RNA-seq with TCM-seq were performed simultaneously on all 110 samples.After demultiplexing the gene expression matrix for every replicate of each cell line, the Pearson correlations of the gene expression matrix of 10 replicates in each of the 11 cell lines were calculated to validate the robustness of TCMseq for characterizing the transcriptomes of multiplexed samples in parallel.

2.7.Extraction and quantitative analysis of PNS and PGS

Lyophilized PGS powder was prepared as follows:PG slices were successively decocted twice with water for 2 h, and then for 1.5 h.The decoction was filtered and pooled, and the filtrate was eluted with water to colorless using a D101 macroporous adsorbent resin column, and then eluted with 60% ethanol.The eluate was collected, and the filtrate was concentrated to an extract with a relative density of 1.06-1.08 (80°C), freeze-dried, and ground to obtain PGS powder.

Lyophilized PNS powder was prepared as follows: PN was crushed into coarse powder with DFT-50 A portable high-speed universal crusher, extracted twice with 70% ethanol for 7 h each time, then filtered; the filtrate was then concentrated under reduced pressure,and filtered again.The pooled filtrate was passed through a column of D101 adsorption resin.Once the adsorption was completed,it was washed with water that was discarded.Then the filtrate was eluted with 80% ethanol and the eluate was concentrated under reduced pressure, decolorized, refined, and concentrated again under reduced pressure to obtain PNS extract that was freeze-dried.

The primary saponins in PGS and PNS were determined employing a previously described method [34] using an Agilent 1100 HPLC system with an Agilent ZORBAX SB-C18column(250 mm × 4.6 mm, 5 μm).The approach was confirmed with strong correlation coefficient (R ＞0.9995), average recovery of 96.7%-102.3%, and relative standard deviation of 2.3%-4.4%.The injection volume for each sample was 20 μL.The flow rate of the mobile phase was 1 mL/min, and the column temperature was 26°C.The detection wavelength was adjusted to 203 nm.The mobile phases A and B were acetonitrile and water-phosphoric acid(100:0.05,V/V).The solvent gradient used was as follows:18%-20% A from 0 to 35 min,20%-21% A from 35 to 40 min,21%-26% A from 40 to 42 min, 26% A from 42 to 70 min, 26%-30% A from 70 to 85 min, 30%-34% A from 85 to 100 min, 34%-100% A from 100 to 101 min,and 100% A from 101 to 110 min.The re-equilibrium time was 10 min.The overall running time was 120 min.To establish calibration curves, mixed standard solutions comprising required concentrations of reference constituents were prepared in 50% (V/V) aqueous methanol.Methanol was used to dissolve lyophilized PGS/PNS powder.

2.8.Cell viability assay

Cells were trypsinized at the logarithmic growth stage and counted with 0.4% trypan blue stain(Gibco;15250061),then plated into 384-well plates at a cell density of 5×103/well and incubated at 37°C for 24 h.Cells treated with PGS or PNS(lyophilized powder dissolved in the corresponding medium) at doses of 0.2, 20, 200,and 400 μg/mL,and control cells were incubated at 37°C for 12 h.The supernatant was discarded from each well,and a drug toxicity test was performed using an 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) cell proliferation and cytotoxicity assay kit(Beyotime,shanghai,China).MTT solution(20 μL)was added to each well and incubated for 4 h.Then, 20 μL of the methanolic solution was added to each well and mixed.The absorbance of each well was measured at 580 nm using Infinite®M1000 PRO multimode reader (Tecan, M¨annedorf, Switzerland),and the cell viability(%)was calculated=(absorbance of treatment group/absorbance of control group) × 100%; triplicate measurements were performed.

2.9.Profiling PGS and PNS at multiple doses across cell lines

Ten cell lines (H9c2, A375, A549, MCF7, HepG2, HUVEC, PC3,U251,AC16,and HT29)were seeded into 384-well plates for 24 h at 5×103cells per well before PGS/PNS treatment.Each cell line was treated with four different doses(0.2,20,200,and 400 μg/mL)PGS and PNS in triplicate; the medium was used as a control.Detailed sample metadata were listed in Table S2.After treatment for 12 h,360 samples from the above cell lines were profiled by TCM-seq using the procedure described in Section 2.3.

2.10.Effects of PGS and PNS on the regulation of genes related to cardiovascular and cerebrovascular diseases

The cardiovascular disease database (CHD@ZJU, http://tcm.zju.edu.cn/chd) and cerebrovascular disease (CBD) database constructed by our group were used to further analyze the regulatory effects of PNS and PGS on cardiovascular and CBD networks,respectively.Cardiovascular disease (CDD)-related genes (1618 genes) from CHD@ZJU were used to investigate the effect of PGS/PNS on human cardiomyocytes (AC16) and HUVEC cell lines.CBDrelated genes (975 genes) of humans from CVD@ZJU were employed to determine the effect of PGS/PNS on human glioma(U251)and HUVEC cell lines.For each of the two cell lines tested for CDD, the gene intersection of 1618 disease genes and the DEGs generated after treatment with four different doses of PGS/PNS compared to the control were utilized to construct a molecular network based on protein-protein interactions (PPI) from STRING database (v11.5, https://cn.string-db.org/) using Cytoscape software.The topological parameters of the molecular network were analyzed using the“Tools-＞Analyze Network”option in Cytoscape,then the number of nodes, and the degree and betweenness centrality of nodes were obtained, respectively, as indicators to compare the regulatory strength of PGS and PNS.In addition, the network perturbation index (NPI) was designed and calculated as shown in the following equation to comprehensively compare the regulation of CDD by PGS and PNS.The analysis process for the regulatory strength of PGS and PNS in CBD was the same as that in CDD.

Where n denotes the number of nodes in the network, D denotes the degree of each node,and B denotes the betweenness centrality of each node.

2.11.Statistical analysis

Statistical analyses were conducted using GraphPad Prism 8.0,and R software.Unless otherwise specified,numeric data displayed in cell viability histograms were indicative of three repeated experiments and were presented as the mean ± standard error of mean.Difference between groups was considered statistically significant if the P ＜0.05.

3.Results

3.1.Design of TCM-seq

To fulfill the demand for parallel transcriptome sequencing of multiplexed samples, we developed TCM-seq, a sample-barcoding technique based on cDNA-specific labeling.TCM-seq primarily uses sample-specific oligonucleotide primers(barcoded primers)to label cDNAs of distinct samples before generating a cDNA library,and then pools cDNAs of all samples to simultaneously construct a mixed cDNA library and sequencing.The TCM-seq experimental technique was presented in Fig.1A.Normal or treated cells from a 384-well plate were first lysed in situ.Then, the cell lysate from each well was transferred to an RNA capture plate with oligo (dT)coated on the inner wall of each well.After hybridization with oligo(dT),the mRNA on the RNA capture plate was purified,eluted,and transferred to a 384-well PCR plate for RT reaction.A well-specific DNA sequence,also known as a barcode,was incorporated into the RT primer, and each well of the PCR plate was assigned a unique barcoded RT primer.The transcripts in each well were subjected to RT simultaneously; lysates from all wells were pooled for column purification, enzymatic cleavage, PCR enrichment, cDNA library generation, and sequencing.The barcoded oligo (dT) primers included a PCR handle, barcode, unique molecular identity (UMI),and oligo(dT)(Fig.1B).Oligo(dT)complements poly-A at the 3′end of mRNA.The barcode is a 10-nucleotide sequence used to uniquely identify the sample or well position of each cDNA.UMI is a 10-nucleotide random base sequence used to remove bias from PCR amplification and sequencing to precisely measure the number of transcripts.After complementary hybridization of the barcoded RT primer with mRNA,the cDNA first strand was synthesized at the 3′end of the RT primer,three C bases were also synthesized at the 3′end with reverse transcriptase,and using a TSO containing three G bases at the 5′end as a template, a sequence was synthesized to pair with the downstream PCR primers (Fig.1B).The barcoded cDNA from all wells was pooled into a mixed cDNA sample and preamplified to create a mixed cDNA library containing the transcript sequence information at the 5′end(Read2)and the sample barcode at the 3′end(Read1).Then,the mixed cDNA library was sequenced at an average depth of 2 million reads per sample.

Fig.1 .Schematic workflow for profiling transcriptome changes for multiplexed samples using RNA sequencing(TCM-seq).(A)The diagram of TCM-seq.Various cell lines are seeded into a 384-well plate and treated with multiple doses of various Chinese medicines (CMs).Then, the mRNA is extracted and followed by barcoding and the generation of complementary DNA(cDNA)library.(B)The molecular schematic diagram for multiplexed RNA-seq in TCM-seq.RT:reverse transcription;Exo I:exonuclease I;PCR:polymerase chain reaction; BC: barcode; TSO: template switch oligo; UMI: unique molecular identity.

3.2.Accuracy and reproducibility of TCM-seq for sample multiplexing

To test whether TCM-seq can accurately retrace the true sample identity of each cDNA in the pooled cDNA library based on the specific sample barcode, a rat and human species-mixing experiment was conducted.The sequencing data of the mixed H9c2 and AC16 cell lines,each of which was contained in 84 wells,were split according to barcodes.The results are shown in Fig.2A.In this figure, the orange and blue represented rat and human samples,respectively, identified by known well positions and barcode sequences of explicit species;the horizontal and vertical coordinates indicate the number of reads in each barcode-identified sample(small circles in the figure) that were matched to each of the two species.The results showed that the species classes of the samples,identified by barcodes and the corresponding Read2 matched to the species classes,were identical,demonstrating that the multiplexed transcriptome sequencing technology developed in this study can accurately tag cDNAs in each well using pre-incorporated barcode sequences, allowing the different tagged samples to be pooled for library construction.To evaluate the accuracy of the TCM-seq for labeling samples, the gene expression matrix obtained after aligning Read2 in each of the barcode-identified samples was subjected to Pearson correlation analysis and clustered for display.As shown in Fig.2B, all samples identified by barcodes were classified into two categories, with samples of the same species clustered together,a strong correlation between samples of the same species,and a poor correlation between samples of different species.These results further illustrate that sample barcoding in TCM-seq can accurately track the origin of individual samples within mixed counterparts.To validate the robustness of TCM-seq, 11 cell lines from three different species (human, mouse, and rat) were sequenced in mixed samples, and the expression profiles of 10 replicate wells from each cell line were correlated (Fig.2C).The correlation of gene expression profiles among the 10 wells of rat cardiomyocytes (H9c2) was extremely high (R ＞0.96), indicating that TCM-seq was highly stable in characterizing the transcriptional profiles of H9c2(Fig.2C).The correlation of gene expression profiles among the 10 wells of mouse cardiomyocytes (HT22) was also extremely high(R ＞0.92)(Fig.2D).In addition,the correlation coefficients of the expression profiles of the 10 wells of each of the nine human cell lines were also evaluated(Figs.2E and S1),and all the correlations were strong.These data reveal the robustness of TCM-seq for transcriptome profiling in cells derived from rats,mice, and humans.

Fig.2 .Species-mixing experiment and reproducibility of profiling transcriptome changes for multiplexed samples using RNA sequencing (TCM-seq) across various cell lines originating from different species.(A)Scatter plot depicting the number of reads from transcriptomes derived from a mixture of barcoded human(AC16)and rat(H9c2)cells.Points are colored based on the barcoded primer assignment.(B) Heatmap of Pearson's correlation coefficient of expression profiles between samples assigned to different species.Heatmap of Pearson's correlation coefficient of expression profiles across every 10 replicates of (C) H9c2, (D) HT22, and (E) AC16 cell lines.

3.3.TCM-seq identifies transcriptional responses to different doses of PGS and PNS

We extracted PGS and PNS, the primary active components of PG and PN(Fig.3A);we also determined the content of 10 saponins in the lyophilized extracts using our previously established method[34] to understand why PG and PN exhibited similar chemical compositions, their efficacy and clinical applications differed greatly.The results showed that the total content of the 10 saponins in PNS and PGS was 91.82%, and 68.97%, respectively.Whereas notoginsenoside R1 only existed in PNS, ginsenoside Rf, ginsenoside Rb2,and ginsenoside Rb3 were only detected in PGS(Table S3 and Fig.S2).To determine the safe dose range of PGS and PNS,four doses (0.2, 20, 200, and 400 μg/mL) were separately applied to 10 different cell lines, including A375 cells, for 12 h.The results showed that the viability of A375 cells treated with the four different doses of PGS and PNS exceeded 90% (Fig.3B).All other nine cell lines also maintained high survival rates (Fig.S3).After determining the safety of the four safe doses,we used TCM-seq to construct a transcriptional profile atlas of the 10 cell lines perturbed by different doses (0.2, 20, 200, and 400 μg/mL) of PGS and PNS(Fig.3C).All the cell lines were plated in 384-well plates for 24 h and treated with four different doses of PGS and PNS.The treated cells were then lysed in situ, the mRNA from each well was captured with an oligo (dT)-coated capture plate, and the primers with specific sample barcodes were subjected to RT reactions.cDNA from all wells was pooled for multiplexed cDNA library construction and sequencing, and the cDNA from the pooled samples was assigned back to their original well position and corresponding specific conditions using the demultiplexing deconvolution algorithm.

Fig.3 .Establishment of transcriptional atlas for cells perturbed by Panax ginseng saponins(PGS)and Panax notoginseng saponins(PNS).(A)Quantitative analysis of the ginsenosides Rg1, Re, Rf, Rb1, Rc, Ro, Rb2, Rb3, Rd, and notoginsenoside R1 in PGS and PNS using high performance liquid chromatography.(B) Cell viability of A375 control cells and those perturbed with different doses(0.2, 20,200, and 400 μg/mL)of PGS and PNS.Data were shown as the mean± standard error of mean,n=3 for each group.(C)Flowchart for the establishment of transcriptional atlas of cells perturbed by different doses of PGS and PNS.HepG2:human hepatocellular carcinoma cells;HUVEC:human umbilical vein endothelial cells; MCF7: Michigan Cancer Foundation-7; PCR: polymerase chain reaction.

3.4.Heterogeneous gene expression patterns in various cell lines exposed to PGS and PNS

TCM-seq was utilized to obtain the gene expression profiles of the control groups and the combination treatment groups.The latter groups included the treatment of cell lines with four doses and two CMs.We first analyzed the overlap of DEGs in the different combination treatment groups compared to the control counterpart in each cell line(Figs.4A and S4).Surprisingly,very few DEGs were common in the four doses of PGS and PNS treatment across all cell lines, whereas there were significantly more DEGs specific for each dose (Figs.4 and S4).GO pathway enrichment analysis was performed for DEGs in different combination treatment groups,and overlapping statistics were performed for the number of GO pathways that were enriched in the different dose treatment groups of each cell line(Figs.4B and S5).The results indicated that there was a greater number of GO pathways shared by both PGS and PNS treatment groups in each cell line.To determine whether tissue or species heterogeneity existed in the transcriptome changes of 10 cell lines treated with PGS and PNS.The expression matrix of DEGs from all combination treatments of the 10 cell lines was downscaled and unsupervised clustered (Fig.4C).The results showed that all samples (points on the graph)of each specific cell line were clustered into the same cluster, that is, samples from treatment groups of the same cell line were clustered into one cluster regardless of whether they were under the influence of PGS or PNS, and the differentiation pattern was consistent.Moreover,recoloring these clusters with the three treatments (control, PGSand PNS-treated) revealed that each cluster contained samples from the control, PGS- and PNS-treated groups (Fig.4C, right),indicating that the gene expression patterns of all samples were not clustered by the three treatments.To further investigate the biological functions of cells affected by PGS and PNS,we performed GO pathway enrichment analysis on all DEGs in each cell line treated with PGS or PNS(Fig.4D).PGS and PNS were found to significantly affect endoplasmic reticulum stress, proteasome-related proteolytic metabolic processes, protein complex catabolism, and transcriptional splicing.All of which are essential pathways for cellular metabolism, indicating that both saponins may have regulatory effects on protein metabolism and transcriptional splicing in cells,regardless of the cell line type, and the regulatory effect of PGS surpasses that of PNS(Fig.4D,first bar chart).However,the effects of PNS and PGS on energy metabolic pathways in the mitochondria,ATP metabolism, oxidative phosphorylation, and electron transfer varied significantly between cell lines.In H9c2 cells, HepG2 cells,HUVEC cells, and U251 cells, both PGS and PNS affected mitochondrial energy metabolism-related pathways, including mitochondrial respiratory chain complex assembly, electron transfer,organic compound oxidation, mitochondrial ATP synthesis, oxidative phosphorylation, and ATP metabolism; whereas in A375 cell,only PNS regulated energy metabolism in mitochondria.In contrast, only PGS perturbed mitochondrial energy metabolism in A549 and AC16 cell lines,and its effect was greater than that of PNS,regardless of cell line type(Fig.4D,second bar chart).The pathways associated with protein localization in the endoplasmic reticulum or cell membrane were significantly different between cell lines treated with PGS and PNS.These pathways were enriched in HepG2, HUVEC, and PC3 cell lines after perturbation by both PNS and PGS.The protein membrane-localization-related pathways were mainly enriched in MCF7 cell line treated with PNS.PGS more strongly regulated protein localization in the membrane surface than PNS, regardless of the cell line type (Fig.4D, third bar chart).Moreover, we discovered that both PGS and PNS, specifically in H9c2 cells,acted on autophagy,nucleotide metabolism,and protein localization pathways.

Fig.4 .Different cell lines vary in transcriptional response to different doses of Panax ginseng saponins(PGS)and Panax notoginseng saponins(PNS).(A)Intersection statistics of the number of differentially expressed genes(DEGs)in AC16 cells perturbed by different doses(0.2,20,200,and 400 μg/mL)of PGS and PNS compared to the control.(B)Intersection statistics of the number of enriched gene ontology(GO)terms by the DEGs in AC16 cells perturbed by different doses(0.2,20,200,and 400 μg/mL)of PGS and PNS compared to the control.(C)t-distributed stochastic neighbor embedding(t-SNE)analysis of the samples perturbed by PGS and PNS and control samples across 10 cell lines.In the left t-SNE map,10 cell line clusters are labeled with different colors;in the right t-SNE map,treatments of PGS,PNS,and control are labeled with different colors.(D)Top 10 GO terms enriched by DEGs in different cell lines after PGS and PNS perturbation compared to control.The bars on the right represent the sum of the number of GO pathways in cell lines treated with PGS and PNS.(E)Heatmap of quantitative module expression(ME)of 12 gene modules from the PNS dataset(left)and 17 gene modules from the PGS dataset(right).ME value indicates the expression feature value of all genes in the module.The MEgrey(grey module)is reserved for genes outside of all modules.(F)Heatmap of the number of shared genes between each gene module originated from the PNS and PGS datasets.The numbers circled in purple and the corresponding gene modules are primarily discussed in the main text.HepG2:human hepatocellular carcinoma cells;HUVEC:human umbilical vein endothelial cells;MCF7:Michigan Cancer Foundation-7;ATP:adenosine triphosphate;SRP:signal recognition particle; ER: endoplasmic reticulum.

In addition, we utilized WGCNA analysis to investigate the coexpressed gene modules in both PGS- and PNS-treated transcriptomic datasets.After filtering out genes with module membership of less than 0.6, 12 gene modules were identified in the PNS-treated dataset (Fig.S6A) and 17 gene modules were identified in the PGS-treated dataset (Fig.S6B).Next, we investigated whether there were cell line-specific gene modules regulated by PNS or PGS by quantifying the overall expression characteristics of the gene modules in different cell lines for both datasets (Fig.4E).Except for HUVEC, which did not form cell line-specific gene modules after PGS treatment, all other cells treated with PGS and PNS formed cell line-specific gene modules.To compare the differences between the gene modules,the number of genes shared in each module in both datasets was counted.We noticed that there were few common genes in the cell line-specific gene modules affected by both PNS and PGS,whereas the majority of genes were shared in the PC3 and H9c2 cell-specific gene modules(Fig.4F).

3.5.Stronger regulation of CDD network by PGS compared to PNS

To compare the regulatory effects of PGS and PNS on CDD and CBD, based on CHD@ZJU (http://tcm.zju.edu.cn/chd/) and CBD database previously established by our group,we used the network pharmacology approach to further compare the differences in pharmacological effects between PGS and PNS.Since cardiomyocytes and vascular endothelial cells play a role in CDDs,we selected AC16 and HUVEC cell lines to investigate the regulatory effects of PGS and PNS on CDD-related genes, whereas U251 and HUVEC were selected to investigate the regulatory effects of both saponins on CBD-related genes.The DEGs in AC16 and HUVEC cell lines treated with PGS and PNS were intersected with CDD-related genes, and then the intersected genes were used to separately construct perturbed molecular networks based on PPI from the STRING database(Fig.5A).The intersected genes from intersecting DEGs,in U251 and HUVEC cell lines treated with PGS and PNS,with CBD-related genes were also employed to separately construct perturbed molecular networks.(Fig.5B).The number of nodes,degree of nodes, and NPI in these molecular networks were counted, respectively, to quantitatively compare the regulatory effects of PGS and PNS on CDD and CBD genes.The results showed that the number of nodes,degree of nodes,and NPI of PGS acting on CDD genes in AC16 and HUVEC were greater than that of PNS(Fig.5C).In contrast,the number of nodes,degree of nodes,and NPI of PNS acting on CBD genes in U251 cells were greater than that of PGS,whereas that of PGS acting on CBD genes in HUVEC cells were greater(Fig.5D).In addition,GO pathway enrichment analysis was performed for CDD genes regulated by PGS and PNS in AC16 and HUVEC(Fig.5E).The figure indicated the top 10 GO pathways were significantly enriched by CDD genes that were regulated by PGS and PNS,indicating that both saponins affected cell aging,coagulation,response to low blood oxygen levels, negative regulation of apoptosis, and oxidative stress in both cell lines.Unlike PGS, PNS positively regulated apoptotic signaling pathways and protein localization.In addition, both PGS and PNS affected intracellular protein transport in AC16 cells, but not in HUVEC cell line.GO pathway enrichment analysis was also performed for CBD-related genes regulated by PGS and PNS in U251 and HUVEC (Fig.5F).This demonstrated that both saponins affected the processes of response to hypoxia and oxidative stress in both cell lines;both PGS and PNS affected the coagulation regulatory process in HUVEC.The latter effect was not observed in U251 glioma cells, indicating that PGS and PNS regulate coagulation via possible cellular specificity,and the effect of PNS on coagulation was significantly stronger than that of PGS.In addition,PGS,but not PNS,was found to affect glial cell apoptosis in HUVEC and U251 cells.

Fig.5 .Network regulatory analysis predicts and compares the effects of Panax ginseng saponins (PGS) and Panax notoginseng saponins (PNS) on cardiovascular disease (CDD) and cerebrovascular disease(CBD).(A)Protein-protein interaction(PPI)maps consist of CDD-related genes affected by PGS and PNS.(B)PPI maps consist of CBD-related genes affected by PGS and PNS.(C)Statistical analysis of the number of nodes(left),degree(middle),and network perturbation index(NPI,right)of the molecular networks in(A).(D)Statistical analysis of the number of nodes(left),degree(middle),and NPI(right)of the molecular networks in(B).Top 10 gene ontology(GO)terms enriched by(E)CDD-related and(F)CBDrelated genes affected by PGS and PNS.The GO pathways with a purple background were discussed in detail in the main text.HUVEC: human umbilical vein endothelial cells.

4.Discussion

PG and PN are valuable Chinese herbs and widely used in China.Modern pharmacological studies have demonstrated that the main active ingredients of PG and PN are dammarane-type ginsenosides,but they have different pharmacological effects and distinct clinical indications.In this study, we first designed TCM-seq based on a DNA barcoding strategy to simultaneously perform RT reactions,PCR amplification, and cDNA library construction on up to 468 samples to characterize the transcriptional response to multiple doses of CMs, and verified the accuracy of TCM-seq in labeling cDNAs from 168 different samples employing species-mixing experiments.TCM-seq could accurately classify each sample into its real species using pre-defined well/species-specific oligo barcodes.We also demonstrated the excellent robustness of TCM-seq by comparing the similarity of transcriptional profiles of 10 replicate samples from each of the 11 cell lines across the three species of humans, mice, and rats.We then used TCM-seq to profile the treatment effect, using four safe doses of PGS and PNS, which are the main active components of PG and PN,respectively,on each of the 10 cell lines;MCF7,PC3,A549,A375,HepG2,and HT29 cell lines were used as references from the CMap database[35,36],and H9c2,AC16, HUVEC, and U251 cell lines were used as model cells for cardiovascular and cerebrovascular studies to compare the differences in efficacy between PGS and PNS at the molecular level.We found that the DEGs in each cell line,treated with the four different doses,were not dose-dependent,while there was a large overlap in the enriched GO pathways, and the unsupervised clustering of differential expression profiles of all cell lines treated with PGS and PNS exhibited cell line-heterogeneous expression pattern.In terms of regulated biological functions,both PGS and PNS had regulatory effects on protein metabolism and transcript splicing in cells, and the regulatory effects of PGS were stronger.There was a significant cell line heterogeneity in the regulatory effects of PGS and PNS on mitochondrial energy metabolism-related pathways and protein membrane-localization-related pathways.Overall, the regulatory effects of PGS were stronger.Moreover, PGS and PNS specifically affected pathways related to autophagy, nucleotide metabolism,and protein localization in H9c2 cells.Using WGCNA analysis of DEGs perturbed by PGS and PNS, we identified cell line-specific gene modules affected by both saponins.PC3 and H9c2 cellspecific gene modules regulated by PGS and PNS shared the majority of genes,while other cell line-specific gene modules had few genes in common.Finally, this study focused on the differences in the regulatory effects of PGS and PNS on CDD- and CBD-related genes, and found that PGS showed stronger regulatory effects on cardiovascular disease.In CBD, the regulatory strength of PGS and PNS on HUVEC and U251 cells differed, and the pathway enrichment results also showed that PNS had stronger coagulation effects on HUVEC.

It is worth mentioning that the few overlapping genes from DEGs spanning four doses and the majority of biological GO pathways enriched by DEGs were shared across the four doses in each cell line, suggesting that the mode of action of CMs in cells should be investigated systematically.The differential expression profiles of each cell line with PGS or PNS perturbation showed evident cell specificity,which is consistent with both holistic and tissue-specific actions of CM and suggests that it is crucial to select suitable model cells for the study of specific pharmacological effects of CM.In this study, we found that the effects of PGS and PNS on rat cardiomyocytes (H9c2) differed from other cell lines in their specific effects on biological pathways, including autophagy, nucleotide metabolism, and protein localization in organelle-related pathways.A previous study demonstrated that PNS exerts anti-acute myocardial infarction effects by regulating autophagy in a mouse model [4], and has a protective effect on mitochondrial damage in rat cardiomyocytes due to autophagy inhibition.Moreover, PG has been reported to promote autophagy in aged mice [37].These results are consistent with the pathways specifically affected in H9c2 cells in our study.This also suggests that characterization of the gene expression profile of cells in vitro after CM treatment,without cell modeling, may predict the effect of CM and its underlying molecular mechanisms.WGCNA analysis revealed that two cell line-specific gene modules in H9c2 and PC3 cells treated with PGS and PNS had similar genes.Finally, our study also found that PGS exerted a stronger regulatory effect on cardiovascular disease,which is consistent with the clinical application of ginseng for heart failure.The analysis of pathway enrichment for DEGs in CDD or CBD showed that both PGS and PNS affected processes such as cellular senescence, coagulation, response to hypoxia, negative regulation of apoptosis, and oxidative stress related to CDD, which is consistent with that ginsenoside Rg3 prevented myocardial hypertrophy by inhibiting oxidative stress through modulation of the SIRT1/NFκB pathway in an AC16 cell model [38], and Korean red ginseng extract could inhibit HUVEC cell senescence[39].Both PGS and PNS regulate CBD-related genes involved in response to hypoxia and oxidative stress.Both PGS and PNS regulated coagulation processes in HUVEC endothelial cells,but not in U251 glioma cells,suggesting that the saponins regulate coagulation processes with some cell specificity.This is also consistent with the involvement of vascular endothelial cells in the coagulation processes;PNS regulated DEGs in CBD more strongly than PGS, which is consistent with current research and clinical applications.It has been reported that ginsenoside Rd can reverse the nicotine-induced decrease in NO and reduce angiotensin II production in HUVEC cells,a process that may be related to coagulation regulation[40].Furthermore,unlike PNS,we found that PGS affected the glial cell apoptotic process in HUVEC and U251 cells,which is consistent with the neuroprotective effects of ginsenosides Rb1, Rg3, and Rg1 through the activation of microglia and inhibition of apoptosis[41,42].

Traditional approaches for investigating the efficacy of PG and PN consist mostly of chemical activity assays,in vitro cellular model evaluation, and in vivo animal models that are constrained by a single evaluation indicator, or are time-consuming and have low throughput.To the best of our knowledge,this is the first study to employ comprehensive parallel multiplexed transcriptome sequencing technology to compare the efficacies of PGS and PNS.In addition, TCM-seq is theoretically applicable to organoids, zebrafish,and even in vivo tissue cells,and not just in vitro cell lines.Due to time limitations, this study was unable to examine the pharmacological effects of PG and PN in other disease models,nor could it evaluate other active components in PN and PG, including polysaccharides.Furthermore, the geographical differences in CM materials that may result in different contents of saponins can be considered for a more general analysis of constituent differences.Subsequent studies can also compare the efficacy of PG, PN, and their related preparations in ex vivo evaluation models based on the relationship between the components and the compatibility principle to further refine the scientific meaning of the difference in efficacy between PG and PN.

5.Conclusions

In this study, we developed a multiplexed transcriptome sequencing technology for up to 468 samples treated with CMs based on a DNA barcoding strategy called TCM-seq.We first verified the accuracy of sample multiplexing in TCM-seq by species-mixing experiments, and demonstrated that TCM-seq was robust across 11 cell lines from human,mouse,and rat species.A comprehensive comparative study was conducted using TCM-seq to profile transcriptional changes in 10 cell lines from multiple tissues treated with four different doses of PGS and PNS.DEGs analysis, unsupervised clustering of expression profiles, WGCNA analysis, and biological functional analysis were performed to systematically compare the perturbing effects of PNS and PGS on the transcriptome profiles of the 10 cell lines studied.Finally,we adopted a network pharmacology approach to predict and compare the similarities and differences between PGS and PNS in the regulation of CDD and CBD networks and elucidated the regulatory effects of PGS and PNS in terms of biological functions.In summary, this article provides a high-throughput and low-cost technique and paradigm for the systematic study and comparison of the pharmacological effects and molecular mechanisms of PG and PN.This is the first comprehensive multiplexed RNA-seq study to simultaneously evaluate hundreds of PG and PN sample conditions using transcriptomic changes.In addition, TCM-seq is theoretically applicable to organoids, zebrafish, and even in vivo tissue cells,allowing for additional in vivo research on PG and PN samples.However,this study did not analyze the pharmacological effects of PG and PN in disease models, evaluate other active components such as polysaccharides, or take into account the geographical variations of CM materials.Subsequent investigations can further analyze the interaction of constituents and the efficacy of PG, PN,and their related preparations based on the compatibility principle and the interplay between components to further refine the scientific significance of the difference in efficacy between PG and PN.

CRediT author statement

Junyun Cheng:Conceptualization, Writing - Original draft preparation;Jie Chen:Methodology, Writing - Reviewing and Editing;Jie Liao:Visualization;Tianhao Wang:Investigation;Xin Shao:Formal analysis;Jinbo LongandPenghui Yang:Validation;Anyao Li:Software;Zheng Wang:Funding acquisition;Xiaoyan Lu:Supervision, Project administration;Xiaohui Fan:Conceptualization, Supervision, Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos.: 81973701 and 81903767), the Innovation Team and Talents Cultivation Program of National Administration of Traditional Chinese Medicine (Grant No.:ZYYCXTD-D-202002), and the Natural Science Foundation of Zhejiang Province(Grant No.:LZ20H290002).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jpha.2023.02.009.