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

DnaK蛋白扭鏈殘基突變體影響其ATPase活性的分子動力學模擬研究

2016-11-24 08:28:24張橋石李灌澍竇薛楷薛友林宋有濤
生物信息學 2016年3期
關鍵詞:模型

張橋石,李灌澍,竇薛楷,薛友林,宋有濤,*

(1.遼寧大學生命科學院, 沈陽 110036;2.遼寧大學環境學院, 沈陽 110036;3.遼寧大學輕型產業學院, 沈陽 110036)

?

DnaK蛋白扭鏈殘基突變體影響其ATPase活性的分子動力學模擬研究

張橋石1,李灌澍2,竇薛楷2,薛友林3,宋有濤1,2*

(1.遼寧大學生命科學院, 沈陽 110036;2.遼寧大學環境學院, 沈陽 110036;3.遼寧大學輕型產業學院, 沈陽 110036)

大腸桿菌分子伴侶蛋白DnaK氮端核苷酸結合域(NBD, nucleotide-binding domain)的II-A和II-B子域之間的一些高度保守的扭鏈殘基突變后(I202A, S203A, G223A, L227A, G228A),其ATPase活性也發生變化原因不清楚。我們通過同源建模的方法構建NBD與小分子ATP相互作用的各蛋白模型,使用分子動力學模擬方法研各模型的結構變化并嘗試找出其與ATPase活性變化的關系。結果表明,除L227A外,所有突變模型T11烴基與ATP-γ磷酸基團間的距離與活性變化間具有明顯規律;但是所有模型中,能影響與DnaJ結合,從而影響ATPase活性的β220(214-221)部分的緊致性變化符合規律,進一步的蛋白對接實驗證實了這一點,所以這些扭鏈殘基突變體可能主要是通過這兩個部分的變化,引起ATPase活性的改變。

DnaK;扭鏈殘基;同源建模;分子動力學模擬;蛋白對接

大腸桿菌DnaK蛋白是熱休克蛋白70(Hsp70, Heat shock protein 70)家族的重要成員,以分子伴侶形式發揮作用,幫助新蛋白質的正確折疊,錯誤折疊和聚集蛋白的重新折疊,調節蛋白的活性控制,細胞器和分泌蛋白的跨膜易位等[1-2],而因為DnaK在蛋白質內穩態網絡中具有的中心作用,對DnaK的結構和功能的研究是很有意義和價值的[3]。

DnaK蛋白與其他Hsp70s相同,都由一個N-端高度保守的核苷酸結合功能域和一個C-端的底物結合功能域(SBD, substrate-binding domain)組成,兩者之間可以通過變構作用進行相互調節,相互影響[4]。其中,NBD的結構類似于凝集素和己糖激酶,由4個子域組成,分別是I-A(殘基3-38;112-184), II-A(殘基185-228;310-388), I-B(殘基39-111), II-B(殘基229-309)。這些子域通過兩個交叉的α-螺旋相連接,在子域的中心形成一個包圍核苷酸與金屬粒子的特定結合口袋,對ATP進行水解,而在NBD對ATP水解時,能夠刺激增強C-端另一功能域SBD部位的底物親和能力,從而降低SBD的底物交換率,影響Hsp70發揮伴侶活性與幫助蛋白重折疊等功能,所以研究核苷酸結合功能域對ATP的水解能力,即ATPase活性,對研究Hsp70的各種功能有重要意義[5-6]。

