999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

eEF2K蛋白同源模建及其抑制劑小分子的虛擬篩選研究

2019-09-10 07:22:44黎玉梅孔研于大永宋昱唐川史麗穎
中國藥房 2019年16期

黎玉梅 孔研 于大永 宋昱 唐川 史麗穎

中圖分類號 R914.2 文獻標志碼 A 文章編號 1001-0408(2019)16-2199-07

DOI 10.6039/j.issn.1001-0408.2019.16.08

摘 要 目的:篩選潛在的真核生物延伸因子2激酶(eEF2K)抑制劑小分子,為eEF2K抑制劑的設計和研發提供參考。方法:采用同源模建技術構建eEF2K蛋白晶體結構模型,并進行Loop優化和分子動力學優化,借助SAVES在線服務器從Verify_3D、EERAT和拉氏圖等3個方面對上述模型進行評估。收集55個eEF2K抑制劑小分子,使用Insight Ⅱ軟件以其中的28個(編號為奇數,設為訓練集)為基礎構建具有活性預測能力的Hypogen藥效團模型,以另外27個(編號為偶數,設為測試集)進行驗證,通過擬合活性[即半數抑制濃度的負對數(pIC50)]預測值與真實值并借助Ligand profiler熱圖篩選最優藥效團模型。結合上述藥效團模型和Lipinski五規則、分子對接方法進行eEF2K抑制劑小分子的虛擬篩選。結果與結論:所建eEF2K蛋白晶體結構模型的整體質量因素得分為93.697,其中83.33%的氨基酸Verify_3D得分≥0.2,且位于不允許區的氨基酸占氨基酸總數的1.7%,其氨基酸構象及骨架結構合理,模型可靠性高。共構建了9個具有活性預測功能的Hypogen藥效團模型(02~10號),其中03號藥效團模型包含2個氫鍵受體和2個共軛芳香環,可更好地區分活性及非活性分子,其pIC50預測值與真實值擬合最好(相關系數為0.665 3),具有較好的預測能力和較高的可靠性。通過虛擬篩選最終獲得9個潛在的eEF2K抑制劑小分子(pIC50預測值為1.074~1.185,分子與蛋白相互作用的Dcoking-score得分為-9.730~-7.467),其中Pro268、Asp267、Gln171、Phe121、Glu212可能是eEF2K抑制劑與靶點蛋白相互作用的關鍵氨基酸,作用方式包括氫鍵、鹽橋、疏水等。上述分子有望成為eEF2K抑制劑研發的先導化合物。

關鍵詞 真核延伸因子2激酶;抑制劑小分子;同源模建;Hypogen藥效團;分子對接;虛擬篩選

Study on Homology Modeling of eEF2K Protein and Virtual Screening of Its Inhibitors Molecules

LI Yumei1,KONG Yan2,YU Dayong1,SONG Yu1,TANG Chuan1,SHI Liying3(1. School of Life Sciences and Technology, Dalian University, Liaoning Dalian 116622, China; 2. Dept. of Oncology, Dalian Third People’s Hospital, Liaoning Dalian 116033, China; 3. Institute of Materia Medica, Dalian University, Liaoning Dalian 116622, China)

