• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A proteomic landscape of pharmacologic perturbations for functional relevance

    2024-03-21 05:51:14ZhiweiLiuShngwenJingBingingHoShuyuXieYingluoLiuYuqiHungHengXuChengLuoMinHungMinjiTnJunYuXu
    Journal of Pharmaceutical Analysis 2024年1期

    Zhiwei Liu , Shngwen Jing , Binging Ho , Shuyu Xie , Yingluo Liu , Yuqi Hung ,Heng Xu , Cheng Luo ,, Min Hung , Minji Tn ,,c,d, Jun-Yu Xu ,,c,*

    a State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai, 201203, China

    b School of Chinese Materia Medica, Nanjing University of Chinese Medicine, Nanjing, 210023, China

    c Zhongshan Institute for Drug Discovery, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Zhongshan, Guangdong, 528400, China

    d State Key Laboratory of Pharmaceutical Biotechnology, Nanjing University, Nanjing, 210023, China

    Keywords:Proteomics Drug Perturbation Drug target Drug combination

    A B S T R A C T Pharmacological perturbation studies based on protein-level signatures are fundamental for drug discovery.In the present study,we used a mass spectrometry(MS)-based proteomic platform to profile the whole proteome of the breast cancer MCF7 cell line under stress induced by 78 bioactive compounds.The integrated analysis of perturbed signal abundance revealed the connectivity between phenotypic behaviors and molecular features in cancer cells.Our data showed functional relevance in exploring the novel pharmacological activity of phenolic xanthohumol,as well as the noncanonical targets of clinically approved tamoxifen, lovastatin, and their derivatives.Furthermore, the rational design of synergistic inhibition using a combination of histone methyltransferase and topoisomerase was identified based on their complementary drug fingerprints.This study provides rich resources for the proteomic landscape of drug responses for precision therapeutic medicine.

    1.Introduction

    Molecular profiling of drug-response phenotypes expands our knowledge of cancer pharmacology and provides guidance for precision medicine [1].Systematic molecular profiling is a useful approach to understand the mechanisms of action(MOA)of drugs in cancer cell lines[2].Currently,different layers of omics data are considered to reflect cell-line perturbation after compound treatment.These layers include genomic data at the DNA or messenger RNA (mRNA) level, as well as metabolomics and proteomics data.As end products encoded by genes,targets of small molecules,and direct executors of all biological functions,proteins are superior in acting as functional and proximal readouts that directly reveal drug actions [3].Genomic sequences are not disturbed during shortterm drug stimulation, which restricts their application in deciphering molecular mechanisms.Metabolites and transcripts do not provide a direct readout of the perturbagen response,limiting their ability to capture slight signal changes.In recent years,breakthroughs in mass spectrometry (MS)-based proteomic technology have promoted the technical feasibility of exploring unbiased proteome perturbations in drug responses.To date,large-scale proteomic drug responses as signature profiles and the integration of drug sensitivity in perturbed cancer cell lines have been reported[4-7].These preliminary studies explored the current understanding of potential biomarker and drug target discovery.However, some important drug types such as the immunomodulatory drugs (IMiDs) were not included in these studies.Recently, a proteome-wide atlas of drug MOA in the colon cancer cell line HCT116 has been reported[8];however,the cancer type of the cell lines dominates the molecular classification[9].Additionally,these databases have not been exhaustively utilized for their potential applications in pharmacological exploration, particularly for exploring novel drug MOA and predicting drug combinations.

    To fill this gap, we used a Tandem Mass Tag (TMT)-based proteomic strategy to establish proteomic signature data for 78 distinct bioactive small molecules in the breast cancer epithelial MCF7 cell line.We then compared and analyzed the signature of the published proteomic data underlying compound treatment with our dataset.Based on the similarity of protein expression profiles, we uncovered novel mechanism or “off-target” effects of bioactive compounds.Furthermore, we revealed and validated a rational design for a drug combination strategy.Our results provide a resource for the systematic characterization of the molecular responses to drug treatments.These datasets can serve as functional proteomic platforms to facilitate drug discovery and inform new therapeutic approaches.

    2.Experimental

    2.1.Compounds

    The compounds were obtained commercially or from the compound library of the Institute of Medicine, Chinese Academy of Sciences.The details of these compounds are listed in Table S1.

    2.2.Cell lines

    MCF7 cells were cultured in Dulbecco's Modified Eagle's Medium (DMEM; Gibco, Waltham, MA, USA) with 2 mM L-glutamine(Amresco,Solon,OH,USA),10%fetal bovine serum(FBS;Gibco),and 100 units/mL of penicillin/streptomycin (Meilunbio,Dalian,China)and incubated at 37°C in 5% CO2.Cells were collected after drug treatment for 12 h.Three experimental periods were conducted,the details of which are listed in Table S1.

    2.3.Quantitative proteomics

    The cells were cultured in a 15-cm dish at approximately 90%confluency, washed with 1× phosphate-buffered saline (PBS,Meilunbio)and lysed using freshly prepared lysis buffer(8 M urea,100 mM NH4HCO3,and 2×protease inhibitor cocktail(Roche,Basel,Switzerland); pH 8.0).The cell lysate was incubated at 4°C for 30 min and then sonicated for 2 min with 30% energy (ultrasonic cell disruptor JY92-II, Ningbo Scientz Biotechnology Co., Ltd.,Ningbo, China).

    For each sample,50 μg of proteins(determined by bicinchoninic acid assay) was reduced with 5 mM dithiothreitol (Sigma-Aldrich,St.Louis, MO, USA) for 30 min at 56°C.Cysteine was alkylated by incubating with 15 mM iodoacetamide (Sigma-Aldrich) at room temperature for 30 min.The reaction was terminated by incubation with 30 mM L-cysteine (Sigma-Aldrich) for 30 min at room temperature.Protein lysis was diluted to 2 M with 100 mM NH4HCO3and digested with trypsin (Hualishi Scientific, Beijing, China) at a ratio of 50:1(protein:enzyme,m/m)for 16 h at 37°C.The samples were treated with trypsin(Hualishi Scientific)at a protein:enzyme ratio of 100:1 (m/m) for 4 h.The digested lysates were desalted using tC18 cartridges (Sep-Pak, Waters, Milford, MA, USA) and dried using a SpeedVac.

    TMT(Thermo Fisher Scientific Inc.,Waltham,MA,USA)reagents(6-plex or 10-plex)were used for peptides labeling according to the manufacturer's protocol.Simply,0.4 mg of TMT reagent was added to each sample and then incubated at room temperature with mild shake for 60 min.The reaction was terminated by 4 μL of hydroxylamine(5%,m/m)and mixed.The peptides were desalted by tC18 cartridges (Sep-Pak, Waters).The same control lysate (MCF7 cell lines treated with dimethyl sulfoxide (DMSO; WAK-Chemie Medical GmbH, Jena, Germany) was used across all TMT-6plex groups with two biological replicates.

    The mixed samples were then re-suspended in 0.1% (V/V) trifluoroacetic acid (Sigma-Aldrich), and separated into 80 fractions using high pH reverse-phase high performance liquid chromatography (HPLC) with XBridge Peptide BEH C18column (5 μm,4.6 mm × 250 mm,130 ?, Waters).The flow rate was set at 1 mL/min.The gradient was set to 90 min from 0%to 95%of buffer(98%(V/V) acetonitrile (ACN; pH 8.5; Thermo Fisher Scientific Inc.).The 80 fractions were concatenated into 20 fractions and dried using a SpeedVac.

    2.4.Liquid chromatography-mass spectrometry (LC-MS) analysis

    Each fraction of the dried peptides was re-suspended in formic acid (0.1%, V/V) and separated using an EASY-nLC 1000 LC system(Thermo Fisher Scientific Inc.).The gradient was set to 70 min from 8%to 46%ACN.A home-made capillary column(75 μm i.d.×20 cm length) packed with C18beads (size 1.8 μm, 100 ?, Dikma Technologies, Beijing, China) was used.

    The peptides were analyzed using an Orbitrap Fusion mass spectrometer (Thermo Fisher Scientific Inc.).The full scan range was m/z 450-1500.The maximum injection times were set to 50 and 80 ms for the full and MS/MS scans,respectively.The resolution at m/z 200 was set to 60,000 and 15,000 for the full and MS/MS scans, respectively.The automatic gain control targets was set to 500,000.Fifteen top ions were used for analysis.The high collision dissociation normalized collision energy was set to 40%.

    2.5.MS data pre-processing

    MaxQuant software(1.5.3.8)with human database from UniProt(96,447 sequences,downloaded on June 6,2019)was used for raw data parsing.The parameters were set with an false discovery rate(FDR) of 0.01 at the level of protein, peptide, and modification.Enzyme specificity was set to trypsin/P with tolerance of 2 missed cleavages.The fixed modification was carbamidomethyl (cysteine residues) with acetyl (protein N-term) and oxidation (M) set as variable modifications.TMT label was set to 6-plex TMT for proteomic data.The MaxQuant settings for the 10-plex labelled samples were the same with that of 6-plex labelled samples.The MaxQuant setting of 6-plex labelled samples was used for our drugresponse proteomic landscape.In addition,the MaxQuant setting of 10-plex labelled samples was used for the drug combination assays(EZH2i and topotecan).

    The reverse or potential contaminant proteins were deleted.The intensities from the same gene were summed and then normalized by the median intensity for each sample.Thirty-two expression matrices were combined by gene symbols.The ratio of sample(drug treatment) abundance to control (DMSO) abundance was calculated as the relative abundance and then log2 transformed.

    To balance the confidence and reservation of the scope of data set to low abundant proteins,the data set was filtered for proteins measured in 50%of all samples.The missing values was imputed by the minimum value of the data set.The batch effect was removed by the R tool,ComBat[10,11].We used time point information(P1,P2, and P3) as batch covariate while other parameters are default(par.prior = TRUE, mean.only = FALSE, BPPARAM = bpparam(“SerialParam”), mod = NULL, ref.batch = NULL) for “Combat”function.

    2.6.Protein-protein interaction (PPI) networks

    The PPI wasanalyzed using the STRING website(https://string-db.org/)using all of the commonly downregulated/upregulated proteins(Student's t-test,P <0.05 and fold change >1.5 in at least 60%drugs for each class type).In total,201 downregulated proteins(59 in the adenosine monophosphate (AMP)-activated protein kinase (AMPK)inhibitor,37 in the bromodomain and extraterminal(BET)inhibitor,one in the cyclin-dependent kinase (CDK) inhibitor, 22 in the DNA alkylating drug,6 in the DNA/RNA synthesis inhibitor,15 in the IMiDs,38 in the phosphatidylinositol-3-kinase (PI3K) inhibitor, 11 in the proteasome inhibitor,10 in the topoisomerase inhibitor,and 2 in the tubulin inhibitor) were identified, resulting in 160 unique proteins.Similarly,229 upregulated proteins(42 in the AMPK inhibitor,38 in the BET inhibitor, 3 in the CDK inhibitor,15 in the DNA alkylating drug, 5 in the DNA/RNA synthesis inhibitor, 2 in the histone deacetylase(HDAC)inhibitor,48 in the IMiDs,26 in the PI3K inhibitor, 28 in the proteasome inhibitor, 16 in the topoisomerase inhibitor,and 6 in the tubulin inhibitor)were identified,resulting in 187 unique proteins.The PPI networks were plotted by the R tools“igraph”.

    2.7.The “querying” process

    The querying signatures were obtained from the proteome data of public articles.Proteomic data were downloaded and mapped to our database using gene names.Significantly regulated genes(Student's t-test,P <0.05,fold-change >1.5)were used as querying signatures and were compared with our database using Spearman's correlation.The similarity between the querying signatures and our database was defined using the following criteria:1)the P value of Spearman's correlation should be less than 0.05; 2) the absolute correlation coefficient of Spearman's correlation should be larger than 0.2;and 3)we mainly considered the top five correlated drugs in our database.

    2.8.Immunofluorescence and assessment

    MCF7 cells were plated in 12-well dishes on glass coverslips and treated with colchicine (10 or 20 nM) or xanthohumol (5, 10, or 15 μM)for 24 h before staining with Tubulin-Tracker Red.Then,the cells were washed in PBS and fixed with 4% paraformaldehyde for 30 min.Next,the fixed cells were washed with 0.1%(V/V)Triton X-100 in PBS for 5 min; this was repeated twice.Subsequently, cells were incubated with Tubulin-Tracker Red for 1 h in dark and then washed thrice with 0.1% (V/V) Triton X-100 in PBS.Immunofluorescence was examined using confocal microscopy (Olympus FV1000) and photographed at ×40 or ×100.Cells were counted from six images per well (×40) with a minimum of 100 cells in total.

    2.9.Determination of 20S proteasome activity

    In vitro assay was performed to estimate 20S proteasome activity after treatment with different compounds.In Tris-HCl buffer(50 mM, pH 7.5), 250 ng commercial human 20S proteasome was reacted with 100 μM Suc-LLVY-AMC at 37°C for about 30 min with the treatment of 10 μM tamoxifen or 4-hydroxy tamoxifen using 0.5 μM carfilzomib as control.Hydrolyzed 7-amino-4-methylcoumarin (AMC) was measured using a microplate reader(Spark Cyto, Tecan, Zürich, Switzerland) with a filter of 360 nm excitation and 460 nm emission.

    2.10.Cell viability assay

    Different cell lines (CRBN+/+HCT-116 or SMMC-7721) were cultured in 96-well plates for 24 h with different compounds(including CC-885, benzethonium chloride, MLN4924, and statin drugs).After 48 h of treatment, the absorbance (450 nm) was measured for Cell Counting Kit-8 (Meilunbio) on a plate reader(Infinite 200pro,Tecan).Each assay was repeated no less than three times.The combination index (CI) was estimated using the CompuSyn (version 1.0) to define the potential synergistic effect.

    2.11.Gene set enrichment analysis

    Proteomic profiling data of SMCC-7721 cells treated under different conditions (DMSO, topotecan, and topotecan + GSK126)were used for gene set enrichment analysis (GSEA).The GSEA was performed using the oncogenic gene set(C6) from MsigDB(v 6.2).The permutation type was selected as “gene-set”.

    2.12.Docking procedure

    Docking was performed using Schr¨odinger (version 2.4, New York City,NY,USA,2010).The LigPrep module was used to prepare small molecules.The Epik was used to determine the protonation states of the molecules [12].Protein preparation and grid generation was conducted using crystal structure of DDB1-CRBN E3 ubiquitin ligase bound to lenalidomide (Protein Data Bank (PDB)ID: 4CI2) before docking.Glide module (version 6) was used to perform docking procedure with standard precision mode.The binding interactions were analyzed and displayed using PyMOL(version 1.8).

    2.13.Data and materials availability

    All MS raw data have been deposited to the iProX Consortium with Project ID:IPX0003183000(URL:https://www.iprox.cn/page/project.html?id=IPX0003183000).

    3.Results

    3.1.Construction of the dynamic proteomic landscape of drug response in MCF7 cell line

    To establish a systematic study of perturbed proteomes induced by the 78 distinct bioactive compounds, we used a 6-plex TMT labelled proteomic approach with biological replicates in the MCF7 cell line(Fig.1A).An overview of these small molecules is provided in Table S1, which covers almost all types of anticancer drugs.The targets of the selected compounds were involved in a broad range of signaling pathways, tyrosine kinases, proteasomes, DNA synthesis/repair enzymes and epigenetic enzymes (Fig.1B and Table S1).In addition, the proteome profiles of IMiDs such as lenalidomide, thalidomide, and their derivatives, which target protein degradation, were collected in our dataset (Fig.1B).The concentrations of the compounds were selected based on their 10-fold half-maximal inhibitory concentration (IC50) values.For compounds with no apparent effect on cell viability (IC50≥1 μM),dosing concentrations were controlled at 10 μM for proteomic analysis(Table S1).The duration of the compounds was set to 12 h to avoid non-apparent signals or extensive cell death effects when the profiles were obtained too early or late.

    A total of 32 batches of TMT-labelled samples were analyzed in our study.The resulting dataset contained quantitative proteomic data for >15,000 proteins(8,889 genes <1%false discovery rate,and≥1 unique peptide) (Table S2).The average number of proteins identified in each TMT group was approximately 6,400 (Figs.S1A and B).All the proteomic data presented a unimodal distribution(Fig.S1C) and passed the quality control process (Figs.S1A-E).Because the different time points for the experiments (P1, P2, and P3) that could result in the batch effect when data from multiple batches were integrated, we removed it by using a widely used Combat method (Figs.S1D and E) according to the previous study[13-15].An overview of proteomic profiling with different drug treatments revealed the dynamics of regulated proteins especially down-regulated proteins, as well as the highest and lowest abundant proteins identified by the median protein regulation across all compounds (Fig.1C).The Spearman's correlation coefficient of protein expression was found to be much higher in corresponding compounds with similar pharmacological mechanisms than in those with distinct pharmacological mechanism (Fig.1D).In addition,principal component analysis(PCA)was performed on the signatures of all the 78 compounds.As expected,most compounds with similar mechanisms clustered closer than the others(Fig.1E),such as cytotoxic drugs, epigenetic inhibitors, and inhibitors targeting the ubiquitin-proteasome system.Compared with these small molecules,compounds targeting the cell cycle were far from each other in the PCA analysis(Fig.1E),which may indicate a wide variation in the diverse mechanisms and biological effects of CDK inhibitors.Known IMiDs-dependent targets [16] were downregulated after treatment of IMiDs; however, some of them were not quantified which may be due to cell-type specificity(Fig.S1F).

    Fig.1.Proteomic landscape of drug treatment.(A) Workflow of the proteomic profiling and bioinformatic analysis.(B) Summary of the information for compounds we used.(C)Overview of the proteomic profiling with different drug treatment.Shown are the dynamics of protein abundances (log2 ratio).The highest- and lowest-abundance proteins,identified by the median of protein regulation across all compounds,are shown in the upper right and bottom-left box respectively.(D)The spearman correlation between samples with different drug treatment.The color represents different class of drugs as shown in Fig.1B.(E)The results of principal components analysis.The color represents different class of drugs as shown in Fig.1B.TMT:tandem mass tag;AMPK:adenosine monophosphate(AMP)-activated protein kinase;BET:bromodomain and extra-terminal domain;CDK:cyclin dependent kinase; HDAC: histone deacetylase; IMiDs: immunomodulatory drugs; PI3K: phosphatidylinositol-3-kinase; Dim: dimension.

    3.2.The application scope of our proteomic signature resource

    We further estimated the application potential of our database to understand the underlying pharmacological mechanisms using proteomic data from different cell lines,different processes of drug treatment, and different previously reported proteomic platforms.For these analyses, proteomic signatures from other groups were collected and used as input data to query the database.In this process, different therapeutic classes of drugs or bioactive “tool”compounds were included (Table S3), and the similarity was compared using proteome data in our and other groups.

    A previous deep proteome dataset showed the proteomic response to treatment with nine compounds derived from different drug types in MCF7 and A549 cell lines in a TMT-based quantitated manner[4].All signatures of the nine drugs in both cell lines were used as input data for the “query” process, and detailed results of the top five drugs with high similarity (P <0.05, rho >0.2) in our dataset were obtained (Fig.2A) [4].Most of the perturbations caused by drug treatment on the proteome were recapitulated in our database in both the MCF7 and A549 cell lines.For example,only one Mdm2 inhibitor (nutlin-3a) was used to establish our database,and our data showed that its signature ranked at the top(rho = 0.45) after the query of nutlin-3a signatures in both cancer cell lines (Fig.2A).Moreover, notably higher similarity was observed between the signature of bortezomib and that of our three proteasome inhibitors (bortezomib, rho = 0.42; carfilzomib,rho = 0.4; and MG-132, rho = 0.4) in the database (Fig.2A).In addition, the signatures of the tubulin inhibitors vincristine(rho = 0.6) and colchicine (rho = 0.58) in our database were captured when the proteomic signature of vincristine treatment was used as input data (Fig.2A).A relatively low similarity was found for cytotoxic drugs, which may be due to their similar pharmacological effects after treatment with DNA synthesis inhibitors, DNA alkylating drugs, and topoisomerase inhibitors(Fig.2A).No drug in our database was similar to raltitrexed(a DNA synthesis inhibitor).For querying the proteomic signature of two tyrosine kinase inhibitors (dasatinib and gefitinib), we only found that the profiles of the BCR/ABL kinase inhibitor imatinib or Src/Abl inhibitor bosutinib were in the top five ranks but other tyrosine kinase inhibitors were out of the top five ranks (Fig.2A), which indicated that the signature, at the phosphorylation level but not the protein level, may be more representative of kinase inhibitors.Our results indicate that our database can be utilized for the MCF7 cell line as well as other cell lines.

    Fig.2.Proteomic alteration revealed the mechanism of drug effect.(A) The querying results of our database using signatures from published proteome data [4] based on Tandem Mass Tag(TMT)labeling method.(B)The querying results of our database using signatures from published proteome data[17-22]based on stable isotope labeling by amino acids in cell culture(SILAC) labeling and label-free method.BET:bromodomain and extra-terminal domain; CDK:cyclin dependent kinase; CHEK: checkpoint;HDAC:histone deacetylase;MAPK:mitogen-activated protein kinase;Mdm2:murine double minute 2;IMiDs:immunomodulatory drugs;NAE:NEDD8-activating enzyme;PARP:poly(adenosine diphosphateribose)polymerase;PDGFR:Platelet-derived growth factor receptor;TGF:transforming growth factor;SAHA:suberoylanilide hydroxamic acid;PI3K:phosphatidylinositol-3-kinase.

    We next explored whether our TMT-based drug signature database could also be used to explore pharmacological mechanisms using proteomic signatures produced from other quantification methods, such as Stable isotope labeling by amino acids in cell culture (SILAC) or label-free quantitative proteomic methods.Our results showed that the proteomic signatures from different types of drugs acquired using both SILAC and label-free proteomic strategies [17-22] correlated well with the proteomic profiles of similar drugs identified in our database (Fig.2B).

    These results showed the broad application potential of our dataset in exploring pharmacological mechanisms using the proteomic data produced by our group and others, although the differences were unavoidable in fields such as cancer cell types,proteomic strategies(TMT-based,SILAC-based,and label-free-based methods),compound concentrations, and duration.This result is consistent with those obtained at the transcriptomic level [1].However, we could not ignore the“false positive”or“false negative”query results,which may be due to the initial experimental designs of our and other groups.

    3.3.The characterization of commonly changed proteins under the treatment of different drugs

    Next, to explore the potential molecular mechanism of drug activity, we determined the commonly regulated proteins in each drug type(Student's t-test,P <0.05 and fold change >1.5,in at least 60%of drugs for each drug type).Using this criterion,we identified 187 commonly upregulated proteins and 160 commonly downregulated proteins.To explore the correlation between each type of drug perturbation,we constructed PPI networks for these proteins using the String database (Fig.3A), in which 78 upregulated proteins with 128 interactions and 94 downregulated proteins with 432 interactions were acquired (Fig.3A).Interestingly, more interactions were observed among commonly downregulated proteins than among commonly upregulated proteins.Interactions among commonly regulated proteins may indicate complementary possibilities of different drug perturbations.

    Fig.3.Proteomic alteration by the treatment of each drug type.(A) Protein-protein interaction (PPI) networks of upregulated/downregulated proteins from nine class of compounds.(B) The representative downregulated proteins with the CRISPR-Cas9 gene knockout scores (CERES) <-0.6 of MCF7 cell line and the mean CERES <-0.6 of 12 estrogen receptor (ER) positive breast cancer cell lines for each type of drugs.The number shows mean of fold-change for protein in each type of drugs.(C) The prognosis power of the potential key proteins.The color of circle represents the protein expression positively(red)or negatively(blue)correlated with the poor prognosis in breast cancer.Log-rank P-value was used for the significance.*P <0.05 and **P <0.01;n.s.:no significance.NA:not available.AMPK:adenosine monophosphate(AMP)-activated protein kinase;BET:bromodomain and extra-terminal domain; IMiDs: immunomodulatory drugs; PI3K: phosphatidylinositol-3-kinase.

    To identify more important proteins for the commonly downregulated proteins of each type of drug perturbations,we estimated the gene dependence of all the commonly downregulated proteins in each drug class using the DepMap database(https://depmap.org/portal/download/).The average gene essentiality scores (CRISPRCas9 gene knockout scores (CERES)), which represent gene dependence, were calculated for the 12 estrogen receptor (ER)-positive breast cancer cell lines (similar to MCF7 in genetic background; Table S4) for the downregulated proteins (Fig.3B).The results showed downregulated proteins with both CERES <-0.6 in MCF7 cell line and average CERES <-0.6(defined as essential genes[23]) in 12 ER-positive breast cancer cell lines as the potential key proteins for pharmacologic activity(Fig.3B).Interestingly,four U.S.Food and Drug Administration (FDA)-approved drug targets defined in the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/search/protein_class:FDA+approved+drug+targets) were identified as potential key proteins for pharmacological activities, including TOP2A, RRM2, DHFR, and ESR1.Furthermore,we explored the prognostic power of the above potential key proteins using The Cancer Genome Atlas breast cancer data analyzed in the HPA database (https://www.proteinatlas.org/humanproteome/pathology).Our results also indicated that for each drug type,higher expression of some key proteins was positively associated with poor prognosis (Fig.3C).The co-inhibition of these key proteins may have positive effects on breast cancer.We provided resources for the discovery of new therapies.

    3.4.Our dataset revealed pharmacological mechanism of natural products

    Fig.4.Drug mechanism querying of natural products.(A)Workflow of the drug mechanism prediction of natural product.(B)MCF7 cells treated with dimethyl sulfoxide(DMSO),colchicine, or xanthohumol (XN) for 24 h.Then the α-tubulin morphology was detected by immunofluorescence with anti-α-tubulin and 4',6-diamidino-2-phenylindole (DAPI)staining.Shown was the representative images which visualize alpha-tubulin (red) and nuclear (blue).(C) Average immunofluorescence intensity was count from 6 different regions.One-way analysis of variance was used in the analysis.Compared with DMSO, ***P <0.001 and ****P <0.0001.AMPK:adenosine monophosphate(AMP)-activated protein kinase; IMiDs: immunomodulatory drugs; BET: bromodomain and extra-terminal domain; CDK: cyclin dependent kinase; TGF: transforming growth factor; HDAC: histone deacetylase; Mdm2: murine double minute 2.

    Natural products, originating mainly from plants and microorganisms, have diverse chemical skeletons with broad biological activities [24].Nearly 40% of FDA-approved compounds were natural products in the past 30 years[24-28].Because natural products are not rationally designed based on known targets,they have been reported to have multiple targets.Our dataset then provided an opportunity for exploring new mechanism based on the “query assay” of proteomic signatures of natural products (Fig.4A).The proteomic data with statistical significance (Student's t-test,P <0.05,fold-change >1.5)from previous study[29-32]were used as signatures (Table S5) queried in our database, and drugs with similar profiles tended to appear at the top of the list with a highly positive Spearman's correlation(P <0.05,rho >0.2)(Fig.4A).

    We queried the proteomic signature of Rosemary extract from a previous study in our database [29] (Table S5).Our results showed that the three different proteasome inhibitors presented a high positive correlation (rho = 0.54 for carfilzomib, rho = 0.52 for MG132,and rho=0.45 for bortezomib),indicating that some major components in Rosemary may inhibit proteasome activity.This hypothesis was ultimately confirmed by a previous study showing that carnosol, a dietary diterpene and major polyphenol in Rosemary extract,inhibited 20S proteasome proteolytic activity[33].Similarly,signature of extract of Lippia origanoides[30](Table S5)was found to be similar to proteasome inhibitors in function (rho = 0.24 for MG132, rho = 0.22 for carfilzomib), which was confirmed by a previous study, as the polyphenol quercetin and apigenin were reported to be abundant in it and showed inhibition of the 20S proteasome[34].In addition,the cyclin dependent kinase inhibitor Ro-3306 ranked at the top of the list after the signature query of odoroside A (rho = 0.27), and its potential activity in causing cell cycle arrest was supported by previous studies[31].

    Chalcones exhibit anticancer activities, induce cell apoptosis,disrupt angiogenesis,and inhibit tubulin assembly[32].We queried the proteomic signature of xanthohumol (a prenylated chalcone found in hops)[32]in our database(Table S5).The results showed that its proteomic pattern was similar to that of the tubulin inhibitors colchicine (rho = 0.39) and vincristine (rho = 0.37).Thus,we explored whether xanthohumol inhibited tubulin assembly using an immunofluorescence assay (Fig.4B).Our results clearly showed that the fluorescence shape size of the cytoskeleton was relatively decreased in a dose-dependent manner when cells were treated with xanthohumol (Fig.S2).The results also showed decreased fluorescence intensity with an increased concentration of xanthohumol in MCF7 cells (one-way analysis of variance,P <0.05),with the tubulin inhibitor colchicine as a positive control(Fig.4C).This result indicates the predictive power of our database and suggests a new strategy for its application in natural product discovery and development.

    3.5.Proteomic signatures from clinically approved drugs identified“off-target” pharmacological mechanisms

    The polypharmacology of drugs presents multiple-target activity,leading to adverse effects on human safety.Meantime, the unprecedented efficacy of polypharmacological drugs for the treatment of multigenic diseases presents opportunities for modern drug discovery and drug repurposing[35].Therefore,effective inspection of unknown mechanisms for FDA-approved drugs could not only compromise drug safety, but also confer superior therapeutic efficacy.By comparing the proteomic perturbation profiles of clinical drugs with the signatures derived from our database, we could predict potential “off-target” pharmacological effects based on the positive Spearman's correlation result(P <0.05,rho >0.2)(Fig.5A).

    Fig.5.Off-target effect prediction of U.S.Food and Drug Administration(FDA)-approved drugs.(A)Workflow of the off-target effect prediction using published proteome profiling treated with US FDA-approved drug.(B) The in-vitro validation of the 20S proteasome inhibition.The prediction used signatures from the published proteome data [37].(C) The pharmacological activities of CC-885(25 nM)(left)and statin drugs(right)in HCT116 cell line.One-way analysis of variance was used in the analysis(n=6 for each group).(D)Cell viability of MCF7 or cereblon(CRBN)knockout(KO)MCF7 cells treated with lovastatin(LOV,25 μM)at different dose.Student's t-test was used in the analysis(n=6 for each group).(E) Molecular docking result of LOV with CRBN protein.*P <0.05, **P <0.01, ***P <0.001, and ****P <0.0001.DMSO: dimethyl sulfoxide; MLN: MLN4924; WT: wide type.

    Tamoxifen and 4-hydroxytamoxifen, which are selective estrogen receptor modulators,have been used clinically to treat patients with ER-positive breast cancer [36].We extracted the proteomic signature of 4-hydroxytamoxifen(Table S5)in a previous study[37]and used this signature for database “query” process.Surprisingly,we found that the proteomic fingerprints of 4-hydroxytamoxifen presented potentially diverse pharmacological activities, such as an AMPK activator (metformin), a proteasome inhibitor(carfilzomib), an Mdm2 inhibitor (Nutlin-3a), and an epigenetic inhibitor (valproic acid and JQ-1).A previous study reported that the tamoxifen derivative ridaifen-F(Fig.S3A)directly inhibited the human 20S proteasome [38].Therefore, we tested whether tamoxifen and 4-hydroxytamoxifen could function as potential proteasome inhibitors.Our results clearly showed that tamoxifen and 4-hydroxytamoxifen could inhibit chymotrypsin-like activity in 20S purified proteasome with Suc-LLVY-AMC labelled peptide as a substrate and carfilzomib as a positive control (Fig.5B).These results indicate the potential application of our database for exploring pharmacological effects.

    Using the same “query” strategy, we found lovastatin [39] and benzethonium chloride [40] (BZN) may have similar pharmacological mechanism with IMiDs (Table S5).IMiDs such as thalidomide and its analogs can bind to the cereblon(CRBN)subunit of the CRL4 CRBN E3 ubiquitin ligase, which confers ubiquitination and proteasomal degradation of specific substrates.This pharmacological mechanism has been clinically approved for the treatment of multiple myeloma [41,42].In addition, neddylation inhibitors such as MLN4924 can inhibit the activity of cullin-RING E3 ligases and reverse the anticancer activity of thalidomide and its derivatives [43].We then conducted a cell viability assay using lovastatin and its derivative atorvastatin in a MLN4924-dependent manner with CC-885 (a novel thalidomide derivative with IMiDs activity) as a positive control.The results showed that lovastatin inhibited the cellular survival rate, and the addition of MLN4924 prevented the cytotoxic activity of these two statin drugs, which was consistent with the pharmacological results for CC-885(Fig.5C).In addition, two other statin drugs (atorvastatin and revastatin) showed similar pharmacological activities (Fig.5C),indicating their potential “off-target” effects.The combination of benzethonium chloride and MLN4924 did not have a similar effect(Fig.S3B).Next, we compared the viability of CRBN+/+and the CRBN-/-cell lines treated with lovastatin to explore the relationship of statin drugs with CRBN subunit,and the results showed that significantly stronger cytotoxic activity of lovastatin was observed in CRBN+/+cell line in increasing drug dose(Fig.5D).These results indicated a possible dependency of CRBN on the pharmacological activity of lovastatin.This hypothesis was further supported by molecular docking results, which showed a potential interaction between lovastatin and the CRBN subunit (Fig.5E).

    Polypharmacology of small molecule was mostly related to their“off-target”effects and was a puzzling question in their clinical use or rational design of synergistic combinations [44-46].Our proteomic pattern-based study contributes to a key element in drug discovery in view of this problem.

    3.6.Rational design of synergistic drug combinations by the complementary proteomic profiles

    Fig.6.Combination therapy prediction for drug resistant exploration.(A)Workflow of the combination therapy prediction using published proteome profiling treated with drugs in resistant versus sensitive cells.(B)The combination index for topotecan+GSK126 combination at different dose in SMMC-7721 cells.The prediction used signatures from published proteome data[48].(C) The oncogenic pathway analysis for proteomic profiling of SMMC-7721 cells treated with dimethyl sulfoxide (DMSO), topotecan,and topotecan+GSK126 combination.(D)The cell viability of SMMC-7721 cells treated with different drugs for triple-combination therapy validation.One-way analysis of variance was used in the analysis.*P <0.05 and ****P <0.0001.GSK126:10 μM;topotecan:500 nM;AZD:AZD8055(50 nM).AMPK:adenosine monophosphate(AMP)-activated protein kinase;NF-κB:nuclear factor kappaB; CDK: cyclin dependent kinase; CHEK: checkpoint; HDAC: histone deacetylase; JAK: janus kinase; NA: missing value; FDR: false discovery rate.

    Developing highly selective molecules targeting a definite“driver”molecular event was as exciting as“targeted therapy”,but its major limitation was the drug resistance through rewiring of cellular signaling network[47].Our proteomic platform can also be used to elucidate mechanistic biomarker candidates and identify drug combination strategies(Fig.6A).For primary drug resistance,the proteomic profiles of sensitive and resistant cell lines can be acquired with or without treatment with the corresponding drugs.The difference between the proteomic signature under drug stimulation could be used as input data for “query” process, and compounds with negative Spearman's correlation coefficient(P <0.05,rho <-0.2) may be combined drugs utilized in the resistant cell line.Similarly, for acquired drug resistance, the proteomic profile difference between the sensitive and acquired resistant cell lines was utilized as input data, and drugs in our database were considered as potential combined drugs when the corresponding Spearman's correlation coefficients were negative (P <0.05, rho<-0.2) (Fig.6A).We then used our dataset to explore potential drug combinations.

    Our previous study reported one of the key drug resistance mechanisms of the methyltransferase EZH2 inhibitor (EZH2i) in solid tumors and showed that combined inhibition with BRD4 could increase the sensitivity of EZH2i in solid tumors[48].To validate the power of our database in exploring potential drug combination,we queried the proteomic signature of EZH2i in primary resistant and sensitive cell lines acquired previously [48] (Table S6).Our results showed that the BRD4 inhibitor, JQ1, had a negative Spearman's correlation coefficient (rho = -0.23) with EZH2i resistance.This drug combination of EZH2i and JQ1 was then validated in an EZH2i resistant cell line in our previous study, which indicated the utilization potential of our database for predicting drug combination[48].In addition, we determined that the topoisomerase inhibitor topotecan presented a notably higher negative Spearman's correlation coefficient(rho=-0.46),suggesting its synergistic effect with the EZH2 inhibitor.Cell viability was significantly different between the individual drugs and the combination of the two drugs(Fig.S4A).We then calculated the CI of EZH2i and topotecan in EZH2i resistant SMMC-7721 cell line, and the results showed a synergistic effect (combination index <1) between the EZH2 inhibitor GSK126 and topotecan at different doses (Fig.6B and Table S6).To determine the molecular mechanism of this synergistic effect, we conducted a TMT-based proteomic strategy to compare the dynamic proteome changes in topotecan treatment with a combination of the two drugs, using DMSO as a control in three biological replicates (Fig.S4B).The GSEA pathway analysis (C6:oncogenic pathway) showed several oncogenic pathways were downregulated, including VEGF, TGFB, MEK, etc., by the combined inhibition of EZH2 and topoisomerase compared with the single treatment with topotecan and DMSO (Fig.6C).It may be the potential reason for the synergistic effect of these two drugs.However,we also found that several oncogenic pathways were upregulated in the combination group, in which the mammalian target of rapamycin (mTOR) signaling kinase pathway was highly enriched(Fig.6C).This result suggests that the triple-combination inhibition of EZH2,topoisomerase, and mTOR signaling pathway may display better efficacy, which was further proven by cell viability experiments with the treatment of these three drugs (Fig.6D).

    Next,we explored potential combination strategies for acquired drug resistance by querying the proteomic differences between acquired resistant and sensitive cell lines[49-54](Table S6).Based on the negative correlation,corresponding drug combinations were estimated to potentially overcome different types of drug resistance (Table S6).

    Therefore, a system-wide proteomic approach at the protein level integrated with pharmacological analysis could raise rational design of potential drug combinations and facilitate new therapeutic options.

    4.Discussion

    Studying intracellular dynamic molecular networks perturbed by drug stimulation can provide useful information for revealing MoAs and “off-target” effects of compounds.Perturbation experiments on cancer cell lines have advanced our understanding of the biological consequences of phenotypic changes under stress [55].Recently, pharmacologic and genomic perturbation studies have been conducted on cell lines at the genomic level, in which genome-wide “cancer dependency maps” [23,56,57] and transcriptome-wide “Connectivity Map” [58] have pioneered this concept and provided a system-level view of phenotypic and cellular effects in cancers.Large-scale drug signature profiles at the transcriptional level have been successfully utilized to identify connections between different compounds, characterize novel pharmacological mechanisms, and explore new drug uses [1,58].However, the relatively low correlation coefficient between transcriptomic and proteomic data,as reported in both cell lines[59,60]and clinical tumor samples [61], reveals the limitations of using mRNA data as a reliable readout.Proteomes, being directly influenced by chemical substances, exhibit superior specificity for compound interference compared to transcriptomes.Therefore,reverse-phase protein array or MS-based protein-level readout assays have been used to quantitatively predict how cancer cells employ their adaptive protein response under the stress of drug treatment [4-7].These two proteome-wide methods are complementary due to their identification scales and study purposes,which highlight the value of the systematic proteomic characterization of behaviors or adaptations in cancers.

    In this study, we present a TMT-based proteomic platform for the global proteome perturbation profiling of 78 compounds in the MCF7 cell line, which covers almost all types of anticancer drugs.We explored the functional relevance of our dataset using a“querying”process based on the similarity of proteomic signatures acquired in other studies.We uncovered the microtubuleinhibitory activity of the chalcone xanthohumol, and the potential“off-target” effects of tamoxifen and lovastatin.In addition, the protein responses of different compounds were used for the rational design of drug combination therapies, and we found a synergistic effect between the EZH2 inhibitor GSK126 and the DNA synthesis inhibitor topotecan.These results indicate the successful utility of our protein response dataset,which facilitates therapeutic opportunities for new uses of conventional drugs.

    We also recognize the current limitations of the rational design of our project.First,we controlled concentration at 10 μM for most of the compounds we used, and the duration was set as 12 h.However, this design is somewhat unsuitable for compounds such as epigenetic inhibitors.Second,we only used the MCF7 cell line in our database.As an important cell line for drug screening,MCF7 has been used to identify transcriptomic drug responses [1].Although our dataset was suitable for query processing of proteome signatures acquired in other cancer cell lines, some false negative or positive results were unavoidable,such as the BZN signature(from the A549 cell line).This may have resulted from the difference in the origin of the tumor tissues, as we only used the MCF7 cell line for dataset establishment.Mutation information(Dep-Map Public 19Q3[62], https://depmap.org/portal/download/), such as PIK3CA E545K and GATA3 frameshift mutations, may contribute to proteomic signatures.We also found a significant difference (Fisher's exact test,P = 0.0073) in the number of non-silent mutations between the MCF7 and A549 cell lines.The power(n=2 for each drug treatment)was unsuitable for multiple testing.Third, proteomic signatures at the post-translational modifications(PTM)level are also important,in which dynamic PTM changes in phosphorylation, epigenetic modification, acetylation and ubiquitination directly reflect molecular perturbations in signaling pathways, gene expression, energy metabolism, and protein turnover events, respectively.Therefore,the integration analysis of proteomic and diverse PTM-omics with phenotypic drug response or pharmacological data provides new insights into drug combinations, especially for kinase inhibitors,ubiquitination system modulators, epigenetic inhibitors, and metabolism inhibitors.Additionally, we noted that the first hit of some signatures was not the drug type itself, such as for 8-azagunaine, and no drug in our database was similar to raltitrexed(a DNA synthesis inhibitor).This discrepancy may be due to the loss of drugs with a signature similar to that of 8-azaguanine and raltitrexed,indicating the incompleteness of our database.

    5.Conclusions

    In conclusion, our data portal raised functional connectivity by comparing proteome profiles to identify novel biological effects,which promoted the integrated analysis of various data types and facilitated the investigation of drug behaviors in cancer therapy.

    CRediT author statement

    Zhiwei Liu:Data curation,Formal analysis,Investigation,Methodology, Visualization, Writing - Original draft preparation,Reviewing and Editing;Shangwen Jiang: Investigation, Methodology;Bingbing Hao,Shuyu Xie, andYingluo Liu: Validation;Yuqi HuangandHeng Xu:Formal analysis;Cheng Luo:Writing-Original draft preparation;Min HuangandMinjia Tan:Writing-Reviewing and Editing;Jun-Yu Xu: Conceptualization, Investigation, Methodology,Writing-Original draft preparation,Reviewing and Editing.

    Declaration of competing interest

    The authors declare that there are no conflicts of interest.

    Acknowledgments

    This work was supported by the Natural Science Foundation of China(Grant Nos.:22225702 and 32322048),the National Key R&D Program of China (Grant No.: 2020YFE0202200), the Shanghai Academic/Technology Research Leader Program, China(Grant No.:22XD1420900), Guangdong High-level New R&D Institute, China(Grant No.:2019B090904008), Guangdong High-level Innovative Research Institute, China (Grant No.:2021B0909050003), the Shanghai Rising-Star Program,China(Grant No.:22QA1411100),the Youth Innovation Promotion Association of the Chinese Academy of Sciences (Grant No.: 2021276), the Young Elite Scientists Sponsorship Program by China Association for Science and Technology,China (Grant No.: 2022QNRC001), and the open fund of State Key Laboratory of Pharmaceutical Biotechnology, Nanjing University,China (Grant No.: KF-202201).We also thank the support of the Innovative Research Team of High-Level Local Universities in Shanghai, China (Grant No.: SHSMU-ZDCX20212700) and Sanofi scholarship program.

    Appendix A.Supplementary data

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

    欧美日韩精品成人综合77777| 99久国产av精品国产电影| 日韩一区二区三区影片| 男人添女人高潮全过程视频| 三级国产精品欧美在线观看| 草草在线视频免费看| 99九九线精品视频在线观看视频| 久久久久国产精品人妻一区二区| 大又大粗又爽又黄少妇毛片口| 日韩熟女老妇一区二区性免费视频| 久久99蜜桃精品久久| 亚洲国产欧美在线一区| 免费观看的影片在线观看| 久久久久久久国产电影| 国产在线免费精品| 观看av在线不卡| 偷拍熟女少妇极品色| 午夜激情福利司机影院| 国产在线男女| 亚洲真实伦在线观看| 国产成人精品婷婷| 老司机影院成人| 免费黄网站久久成人精品| 中文资源天堂在线| 亚洲熟女精品中文字幕| 久久精品久久精品一区二区三区| 久热久热在线精品观看| 久久精品国产亚洲av涩爱| a级毛片在线看网站| 国产成人精品婷婷| 国产淫片久久久久久久久| 国产在线一区二区三区精| 少妇高潮的动态图| 精品一区二区免费观看| 久久午夜福利片| 成年av动漫网址| 国产毛片在线视频| 69精品国产乱码久久久| 免费观看性生交大片5| 麻豆乱淫一区二区| 日日摸夜夜添夜夜添av毛片| 国产亚洲5aaaaa淫片| 久久人人爽人人爽人人片va| 午夜老司机福利剧场| 99久久人妻综合| 成年人免费黄色播放视频 | 精品少妇久久久久久888优播| 三级国产精品片| 精品一品国产午夜福利视频| freevideosex欧美| 人人澡人人妻人| 色视频www国产| 中文资源天堂在线| 精华霜和精华液先用哪个| 精品少妇久久久久久888优播| 国产精品久久久久久精品古装| 熟女av电影| 免费人妻精品一区二区三区视频| 国产高清有码在线观看视频| 久久精品国产a三级三级三级| 成人亚洲精品一区在线观看| 亚洲欧美成人综合另类久久久| 一二三四中文在线观看免费高清| 晚上一个人看的免费电影| a级片在线免费高清观看视频| 亚洲av中文av极速乱| 大片电影免费在线观看免费| 黄色视频在线播放观看不卡| 亚洲,欧美,日韩| 91aial.com中文字幕在线观看| 欧美另类一区| 久久久久久久大尺度免费视频| 成年人免费黄色播放视频 | 免费在线观看成人毛片| 日韩,欧美,国产一区二区三区| 日韩中字成人| 夜夜骑夜夜射夜夜干| 国产亚洲av片在线观看秒播厂| 九九在线视频观看精品| 亚洲内射少妇av| 又粗又硬又长又爽又黄的视频| 国产高清有码在线观看视频| 又黄又爽又刺激的免费视频.| 欧美xxⅹ黑人| 国产精品免费大片| 国产日韩欧美视频二区| 性高湖久久久久久久久免费观看| 日日撸夜夜添| 老司机影院成人| 在线观看一区二区三区激情| 亚洲av国产av综合av卡| 观看美女的网站| 欧美日韩亚洲高清精品| 欧美日韩一区二区视频在线观看视频在线| av又黄又爽大尺度在线免费看| 国产淫片久久久久久久久| 中文字幕久久专区| 国产精品成人在线| 久久午夜福利片| 久久久a久久爽久久v久久| 九九在线视频观看精品| 97在线视频观看| 51国产日韩欧美| 午夜视频国产福利| 美女主播在线视频| 午夜91福利影院| av在线老鸭窝| 亚洲av中文av极速乱| 国产高清国产精品国产三级| 狠狠精品人妻久久久久久综合| 欧美最新免费一区二区三区| 免费不卡的大黄色大毛片视频在线观看| 欧美日韩国产mv在线观看视频| 亚洲av成人精品一区久久| 交换朋友夫妻互换小说| 免费大片18禁| 中文在线观看免费www的网站| 午夜福利在线观看免费完整高清在| 各种免费的搞黄视频| 欧美日韩av久久| 尾随美女入室| 亚洲欧洲国产日韩| 亚洲欧美清纯卡通| 久久av网站| av福利片在线| 久久久久精品性色| 国产成人精品福利久久| 国产男人的电影天堂91| 国产精品嫩草影院av在线观看| 成人综合一区亚洲| 在线观看免费日韩欧美大片 | 中文字幕亚洲精品专区| 国产精品.久久久| 亚洲伊人久久精品综合| 久久久久网色| 国产伦理片在线播放av一区| 亚洲欧洲国产日韩| 亚洲国产色片| 久久久久久久久久成人| 国产成人91sexporn| 亚洲va在线va天堂va国产| tube8黄色片| 精品人妻偷拍中文字幕| 青青草视频在线视频观看| 特大巨黑吊av在线直播| 永久免费av网站大全| 亚洲欧洲日产国产| 欧美成人精品欧美一级黄| 少妇高潮的动态图| 色哟哟·www| 亚洲精品中文字幕在线视频 | 青青草视频在线视频观看| 女的被弄到高潮叫床怎么办| 99精国产麻豆久久婷婷| av天堂久久9| 中文字幕制服av| 大陆偷拍与自拍| 大片电影免费在线观看免费| 久久久久久久国产电影| 寂寞人妻少妇视频99o| 97在线视频观看| 欧美3d第一页| 在线天堂最新版资源| 丰满乱子伦码专区| 乱人伦中国视频| 中文天堂在线官网| 久久青草综合色| 乱码一卡2卡4卡精品| 亚洲av欧美aⅴ国产| 亚洲情色 制服丝袜| 亚洲精品乱码久久久v下载方式| 久久人妻熟女aⅴ| av一本久久久久| 美女视频免费永久观看网站| 成人毛片60女人毛片免费| 亚洲欧洲精品一区二区精品久久久 | 性色avwww在线观看| 日韩一区二区视频免费看| 这个男人来自地球电影免费观看 | 亚洲不卡免费看| 寂寞人妻少妇视频99o| 亚洲精品国产色婷婷电影| 欧美 日韩 精品 国产| 久久女婷五月综合色啪小说| 国产精品一区二区在线不卡| 成人美女网站在线观看视频| 看非洲黑人一级黄片| 国语对白做爰xxxⅹ性视频网站| 国产色婷婷99| 3wmmmm亚洲av在线观看| 在线免费观看不下载黄p国产| 精品少妇黑人巨大在线播放| 最近最新中文字幕免费大全7| 亚洲国产精品成人久久小说| 免费大片黄手机在线观看| 国产男人的电影天堂91| 纯流量卡能插随身wifi吗| 欧美三级亚洲精品| 99九九线精品视频在线观看视频| 国产精品麻豆人妻色哟哟久久| 国产极品天堂在线| 国产成人a∨麻豆精品| 精品亚洲成国产av| 国产一区二区在线观看日韩| 国产乱来视频区| 精品酒店卫生间| 国产精品国产三级国产av玫瑰| a 毛片基地| 伦理电影大哥的女人| 精品久久国产蜜桃| 日日啪夜夜爽| 亚洲av国产av综合av卡| 女性被躁到高潮视频| 国产在线视频一区二区| 亚洲精品国产色婷婷电影| 久久精品国产鲁丝片午夜精品| 成人毛片60女人毛片免费| 国产综合精华液| 国产高清三级在线| 又大又黄又爽视频免费| 天堂俺去俺来也www色官网| 99精国产麻豆久久婷婷| 啦啦啦视频在线资源免费观看| 中文在线观看免费www的网站| 男人舔奶头视频| 大片免费播放器 马上看| 午夜视频国产福利| 777米奇影视久久| 中文字幕制服av| 在线看a的网站| 日本爱情动作片www.在线观看| 国产高清三级在线| 日韩制服骚丝袜av| 韩国高清视频一区二区三区| 亚洲欧美一区二区三区黑人 | 夫妻性生交免费视频一级片| 成人18禁高潮啪啪吃奶动态图 | 尾随美女入室| 一级片'在线观看视频| 欧美 日韩 精品 国产| 成人毛片a级毛片在线播放| 中文字幕av电影在线播放| 青春草视频在线免费观看| 国产真实伦视频高清在线观看| 亚洲精品国产av成人精品| 久久毛片免费看一区二区三区| www.色视频.com| 国产精品国产三级国产专区5o| av.在线天堂| 精品一区在线观看国产| 久久国产精品男人的天堂亚洲 | freevideosex欧美| 亚洲欧美中文字幕日韩二区| 大香蕉久久网| 国产视频首页在线观看| av在线app专区| 午夜av观看不卡| 亚洲国产精品成人久久小说| 国产探花极品一区二区| 久久久国产一区二区| 久久久久久久国产电影| 免费人妻精品一区二区三区视频| 狂野欧美激情性bbbbbb| 国产免费福利视频在线观看| 亚洲色图综合在线观看| 亚洲国产精品999| 777米奇影视久久| 欧美成人午夜免费资源| 亚洲一区二区三区欧美精品| 亚洲精品日本国产第一区| 人人妻人人爽人人添夜夜欢视频 | 国产在线一区二区三区精| 22中文网久久字幕| 91精品伊人久久大香线蕉| 秋霞伦理黄片| 一级毛片 在线播放| 特大巨黑吊av在线直播| 少妇人妻一区二区三区视频| 丰满饥渴人妻一区二区三| 寂寞人妻少妇视频99o| 国产高清有码在线观看视频| 亚洲久久久国产精品| 亚洲国产欧美在线一区| 国产91av在线免费观看| 国产熟女欧美一区二区| 青青草视频在线视频观看| 国产精品麻豆人妻色哟哟久久| 26uuu在线亚洲综合色| 精品一区二区三区视频在线| 寂寞人妻少妇视频99o| 午夜免费鲁丝| 麻豆成人av视频| 五月玫瑰六月丁香| 又黄又爽又刺激的免费视频.| 又大又黄又爽视频免费| 国产欧美日韩一区二区三区在线 | 美女xxoo啪啪120秒动态图| 久久久a久久爽久久v久久| 观看av在线不卡| 久久精品国产亚洲网站| 成人毛片a级毛片在线播放| 亚洲综合精品二区| 精品少妇黑人巨大在线播放| 观看美女的网站| 人人妻人人看人人澡| 五月玫瑰六月丁香| 欧美一级a爱片免费观看看| 99热网站在线观看| 精品99又大又爽又粗少妇毛片| 亚洲精品自拍成人| 女性被躁到高潮视频| 色网站视频免费| 天美传媒精品一区二区| 精品少妇黑人巨大在线播放| 久久久久久久精品精品| 99re6热这里在线精品视频| 欧美另类一区| 人人妻人人看人人澡| 免费久久久久久久精品成人欧美视频 | 最近最新中文字幕免费大全7| 美女视频免费永久观看网站| 哪个播放器可以免费观看大片| 国产精品久久久久久精品电影小说| 久久国产精品大桥未久av | 全区人妻精品视频| 最后的刺客免费高清国语| 国产在线免费精品| 91精品伊人久久大香线蕉| 尾随美女入室| 国产精品一区二区在线不卡| 国产精品不卡视频一区二区| 色吧在线观看| 亚洲国产毛片av蜜桃av| 在线观看三级黄色| 日本欧美视频一区| 亚洲内射少妇av| 成人美女网站在线观看视频| 又粗又硬又长又爽又黄的视频| 伊人久久国产一区二区| 日韩中字成人| 中文字幕亚洲精品专区| 久久99蜜桃精品久久| 水蜜桃什么品种好| 亚洲av日韩在线播放| 嫩草影院入口| 欧美区成人在线视频| 王馨瑶露胸无遮挡在线观看| 夫妻性生交免费视频一级片| 久久 成人 亚洲| 久久久久久久久久成人| 久久精品久久精品一区二区三区| 亚洲国产精品999| 男女国产视频网站| 欧美变态另类bdsm刘玥| 国产精品秋霞免费鲁丝片| 国产精品久久久久久久久免| 国产精品偷伦视频观看了| 久久精品夜色国产| 91久久精品电影网| 在线免费观看不下载黄p国产| 久久99一区二区三区| 夜夜骑夜夜射夜夜干| 在线天堂最新版资源| 国语对白做爰xxxⅹ性视频网站| 伦理电影大哥的女人| 日日摸夜夜添夜夜爱| 中文资源天堂在线| 久久久久久久精品精品| 国产精品.久久久| 免费观看性生交大片5| 中文资源天堂在线| 久久久久久久精品精品| 国产白丝娇喘喷水9色精品| 一二三四中文在线观看免费高清| av有码第一页| 午夜老司机福利剧场| 国产av精品麻豆| videossex国产| 日韩熟女老妇一区二区性免费视频| 在线观看免费高清a一片| 99re6热这里在线精品视频| 黄色怎么调成土黄色| 伊人久久精品亚洲午夜| 一二三四中文在线观看免费高清| 人妻一区二区av| 99九九线精品视频在线观看视频| h视频一区二区三区| 一边亲一边摸免费视频| 国产一区有黄有色的免费视频| 18禁裸乳无遮挡动漫免费视频| 99久久人妻综合| 少妇被粗大猛烈的视频| 国产深夜福利视频在线观看| 国产成人免费观看mmmm| 丝袜喷水一区| 夫妻午夜视频| 尾随美女入室| 偷拍熟女少妇极品色| 在线精品无人区一区二区三| 国产av国产精品国产| 九草在线视频观看| 久久 成人 亚洲| 黄色日韩在线| 久久韩国三级中文字幕| 国产高清不卡午夜福利| 精品少妇内射三级| 最近中文字幕2019免费版| 黄色毛片三级朝国网站 | 噜噜噜噜噜久久久久久91| 两个人免费观看高清视频 | 性高湖久久久久久久久免费观看| 人妻人人澡人人爽人人| 国产免费一级a男人的天堂| 麻豆成人午夜福利视频| 亚洲真实伦在线观看| 久久精品国产亚洲网站| .国产精品久久| 最近中文字幕高清免费大全6| 另类亚洲欧美激情| 亚洲av福利一区| 精品99又大又爽又粗少妇毛片| 天堂中文最新版在线下载| 午夜免费鲁丝| 免费少妇av软件| 五月伊人婷婷丁香| 国产中年淑女户外野战色| 我的老师免费观看完整版| 91精品伊人久久大香线蕉| 不卡视频在线观看欧美| 偷拍熟女少妇极品色| 两个人免费观看高清视频 | 国产成人精品福利久久| 高清欧美精品videossex| 日韩欧美 国产精品| 亚洲欧美日韩另类电影网站| 精品久久久久久久久亚洲| 麻豆乱淫一区二区| 国产黄频视频在线观看| 国产在视频线精品| 成人毛片60女人毛片免费| 欧美日韩av久久| freevideosex欧美| 久久久久久久国产电影| 99久久精品热视频| 男女无遮挡免费网站观看| 国产高清国产精品国产三级| av有码第一页| 在现免费观看毛片| av在线观看视频网站免费| 中文字幕亚洲精品专区| 亚洲av男天堂| 国产日韩欧美在线精品| 国产男女内射视频| 99re6热这里在线精品视频| 国精品久久久久久国模美| 人体艺术视频欧美日本| 热99国产精品久久久久久7| 丝袜在线中文字幕| 日韩制服骚丝袜av| 国产成人精品一,二区| 中文精品一卡2卡3卡4更新| 秋霞伦理黄片| 国产精品女同一区二区软件| 女性生殖器流出的白浆| 国产精品熟女久久久久浪| 国产成人精品无人区| 青春草国产在线视频| av免费观看日本| 国产深夜福利视频在线观看| 亚洲内射少妇av| 伦精品一区二区三区| 免费观看性生交大片5| 亚洲成色77777| 一级毛片我不卡| 亚洲国产精品一区二区三区在线| 韩国av在线不卡| 国产午夜精品久久久久久一区二区三区| 极品教师在线视频| 久久久久视频综合| kizo精华| 啦啦啦视频在线资源免费观看| 亚洲电影在线观看av| 亚洲人与动物交配视频| 又大又黄又爽视频免费| 免费人成在线观看视频色| 国产日韩一区二区三区精品不卡 | 你懂的网址亚洲精品在线观看| 狂野欧美激情性bbbbbb| 久久精品久久久久久噜噜老黄| 国产成人精品一,二区| 啦啦啦在线观看免费高清www| 午夜福利,免费看| 99九九在线精品视频 | 久热久热在线精品观看| 看非洲黑人一级黄片| 乱码一卡2卡4卡精品| 成人特级av手机在线观看| 久久久久国产网址| 中文字幕精品免费在线观看视频 | 天天操日日干夜夜撸| 啦啦啦啦在线视频资源| 精品久久久久久电影网| 国产精品国产av在线观看| 久久 成人 亚洲| 我的老师免费观看完整版| 超碰97精品在线观看| 99久久精品热视频| 国产精品国产av在线观看| 在线观看免费高清a一片| 边亲边吃奶的免费视频| 国产黄色视频一区二区在线观看| 一级,二级,三级黄色视频| 国产黄色视频一区二区在线观看| 亚洲真实伦在线观看| 国产视频内射| 黄片无遮挡物在线观看| 午夜激情福利司机影院| 五月伊人婷婷丁香| 男人添女人高潮全过程视频| 91在线精品国自产拍蜜月| 欧美三级亚洲精品| 成人二区视频| 日本午夜av视频| 少妇人妻精品综合一区二区| 色网站视频免费| 日韩av在线免费看完整版不卡| 色哟哟·www| 亚洲三级黄色毛片| 2018国产大陆天天弄谢| 黄色欧美视频在线观看| 国产精品秋霞免费鲁丝片| 狂野欧美激情性bbbbbb| 亚洲av中文av极速乱| 国产精品久久久久久精品电影小说| 国产av精品麻豆| 国产在视频线精品| 久久国产精品大桥未久av | 丝袜脚勾引网站| 97超碰精品成人国产| 免费av中文字幕在线| 麻豆成人av视频| 国精品久久久久久国模美| 女人久久www免费人成看片| 久久国内精品自在自线图片| 永久免费av网站大全| 国产老妇伦熟女老妇高清| 欧美日韩一区二区视频在线观看视频在线| 自拍欧美九色日韩亚洲蝌蚪91 | 日本91视频免费播放| 十分钟在线观看高清视频www | 九草在线视频观看| 天美传媒精品一区二区| www.av在线官网国产| 六月丁香七月| 最近中文字幕2019免费版| 大香蕉97超碰在线| 我的老师免费观看完整版| 国产片特级美女逼逼视频| 男女免费视频国产| 丝袜脚勾引网站| 亚洲精品,欧美精品| 麻豆精品久久久久久蜜桃| 国产av码专区亚洲av| 最近手机中文字幕大全| 精品国产露脸久久av麻豆| 男人和女人高潮做爰伦理| 日韩欧美 国产精品| 久久99热6这里只有精品| 成人国产麻豆网| 亚洲国产精品国产精品| 久久久久视频综合| 欧美老熟妇乱子伦牲交| 黄色毛片三级朝国网站 | 美女cb高潮喷水在线观看| 我的老师免费观看完整版| 男人和女人高潮做爰伦理| 交换朋友夫妻互换小说| 99热6这里只有精品| 欧美日韩在线观看h| 下体分泌物呈黄色| 亚洲欧美中文字幕日韩二区| 国产免费福利视频在线观看| 在线 av 中文字幕| 日韩av在线免费看完整版不卡| 熟女人妻精品中文字幕| 免费高清在线观看视频在线观看| 欧美成人精品欧美一级黄| 一区在线观看完整版| av线在线观看网站| 日韩精品有码人妻一区| videos熟女内射| 99精国产麻豆久久婷婷| 高清欧美精品videossex| 纵有疾风起免费观看全集完整版| 51国产日韩欧美| 69精品国产乱码久久久| 成年av动漫网址| 日韩精品有码人妻一区| 99热这里只有是精品50| 丰满迷人的少妇在线观看| 久久久久久久久久成人| 热re99久久国产66热| 欧美xxⅹ黑人| 亚洲精品第二区| 一区在线观看完整版| 少妇精品久久久久久久| 偷拍熟女少妇极品色| 国产精品不卡视频一区二区| 青春草亚洲视频在线观看| 国产成人一区二区在线| 免费观看性生交大片5| 色94色欧美一区二区| 我要看黄色一级片免费的| 国产精品成人在线| 国产亚洲av片在线观看秒播厂|