——兼論二代測序技術在雜種鑒定中的應用"/>
999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

NBA1/MERIT40在3種沙棘中的適應性分化研究
——兼論二代測序技術在雜種鑒定中的應用

2020-07-14 01:00:12丁雪洋李雪麗劉本立
植物研究 2020年5期

孫 坤 丁雪洋 張 輝 李雪麗 汪 穎 王 娟 劉本立

(西北師范大學生命科學學院,蘭州 730070)

沙棘屬(HippophaeL.)廣泛分布于歐亞大陸溫帶地區,全球共有7種11亞種,其中我國具7種7亞種,分為有皮組和無皮組[1~2]。無皮組代表中國沙棘(Hippophaerhamnoidesssp.sinensis)為較原始的種群,一般分布在海拔3 000米以下的地區[3]。肋果沙棘(H.neurocarpa)為青藏高原隆起過程中形成的適應更寒冷生態環境的有皮組一員,其生長于海拔3 400~4 300 m的地區[4~5]。棱果沙棘(H.goniocarpa)總是分布在肋果沙棘與中國沙棘同域分布的地區。形態研究及分子證據表明棱果沙棘為中國沙棘與肋果沙棘同倍化雜交形成,且是雙向雜交[6~7]。海拔導致的不同生態環境,使得生存于此的中國沙棘與肋果沙棘在表型上已有較大差異,如肋果沙棘已經進化出適應高海拔生境、較中國沙棘植株更矮小,葉片更小更厚等特征[8~9]。棱果沙棘表型大多呈中國沙棘與肋果沙棘的居間性質。但目前對三者的生態適應性研究大多數僅停留在表型層面,鮮少有人在基因層面研究不同海拔下三者適應性分化的分子機制。

BRCA1(breast cancer suppressor gene 1)是位于17號染色體上與乳腺癌和卵巢癌有關的抑癌基因,主要通過修復受損DNA及控制細胞周期檢查點等來維持基因組的穩定。BRCA1突變會提高患遺傳性乳腺癌與卵巢癌的幾率[10~11]。BRCA1的C末端包含兩個保守的BRCT結構域,該結構域結合Abraxas/CCDC98(coiled-coil domain-containing protein)[12~13]、RAP80(receptor associated protein80/ubiquitin interaction motif containing 1)[14]、BRCC36(BRCA1-BRCA2-containing complex subunit 3)、BRCC45(BRCA1-BRCA2-containing complex subunit 45)/BRE(brain and reproductive organs)[15~16]、NBA1(new component of the BRCA1-A complex)/MERIT40(mediator of RAP80 interactions and targeting 40 kDa)[17~18]等亞基組成BRCA1-A復合體(breast cancer type 1 susceptibility pro-tein);該復合體位于細胞核中,主要功能是募集BRCA1至雙鏈斷裂處、控制G2/M檢查點、去泛素化Lys-63泛素鏈信號、防止過度切除DSB(double-stranded DNA breaks)末端[19],這些功能由每個亞基之間相互協作完成。Abraxas亞基通過其磷酸化的SXXF基序直接與BRCT域結合,也招募具UIM(ubiquitin-interacting motif)域的RAP80亞基結合至BRCT域[20]。BRCC36為鋅離子金屬蛋白質酶體MPN+/JAMM家族的一個成員,能特異性地剪切K-63鏈接的多聚泛素鏈[21~22],其與細胞有絲分裂G2/M的檢查點有關[21];也在減數分裂時期的姐妹染色單體互作中扮演著重要的角色[23]。NBA1為該復合體中重要的接頭蛋白;與BRCC45亞基組成一個分子支架來維持該復合體的穩定[24~25];NBA亞基也是抵抗電離輻射所必需的,其與RAP80亞基相互作用募集BRCA1/BARD1 E3連接酶至DNA損傷位點,敲低NBA1導致細胞對化學試劑導致的DNA損傷、電離輻射更加敏感[23,26];NBA1具有介導蛋白酶體加工蛋白質的VWA結構域[27];與其同源的19S調節亞基Rpn10也具有VWA結構域,為募集修復因子至DNA損傷點所必需[28]。BRCC36、NBA1、BRE等亞基也參與形成位于細胞核及細胞質的BRISC復合體(BRCC36 isopeptidase complex)[20~21];該復合體主要在細胞核外行使功能,與有絲分裂紡錘體的組裝、泛素化信號及病毒的蛋白降解有關[21,29];BRISC復合體在細胞質中具有多個底物,如細胞質中的代謝酶SHMT2(serine hydroxy methyltransferase 2)與其結合,導致二者的酶活性受到抑制。NBA1亞基在BRISC復合體與BRCA1-A復合體中的功能稍有不同,前者主要與泛素結合有關,后者不僅參與泛素化結合,也與RAP80亞基的整合有關[29]。目前對BRCA1-A&BRISC復合體的研究,主要集中在哺乳動物上。植物中對BRCA1-A&BRISC復合體研究較少。Block等人在擬南芥中發現了該復合體亞基BRCC36的同源蛋白質AtBRCC36A(At1g80210)及AtBRCC36B(At3g06820),二者與DNA的交叉連接修復有關[25]。