ABSTRACT? ?OBJECTIVE: To screen potential eEF2K inhibitor molecules, and to provide reference for the design and R&D of eEF2K inhibitor. METHODS: The eEF2K crystal structure model was constructed by homology modeling technique. The model was optimized by Loop optimization and molecular dynamics. With the help of SAVES online server, the above models were evaluated from three aspects such as Verify_3D, EERAT and Laplace diagram. Totally 55 eEF2K inhibitor molecules were collected. Hypogen pharmacophore model with activity prediction ability was constructed based on 28 of them (odd number, as training set) by Insight Ⅱ software and validated by other 27 (even number, as test set). The optimal pharmacophore model was screened by fitting the predicted and experimental values of activity [i.e. negative logarithm of half inhibitory concentration (pIC50)] and using Ligand profiler thermogram. The virtual screening of small molecules of eEF2K inhibitors was carried out by combining the above pharmacophore model, Lipinski’s five rules and molecular docking method. RESULTS & CONCLUSIONS: The overall quality factor score of the crystal structure model of eEF2K protein was 93.697. Among them, 83.33% of the amino acid Verify_3D score was more than or equal to 0.2, and 1.7% of the total amino acids were located in the non-permissible region. The amino acid conformation and skeleton structure of the model were reasonable and the reliability of the model was high. Totally 9 Hypogen pharmacophore models (No. 02-10) with active predictive function were constructed, among which No. 03 pharmacophore model included 2 hydrogen bond receptors and 2 conjugated aromatic rings, which could better distinguish active and inactive molecules. The predicted value of pIC50 fitted the experimental value best (the correlation coefficient was 0.665 3), and it had good predictive ability and high reliability. Finally, 9 potential eEF2K inhibitor molecules were obtained through virtual screening (pIC50 ranged from 1.074 to 1.185, and Dcoking-score of protein-molecule interaction ranged from -9.730 to -7.467). Pro268, Asp267, Gln171, Phe121 and Glu212 may be the key amino acids for the interaction between eEF2K inhibitors and target proteins, including hydrogen bonds, salt bridges and hydrophobicity. These 9 molecules are expected to be the lead compounds for the development of eEF2K inhibitors.

KEYWORDS? ?eEF2K; Inhibitor molecule; Homology modeling; Hypogen pharmacophore; Molecular docking; Virtual screening

真核生物延伸因子2激酶(Eukaryotic elongation factor 2 kinase,eEF2K)又稱鈣調素依賴性真核生物延伸因子2激酶(Calmodulin-dependent protein kinase Ⅱ,CaMKⅡ),是由人類基因組eEF2K基因編碼的一種激酶[1]。該酶在多種惡性腫瘤組織中均呈高表達,如惡性膠質瘤、髓母細胞瘤、乳腺癌等[2-4]。Parmer TG等[5]檢測了不同細胞周期人乳腺癌MCF-7細胞中eEF2K的活性,發現S期細胞中eEF2K活性明顯增強,提示其活性高低與腫瘤細胞周期密切相關。此外,eEF2K還可參與腫瘤細胞自噬[6]、凋亡[7]、血管新生[8]、侵襲與轉移[9]等過程的調控。由此可見,作為潛在的抗癌藥物作用靶點,eEF2K在腫瘤細胞生長、增殖過程中發揮著極其重要的作用,研究與開發新型的eEF2K抑制劑具有重要的臨床意義,現已成為當前抗腫瘤藥物研發的熱點之一。

計算機輔助藥物設計(Computer aided drug design,CADD)是以計算化學為基礎,通過計算機進行模擬、計算和預測藥物與受體之間的相互關系,設計并優化藥物先導化合物的藥物設計方法。該法可大大提高藥物研發的成功率,減少研究的盲目性,具有成本低、周期短等優點,是藥物研發的重要手段之一[10]?;诮Y構的藥物設計、基于片段的藥物設計和基于配體的藥物設計是CADD的主要應用技術,其中同源模建(Homology modeling)、分子對接(Molecular docking)和虛擬篩選(Vir- tual screening)等方法對藥物的研發意義重大。鑒于此,本研究以CADD為手段,首先通過同源模建的方法構建了eEF2K蛋白的晶體結構模型,以文獻報道的eEF2K抑制劑分子[11-12]為基礎構建了具有活性預測能力的Hypogen藥效團模型,隨后結合該藥效團模型和類藥性五規則(Lipinski 5規則)[13]、分子對接方法篩選出潛在的eEF2K抑制劑小分子,旨在為eEF2K抑制劑的優化設計和進一步研發提供理論依據。

1 資料來源

