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

阿爾茨海默癥miRNA靶基因預測研究

2015-05-11 06:13:46牟曉陽
福州大學學報(自然科學版) 2015年6期

孔 薇, 張 昕, 牟曉陽

(1. 上海海事大學信息工程學院, 上海 201306; 2. 美國羅文大學生物化學系, 新澤西 08028)

阿爾茨海默癥miRNA靶基因預測研究

孔 薇1, 張 昕1, 牟曉陽2

(1. 上海海事大學信息工程學院, 上海 201306; 2. 美國羅文大學生物化學系, 新澤西 08028)

為深入了解阿爾茨海默癥的發病機制, 利用匹配的基因表達數據和miRNA表達數據, 對miRNA的靶基因進行預測分析. 首先對選擇出的差異表達miRNA進行靶基因預測; 然后利用HCTarget算法驗證預測的靶基因; 最后對miRNA及其靶基因構建調控網絡. 調控網絡中包含11個已確定的與AD有關的miRNA及6個AD致病基因. 通過分析網絡發現miRNA及其靶基因在正常與患病兩種情況下的活性變化趨勢, 生物學分析也證實了它們在AD中的重要作用, 為AD發病機制研究提供了依據.

微小RNA; 阿爾茨海默癥; 靶基因; 預測; HCTarget算法

0 引言

阿爾茨海默癥(Alzheimer’s disease, AD)是一種典型的中樞神經退行性疾病, 以進行性癡呆為主要臨床表現. AD的特征性病理改變為大腦皮層及腦區的纖維蛋白沉積, 即細胞外間隙的β淀粉樣蛋白(β-amyloid,β-A)和細胞內多聚Tau蛋白的沉積, 病理形態學上分別表現為老年斑(SP)和神經纖維纏結(NFT)[1-2]. 近年來對AD的生物信息學研究主要集中在基因表達數據的處理、 基因調控網絡的構建等方面. 隨著miRNA研究的深入, 其在AD發生發展中的作用也漸漸為人們重視. miRNA對AD發生發展的影響涉及到對APP、 BACE1、 神經元的凋亡等多環節的調控作用[3-5]. 融合基因表達數據和匹配的miRNA數據, 進行miRNA靶基因預測并構建miRNA—基因調控網絡.

miRNAs(micro RNAs)是一類能在轉錄水平或轉錄后水平調控基因表達的長度為19~27個核苷酸的內源性非編碼的小分子RNA, 廣泛分布于中樞神經系統中, 對神經發育、 分化和成熟起著重要作用[6]. 目前, 對miRNAs的研究主要集中在其mRNA靶基因的識別上, 研究人員從實驗和計算機方法兩方面對miRNAs和mRNAs靶基因進行識別. 計算機方法成為miRNAs和mRNAs靶基因識別的主要方法. 近年來, 已經有GenMiR++, Talasso等用表達數據來預測miRNA-mRNA相互作用關系的算法. 它們結合序列基本信息和表達數據來獲得更可靠的miRNA靶基因. 使用由蘇乃芳等人在GenMir++基礎上提出的HCtarget算法預測靶基因, 并使用了馬爾科夫鏈(MCMC)對GenMir++進行了改進, 其具備更好的穩健性和準確性[7].

現有的miRNA研究大部分是基于單個miRNA在疾病中的功能及受其調控的基因的研究, 對于多個與疾病相關的miRNA綜合研究及其靶基因預測方面的研究相對較少[8-10]. 使用miRNA與靶基因匹配數據, 對與AD相關的miRNAs進行靶基因預測、 驗證及調控網絡的構建. 首先使用TargetScan、 Pictar、 miRanda等在線分析工具對篩選出的miRNA進行靶基因預測, 取三者的交集, 提高預測靶基因的準確性; 其次用HCtarget算法, 利用基因表達數據和miRNA表達數據, 結合miRNA—靶基因調控關系, 對預測的靶基因進行驗證, 找出相互關系概率高的miRNA—靶基因對; 最后對篩選出的miRNA—靶基因對繪制調控網絡.

1 算法

1.1 miRNA靶基因預測算法