之前,Peter等人的研究指出,DnaK蛋白中位于NBD的II-A和II-B之間的一些在域間通信中起扭鏈作用的殘基,在變構調節中發揮重要作用[7],他們選出一些重要的扭鏈殘基(見圖1),對其進行突變后,發現其ATPase活性發生了變化。其中,突變體S203A、G228A的ATPase活性增強,突變體I202A、L227A的ATPase活性減弱,突變體G223A和WT的ATPase活性基本一致。但是,Peter等人并沒有從分子水平去探討扭鏈殘基突變體ATPase活性變化的原因,對DnaK蛋白中扭鏈殘基突變與ATPase活性變化之間的關系研究仍不太清楚,尚需進一步研究。另外,McKay等人研究發現牛Hsc70(Heat shock cognate 70 protein)的T13位點(同源于DnaK蛋白T11)的自身磷酸化與ATPase活性有重要關系[8],并且T13側鏈的烴基對ATPase活性有重要作用[9]。此外,Chiappori 等人研究Hsp70與Hsp40相互作用時發現,DnaK蛋白中β220(214-221)部位的緊致性對ATPase活性也有重要影響,這部位結構越緊致,越容易與DnaJ蛋白的J結構域相結合,從而促進ATP的水解[15],并且β220與這些突變的扭鏈殘基位于同一部位上,這些扭鏈殘基突變后,可能使相連的β220部位的緊致性產生了改變,從而使其與DnaJ蛋白J結構域相互作用情況發生改變,影響DnaJ對ATP的水解刺激作用,影響ATPase活性,而我們的模擬結果也證實了這一點。

本研究不僅是對Peter等人生化試驗數據進行分子動力學上的解釋,也是對Chiappori 等人研究結果的進一步驗證,對后續其他Hsp70蛋白家族中氮端NBD,特別是扭鏈殘基的研究具有重要的借鑒與指導意義。

圖1 WT模型平衡后構象展示及各扭鏈殘基與重要部分Fig.1 The show of WT models and the hinge residues and important parts after equilibrium

1 實驗過程

1.1 蛋白質來源及同源建模

本研究采用RCSB蛋白數據庫中大腸桿菌DnaK蛋白經典核苷酸結合域模型(PDB編號1DKG: D)[10],并參考與DnaK蛋白同族的包含有小分子ATP的Hsp70與ATP結合狀態模型(PDB code:4B9Q: A)為模板[11],通過同源建模軟件Modeller9v8構建DnaK蛋白NBD-ATP 結合狀態的WT與突變體I202A、S203A、G223A、L227A、G228A三維結構模型。所有結構模型均使用Procheck對其合理性進行了評估。

1.2 動力學模擬過程

分子動力學模擬均在GROMACS 4.5.5[11]軟件包下完成, 模擬體系DnaK蛋白核苷酸結合域直接使用GROMOS96 43a1力場[12],小分子ATP使用PRODRG2服務器(http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg)添加基于GROMACS的小分子ATP的力場[13]。然后將蛋白配體模型溶于約包含26 186個SPC216水分子的立方體盒子中,蛋白邊緣與盒子邊緣的最小距離為1.0 nm。采用LINCS(linear constraint solver)算法對體系中所有鍵長進行約束。范德華力通過LJ勢(Lennard-Jones potential)方法進行估算,截斷半徑為1.0 nm。靜電相互作用采用PME(Particle mesh Ewald)算法進行估算,截斷半徑為1.0 nm。體系采用NPT系綜,溫度和壓力通過采用V-rescale 和Parrinello-Rahman 算法分別維持在300 K和1 bar,pH為7.0。體系中加入Na+或Cl-中和多余電荷,中和后的中性體系先進行2 000步的能量最小化,然后進行80 ps的限制性模擬。最后各體系在恒定的溫度和壓力下進行10 ns分子動力學模擬。分子的運動軌跡每4 ps保持一次,用于后續的數據分析。

1.3 蛋白對接

本模擬利用PatchDock網站服務器進行蛋白-蛋白對接實驗[14],DnaJ模型使用(PDB編號1XBL: A),相應結合位點根據Chiappori等人文獻中確定[15],其他設置默認,然后與平衡后各模型分別進行對接,根據打分函數對結果進行相應排序,用網站上后續工具FireDock評測其中最優的10個構象結合所需的結合能[17]。

2 結果與討論