eEF2K氨基酸序列(含725個氨基酸殘基)來源于Uniprot數據庫(https://www.uniprot.org/);用于同源模建的模板蛋白來源于PDB數據庫(http://www.rcsb.org/)。本研究所需的55個eEF2K抑制劑小分子(編號1~55,結構見圖1)均來源于相關文獻[11-12],將其結構導入藥物設計工具Schrodinger Suite 2017軟件,利用其中的“Ligprep”工具以“OPLC_2005”力場進行結構優化及3D結構轉化,將結果保存為“mol2”文件格式,用于后續藥效團模型和分子對接的構建與篩選。

2 方法

2.1 同源模建

2.1.1 模型構建與優化 利用美國國立生物技術信息中心(NCBI)的“BLAST”程序(https://blast.ncbi.nlm.nih.gov/blast.cgi)進行模板搜索,采用生物分子模擬軟件Insight Ⅱ的“Align sequence to templates”模塊進行序列比對,根據比對結果進行適當的篩除修飾,以提高目標蛋白與模板蛋白的同源性,選擇比對結果最好(即同源性和相似性最高)[14]的作為模板蛋白。以Insight Ⅱ軟件中的“Modeller 9.2”程序進行同源模建,構建eEF2K蛋白的三維結構,將運算參數中的“Number of models”設置為“10”(即構建10個蛋白模型,下同),“Optimization level”選擇“High”(即高優化水平,下同),其余參數均選擇“缺省值”。以輸出模型的“PDF total energy”和“DOPE score”選擇最優模型,兩者得分越低,表明模型在同源約束條件下優化得越好、模型質量越可靠[13]。利用Insight Ⅱ軟件“Loop refinement”模塊和“Standard dynamics cascade”模塊分別進行Loop優化和分子動力學優化,以消除原子間的不合理接觸。結合前期預試驗,將Loop優化參數“Number of models”設為“10”,“Optimization level”設為“High”。分子動力學優化時,首先賦予eEF2K蛋白“CHARMM”力場,溶劑模型選擇“隱形溶劑模型(GBSW)”;隨后使用“最陡下降法(Steepest descent method)”優化,設置“收斂RMS參數(RMS gradient)”為“0.05(kcal/mol)·?”(1 kcal=4.186 kJ,1?=0.1 nm)、“最大迭代步數(Max iterations)”為“10 000”;最后使用“共軛梯度法”,設置“收斂RMS參數”為“0.001(kcal/mol)·?”、最大迭代步數為“100 000”,對模型的分子構象進行進一步優化。

2.1.2 模型評估 借助SAVES在線服務器(http://servicesn.mbi.ucla.edu/SAVES/),從Verify_3D、ERRAT和拉氏圖(Ramachandran Plot)等3個方面對eEF2K蛋白模型進行評估。其中,Verify_3D用以評估模型和氨基酸一級結構的相關性,當80%以上的氨基酸得分≥0.2即表示氨基酸側鏈構象合理[15]。ERRAT用以計算在空間距離0.35 nm內,不同原子類型間非鍵相互作用側鏈(包括C—C、C—N、C—O、N—N、N—O、O—O)的數量,并計算整體質量因素(Overall quality factor)得分,若該得分>85則表明模型可靠性高[16]。拉氏圖是通過蛋白中非鍵合原子間最小接觸距離來評價兩個相鄰肽單位骨架構象的合理性,并以此表示蛋白中允許和不允許的構象,位于不允許區的氨基酸數應小于氨基酸總數的5%[16]。

2.2 Hypogen藥效團模型構建

將圖1中編號為奇數的eEF2K抑制劑小分子作為訓練集(共28個),編號為偶數的作為測試集(共27個)。利用Insight Ⅱ軟件中“3D QSAR pharmacphore generation”模塊以訓練集分子為基礎構建具有活性預測能力的Hypogen藥效團模型,擬合活性以半數抑制濃度(IC50)的負對數(pIC50)表示。設置“Minimum interfeature distance”為“1.5”、“Conformation generation”為“Best”、“Energy threshold”為“10”,其余參數為默認值,表征各分子的構象空間,保留能量值<10 kcal/mol的分子構象[17],構建Hypogen藥效團模型。所建模型通過測試集分子進行驗證,以篩選出可靠性及預測能力較好(即預測值與真實值擬合度較好,真實值來源于文獻[11-12])的藥效團模型,模型驗證利用Insight Ⅱ軟件中的“Ligand profiler”程序進行,“Conformation generation”設為“Best”、“Energy threshold”設為“10”、“Maximum omitted features”設為“-1”以允許小分子與Hypogen藥效團特征元素部分匹配,其余參數為默認值,以Ligand profiler熱圖確定最優藥效團模型(即可區分活性分子和非活性分子的模型)用于eEF2K抑制劑的虛擬篩選。

2.3 eEF2K抑制劑虛擬篩選

2.3.1 Hypogen藥效團篩選 從ZINC數據庫(http://zinc.docking.org/)搜索用于虛擬篩選的小分子,并利用Lipinski 5規則進行初步篩選。篩選標準[13]:①化合物分子量<500 Da;②化合物結構中氫鍵給體(包括羥基、氨基等)的數量不超過5個;③化合物結構中氫鍵受體的數量不超過10個;④化合物脂水分配系數的對數值(lgP)為-5~-2;⑤化合物中可旋轉的化學鍵數不超過10個。隨后使用Insight Ⅱ軟件中“Ligand pharmacophore mapping”程序對“2.2”項下所獲Hypogen藥效團模型進行進一步篩選,并對篩選所得小分子進行活性預測。程序中“Maximum omitted features”設置為“0”、“Conformation generation”設為“Best”、“Energy threshold”設為“10”,其余參數為默認值,保留小分子所有特征與藥效團匹配其能量值均小于10 kcal/mol的構象。

2.3.2 分子對接篩選 用于分子對接目標蛋白為同源模建所構建的eEF2K蛋白,使用Schrodinger Suite 2017軟件中“Ligand docking”模塊對由Hypogen藥效團初步篩選所得小分子進一步篩選。通過該軟件的“Sitmap”程序對目標蛋白最具潛力的活性位點進行精確對接,以“Docking-score”得分函數值評價小分子與eEF2K蛋白的相互作用,該函數綜合考慮了氫鍵、疏水、范德華力等相互作用,其絕對值越大表明小分子與目標蛋白的對接復合物越穩定、匹配結合作用越好[18]。

3 結果

3.1 同源模建和模型評估

3.1.1 同源模建結果 通過BLAST程序搜索得與eEF2K蛋白同源性較高的蛋白晶體5KS5、3PDT、3LKM、3LMH,其相關參數見表1(蛋白名稱后的“A”表示該晶體有一條A鏈,下同)。其中,“最大分”和“總分”均為蛋白晶體結構評價值,兩者值越高表明晶體結構越優;“覆蓋率”為序列覆蓋率,其值越大表明晶體結構越優;“E”值為相似性分值遠離理論值可能性的概率,該值越低表明序列相似性越高;“I”值表示序列經BLAST比對的同源性,其值越高表明晶體同源性越好[19]。以上述4個蛋白作為預選模板蛋白的序列比對結果見表2。經篩除修飾后,3PDT的同源性和相似性最高,故將其作為模板蛋白進行同源模建。選取構建的10個蛋白模型中的最優模型(PDF total energy為1 316.46、DOPE score為-688.31)進行后續研究。

3.1.2 模型評估 Verify_3D評估結果顯示,eEF2K蛋白模型中,有83.33%的氨基酸得分≥0.2,表明模型中氨基酸側鏈構象合理;ERRAT評估結果顯示,本模型的整體質量因素得分為93.697,表明模型可靠性好;拉氏圖評估結果如圖2(圖中,“ ”表示位于不允許區的氨基酸)所示,本模型共分為4個區域,其中位于不允許區的氨基酸占氨基酸總數的1.7%(12/689),小于5%,表明模型骨架結構合理。

3.2 Hypogen藥效團模型構建結果

本研究共構建了9個Hypogen藥效團模型(序號:02~10),進一步通過測試集分子進行模型驗證,其Ligand profiler熱圖見圖3(圖中,熱值數值越大表明對應小分子的預測活性越高。其中,“-”表示非活性分子,“*”表示中等活性分子,“+”表示活性較高分子)。

由圖3可見,03號Hypogen藥效團模型能更好地區分活性分子和非活性分子,其預測值與真實值見表3,線性關系見圖4。由表3、圖4可見,該藥效團模型預測值與真實值相差較小,兩者相關系數(R2)為0.665 3,提示該藥效團模型具有較好的預測能力和可靠性。

03號Hypogen藥效團模型由4個部分構成,包括2個氫鍵受體(A1、A2)和2個共軛芳香環(B1、B2),詳見圖5。

3.3 eEF2K抑制劑虛擬篩選結果

本研究從ZINC數據庫中共篩選到符合Lipinski 5規則的小分子共10萬個,選取經03號Hypogen藥效團模型篩選出的pIC50預測值由大到小排序前500位的小分子進行分子對接,其中Docking-score得分排序前9位的小分子化合物(序號1~9)見表4(表中,“-”表示無相應的作用方式)。由表4可見,Pro268、Asp267、Gln171、Phe121、Glu212(同一作用方式出現頻次較高)可能為eEF2K抑制劑小分子與靶點蛋白相互作用的關鍵氨基酸。

以1號小分子為例,其與eEF2K的相互作用方式見圖6。由圖6可見,該小分子可與氨基酸Lys153、Thr266、Glu212、Ile215形成氫鍵,與氨基酸Tyr219、Tyr214、Met95、Ile215具有疏水作用,與氨基酸Tyr219具有Pi-cation作用。

4 討論

由于腫瘤細胞具有無限增殖的特點,需要消耗大量的營養物質,因此腫瘤細胞常處于缺乏營養和缺氧等條件下,此時其細胞中eEF2K的活性會增強,使得真核生物延伸因子2(Eukaryotic elongation factor 2,eEF2)的活性受到抑制,從而阻止了細胞蛋白合成過程中多肽鏈的延伸,降低了細胞對營養物質及能量的消耗,協助腫瘤細胞抵抗不良環境,有助于其增殖[20]。由此可見,eEF2K可作為腫瘤治療的潛在靶點之一,但針對該靶點的抗腫瘤藥物研發尚處于初期。相較于傳統藥物研發方法,CADD技術能夠顯著提高藥物設計速度,降低藥物研發成本,縮短研發周期[10]。其中,虛擬篩選技術可在生物活性篩選前借助計算機技術對化合物分子進行提前篩選,以降低實際篩選化合物的數量。目前,常見的藥效團篩選方法包括基本分子共同特征藥效團模型和基于受體-配體復合物的藥效團模型等。筆者認為相較于上述兩種藥效團模型,Hypogen藥效團模型可針對具有明確活性的特定化合物分子構建出具有活性預測能力的模型,對優化小分子構象具有重要的指導意義。

本研究首次針對eEF2K蛋白,利用同源模建技術構建蛋白模型,并從Verify_3D、ERRAT和拉式圖等3個方面對模型進行評估。本研究結果顯示,模型中有83.33%的氨基酸Verify_3D得分≥0.2,整體質量因素評分為93.697,位于不允許區的氨基酸數占1.7%,提示本研究所建模型的氨基酸構象及骨架結構合理,且可靠性高。同時,本研究以eEF2K抑制劑小分子(訓練集)為基礎,構建了9個(02~10號)具有活性預測能力的Hypogen藥效團模型,并以測試集小分子進行模型驗證。結果顯示,03號藥效團模型包括2個氫鍵受體(A1、A2)和2個共軛芳香環(B1、B2),其pIC50的預測值與真實值擬合較好,R 2為0.665 3,具有較好的預測能力和可靠性。本研究進一步通過03號藥效團和分子對接方法對符合Lipinski 5規則的10萬個小分子進行虛擬篩選,共得到9個抑制劑小分子,其中Pro268、Asp267、Gln171、Phe121、Glu212可能是eEF2K抑制劑與靶點蛋白相互作用的關鍵氨基酸;這9個小分子可作為eEF2K抑制劑研發的先導化合物。

綜上所述,本研究構建了可靠的eEF2K蛋白模型,該模型可用于eEF2K蛋白分子對接、虛擬篩選等分子模擬方面的研究。本課題組后續將對本研究所得的9個小分子進行結構修飾,并對其eEF2K活性抑制作用進行驗證和量效關系研究,以期獲得高活性的eEF2K抑制劑。

參考文獻

[ 1 ] XIE JL,VAN DAMME PETRA,FANG D,et al. Ablation of elongation factor 2 kinase enhances heat-shock protein 90 chaperone expression and protects cells under proteotoxic stress[J]. J Biol Chem,2019. DOI:10.1074/jbc.AC119.008036.

[ 2 ] ZHU H,YANG X,LIU J,et al. Eukaryotic elongation factor 2 kinase confers tolerance to stress conditions in cancer cells[J]. Cell Stress Chaperones,2015,20(2):217- 220.

[ 3 ] PROUD CG. Regulation and roles of elongation factor 2 kinase[J]. Biochem Soc Trans,2015,43(3):328-332.

[ 4 ] WANG X,REGUFE DA MOTA S,LIU R,et al. Eukaryo- tic elongation factor 2 kinase activity is controlled by? multiple inputs from oncogenic signaling[J]. Mol Cell Biol,2014,34(22):4088-4103.

[ 5 ] PARMER TG,WARD MD,YURKOW EJ,et al. Activity and regulation by growth factors of calmodulin-dependent protein kinase Ⅲ(elongation factor 2-kinase)in human breast cancer[J]. Br J Cancer,1999,79(1):59-64.

[ 6 ] KENNEY JW,MOORE CE,WANG X,et al. Eukaryotic elongation factor 2 kinase,an unusual enzyme with multiple roles[J]. Adv Biol Regul,2014. DOI:10.1016/j.jbior.2014.

04.003.

[ 7 ] LIU XY,ZHANG L,ZHANG Y,et al. Roles of eEF-2 kinase in cancer[J]. Chin Med J:Engl,2012,125(16):2908-2913.

[ 8 ] KHAN AA,DACE DS,RYAZANOV AG,et al. Resveratrol regulates pathologic angiogenesis by a eukaryotic elongation factor-2 kinase-regulated pathway[J]. Am J Pathol,2010,177(1):481-492.

[ 9 ] MANNING BD. Adaptation to starvation:translating a ma- tter of life or death[J]. Cancer Cell,2013,23(6):713-715.

[10] 劉維國.抗藥物依賴功效藥效團模型的構建及其虛擬篩選[D].北京:中央民族大學,2013.

[11] PAN Z,CHEN Y,LIU J,et al. Design,synthesis,and biological evaluation of polo-like kinase 1/eukaryotic elongation factor 2 kinase(PLK1/EEF2K)dual inhibitors for regulating breast cancer cells apoptosis and autophagy[J]. Eur J Med Chem,2018. DOI:10.1016/j.ejmech.2017.12.

046.

[12] GUO Y,ZHAO Y,WANG G,et al. Design,synthesis and structure-activity relationship of a focused library of β-phenylalanine derivatives as novel eEF2K inhibitors with apoptosis-inducing mechanisms in breast cancer[J]. Eur J Med Chem,2018. DOI:10.1016/j.ejmech.2017.11.

065.

[13] 張潔,譚初兵,徐為人. Lipinski五規則的研究進展[J].藥物評價研究,2011,34(6):451-455.

[14] 李微.豬繁殖與呼吸綜合癥糖基化蛋白GP2/3與CD163/9同源模建及相互作用模擬研究[D].長春:吉林大學,2016.

[15] 黃佳俊,林淑玲,歐陽萍蘭,等. DBAT酶的同源建模及與多樣性底物對接[J].化學研究與應用,2018,30(8):1246-1251.

[16] 張瑩,楊洪干,李娟,等.人類CCR5三維結構的同源模建[J].藥物生物技術,2012,19(5):432-436.

[17] 宋昱,任聰,楊彭真,等. GPR40受體苯丙酸類激動劑三維定量構效關系研究[J].化學通報,2019,82(5):446- 451.

[18] CAVASOTTO CN,ABAGYAN RA. Protein flexibility in ligand docking and virtual screening to protein kinases[J]. J Mol Biol,2004,337(1):209-225.

[19] 曹洪玉,吳艷華,周興智,等.虛擬篩選競爭性抑制NDRG3與L-Lactate結合的小分子研究[J].華南師范大學學報(自然科學版),2018,50(3):58-64.

[20] XIE J,DE SOUZA ALVES V,VON DER HAAR T,et al. Regulation of the elongation phase of protein synthesis enhances translation accuracy and modulates lifespan[J]. Curr Biol,2019,29(5):737-749.

(收稿日期:2018-12-27 修回日期:2019-06-11)

(編輯:張元媛)

主站蜘蛛池模板: 欧美一区二区啪啪| 日韩精品一区二区三区swag| 一边摸一边做爽的视频17国产| 一级爱做片免费观看久久| 欧美区一区| 国产福利免费观看| 国产精品视频白浆免费视频| 国产毛片网站| 国产午夜无码片在线观看网站| 久久亚洲中文字幕精品一区| 四虎在线高清无码| 国产呦视频免费视频在线观看| 欧美综合区自拍亚洲综合绿色| 亚洲国产91人成在线| 亚洲精品色AV无码看| 亚洲码一区二区三区| 国产精品青青| 2020极品精品国产| 乱系列中文字幕在线视频| 欧美.成人.综合在线| 日本黄色不卡视频| 99成人在线观看| 在线欧美日韩国产| 制服丝袜在线视频香蕉| 凹凸国产熟女精品视频| 久久 午夜福利 张柏芝| 女人18毛片一级毛片在线| 久久精品国产亚洲麻豆| 91破解版在线亚洲| 99久久人妻精品免费二区| 色婷婷丁香| 国产一区二区三区在线观看免费| 亚洲精品动漫在线观看| 欧美精品高清| 日本不卡在线| 亚洲精品va| 超薄丝袜足j国产在线视频| 99re热精品视频中文字幕不卡| 在线视频一区二区三区不卡| 99国产精品免费观看视频| 国产网站免费看| 91亚洲视频下载| 日本欧美一二三区色视频| 99r在线精品视频在线播放| 国模在线视频一区二区三区| 99久久国产综合精品女同| 国产XXXX做受性欧美88| 亚洲人成网7777777国产| 国产无码制服丝袜| 国产成人精品一区二区秒拍1o| 日韩成人午夜| 全部无卡免费的毛片在线看| 国产导航在线| 亚洲国产成人无码AV在线影院L| 久久无码高潮喷水| a级毛片在线免费| 无码日韩人妻精品久久蜜桃| 欧美午夜理伦三级在线观看| 97在线碰| 永久在线播放| 欧美日韩中文字幕二区三区| 欧美日本在线播放| 午夜免费视频网站| 天天色天天操综合网| 欧美日韩国产综合视频在线观看 | 91精品人妻互换| 免费观看欧美性一级| 国产精品美女在线| 久久久久国色AV免费观看性色| 伊人AV天堂| 精品免费在线视频| 欧美一区二区三区国产精品| 美女无遮挡免费网站| 亚洲 欧美 中文 AⅤ在线视频| 久久综合色天堂av| 亚洲午夜天堂| 亚洲视频免| 亚洲国产精品一区二区第一页免 | 手机永久AV在线播放| 在线国产欧美| 不卡的在线视频免费观看| 福利在线不卡一区|