近年來, 研究者提出了很多miRNA的靶基因預測方法, 其中比較公認的是TargetScan[11]、 Pictar[12]和miRanda[13]三種方法. TargetScan(http://www.targetscan. org/index.html)算法是基于mRNA的3’非編碼區搜尋與miRNA的5’端第2~8位核苷酸完全互補的種子序列, 并綜合考慮RNAFold軟件計算得到的結合位點的熱力學穩定性值, 最后選取評分最高的mRNA序列. miRanda(http://www.miRNA.org/miRNA)算法原理依據主要是通過得分矩陣計算出互補程度大小, 尋找互補性最高的3’UTRs, 強調miRNA 與靶基因結合位點的保守性, 同時以miRNA 序5’端搜索靶基因. miRanda 利用Vien2naRNA計算miRNA 靶標復合體熱力學穩定性, 并淘汰不能形成雙連體的假陽性靶標. Pictar(http://pictar.mdc-berlin.de/)認為基因3’UTR序列是由miRNA綁定點及背景序列組成. PicTar利用Baum-welch算法來計算3’UTR序列是由此隱馬爾科夫模型產生的最大似然概率. 但是以上算法由于對完全種子匹配的要求過于嚴格, 訓練數據集不夠等導致了較高的假陽性率, 因此使用HCtarget算法結合序列信息和表達數據預測miRNA靶基因.

1.2 HCTarget算法原理

HCtarget用一個線性模型來描述miRNA和mRNA的關系, 并使用馬爾科夫鏈-蒙特卡羅(Markov Chain Monto Carlo, MCMC)算法計算靶標的概率.

1.2.1 模型

假設miRNA-mRNA對中, miRNA個數為M, 基因個數為N, 樣本數為T. 輸入矩陣Z代表M個miRNA在T個樣本中的表達值; 輸入C矩陣代表基因與miRNA的調控關系, 受miRNA調控則對應寫1(Cij=1), 否則寫0(Cij=0)[14]; 輸出yit代表N個基因在T個樣本中的表達值, 組成表達值矩陣Y. 這里i=1, …,N;j=1, …,M.

mRNA及miRNA表達譜間的關系用公式表示如下:

這里的yit和zjt分別表示mRNA及miRNA在T個樣本中表達值.βjt代表miRNA在T個樣本中調控強度,βot是樣本t的背景影響.

使用R這個潛在的二進制數來表示miRNA與受其調控的靶基因之間的關系. 此算法模型的目標就是求出R. 在矩陣C中, 假設R服從伯努利分布(二項分布), 也就是說, 在Cij=1和rij=0(Cij=0)的情況下,rij~bernoulli(π), 假設R:

這里的π可以作為預測序列的準確性. 此假設可以減少先前預測的假陽性率.

1.2.2 統計推斷算法MCMC

基于上面的模型, 假設觀察數據概率分布為:

為了估計參數θ= (β,σ2,π)及潛在變量R, 應用貝葉斯方法和MCMC算法. 結合適合的先驗假設,R和θ可通過如下迭代使用MCMC算法直接計算得出:

1) 基于更新的潛在變量取參數θ;

2) 基于更新的參數取變量R.

校正參數:

這里v=N-M-1, 且

對于π, 結合先驗的π~Beta(a0,b0),π的后驗分布為:

校正潛在變量:

HCTarget算法:

基于上述討論, 使用一個傳統的MCMC方法來反復評估參數和變量:

③ 由校正的參數來取樣潛在變量rij;

④ 重復上面兩個步驟直至收斂.

輸出的pij, 作為預測的miRNA調控mRNA的后驗概率.miRNA-mRNA對的p值大于某一閾值時, 此miRNA-mRNA對就認為是推斷正確的miRNA-mRNA對.

2 實驗與結果

2.1 數據來源

所使用的阿爾茨海默癥mRNA表達數據及miRNA表達數據來自美國國立生物技術信息中心(NationalCenterforBiotechnologyInformation,NCBI)網站的基因表達綜合數據庫(geneexpressionomnibus,GEO)中的數據集GSE16759. 基于AffymetrixHumanGenomeU133Plus2.0Array平臺的數據. 此數據來自AD患者及年齡相仿的正常實驗者大腦皮層的匹配的miRNA和mRNA表達數據. 包括16組AD樣本, 其中8個樣本為mRNA表達數據, 共54 675個基因, 8個樣本為miRNA表達數據, 共940個miRNA值. 兩類樣本都是按照4組正常對照樣本和4組患病樣本分類. 阿爾茨海默癥的miRNA來自人類小RNA疾病數據庫. 該數據庫提供實驗驗證的與人類疾病相關的miRNA, 及其與對應疾病的關系.