2.1 結構穩定性

RMSD ( Root mean square deviation,均方根偏差) 是評價蛋白質穩定性的一個重要參數,通過比較RMSD曲線我們發現各模型在7ns后均達到結構保持穩定的平衡狀態(見圖2)。

圖2 NBD主鏈的RMSD*Fig.2 The RMSD values of the backbone of NBD

注:*彩圖見電子版(http://swxxx.alljournal.cn/ch/index.aspx)(2016年第3期doi:10.3969/j.issn.1672-5565.2016.03.03)。

其中除S203A外,所有突變體的RMSD波動都比較大,但活性增強突變體S203A其實在5-7ns時波動也較野生型大,平衡后構象也與野生型也有很大差別(見圖3)。這說明所有突變體與野生型相比,其構象都發生了一定變化。平衡過程中,各模型的整體RMSD值,分別如下:WT(0.24 nm),I202A(0.27 nm),S203A(0.26 nm),G223A(0.30 nm),L227A(0.32 nm),G228A(0.35 nm),并且RMSD值明顯隨扭鏈殘基所在部位不同,呈現分級趨勢,202位點與203位點,雖然ATPase活性變化相反,但兩位點卻相鄰在一起,而223、227、228這三個位點RMSD值比其他都大,它們也都在另一相連部位(見圖3)。從這可以看出,突變不同部位的扭鏈殘基,其引起的波動程度也是有差別的。而為什么活性基本不變的突變體G223A的RMSD波動也非常大,可能是其引起波動的部位與ATPase活性無關或者影響作用被抵消的原因,這需要進一步研究。

圖3 野生型WT與活性增強突變體S203A平衡后三維構 象對比Fig.3 The comparison of three-dimensional conformation between WT and the positive mutant S203A after equilibrium

2.2 平衡過程中NBD與ATP相互作用分析

ATP與NBD間的相互作用是影響DnaK蛋白ATPase活性變化的直接原因之一,主要有氫鍵、鹽橋、疏水作用,結果如下(見表1),其中G223A突變體活性雖然與野生型差不多,但與ATP相互作用是增強的,尤其是G223A的氫鍵作用是遠遠強于其他突變體的,NBD與ATP間相互作用結果并沒有什么明顯的規律。

不過進一步查看NBD與ATP間氫鍵作用存活時間發現(見圖4),影響ATPase活性的重要殘基T11與ATP之間的氫鍵作用具有明顯規律?;钚圆蛔兺蛔凅wG223A雖然整體與ATP氫鍵較多,但其T11:ATP只有一條(100%),與野生型的類似T11:ATP(98%),活性減弱突變體I202A沒有T11:ATP的氫鍵,而L227A只有一條T11:ATP(92%),但存活時間比野生型弱,這可能與其活性減弱有一定關系。并且活性增強突變體S203A、G228A都具有兩條強T11:ATP的氫鍵,S203A的兩條還強于G223A突變體的,ATPase活性規律符合,這很可能是影響扭鏈殘基突變體活性變化的重要原因。

表1 模擬過程中各模型NBD與ATP間相互作用結果

注:* Peter等人的ATPase活性實驗結果。單位μM Pi/μg DnaK/hr。

圖4 模擬過程中NBD與ATP間氫鍵存活時間Fig.4 The Hydrogen bond existence map of NBD and ATP during MDs

2.2 重要位點T11與ATP距離和相互作用分析

進一步根據McKay等人的研究,提取T11與ATP殘基的距離發現(見圖5a),活性減弱突變體I202A、L227A的T11在整個模擬過程中,與野生型和活性增強、活性不變突變體相比,T11明顯的遠離了ATP,使其更難與ATP相互作用,這可能是其與ATP間氫鍵作用減弱的主要原因。

而繼續提取T11的烴基與ATP的γ-磷酸基團間的距離(見圖5b)發現,扭鏈殘基突變體T11距離的變化與ATPase活性間的關系更加明顯,其平衡過程中的T11:ATP平均距離如下:WT(0.43 nm),活性增強突變體S203A(0.40 nm)、G228A(0.39 nm),活性減弱突變體I202A(0.59 nm)、L227A(0.45 nm),活性基本不變突變體G223A(0.42 nm)?;钚栽鰪娡蛔凅w的T11烴基與ATP的γ-磷酸距離變近,更易接觸,形成更多的氫鍵作用;活性減弱突變體間的T11烴基與ATP的γ-磷酸距離變遠,更難接觸,形成較少或沒有氫鍵作用,雖然L227A的距離增大并不太明顯,氫鍵結果也顯示了這一點,只是稍微減弱,可能還有其他影響L227A的ATPase活性變化的原因;另外,活性基本不變突變體T11烴基與ATP的γ-磷酸距離稍微減小,基本不變,氫鍵作用強弱也不改變,對ATP的水解活性也基本不變,在一定程度上能說明Peter等人的生化實驗結果。

圖5 重要位點T11與ATP距離*

Fig.5 The distance of important loci T11 and ATP

注:*彩圖見電子版(http://swxxx.alljournals.cn/ch/index.aspx)(2016年第3期doi:10.3969/j.issn.1672-5565.2016.03.03)。

2.3 扭鏈殘基周圍與ATPase活性相關重要部位β220緊致性分析

進一步查閱文獻,Chiappori等人的研究結果表明,DnaK的β220(214-221)部分緊致性能影響與DnaJ結合,從而影響ATPase活性[15],而扭鏈殘基都與β220部分相連,推測是不是這些扭鏈殘基突變后,從而影響了相連的與DnaJ結合β220部分發生改變,從而造成ATPase活性發生變化。

而蛋白質結構的緊致性可以通過其Rg(Radius of gyration,回旋半徑)來分析[16],Rg越大,表明其緊致性越低,通過提取各模型在模擬過程中β220的Rg發現(見圖6),活性減弱突變體I202A、L227A的Rg增大很多,I202A(0.51 nm)、L227A(0.51 nm),其β220部位緊致性降低,其中,I202A的Rg在整個模擬過程中都比WT高很多,L227A的Rg在平衡過程中急劇升高,最后體系平衡后,活性減弱突變體的Rg都保持在一個很相似的程度,這也表明,平衡之后的Rg值可能更準確。另外如圖5所示,平衡過程中,活性基本不變突變體G223A與WT的Rg基本保持一致G223A(0.48 nm)、WT(0.48 nm),而活性增強突變體S203A、G228A的Rg減小,S203A(0.47 nm)、G228A(0.47 nm),緊致性增強,可見,與突變點相連的β220部位的緊致性在ATPase活性變化中具有一定的規律,與最后ATPase活性變化結果符合,故扭鏈殘基突變后,可能也通過引起相連的β220部位緊致性發生改變,從而影響與DnaJ的結合,影響了ATPase的活性,對T11烴基與ATP的γ-磷酸距離變化規律形成補充,共同影響造成最后ATPase活性變化結果。

2.4 蛋白-蛋白對接與結合能分析

通過PatchDock網站進行DnaK與DnaJ對接,結果見圖7,以對接后得分最高構象為例,DnaK與DnaJ對接后,野生型與活性基本不變突變體DnaK和DnaJ間的平均距離差別不大,WT(10.6 nm)、G223A(10.7 nm),活性減弱突變體與DnaJ的平均距離增大,I202A(11.4 nm)、L227A(11.7 nm),活性增強突變體與DnaJ的平均距離減小,S203A(9.2 nm)、G228A(9.5 nm),這一結果能在一定程度上反應DnaK與DnaJ蛋白作用情況,但兩蛋白之間相互作用距離只是反映其相互作用強弱的一個方面,進一步對大范圍各模型與DnaJ對接結果進行分析,主要是能直接反應相互作用情況的平均所需結合能的分析(見表2)。

圖6 模擬過程中β220 Rg隨時間的變化*Fig.6Time dependence of the Rg in region β220 during MDs

注:*彩圖見電子版(http://swxxx.alljournals.cn/ch/index.aspx)(2016年第3期doi:10.3969/j.issn.1672-5565.2016.03.03)。

發現其中所有突變體都比野生型形成更多構象,DnaK與DnaJ接觸面積增加,雖然活性減弱突變體反而傾向形成最多的構象,I202A(638個)、L227A(622個),遠遠多于其他突變體。不過進一步的結合所需能量結果發現,活性減弱突變體雖然傾向形成更多構象,但其結合所需能量遠遠高于其他突變體,更難與DnaJ結合,影響DnaJ刺激ATP水解。而活性不變突變體G223A所需結合能顯示與野生型相差不多,活性增強突變體S203A、G228A所需結合能量比野生型減少,更有利于結合。

圖7 各DnaK模型與DnaJ蛋白對接后得分最高構象及平均距離展示Fig.7 The highest score conformation after DnaK model docking with the DnaJ and the average distance表2 蛋白對接結果統計Table 2 The statistical results of protein docking

突變體構象數量/個構象平均得分(前10個)平均接觸面積/(nm2)(前10個)平均所需結合能/(kcal.mol)(前10個)WT2548795127919.52I202A63896261331112.24S203A4049921143610.08G223A3689455135619.69L227A6229961149866.38G228A3548965131518.86

這說明突變體的緊致性變化的確影響會其與DnaJ的結合情況,從而影響ATPase活性變化。至于為什么活性減弱突變體傾向形成最多的突變體,這可能與其β220下面部分形成更多的β片層有關。另外,也不排除β220及其相連部位也是NBD向SBD域間信號傳遞的重要扭鏈區域之一,將核苷酸信號改變通過一系列扭鏈區域作用傳遞到域間linker和SBD中[7],那么反過來,β220及其相連部位扭鏈的結構改變,如緊致性、β片層程度變化,也可能會對ATPase活性產生直接影響,但具體是怎樣影響和傳遞的,本實驗還沒有找到明確規律,尚需要進一步的研究。

3 結 論

通過分子動力學模擬的實驗數據可知,活性增強突變體S203A、G228A突變引起ATPase活性區域重要殘基T11的烴基與ATP的γ-磷酸基團靠近,使ATP更容易與T11發生反應,提高ATP水解速率,從而增強ATPase活性;活性基本不變突變體G223A突變沒有使T11烴基與ATP的γ-磷酸基團距離發生改變,從而使ATP水解速率與野生型時一樣,不改變ATPase活性;活性減弱突變體I202A突變使T11的烴基與ATP的γ-磷酸基團遠離,使ATP更難與T11發生作用,從而降低ATP水解速率,降低ATPase活性,L227A雖然遠離的不明顯,但其與突變殘基相連的β220(214-221)部位的緊致性降低,所需結合能增大,難與DnaJ結合,從而影響了其ATPase活性。另外,活性增強突變體S203A、L227A突變使相連的β220的緊致性增強,更易結合,所需結合能減小?;钚曰静蛔兺蛔凅wG223A緊致性與所需結合能結果也與ATPase活性結果相符,所以,扭鏈殘基突變體很可能是通過突變引起T11部分與β220部分結構改變,從而影響了其ATPase活性的。

本研究為我們從分子水平上解釋Peter等人的生化實驗結果提供了數據支持,同時也從分子水平對McKay等人研究發現的牛Hsc70的T13位點(同源于DnaK蛋白T11)側鏈的烴基對ATPase活性有重要影響在DnaK中進行了驗證,另外,也對Chiappori 等人研究發現的DnaK蛋白中β220(214-221)部位的緊致性對ATPase活性也有重要影響的結果進行了進一步的驗證,對后續的Hsp70家族NBD部分突變及機制的研究有重要的借鑒作用。

References)

[1]MCCLELLAN A J, TAM S, KAGANOVICH D, et al. Protein quality control: chaperones culling corrupt conformations[J]. Nature Cell Biology, 2005, 7(8): 736-741.

[3]PATURY S, MIYATA Y, GESTWICKI J E. Pharmacological targeting of the Hsp70 chaperone[J]. Current Topics in Medicinal Chemistry, 2009, 9(15): 1337.

[5]MAYER M P, BUKAU B. Hsp70 chaperones: cellular functions and molecular mechanism[J].Cellular and Molecular Life Sciences, 2005, 62(6): 670-684.

[6]HARTL F U, HAYER-HARTL M. Converging concepts of protein folding in vitro and in vivo[J].Nature Structural & Molecular Biology, 2009, 16(6): 574-581.

[7]UNG P M U, THOMPSON A D, CHANG L, et al. Identification of key hinge residues important for nucleotide-dependent allostery in E. coli Hsp70/DnaK[J].PLOS Computational Biology,2013, 9(11): e1003297.

[8]CHANG L, THOMPSON A D, UNG P, et al. Mutagenesis reveals the complex relationships between ATPase rate and the chaperone activities of Escherichia coli heat shock protein 70 (Hsp70/DnaK)[J]. Journal of Biological Chemistry, 2010, 285(28): 21282-21291.

[9]SOUSA M C, MCKAY D B. The hydroxyl of threonine 13 of the bovine 70-kDa heat shock cognate protein is essential for transducing the ATP-induced conformational change[J]. Biochemistry, 1998, 37(44):15392-15399.

[10]HARRISON C J, HAYER-HARTL M, DI LIBERTO M, et al. Crystal structure of the nucleotide exchange factor GrpE bound to the ATPase domain of the molecular chaperone DnaK[J].Science, 1997, 276(5311): 431-435.

[11]KITYK R, KOPP J, SINNING I, et al. Structure and dynamics of the ATP-bound open conformation of Hsp70 chaperones[J].Molecular Cell, 2012, 48(6): 863-874.

[12]GUEX N, PEITSCH M C. SWISS-MODEL and the Swiss-Pdb Viewer: an environment for comparative protein modeling[J].Electrophoresis, 1997, 18(15): 2714-2723.

[14]SCHNEIDMAN-DUHOVNY D, INBAR Y, NUSSINOV R, et al. PatchDock and SymmDock: servers for rigid and symmetric docking[J].Nucleic Acids Research, 2005, 33(suppl 2):W363-W367.

[15]AHMAD A, BHATTACHARYA A, MCDONALD R A, et al. Heat shock protein 70 kDa chaperone/DnaJ cochaperone complex employs an unusual dynamic interface[J]. Proceedings of the National Academy of Sciences, 2011, 108(47): 18966-18971.

[17]MASHIACH E, NUSSINOV R, WOLFSON H J. FiberDock: a web server for flexible induced-fit backbone refinement in molecular docking[J].Nucleic Acids Research, 2010, 38(suppl 2): W457-W461.

Using molecular dynamics simulation to study the effects of the ATPase activity in mutants of hinge residues of Dnak

ZHANG Qiaoshi1, LI Guanshu2, DOU Xuekai2, XUE Youlin3, SONG Youtao1,2*

(1.CollegeofLifeScience,LiaoningUniversity,Shenyang110036,China;2.SchoolofEnvironmentalScience,LiaoningUniversity,Shenyang110036,China;3.CollegeofLightIndustry,LiaoningUniversity,Shenyang110036,China)

When the highly conserved hinge residues (I202, S203, G223, L227, G228), which are located in the subdomain II-A and II-B of NBD (Nucleotide binding domain) in theE.coil’sDnak, mutate into alanine. It is not clear that the change reason of ATPase activity. We build all of the wild type and mutant protein model which contain small molecule ATP by using the method of homologous modeling, than using the molecular dynamics simulation (MDs) to study the comformational change of these mutants, and want to find the relationship with the change of ATPase activity. Results show that the distances between the hydroxyl of Tyr11 residue and the γ phosphate of ATP in all models expect L227A have obvious rules with the activity of ATPase; but the change of compactness in β220 (214-221), which can impact the binding of DnaJ and effect the activity of ATPase, conform to the rules. The further experiment of protein docking confirm it, so the mutants of hinge residues may influence the ATPase activity by the changes of these two parts.

DnaK; Hinge residues; Homologous modeling; Molecular dynamics simulation; Protein docking

2015-12-16;

2016-02-19.

國家自然科學基金項目(31570154);國家自然科學基金項目(31201285)。

張橋石,男,碩士研究生。研究方向:生物信息學、熱休克蛋白;E-mail: haxqdkks@hotmail.com.

*通信作者:宋有濤,男,教授、博士生導師;研究方向:朊病毒聚集機制、環境化學;E-mail: ysong@lnu.edu.cn.

10.3969/j.issn.1672-5565.2016.03.03

Q523

A

1672-5565(2016)03-139-07

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 99精品福利视频| 亚洲看片网| 亚洲午夜18| 久久香蕉国产线看观看亚洲片| 久久一级电影| 麻豆精品在线播放| 亚洲国产日韩欧美在线| 婷婷综合在线观看丁香| 日韩在线播放中文字幕| 欧洲在线免费视频| 久久亚洲中文字幕精品一区| 久久婷婷五月综合97色| 女人爽到高潮免费视频大全| av午夜福利一片免费看| 精品国产www| 国产精品无码AV中文| 成人午夜在线播放| 欧美一级高清视频在线播放| 久久6免费视频| 亚洲视频黄| 色综合久久88色综合天天提莫| 另类重口100页在线播放| 婷婷在线网站| 国产成人三级| 欧美日韩国产在线播放| 亚洲人成网站色7799在线播放| 免费毛片网站在线观看| 亚洲福利视频一区二区| 中文字幕日韩欧美| 人妻无码中文字幕一区二区三区| 国产乱肥老妇精品视频| 91福利国产成人精品导航| 日本高清在线看免费观看| 91啪在线| 999国内精品视频免费| 少妇精品久久久一区二区三区| 国产在线观看一区精品| 亚洲中文字幕无码爆乳| 久久99国产综合精品女同| 伊人丁香五月天久久综合| 国产激情无码一区二区免费| 免费在线看黄网址| 香蕉久人久人青草青草| 在线无码九区| 国产在线精品人成导航| 免费一极毛片| 夜精品a一区二区三区| 久久久四虎成人永久免费网站| 午夜啪啪网| 四虎成人精品在永久免费| 色久综合在线| 丝袜久久剧情精品国产| 国产精品视频系列专区| 黄色免费在线网址| 亚洲人在线| 国产丝袜啪啪| 欧美激情,国产精品| 999精品视频在线| 精品少妇人妻av无码久久| 国产高颜值露脸在线观看| 丁香五月亚洲综合在线| 国产拍揄自揄精品视频网站| 99久久99这里只有免费的精品| 国产香蕉一区二区在线网站| 在线不卡免费视频| 宅男噜噜噜66国产在线观看| 精品一區二區久久久久久久網站 | 精品国产中文一级毛片在线看| 国模私拍一区二区| 少妇人妻无码首页| 免费毛片视频| 国产精品亚洲专区一区| 一级爆乳无码av| 婷婷六月综合网| 亚洲天堂.com| 农村乱人伦一区二区| 国产在线精品99一区不卡| 国产精品yjizz视频网一二区| 一级爆乳无码av| 国内精品九九久久久精品| 色播五月婷婷| 国产成人av大片在线播放|