• <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.

    国产精品电影一区二区三区| 国产免费现黄频在线看| 99在线视频只有这里精品首页| 亚洲avbb在线观看| 热99国产精品久久久久久7| 欧美日韩亚洲高清精品| 美女高潮喷水抽搐中文字幕| 波多野结衣一区麻豆| 9热在线视频观看99| 99久久人妻综合| 啦啦啦 在线观看视频| 一区二区三区精品91| 老司机午夜福利在线观看视频| 色播在线永久视频| 99久久人妻综合| 亚洲第一青青草原| 久久午夜亚洲精品久久| 最近最新中文字幕大全电影3 | 一级片免费观看大全| 成人国语在线视频| 女性被躁到高潮视频| 亚洲自拍偷在线| 免费高清在线观看日韩| 后天国语完整版免费观看| 久久久久精品国产欧美久久久| 曰老女人黄片| 欧美一区二区精品小视频在线| 亚洲成av片中文字幕在线观看| 69av精品久久久久久| 国产97色在线日韩免费| 国产精品乱码一区二三区的特点 | 精品国产乱码久久久久久男人| 免费在线观看亚洲国产| 老司机深夜福利视频在线观看| av天堂在线播放| videosex国产| 久久久精品国产亚洲av高清涩受| 看黄色毛片网站| 国产精品综合久久久久久久免费 | 国产av精品麻豆| 国产成人av教育| 国产精品电影一区二区三区| 欧美日本亚洲视频在线播放| 91成年电影在线观看| 在线观看免费午夜福利视频| 黑人操中国人逼视频| 好男人电影高清在线观看| 国产一区二区三区综合在线观看| 精品久久久久久电影网| 两性午夜刺激爽爽歪歪视频在线观看 | 99国产综合亚洲精品| 久久这里只有精品19| 丁香六月欧美| 18禁观看日本| 丝袜美腿诱惑在线| 国产亚洲欧美精品永久| 国产乱人伦免费视频| 国产高清激情床上av| 男女下面插进去视频免费观看| 亚洲第一欧美日韩一区二区三区| 在线观看66精品国产| 国产av一区二区精品久久| 夜夜爽天天搞| 我的亚洲天堂| 韩国精品一区二区三区| 18禁美女被吸乳视频| 中出人妻视频一区二区| 成年版毛片免费区| 欧美黑人欧美精品刺激| 69精品国产乱码久久久| 好看av亚洲va欧美ⅴa在| 人人妻,人人澡人人爽秒播| 一个人免费在线观看的高清视频| 久久久国产成人精品二区 | 免费看十八禁软件| 黄色怎么调成土黄色| 亚洲精品中文字幕在线视频| 久久久精品欧美日韩精品| 亚洲人成伊人成综合网2020| 在线十欧美十亚洲十日本专区| 国产av在哪里看| 叶爱在线成人免费视频播放| 亚洲男人的天堂狠狠| 国产精品一区二区精品视频观看| 99国产精品免费福利视频| 国产精品偷伦视频观看了| 丰满迷人的少妇在线观看| 妹子高潮喷水视频| 啦啦啦 在线观看视频| 人人妻人人澡人人看| 午夜日韩欧美国产| 日韩 欧美 亚洲 中文字幕| 精品福利永久在线观看| 成人18禁高潮啪啪吃奶动态图| 国产亚洲精品久久久久久毛片| 91麻豆精品激情在线观看国产 | 亚洲伊人色综图| 国产av精品麻豆| 精品久久蜜臀av无| 亚洲欧美激情在线| 亚洲精品国产精品久久久不卡| 国产欧美日韩综合在线一区二区| 国产免费男女视频| 精品国产一区二区三区四区第35| 18美女黄网站色大片免费观看| 久久影院123| 国产99白浆流出| 亚洲精品粉嫩美女一区| 久久久久国产一级毛片高清牌| 在线十欧美十亚洲十日本专区| 欧美黑人精品巨大| 成人黄色视频免费在线看| 成人精品一区二区免费| 午夜a级毛片| 又大又爽又粗| 丝袜美足系列| 天天躁狠狠躁夜夜躁狠狠躁| 中国美女看黄片| 50天的宝宝边吃奶边哭怎么回事| 每晚都被弄得嗷嗷叫到高潮| 高清黄色对白视频在线免费看| 久久伊人香网站| 午夜两性在线视频| 大陆偷拍与自拍| 久久久久国产一级毛片高清牌| 国产精品98久久久久久宅男小说| 欧美老熟妇乱子伦牲交| 性少妇av在线| 亚洲人成网站在线播放欧美日韩| 国产成人精品在线电影| www.精华液| 亚洲人成网站在线播放欧美日韩| 黑丝袜美女国产一区| 首页视频小说图片口味搜索| 国产视频一区二区在线看| 亚洲熟妇熟女久久| 性欧美人与动物交配| 美女国产高潮福利片在线看| 久热这里只有精品99| 啦啦啦 在线观看视频| 久99久视频精品免费| 香蕉久久夜色| 1024香蕉在线观看| 精品久久久久久电影网| 日韩精品中文字幕看吧| 深夜精品福利| 国产单亲对白刺激| 精品少妇一区二区三区视频日本电影| 日韩视频一区二区在线观看| 正在播放国产对白刺激| 日本五十路高清| 大型av网站在线播放| av天堂久久9| 亚洲成人久久性| 嫩草影院精品99| 人妻久久中文字幕网| 久久精品91蜜桃| 欧美日韩中文字幕国产精品一区二区三区 | 超色免费av| 国产成人免费无遮挡视频| 91成人精品电影| 精品福利永久在线观看| 久久久久久久久中文| 美女大奶头视频| 国产黄色免费在线视频| 咕卡用的链子| 国产成年人精品一区二区 | 亚洲片人在线观看| 国产成人av教育| 色老头精品视频在线观看| 免费在线观看完整版高清| bbb黄色大片| www日本在线高清视频| 免费在线观看完整版高清| 日本wwww免费看| 757午夜福利合集在线观看| 亚洲激情在线av| 亚洲人成电影免费在线| 变态另类成人亚洲欧美熟女 | 一进一出抽搐动态| 国产精品国产高清国产av| 日韩欧美三级三区| 一级片'在线观看视频| 丰满饥渴人妻一区二区三| 久久久国产精品麻豆| 精品日产1卡2卡| 一区二区三区激情视频| 午夜福利在线免费观看网站| www.熟女人妻精品国产| 精品久久久久久,| 久久久久九九精品影院| 亚洲中文日韩欧美视频| 宅男免费午夜| 婷婷丁香在线五月| 18禁观看日本| 亚洲五月婷婷丁香| 色婷婷久久久亚洲欧美| 女性生殖器流出的白浆| 9色porny在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 久久香蕉激情| 欧美精品亚洲一区二区| 久久人妻福利社区极品人妻图片| 黄色怎么调成土黄色| 成人三级黄色视频| 在线观看免费午夜福利视频| 中文亚洲av片在线观看爽| 亚洲狠狠婷婷综合久久图片| 亚洲一区中文字幕在线| 日韩欧美一区二区三区在线观看| 狂野欧美激情性xxxx| 动漫黄色视频在线观看| 精品福利永久在线观看| 免费一级毛片在线播放高清视频 | av有码第一页| 国产精品 欧美亚洲| 91字幕亚洲| 电影成人av| 精品国内亚洲2022精品成人| 变态另类成人亚洲欧美熟女 | 女性被躁到高潮视频| 精品国产亚洲在线| 国产三级黄色录像| 国产视频一区二区在线看| 国产视频一区二区在线看| 成年人黄色毛片网站| 亚洲黑人精品在线| 久久人妻熟女aⅴ| 亚洲aⅴ乱码一区二区在线播放 | 少妇被粗大的猛进出69影院| 国产一卡二卡三卡精品| 琪琪午夜伦伦电影理论片6080| 俄罗斯特黄特色一大片| 最近最新中文字幕大全免费视频| 亚洲黑人精品在线| 中文欧美无线码| 亚洲一区高清亚洲精品| 亚洲欧洲精品一区二区精品久久久| 人妻丰满熟妇av一区二区三区| 黄色a级毛片大全视频| 日日摸夜夜添夜夜添小说| 国产精品久久久久成人av| 一级黄色大片毛片| 亚洲 欧美 日韩 在线 免费| 伊人久久大香线蕉亚洲五| 在线观看一区二区三区| 免费高清视频大片| 亚洲狠狠婷婷综合久久图片| 视频在线观看一区二区三区| 香蕉国产在线看| 无人区码免费观看不卡| 在线播放国产精品三级| 50天的宝宝边吃奶边哭怎么回事| 欧美精品啪啪一区二区三区| 国产精品亚洲一级av第二区| 免费在线观看影片大全网站| av中文乱码字幕在线| 人人妻人人爽人人添夜夜欢视频| 五月开心婷婷网| 51午夜福利影视在线观看| 欧美成狂野欧美在线观看| 成熟少妇高潮喷水视频| 国产亚洲精品综合一区在线观看 | 一区二区三区国产精品乱码| 国产精品一区二区精品视频观看| 亚洲三区欧美一区| 美女大奶头视频| www.999成人在线观看| 91av网站免费观看| cao死你这个sao货| 欧美日韩福利视频一区二区| 久久久国产一区二区| 精品高清国产在线一区| 国产精品亚洲一级av第二区| 男女午夜视频在线观看| 亚洲国产精品999在线| 国产精品二区激情视频| 成人黄色视频免费在线看| 男女高潮啪啪啪动态图| 久热爱精品视频在线9| 精品一品国产午夜福利视频| av有码第一页| 在线观看免费午夜福利视频| 欧美中文综合在线视频| 久久精品国产综合久久久| 人人妻,人人澡人人爽秒播| 日韩三级视频一区二区三区| 午夜免费激情av| 精品日产1卡2卡| 亚洲国产欧美网| 色综合站精品国产| 高清黄色对白视频在线免费看| 高清av免费在线| 欧美不卡视频在线免费观看 | 国产成+人综合+亚洲专区| 国产视频一区二区在线看| 亚洲美女黄片视频| 欧美黑人欧美精品刺激| 男女下面插进去视频免费观看| 欧美人与性动交α欧美精品济南到| 国产精品 欧美亚洲| 1024香蕉在线观看| 久久久久久久久中文| 欧美人与性动交α欧美软件| 丝袜在线中文字幕| 欧美老熟妇乱子伦牲交| 黄片大片在线免费观看| 99久久综合精品五月天人人| 久久伊人香网站| 色老头精品视频在线观看| 亚洲伊人色综图| 老汉色av国产亚洲站长工具| 国产精品久久视频播放| 久久久国产一区二区| 咕卡用的链子| 亚洲国产欧美日韩在线播放| 亚洲精品一二三| av网站免费在线观看视频| 夜夜夜夜夜久久久久| 可以在线观看毛片的网站| 午夜两性在线视频| 国产主播在线观看一区二区| 国产成人精品久久二区二区91| 老司机午夜十八禁免费视频| 50天的宝宝边吃奶边哭怎么回事| 亚洲精品成人av观看孕妇| avwww免费| 免费女性裸体啪啪无遮挡网站| 亚洲人成电影免费在线| 久久久久久久久免费视频了| 国产精品乱码一区二三区的特点 | 欧美最黄视频在线播放免费 | 91字幕亚洲| 99热只有精品国产| 最近最新免费中文字幕在线| 日本三级黄在线观看| 搡老熟女国产l中国老女人| 村上凉子中文字幕在线| 国产精品二区激情视频| 精品福利永久在线观看| 免费久久久久久久精品成人欧美视频| 免费在线观看视频国产中文字幕亚洲| 亚洲精品中文字幕在线视频| www.999成人在线观看| 香蕉国产在线看| 久久热在线av| 一区二区三区激情视频| 久久精品aⅴ一区二区三区四区| 一个人观看的视频www高清免费观看 | 国产亚洲欧美在线一区二区| 成人手机av| 精品国产乱子伦一区二区三区| 亚洲中文字幕日韩| 国产av一区二区精品久久| 啪啪无遮挡十八禁网站| 久久伊人香网站| 又紧又爽又黄一区二区| 88av欧美| 国产av又大| 亚洲色图综合在线观看| 日本免费一区二区三区高清不卡 | 久久国产精品男人的天堂亚洲| 最近最新中文字幕大全免费视频| 老司机福利观看| 亚洲一区二区三区不卡视频| 国产一区在线观看成人免费| 久久人妻熟女aⅴ| 欧美日韩亚洲国产一区二区在线观看| 又紧又爽又黄一区二区| 亚洲精品一二三| 久久久久久久久久久久大奶| 亚洲精品国产色婷婷电影| 精品国产超薄肉色丝袜足j| 亚洲激情在线av| 久久香蕉激情| netflix在线观看网站| 成人三级黄色视频| 电影成人av| 久久久久久久久免费视频了| 国产成人欧美| videosex国产| 操美女的视频在线观看| 久久精品亚洲熟妇少妇任你| 亚洲精品久久成人aⅴ小说| 午夜福利在线观看吧| 正在播放国产对白刺激| 日韩中文字幕欧美一区二区| 在线天堂中文资源库| 欧美黄色淫秽网站| 亚洲欧美精品综合一区二区三区| 99热只有精品国产| av福利片在线| 91在线观看av| 日韩 欧美 亚洲 中文字幕| av超薄肉色丝袜交足视频| 国产欧美日韩精品亚洲av| e午夜精品久久久久久久| 亚洲一区二区三区色噜噜 | 欧美在线一区亚洲| xxxhd国产人妻xxx| 欧美在线黄色| 国产极品粉嫩免费观看在线| 亚洲五月色婷婷综合| 国产一区二区三区在线臀色熟女 | 日本五十路高清| 午夜免费激情av| 女性被躁到高潮视频| 啦啦啦 在线观看视频| 精品日产1卡2卡| 精品少妇一区二区三区视频日本电影| 女性生殖器流出的白浆| 亚洲精品国产一区二区精华液| 嫩草影视91久久| 国产成人影院久久av| 一级毛片高清免费大全| 国产精品久久久av美女十八| 69精品国产乱码久久久| 变态另类成人亚洲欧美熟女 | 日韩精品免费视频一区二区三区| av视频免费观看在线观看| 欧美在线一区亚洲| 国产亚洲欧美98| 在线国产一区二区在线| 女人高潮潮喷娇喘18禁视频| 亚洲性夜色夜夜综合| 中出人妻视频一区二区| 亚洲国产精品合色在线| av欧美777| 18禁裸乳无遮挡免费网站照片 | 欧美日韩亚洲高清精品| 深夜精品福利| 在线观看免费午夜福利视频| 高清在线国产一区| 黑人巨大精品欧美一区二区mp4| 午夜a级毛片| 久久99一区二区三区| 亚洲欧美精品综合久久99| 手机成人av网站| av有码第一页| www国产在线视频色| 成人特级黄色片久久久久久久| 亚洲狠狠婷婷综合久久图片| 国产高清国产精品国产三级| 村上凉子中文字幕在线| 桃红色精品国产亚洲av| 亚洲午夜理论影院| 国产精品免费一区二区三区在线| 欧美中文日本在线观看视频| 中文亚洲av片在线观看爽| 亚洲av第一区精品v没综合| 91国产中文字幕| 久久久久久亚洲精品国产蜜桃av| 9热在线视频观看99| 99国产精品免费福利视频| 久久久国产成人精品二区 | 天堂√8在线中文| 久久久久久亚洲精品国产蜜桃av| 国产精品1区2区在线观看.| 国产一区二区在线av高清观看| 婷婷丁香在线五月| 老汉色∧v一级毛片| 天堂中文最新版在线下载| 9色porny在线观看| 久久国产乱子伦精品免费另类| 一区二区三区精品91| 久久精品国产综合久久久| 在线永久观看黄色视频| 日韩有码中文字幕| 老熟妇乱子伦视频在线观看| 天堂√8在线中文| 大码成人一级视频| 欧美日韩黄片免| 男人的好看免费观看在线视频 | 精品人妻在线不人妻| www.精华液| 少妇被粗大的猛进出69影院| 久久久久精品国产欧美久久久| 国产高清激情床上av| 婷婷丁香在线五月| 久久狼人影院| 十分钟在线观看高清视频www| 欧美日本亚洲视频在线播放| 日韩有码中文字幕| 12—13女人毛片做爰片一| tocl精华| 丝袜在线中文字幕| 欧美成狂野欧美在线观看| 在线观看免费视频日本深夜| 日韩高清综合在线| 亚洲一区二区三区色噜噜 | av电影中文网址| 91av网站免费观看| 欧美日韩瑟瑟在线播放| 国产成人精品久久二区二区91| 日韩国内少妇激情av| 欧美激情高清一区二区三区| 成人永久免费在线观看视频| 中文字幕另类日韩欧美亚洲嫩草| 久久午夜综合久久蜜桃| 亚洲精品久久成人aⅴ小说| 多毛熟女@视频| 黄片小视频在线播放| 91老司机精品| 久久精品影院6| 搡老乐熟女国产| 日本 av在线| 久久久国产一区二区| 国产xxxxx性猛交| 日本撒尿小便嘘嘘汇集6| 成人特级黄色片久久久久久久| 狠狠狠狠99中文字幕| 美女福利国产在线| 国产男靠女视频免费网站| 黄片大片在线免费观看| 曰老女人黄片| 俄罗斯特黄特色一大片| 亚洲中文日韩欧美视频| 男女下面进入的视频免费午夜 | 国产成人啪精品午夜网站| 成人国产一区最新在线观看| 久久久久国内视频| 亚洲自偷自拍图片 自拍| 国产精品九九99| 久久国产精品人妻蜜桃| 亚洲国产看品久久| 精品一区二区三区四区五区乱码| 在线观看www视频免费| cao死你这个sao货| 91在线观看av| 99国产精品免费福利视频| 色尼玛亚洲综合影院| 涩涩av久久男人的天堂| 精品熟女少妇八av免费久了| 黄网站色视频无遮挡免费观看| av超薄肉色丝袜交足视频| 高清欧美精品videossex| 久久精品人人爽人人爽视色| 淫妇啪啪啪对白视频| 亚洲片人在线观看| 极品人妻少妇av视频| 国产欧美日韩一区二区三区在线| 一区二区三区精品91| 九色亚洲精品在线播放| 性色av乱码一区二区三区2| 看黄色毛片网站| 亚洲成av片中文字幕在线观看| 色尼玛亚洲综合影院| 啦啦啦免费观看视频1| 女警被强在线播放| 超碰97精品在线观看| 国产黄色免费在线视频| 757午夜福利合集在线观看| 午夜福利在线免费观看网站| 一级毛片精品| 中文字幕人妻丝袜一区二区| 天堂影院成人在线观看| 国产一卡二卡三卡精品| 人人妻人人添人人爽欧美一区卜| 天堂动漫精品| 欧美最黄视频在线播放免费 | 国产精品久久视频播放| 亚洲欧美日韩另类电影网站| 一级片免费观看大全| 久久精品亚洲熟妇少妇任你| 国产精品av久久久久免费| 9191精品国产免费久久| 久久久久久大精品| 黄色毛片三级朝国网站| 深夜精品福利| 久久午夜综合久久蜜桃| 日韩欧美国产一区二区入口| 国产精品 国内视频| 亚洲色图av天堂| 大陆偷拍与自拍| 搡老乐熟女国产| 久久精品aⅴ一区二区三区四区| 夫妻午夜视频| 在线观看免费视频网站a站| 人妻久久中文字幕网| 日韩精品免费视频一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲欧美日韩无卡精品| 亚洲第一欧美日韩一区二区三区| 国产在线观看jvid| 老司机福利观看| 啦啦啦免费观看视频1| 自拍欧美九色日韩亚洲蝌蚪91| 久久狼人影院| 欧美大码av| 免费高清视频大片| 中文字幕av电影在线播放| 美女高潮到喷水免费观看| 久久 成人 亚洲| 精品久久久久久久毛片微露脸| 午夜91福利影院| 中文字幕另类日韩欧美亚洲嫩草| 99国产精品一区二区三区| 欧美成人性av电影在线观看| 日韩精品中文字幕看吧| 色精品久久人妻99蜜桃| 国产精品亚洲一级av第二区| 国产精品成人在线| 亚洲国产欧美日韩在线播放| 视频区图区小说| 男女午夜视频在线观看| 精品久久蜜臀av无| 丰满的人妻完整版| 亚洲情色 制服丝袜| 视频区欧美日本亚洲| 免费在线观看视频国产中文字幕亚洲| 国产成人精品久久二区二区91| 亚洲欧美日韩另类电影网站| 精品熟女少妇八av免费久了| 欧美av亚洲av综合av国产av| 在线播放国产精品三级|