2.2 實驗數據處理

為了去除基因表達數據中大量的冗余數據和噪聲, 首先用T統計對數據進行預處理. 將4個正常樣本與4個患病樣本的基因表達數據和miRNA數據進行T統計分析, 最終選取t>1.0(p<0.048)的7 447個基因作為差異表達基因. 選取t>2.455(p<0.05)的169個miRNA作為差異表達miRNA. 將此169個差異表達的miRNA與HMDD中與AD相關的miRNA進行匹配, 得到11個miRNA, 這11個miRNA及其與AD的關系如表1所示, 表中Description為對應miRNA在HMDD中得到的與AD關系的描述. 這些描述來自實驗方法、 遺傳學實驗、 表觀遺傳學實驗及miRNA靶標交互實驗.

表1 11個miRNA及其與AD的關系

2.3 miRNA靶基因預測

首先, 將上述選取的11個miRNA通過軟件TargetScan、 Pictar和miRanda預測它們的靶基因. 將每個miRNA預測的靶基因與上述選取的差異表達基因進行匹配, 得到對應的基因表達數據. 最終得到11個miRNA共調控806個靶基因, 每個miRNA調控的靶基因數目如表2所示.

表2 11個miRNA及其調控的靶基因數目

2.4 miRNA調控網絡構建及結果

上述得到的miRNA及其靶基因對中, 有單個miRNA調控多個靶基因, 也有同一靶基因受多個miRNA調控的情況. 為了提高預測的準確性, 利用HCTarget算法檢驗得到的1 374個成對的miRNA-mRNA, 經過取不同閾值進行試驗, 在提高預測靶基因的準確性及結果可靠性的前提下, 設定閾值為0.8, 最后得到504對miRNA-mRNA, 其中mRNA共292個. 為了更直觀地觀測, 用Cytoscape軟件繪制出它們之間的調控網絡圖, 如圖1所示. 其中靶基因的表達值由源數據中基因在各樣本中的表達值取均值所得, miRNA表達活性由源數據中miRNA在各樣本中表達活性取均值所得. 圖中圓中心的11個菱形節點代表miRNA, 外層兩圈圓形節點代表靶基因; 節點為紅色代表表達水平上升, 綠色代表表達水平下降, 顏色越深表示表達水平越高.

圖1 11個miRNA及其靶基因在正常、 患病樣本下調控圖

2.5 調控網絡分析

從調控圖可以直觀地看出, miR-20a、 miR-29a、 miR-195、 miR-21、 miR-107和miR-137在正常樣本和患病樣本中表達活性都較高. 與正常樣本相比較, 除了miR-107以外的10個miRNA在患病樣本中的表達活性都上升. miR-107在早期AD發病進展中在顳葉皮層灰質區表達活性降低, 鄰近大腦組織中神經元血小板的計數和神經元纖維纏結數的增加與miR-107的表達水平降低有關[15]. 計算分析表明, 3’UTR區域的mRNA-BACE1受miR-107調控. BACE1是一種叫做淀粉樣蛋白前β位分解酶1的酶, 其促進β-淀粉樣蛋白生成和調控另一個細胞過程導致記憶喪失而引起AD的發病. miR-29a調控的靶基因NAV3在AD患者大腦中表達上升, 免疫組織化學研究表明, NAV3的表達在AD患者大腦皮層退化錐體神經元中明顯增強. miR-29a就是通過增強AD大腦中的神經元NAV3的表達來影響神經退行性過程的[16]. miR-146a調節免疫反應并影響AD的發病. miR-146a參與重要的細胞功能, 其調控的靶基因bcl2控制線粒體功能和細胞老化, 與AD的發生發展有重要關系[17].