極端環境下生存的植物一定有其獨特的適應機制[30]。我們推測肋果沙棘和中國沙棘中與低溫、高輻射等相關的基因已經產生了一系列適應性分化,本研究篩選出中國沙棘、肋果沙棘與棱果沙棘RNA-seq數據中的功能基因NBA1,從該基因在三者體內的表達量與功能兩方面來研究其在三者中的適應性分化。此外,隨著測序技術不斷發展,基于轉錄組或基因組序列鑒定雜交物種已得到廣泛應用[31],但目前還未見文獻評價第二代測序是否適用于鑒定雜交物種。本研究以NBA1基因為例,探討了二代測序結果及在雜交物種鑒定中的應用。

1 材料與方法

1.1 材料

1.1.1 實驗材料

采摘中國沙棘、肋果沙棘及棱果沙棘的嫩葉,隨后放入液氮中迅速冷卻。每種沙棘取3個個體重復。采樣情況具體見表1。

1.1.2 數據來源

NCBI數據庫(https://www.ncbi.nlm.nih.gov/);Phytozome數據庫(https://phytozome.jgi.doe.gov/pz/portal.html);Piceaabies基因組官網(http://congenie.org/),TAIR數據庫(https://www.arabidopsis.org/index.jsp)(見表2)。

表1 實驗材料及采集地

注:標本均保存于西北師范大學生命科學學院標本館

Note:The specimens stored in College of Life,Northwest Normal University

表2 15個物種蛋白質序列數據來源及版本

注:除中國沙棘和肋果沙棘外的其他13種植物選自植物系統進化樹上關鍵節點

Note:The other 13 species exceptH.rhamnoidesssp.sinensisandHippophaeneurocarpaare selected from the key nodes of the phylogenetic tree of plants

表3 克隆NBA1基因所選用的個體及樣本

Table 3 Experimental individuals and samples using cloningNBA1 gene

物種SpeciesPCR選用個體Individuals of PCRSanger法測序的克隆樣本Samples of cloning using Sanger sequencing中國沙棘H.rhamnoides ssp.sinensisR6;R24R6-1,3;R24-2,4肋果沙棘H.neurocarpaN20;N21N20-1,3;N21-1,5棱果沙棘H.goniocarpaGc;GdGd’-3,5;Gd-1,5

注:表中的大寫字母表示選用物種;R.中國沙棘;N.肋果沙棘;G.棱果沙棘;下標數字和小寫字母表示選用的個體編號;“Sanger法測序的克隆樣本”一欄中數字編號表示選用的克隆樣本編號;Gd.所加引物為中國沙棘NBA1基因的引物;Gd’.所加引物為肋果沙棘NBA1基因的引物

Note:The capital letters in the table indicate the selected species; R.H.rhamnoidesssp.sinensis; N.H.neurocarpa; G.H.goniocarpa; in which the subscript numbers and lowercase letters indicate selected individuals’ number; The number in the “Samples of cloning using Sanger sequencing” column represent number of selected samples of cloning; Gd indicates that the primer isNBA1 gene ofH.rhamnoidesssp.sinensis; Gd’ indicates that the primer isNBA1 gene ofH.neurocarpa

1.2 方法

1.2.1 克隆NBA1基因

采用Qiagen試劑盒分別提取3種沙棘的總RNA,隨后采用illumina平臺的第二代測序技術對3種沙棘進行轉錄組測序;轉錄組數據存于西北師范大學生物信息學平臺。每種沙棘隨機選取2個個體利用Takara反轉錄試劑盒進行反轉錄后設計引物,用于目的基因NBA1擴增(見表3)。其中,中國沙棘和肋果沙棘的NBA1基因的引物序列分別為HrF:5′-AGAGAGAGAGAGGAGAGAGAG-3′,HrR:5′-ATCTAAAAATATATTCAAAATGTCCAG-3′和HnF:5′-GAGAGAGAGAGAGAAGAAGAGA-3′,HnR:5′-TAGAATCTGATAAGTCAATTAAACCA-3′,而雜交種棱果沙棘則分別利用中國沙棘和肋果沙棘的引物擴增目的基因。

PCR產物經Takara純化試劑盒純化后,使用PMDTM19-T載體進行克隆,選取條帶清晰的克隆(見表3)進行Sanger法測序。

1.2.2 生物信息學分析

篩選正選擇基因:基于OrthoMCL找出的直系同源基因對,以PAML中codeml程序下的“選擇模型”計算dN/dS值(即w)。w>1即為受到正選擇的基因[32~33]。

查找同源序列:利用Chromas軟件分析Sanger法測序結果并導出[34];再利用NCBI網站BLASTX程序進行對比[35]。過濾條件為e值<10-6。

蛋白質結構分析:選擇同源性最高且e值最小的BLAST結果,用DNASTAR軟件中的Editseq及SeqBuilder程序截取目的基因的CDS區并翻譯為氨基酸序列[36]。利用ProtParam工具(http://www.expasy.org/tools/protparam.html)分析蛋白質的分子量、等電點、親水性等基本理化性質[37~38]。利用https://ch.embnet.org/software/TMPRED_form.html網站預測蛋白質的定位。蛋白質二級結構的預測使用Predictprotein(https://predictprotein.org/)[39~40]。將氨基酸序列提交至Zhanglab網站(https://zhanglab.ccmb.med.umich.edu/I-TASSER/)進行蛋白質三級結構預測[41]。用PyMOL軟件分析蛋白質三維結構[42]。

直系同源基因的篩選:在Linux系統下進行本地BLAST來抽提13個物種中的同源蛋白質[43]。過濾條件:e<10-6。

多序列比對及進化樹分析:利用ClustalX軟件進行多重序列比對[44]。利用Prottest軟件基于BIC算法計算14個物種中16條氨基酸序列最佳替代模型[45]。隨后在MEGA中利用最大似然法選擇最佳模型進行建樹[46]。

2 結果與分析

2.1 Sanger測序驗證雜交相關的3種沙棘中第二代測序結果的準確性

分析OrthoMCL軟件的計算結果發現,NBA1基因在親本種中國沙棘與肋果沙棘中為純合的單拷貝基因,在雜種棱果沙棘的個體中以單個雜合位點的共顯性方式存在。在中國沙棘R6號與R24號個體及肋果沙棘N20與N21號個體中均克隆出一條NBA1基因,而在二者的雜種棱果沙棘個體中同時克隆出兩條NBA1基因。如圖1所示,以Gd號棱果沙棘為例,Sanger測序顯示該個體中同時存在中國沙棘與肋果沙棘的NBA1基因;二代測序的Gb及Gc號棱果個體與中國沙棘聚為一支,Ga及Gd號棱果個體與肋果沙棘聚為一支。這表明NBA1基因在親本種中國沙棘與肋果沙棘中為單拷貝的純合基因型,在雜種棱果沙棘中以共顯性的雜合基因型存在。隨后利用Clustal軟件對Sanger法測序序列與二代測序序列進行多重序列比對,以驗證親本種中二代測序序列的準確性。結果表明,對中國沙棘而言,Sanger法測序與二代測序結果無差異。肋果沙棘N21個體利用Sanger法測出在CDS區的465位點堿基為T,二代測序為G;N16及N20個體Sanger法測序結果與二代測序結果均為G。在二代測序序列計算dN/dS值為1.112,dN/dS>1,即該基因在3個物種中受到達爾文正選擇[33];基于Sanger測序序列(GeneBank號為MK952033,MK952034,MK952035,MK952036,MK952037,MK952038,MK952039)計算dN/dS值為2.560,dN/dS>1,表明該基因在3個物種中確實受到了達爾文正選擇。

圖1 NBA1基因克隆結果及Sanger測序與二代測序序列進化樹 A.NBA1基因PCR結果;B.Sanger法測序與二代測序序列進化樹(帶下劃線的為Sanger法測序序列);C.NBA1基因菌液PCR部分結果Fig.1 Cloning of the NBA1 gene and phylogenetic tree using Sanger sequencing and next generation sequencing A.PCR results of NBA1 gene; B.Phylogenetic tree using Sanger sequencing and next sequencing techonoly(Underlined is sequences of Sanger sequencing); C.Bacteria liquid PCR results of NBA1 gene

圖2 NBA1基因在3種沙棘不同海拔下的FPKM值Fig.2 FPKM values of NBA1 gene in three species of Hippophae L. at different altitudes

2.2 NBA1基因在3種沙棘中的表達量分析

基因的表達水平可用轉錄本的豐度來體現。轉錄本豐度越高,即代表基因表達水平越高。在轉錄組數據分析中,可通過轉錄本的reads/fragments數來計算基因的表達水平[47]。本研究統計并分析了該基因在中國沙棘、肋果沙棘、棱果沙棘中的reads數量,并將測序深度與測序長度標準化,以FPKM(Fragments per Kilobase Million)值來表示基因的表達量水平。NBA1基因在中國沙棘(HrNBA1)、肋果沙棘(HnNBA1)、棱果沙棘中(HgNBA1)的FPKM均值分別為35.772、43.542、46.805。HrNBA1的FPKM均值稍低于HnNBA1與HgNBA1。以海拔為橫坐標,FPKM值為縱坐標對HrNBA1與HnNBA1作線性回歸分析,結果顯示,HrNBA1的FPKM值隨著海拔的升高而增加,R2值為0.237。HnNBA1的FPKM值隨海拔的升高而減小,R2值為0.916。

2.3 NBA1基因在親本種中國沙棘與肋果沙棘中的序列差異分析

利用NCBI網站的BLASTX程序查找NBA1基因的同源蛋白質序列,結果發現該序列與Tremaorientale的BRCA1-A&BRISC復合體序列相似性最高,E值為5e-145,序列覆蓋度為72%。利用DNASTAR軟件中的Seqbuilder與Editseq程序獲得HrNBA1及HnNBA1的CDS區。全長均為771 bp,共編碼256個氨基酸。二者在618,652及706的堿基位點的突變導致其編碼的氨基酸序列在218位點及236位點有差異。如圖3所示,將HrNBA1及HnNBA1與人類中NBA1氨基酸序列比對得知,HrNBA1及HnNBA1的VWA結構域位于135~218區域。HrNBA1在位于VWA結構域的218位點處為甲硫氨酸(M),HnNBA1為亮氨酸(L);236位點處HrNBA1為絲氨酸(S),HnNBA1為脯氨酸(P)。這兩處氨基酸位點的變化可能導致HrNBA1及HnNBA1的高級結構產生差異。

2.4 NBA1基因在親本種中國沙棘與肋果沙棘中的蛋白質三級預測結構及差異分析

ProtParam工具分析得出,HrNBA1與HnNBA1的分子式分別為:C1257H1998N352O389S14,C1260H2002N352O388S13;分子量分別為:28 714.66,28 706.66;等電點均為5.910;且為親水性蛋白。分析其在細胞中的定位得知,該蛋白質位于細胞核內,且核定位序列為209~240區域。Zhanglab網站I-TASSER程序預測HrNBA1與HnNBA1蛋白的三級結構如圖4所示,二者的蛋白質三級結構與小鼠中BRCA1-A復合體亞基NBA1的蛋白質三級結構高度相似[48]。也與TFIIH轉錄因子的Tfb4/p44亞基的蛋白質三級結構相似[49]。該蛋白質的結合位點由26,59,63,74,141,165,167,197,198,200,211等十一個位點的氨基酸組成;HrNBA1與HnNBA1在218位點處氨基酸的差異導致二者218位點后的蛋白質結構發生變化;236位點處HnNBA1的脯氨酸導致此處蛋白質的無規則卷曲發生折疊,HrNBA1較之更為舒展。

圖3 HrNBA1及HnNBA1與人類的NBA1氨基酸序列Fig.3 Amino acid sequence of NBA1 in Human and HrNBA1 and HnNBA1

圖4 NBA1的蛋白質三級預測結構 A,C.HrNBA1的三級預測結構(cartoon模式);B,D.HnNBA1的三級預測結構(mesh模式,紅色標記為該蛋白的配體結合位點)Fig.4 Predicted results of the tertiary structure of protein of NBA1 A,C.The tertiary prediction structure of HrNBA1(in which A,B is the cartoon mode); B,D.The tertiary prediction structure of HnNBA1(C,D is the mesh mode,red is labeled as the ligand binding site of the protein)

圖5 14個物種中NBA1蛋白的系統進化樹 該系統發育樹中支持率小于50的均未顯示Fig.5 Phylogenetic tree for NBA1 protein in 14 species The support rate less than 50 is not shown in this phylogenetic tree

2.5 NBA1基因在14個陸生植物中的進化關系分析

BRCA1-A在動物中為重要的抑癌基因[27~28],其亞基NBA1在復合體穩定中扮演著重要的角色[17]。NBA1基因在哺乳動物中有大量研究[19],植物中鮮少有人研究。本地BLAST結果顯示,Chlamydomonasreinhardtii中無該基因的同源蛋白質。將其他12物種中共得到14條同源蛋白質序列,與本試驗中的HrNBA1與HnNBA1做進化樹分析。利用Prottest軟件基于BIC模型計算得出,“JTT+I”為最佳氨基酸模型。利用最大似然法選擇“JTT+I”模型建樹。結果表明,該復合體基因在14個陸生植物中廣泛存在,且序列較保守。如圖4所示,該基因在這14個陸生植物中的進化關系符合APG version14的植物系統發育樹關系[50]。

3 討論

通過的RNA-seq數據計算直系同源基因對發現,NBA1基因在中國沙棘與肋果沙棘中為單拷貝基因,在棱果沙棘中以共顯性的方式存在。

3.1 NBA1基因在3種沙棘中的適應性分化

青藏高原的沙棘屬植物分布于不同的海拔,肋果沙棘較中國沙棘而言生長于更高的海拔,因此其長期暴露在更為寒冷、輻射更高的環境中。而生長于極端環境的植物已經進化出一套獨特的適應機制[30],故推測中國沙棘與肋果沙棘已進化出一系列基因以應對高海拔引起的低溫高輻射等極端環境。本研究通過比較中國沙棘NBA1基因與肋果沙棘NBA1基因的堿基序列發現,有3個堿基(618,652和706位點)發生突變,其中兩處(652與706位點)為非同義突變,這兩處堿基突變引起二者在218位點及236位點處氨基酸發生變化。Shao G等學者研究人類中的NBA1發現,該蛋白質由329個氨基酸組成,其中第95至第298位的氨基酸片段構成VWA結構域[24],該結構域可能與泛素化有關[29]。本研究中HrNBA1及HnNBA1由256個氨基酸組成,其中VWA結構域為第135位至第218位。HrNBA1及HnNBA1的218差異位點正好位于VWA結構域。將HrNBA1及HnNBA1氨基酸序列與13個位于植物系統進化樹關鍵節點的物種的NBA1同源蛋白序列進行比對,其中在Chlamydomonasreinhardtii中未發現與NBA1的同源蛋白序列,在其他12個陸生植物中均找到了NBA1的同源蛋白。對包括中國沙棘與肋果沙棘在內的14種植物中NBA1的氨基酸序列進行比對發現,該蛋白質序列在這些物種中非常保守,這與Shao G等對人類等哺乳動物中NBA1蛋白序列的保守性研究結果一致[24]。進一步分析14種植物的NBA1氨基酸序列得知,位于VWA結構域的218位點在HrNBA1中為甲硫氨酸,在其他13個物種中均為保守的亮氨酸。這表明相比于中國沙棘而言,肋果沙棘的NBA1具有更加保守的VWA結構域,也具更加完備的NBA1蛋白功能,推測生長于更高海拔的肋果沙棘通過更加保守的蛋白質序列來保持NBA1基因的功能。BRISC復合體主要在細胞核外行使功能[21],BRCA1-A復合體位于細胞核中[29],該蛋白的亞細胞定位結果顯示,HrNBA1與HnNBA1均位于細胞核內。但由于細胞核內也存在少量BRISC復合體,因此也不排除本研究中的NBA1蛋白是位于細胞核中少量BRISC復合體的NBA1亞基。目前尚不清楚核內的BRISC復合體行使何種功能,因此后續基于NBA1蛋白為BRCA1-A復合體亞基來分析本研究中的NBA1蛋白,BRCA1-A復合體中的NBA1亞基可與ABRAXAS-BRCC36二聚體及BRE結合,構成該復合體的蛋白支架[29],隨后RAP80亞基結合至該支架之上形成完整的BRCA1-A復合體,裝配完整的復合體結合至DNA雙鏈斷裂位點,從而使BRCA1遠離DNA損傷位點,通過抑制切除來達到修復斷裂的DNA雙鏈的目的[23~24,29]。本研究中計算機蛋白質三維結構建模結果進一步顯示,HrNBA1與HnNBA1的三級結構有較大差異,且二者的配體結合位點在空間排布上發生顯著變化。推測218位點的差異使中國沙棘的NBA1與肋果沙棘的NBA1在蛋白質三級結構上發生了較大變化,導致了二者結合位點發生變化,從而使HrNBA1與HnNBA1結合其他亞基出現差異,進一步導致二者在修復高輻射引起的DNA損傷功能方面有差異。除NBA1基因序列的差異之外,本研究發現該基因在3種沙棘中的表達量也有所不同,其在3種沙棘中均大量表達。Shao G等學者敲低Hela細胞NBA1基因表達量后發現細胞對電離輻射更加敏感[18,24],Bolton等研究發現上皮卵巢癌細胞中NBA1基因的過表達,保護了BRCA1功能失調的細胞中修復雙鏈DNA斷裂的活性,從而使細胞對DNA損傷更加耐受[48]。推測中國沙棘通過升高NBA1基因的表達量來應對隨著海拔升高而增強的高輻射引起的DNA損傷,而肋果沙棘則通過保持該基因序列的保守性來保證其在DNA損傷修復方面的功能,這些差異可能使得中國沙棘與肋果沙棘對不同生境的適應產生差異。

3.2 第二代測序技術在鑒定雜種方面的適用性

目前鑒定雜交物種通常采用形態鑒定、細胞學標記、分子標記、基因序列的系統分析等方法。上述方法各有優缺點,如由于形態特征的可塑性,常常模糊親代與子代的界限。細胞學標記需要明確染色體倍性與親本范圍。分析基因序列是目前鑒定雜交物種最常用的手段;單拷貝核基因為理想的分析序列[31]。本研究的Sanger測序結果表明,NBA1基因的二代測序序列準確。分析二代測序序列與Sanger測序序列進化樹得知,中國沙棘與肋果沙棘明顯聚為兩枝;在二代測序的棱果沙棘個體中,Gb,Gc個體與中國沙棘聚為一枝,Ga,Gd個體與肋果沙棘聚為一枝。以棱果沙棘Gd個體為例,進行Sanger法測序結果表明,Gd個體同時存在HrNBA1與HnNBA1。這是由于棱果沙棘為中國沙棘與肋果沙棘的雙向同倍化雜交物種[2,7],故雙親純合基因型中的等位基因以共顯性的方式存在于棱果沙棘的雜合基因型中,這與二代測序不適用于雜交物種個體測序的傳統認識是一致的。二代測序技術不適用于在個體水平鑒定雜交物種,這是因為二代測序原理與組裝方式的統計學特點,會隨機丟失測序個體體內的部分共顯性基因[51],但若選取居群內多個雜種個體進行測序,則可在一定程度上隨機反應雜種的親本來源。

綜上所述,本文研究了中國沙棘、肋果沙棘與棱果沙棘中受環境選擇的功能基因NBA1在三者中的適應性分化及與三者不同海拔適應的關系。并以NBA1基因為例,探討了二代測序的準確性以及其并不適用于在個體水平上鑒定雜種。以期為選擇二代測序鑒定雜交物種提供借鑒;也為進一步揭示不同海拔下沙棘屬植物的適應性進化提供基因數據基礎。

主站蜘蛛池模板: 国产成人你懂的在线观看| 国产偷国产偷在线高清| 在线精品亚洲一区二区古装| 色综合热无码热国产| 国产主播在线一区| AV老司机AV天堂| 丁香婷婷综合激情| 视频二区中文无码| 黄色网址手机国内免费在线观看| 久久99热66这里只有精品一 | 第一页亚洲| 无码啪啪精品天堂浪潮av| 精品超清无码视频在线观看| 久久 午夜福利 张柏芝| 亚洲综合二区| 精品一区二区三区自慰喷水| 欧美福利在线| 日韩中文精品亚洲第三区| 亚洲欧美日韩另类在线一| 中文字幕中文字字幕码一二区| 午夜国产精品视频| 日韩免费毛片视频| 日本三级精品| 国产成人精品无码一区二| 免费无码AV片在线观看国产| 四虎永久在线视频| 一级一级特黄女人精品毛片| 欧美97色| 高h视频在线| 亚洲天堂日韩av电影| 午夜视频免费试看| 欧美色图久久| 中国黄色一级视频| 免费高清自慰一区二区三区| 亚洲av日韩综合一区尤物| 亚洲伊人久久精品影院| 亚洲人成网站在线观看播放不卡| 青青青草国产| 免费欧美一级| 精品一区二区三区无码视频无码| 夜夜操国产| 国产午夜无码片在线观看网站 | 亚洲码一区二区三区| 91久久精品国产| 亚洲天堂网在线播放| 久草性视频| 成人在线天堂| 四虎AV麻豆| 国产欧美日本在线观看| 伊人欧美在线| 亚洲人成网站日本片| 欧美69视频在线| 免费观看国产小粉嫩喷水| 精品视频一区二区三区在线播| 丁香五月亚洲综合在线| 欧美不卡视频在线| 亚洲综合婷婷激情| 国产精品免费电影| 四虎影视库国产精品一区| 99ri精品视频在线观看播放| 男女性午夜福利网站| 色婷婷久久| 国产网友愉拍精品视频| 伊人久久影视| 亚洲 欧美 中文 AⅤ在线视频| 久久久久中文字幕精品视频| 毛片免费观看视频| 国产91丝袜在线播放动漫 | 国产丝袜丝视频在线观看| 欧美精品在线观看视频| 亚洲欧美另类日本| 色综合天天娱乐综合网| 国产一线在线| 久久精品一品道久久精品| 免费Aⅴ片在线观看蜜芽Tⅴ| av在线无码浏览| 亚洲成人免费看| 亚洲日本韩在线观看| 91成人精品视频| 色哟哟精品无码网站在线播放视频| 首页亚洲国产丝袜长腿综合| 色屁屁一区二区三区视频国产|