為了進一步分析圖1中miRNA及其靶基因與AD的關系, 使用在線分析網站DAVID (http://david.abcc.ncifcrf.gov/home.jsp)對預測的靶基因進行KEGG通路、 生物過程、 分子功能分析. 表3為預測靶基因的KEGG通路分析結果.

由表3可以看出, 預測的靶基因中大量基因參與了癌癥相關的信號通路, MAPK信號通路, TGF-β信號通路, jak-STAT信號通路, 這些都是與AD相關的通路. 在阿爾茨海默癥中, 通過MAPK信號傳導通路, 淀粉樣纖維激活細胞壞死信號通路. MAPK屬于脯氨酸依賴的蛋白激酶的一種, 在病人大腦中誘使tau蛋白過度磷酸化, 繼而導致AD. TGF-β誘導生物鐘基因的失調, 導致神經通路的改變, 這與AD患者的睡眠—覺醒節律異常是有因果聯系的[18]. 與此同時, TGF-β通路在腫瘤的發生發展中有重要的作用, 通過研究其與miRNA機制可以得到新的潛在治療靶標[19]. jak-STAT信號通路是一條由細胞因子刺激的信號轉導通路, 參與細胞的增殖、 分化、 凋亡以及免疫調節等許多重要的生物過程. 雖然沒有確切的解答答案, 但是現有的大量研究表明, 癌癥和阿爾茨海默癥呈現一種負相關的關系. 因此, 在表3中, 預測的靶基因參與了很多癌癥相關通路, 這并不表示這些基因與AD無關. 如p53基因是一種重要的抑癌基因, 它的缺失或突變將明顯增加惡性腫瘤的易感性. 有研究顯示, p53與AD也有著密切的聯系. p53的激活對促進細胞老化具有直接作用, AD患者的皮質神經元中大量表達p53[20].

表3 KEGG通路分析結果

以miR-195為例來詳細討論miRNA在調控網絡中的作用. 研究表明, miR-195表達水平與BACE1表達水平有關, miR-195表達上升, 則BACE1表達下降; 而miR-195表達下降會導致BACE1的表達上升[21]. 其在網絡中調控65個靶基因, 如圖2所示. 通過KEGG通路分析, 得出結果如表4.

miR-195調控的基因有幾個顯著的功能: 蛋白質復合物捆綁(GO:0032403), 蛋白激酶活性(GO:0004672), 轉錄調節因子粘合物(GO:0030528)和轉錄因子粘合物(GO:0008134). 這些都是生物的基本功能, 說明miR-195的重要作用.

圖2 miR-195調控的靶基因

表4 miR-195靶基因KEGG通路分析結果

在此基礎上, 對預測靶基因進行了生物過程、 分子功能分析. 結果如表5~6.

表5 GO注釋—生物過程(BP)注釋

表6 GO注釋—分子功能(MF)注釋

從上述的生物過程、 分子功能和細胞組成分析可知, 參與生物過程靶基因數目最多的是: 管理細胞增殖、 正調控細胞生物合成過程、 酶鏈接受體蛋白信號通路、 堿基、 核苷、 核苷酸和核酸代謝過程的正調控等. 說明這些生物過程在AD發生發展過程中都扮演了重要的角色. 主要的分子功能為DNA粘合物、 過渡金屬離子活性、 鋅離子粘合物、 轉錄調節因子活性、 轉錄因子活性等. 這些基因參與的金屬離子活性和轉錄活性等分子功能, 都與AD產生密切相關.

3 結語

采用TargetScan、 miRanda和Pictar三種軟件預測miRNA的靶基因, 并采用HCTarget算法對預測的靶基因進行檢驗, 融合匹配的miRNA表達數據和基因表達數據, 構建miRNA-mRNA調控網絡, 分析miRNA和靶基因在AD中的表達活性. 同時利用此算法對AD顯著基因進行對應miRNA的預測. 將參與AD的顯著基因, 利用預測算法得到調控它們的miRNA, 再結合基因表達數據分析, 從而得到與AD相關的新的miRNAs.

miRNAs通過轉錄后水平調控細胞蛋白質的表達, 在神經系統的生長發育、 分化及功能執行中發揮重要的作用. 腦組織內miRNA的異常表達可通過多種途徑影響阿爾茨海默癥的發生和發展. 對于miRNA的研究將有助于深入了解阿爾茨海默癥的發病機制. 選取miR-195、 miR-107、 miR-29a等與AD發病密切相關的miRNA, miR-195和miR-107調控的靶基因BACE1通過促進β-淀粉樣蛋白生成和調控另一個細胞過程導致記憶喪失而引起AD的發病. 受miR-29a調控的靶基因NAV3影響神經退行性過程, 與AD的發病密切相關. 選取的miRNA貫穿AD的發生發展過程, 為AD的致病機理研究, 臨床診斷提供新的方法.

[1] 應俠, 吳振, 雷嚴, 等. 阿爾茨海默病的發病機制及治療藥物研究進展[J]. 中國藥房, 2014, 25(33): 3 152-3 155.

[2] 董賢慧, 柴錫慶. 阿爾茨海默病發病機制研究進展[J]. 中國老年學雜志, 2014(20): 148.

[3] Nunez-Iglesias J, Liu C C, Morgan T E,etal. Joint genome-wide profiling of miRNA and mRNA expression in Aizheimer’s disease cortex reveals altered miRNA regulation[J]. The Public Library of Science One, 2010, 5(2): e8898.

[4] Wang W X, Rajeev B W, Stromberg A J,etal. The expression of microRNA miR 107 decreases early in Alzheimer’s disease and may accelerate disease progression through regulation of beta-site amyloid precursor protein-cleaving enzyme l[J]. The Journal of Neuroscience, 2008, 28(5): 1 213-1 223.

[5] Justin L M, Lahiri D K.MieroRNA 101 downregulates Alzheimer’s amyloid-8 precursor protein levels in human cell cultures and is differentially expressed[J]. Biochemical and Biophysical Research Communications, 2011, 404: 889-895.

[6] Sun E, Shi Y. Micrornas: small molecules with big roles in neurodevelopment and diseases[J]. Experimental Neurology, 2015, 268: 46-53.

[7] Su N, Wang Y, Qian M,etal. Predicting microRNA targets by integrating sequence and expression data in cancer[C]// 2011 IEEE International Conference on Systems Biology (ISB). Zhuhai: IEEE, 2011: 219-224.

[8] Tiribuzi R, Crispoltoni L, Porcellati S,etal. miR128 up-regulation correlates with impaired amyloidβ(1-42) degradation in monocytes from patients with sporadic Alzheimer’s disease[J]. Neurobiology of Aging, 2014, 35(2): 345-356.

[9] Long J M, Ray B, Lahiri D K. MicroRNA-339-5p down-regulates protein expression ofβ-site amyloid precursor protein-cleaving enzyme 1 (BACE1) in human primary brain cultures and is reduced in brain tissue specimens of Alzheimer disease subjects[J]. The Journal of Biological Chemistry, 2014, 289(8): 5 184-5 198.

[10] Lee K H, Lin F C, Hsu T I,etal. MicroRNA-296-5p (miR-296-5p) functions as a tumor suppressor in prostate cancer by directly targeting Pin1[J]. Biochimica et Biophysica Acta (BBA)-Molecular Cell Research, 2014, 1 843(9): 2 055-2 066.

[11] Lewis B P, Burge C B, Bartel D P. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets[J]. Cell, 2005, 120(1): 15-20.

[12] Krek A, Grün D, Poy M N,etal. Combinatorial microRNA target predictions[J]. Nature Genetics, 2005, 37(5): 495-500.

[13] Enright A J, John B, Gaul U,etal. MicroRNA targets in Drosophila[J]. Genome Biology, 2003, 5(1): R1.

[14] Tu K, Yu H, Hua Y J,etal. Combinatorial network of primary and secondary microRNA-driven regulatory mechanisms[J]. Nucleic Acids Research, 2009, 37(18): 5 969-5 980.

[15] Nelson P T, Wang W X. MiR-107 is reduced in Alzheimer’s disease brain neocortex: validation study[J]. Journal of Alzheimer’s Disease, 2010, 21(1): 75-79.

[16] Shioya M, Obayashi S, Tabunoki H,etal. Aberrant microRNA expression in the brains of neurodegenerative diseases: miR-29a decreased in Alzheimer disease brains targets neurone navigator 3[J]. Neuropathology and Applied Neurobiology, 2010, 36(4): 320-330.

[17] Rippo M R, Olivieri F, Monsurrò V,etal. MitomiRs in human inflamm-aging: a hypothesis involving miR-181a, miR-34a and miR-146a[J]. Experimental Gerontology, 2014, 56: 154-163.

[18] Gast H, Gordic S, Petrzilka S,etal. Transforming growth factor‐beta inhibits the expression of clock genes[J]. Annals of the New York Academy of Sciences, 2012, 1261(1): 79-87.

[19] Butz H, Rácz K, Hunyady L,etal. Crosstalk between TGF-βsignaling and the microRNA machinery[J]. Trends in Pharmacological Sciences, 2012, 33(7): 382-393.

[20] 石海濱. 阿爾茨海默病與惡性腫瘤相關性的研究進展[J]. 醫藥前沿, 2014 (2): 117-118.

[21] Zhu H C, Wang L M, Wang M,etal. MicroRNA-195 downregulates Alzheimer’s disease amyloid-βproduction by targeting BACE1[J]. Brain Research Bulletin, 2012, 88(6): 596-601.

(責任編輯: 洪江星)

Target genes prediction for miRNA in Alzheimer’s disease

KONG Wei1, ZHANG Xin1, MOU Xiaoyang2

(1. Information Engineering College, Shanghai Maritime University, Shanghai 201306, China;2. Department of Chemistry and Biochemistry, Rowan University, New Jersey 08028, USA)

The research of miRNAs will help to further understand the pathogenesis of Alzheimer’s disease. In this paper, we use matching gene expression data and the miRNA expression data to estimate and analyze the miRNAs’ target genes. First we estimate the target genes of the chosen differentially expressed miRNAs. Then we use HCTarget algorithm to validate the forecast target genes. At last we built regulation network of miRNAs and its target genes. The regulation network contains 11 identified miRNAs associated with AD and 6 AD disease-causing genes. By analyzing the network we found the miRNAs and its target genes’ active trend in two cases of AD: normal and disease, biological analysis also confirmed their important role in the AD, this provided the basis for the research of the pathogenesis of AD.

mircoRNA; Alzheimer’s disease; target gene; prediction; HCTarget algorithm

2014-10-22

孔薇(1977-), 教授, 主要從事生物信息學方面的研究, weikong@shmtu.edu.cn

國家自然科學基金資助項目(61271446; 61003093)

10.7631/issn.1000-2243.2015.06.0851

1000-2243(2015)06-0851-08

Q343.1

A

主站蜘蛛池模板: 日韩欧美国产中文| 欧美日韩午夜| 国产一区二区福利| 99视频在线免费| 尤物精品视频一区二区三区| 日韩在线播放欧美字幕| 久久精品无码专区免费| 91黄视频在线观看| 国产成人91精品免费网址在线| 亚洲国产精品日韩专区AV| 性色一区| 污污网站在线观看| 亚洲中文字幕无码爆乳| 久久精品国产国语对白| 日韩在线视频网| 精品视频第一页| 国产91无毒不卡在线观看| 亚洲欧美激情另类| 国产欧美视频综合二区| 亚洲欧洲日韩久久狠狠爱| 思思99思思久久最新精品| 色综合久久无码网| 亚洲AV电影不卡在线观看| 永久免费精品视频| 2021天堂在线亚洲精品专区| 午夜成人在线视频| 国产精品99一区不卡| 一级毛片免费观看不卡视频| 亚洲色图欧美一区| 国产亚洲视频中文字幕视频| 香蕉网久久| 99久久无色码中文字幕| 欧美α片免费观看| 噜噜噜综合亚洲| 欧美中文一区| 亚洲色图欧美视频| 2021亚洲精品不卡a| 欧美不卡视频在线观看| 欧美中文字幕一区| 亚洲av无码牛牛影视在线二区| 日本道中文字幕久久一区| 伊人久久大香线蕉综合影视| 国产尤物视频在线| 精品国产美女福到在线直播| 亚洲色精品国产一区二区三区| 欧美色伊人| 久久久久88色偷偷| 九色在线视频导航91| 国产aaaaa一级毛片| 四虎永久在线视频| 亚洲欧美自拍视频| 国产打屁股免费区网站| 午夜精品久久久久久久无码软件| 麻豆精品在线播放| 亚洲六月丁香六月婷婷蜜芽| 国产日韩AV高潮在线| 亚洲国产欧美国产综合久久| 久久99精品久久久大学生| 伊人色天堂| 午夜在线不卡| 亚洲男人的天堂久久香蕉 | 亚洲第一成人在线| 国产精品嫩草影院av| 丝袜久久剧情精品国产| 色欲国产一区二区日韩欧美| 久久久久亚洲精品无码网站| 国产成人久久777777| 最新无码专区超级碰碰碰| 国产屁屁影院| 国产xx在线观看| 91精品国产自产在线观看| 中文字幕中文字字幕码一二区| 欧美中文字幕第一页线路一| 永久免费AⅤ无码网站在线观看| 午夜福利在线观看入口| 91视频首页| 91久久国产成人免费观看| 视频二区国产精品职场同事| 二级特黄绝大片免费视频大片| 全午夜免费一级毛片| 日本欧美在线观看| www.91在线播放|