• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于毛竹重測(cè)序的基因可變剪接研究

    2019-02-12 12:12:42王炯亮趙韓生高志民
    世界竹藤通訊 2019年1期
    關(guān)鍵詞:毛竹附表木質(zhì)素

    王炯亮 趙韓生 高志民

    (國(guó)際竹藤中心北京100102)

    竹子是生長(zhǎng)在溫帶和熱帶森林中的多年生草本植物,其纖維素和半纖維素含量與木本植物相媲美。竹子在全基因組層面的研究有限,僅在2013年首次報(bào)道了毛竹基因組草圖 (2.05 Gb),Scaffold N50為328 Kb,并預(yù)測(cè)出31 987個(gè)基因。得益于測(cè)序技術(shù)和分析方法的改進(jìn),獲得高精確度、連續(xù)性更好的染色體水平毛竹基因組將會(huì)促進(jìn)竹子基因功能和進(jìn)化的研究。

    在真核生物中,可變剪接 (Alternative splicing,AS)是由有限的基因增加蛋白質(zhì)復(fù)雜性和多樣性的主要機(jī)制。經(jīng)預(yù)測(cè),人類中超過95%的多外顯子基因可以表達(dá)生成多個(gè)不同的剪接體 (isoform),而且據(jù)報(bào)道,在擬南芥、大豆、二穗短柄草、棉花、玉米、水稻中分別約有 61%、52%、42%、40%、40%、33%的基因存在可變剪接事件。單個(gè)基因可產(chǎn)生不同的剪接體,是基因具有功能可塑性的主要原因,可能在植物生長(zhǎng)發(fā)育、脅迫反應(yīng)、信號(hào)轉(zhuǎn)導(dǎo)和開花時(shí)間等生物過程中扮演著重要角色。在一定程度上,物種特有的可變剪接將使有限蛋白質(zhì)編碼的基因具有功能多樣性。然而,可變剪接調(diào)控植物進(jìn)化的機(jī)制尚不清楚;可變剪接在保守程度不同的基因中的特性同樣仍不清楚。

    基因組質(zhì)量和轉(zhuǎn)錄組測(cè)序深度對(duì)可變剪接的鑒定有著至關(guān)重要的作用。由于毛竹第一版基因組組裝不完整、且片段分散,轉(zhuǎn)錄組的樣本組織數(shù)量少覆蓋率低,所以很難得到完整的可變剪接圖譜,全面鑒定毛竹可變剪接需要高質(zhì)量的基因組和更多的轉(zhuǎn)錄組數(shù)據(jù)。因此,本研究利用從中國(guó)6個(gè)主要竹產(chǎn)區(qū)采集的毛竹 (Phyllostachys edulis)26個(gè)不同組織樣本測(cè)序獲得的轉(zhuǎn)錄組數(shù)據(jù),進(jìn)行全基因組分析,最終獲得毛竹的可變剪接圖譜。轉(zhuǎn)錄組數(shù)據(jù)由Illumina和Pacific Biosciences(PacBio) 平臺(tái)測(cè)序產(chǎn)生。本研究鑒定得到大量可變剪接基因和可變剪接事件,并對(duì)它們進(jìn)行了分類。另外,通過全基因組范圍的研究,確定了氨基酸保守性與可變剪接之間的關(guān)系,并研究木質(zhì)素生物合成相關(guān)基因可變剪接的進(jìn)化規(guī)律。本研究不僅為進(jìn)一步研究竹子基因功能和調(diào)控網(wǎng)絡(luò)提供了全面的可變剪接圖譜,而且從進(jìn)化角度揭示了可變剪接的作用。

    1 材料方法與數(shù)據(jù)分析

    1.1 樣本采集

    為得到完整的可變剪接圖譜,于2015年春從中國(guó)6個(gè)主要竹產(chǎn)區(qū)采集毛竹樣本進(jìn)行測(cè)序。6個(gè)采集地分別是:江蘇宜興 (31°15′08.41″N, 119°43′42.55″E,海拔212 m)、 浙江天目山 (30°19′13.42″N,119°26′55.21″E,海拔 480 m)、湖北咸寧 (29°81′10.02″N, 114°31′21.12″E, 海拔 150 m)、 湖南桃江 (28°28′39.74″N, 112°11′18.62″E, 海拔 320 m)、 廣西桂林(28°28′39.74″N, 112°11′18.62″E, 海拔 216 m)、 貴州赤水 (28°28′15.27″N,105°59′41.43″E, 海拔 120 m)。收集了處于不同生長(zhǎng)階段毛竹的26個(gè)組織,包括地下莖、根、筍、葉、籜片和芽。有關(guān)生物樣本采集的詳細(xì)信息參見原文附表S19。

    1.2 基因組測(cè)序、組裝、注釋

    使用WGS和Hi-C策略對(duì)毛竹基因組進(jìn)行組裝,并使用前人的研究方法進(jìn)行注釋。詳細(xì)描述見Protocols.io。Hi-C庫(kù)是按照前文所描述的流程進(jìn)行準(zhǔn)備,具體說明見原文附件。

    1.3 RNA的分離、Illumina RNA-seq建庫(kù)與測(cè)序、數(shù)據(jù)分析

    使用前人的方法對(duì)RNA進(jìn)行分離、純化、濃度測(cè)定、反轉(zhuǎn)錄和cDNA文庫(kù)的構(gòu)建。所有cDNA文庫(kù)的構(gòu)建和標(biāo)準(zhǔn)化過程信息詳見原文附件。在質(zhì)控之后,對(duì)文庫(kù)進(jìn)行光學(xué)檢驗(yàn),隨后使用HiSeq-4000平臺(tái) (雙端測(cè)序,插入長(zhǎng)度為150 bp)進(jìn)行測(cè)序。然后,計(jì)算了序列的質(zhì)量,并使用FastQC(軟件版本0.11.3)在默認(rèn)參數(shù)下對(duì)低質(zhì)量序列進(jìn)行過濾。最后,使用RNA-SeQC(軟件版本1.1.8)在默認(rèn)參數(shù)下計(jì)算得到RNA-seq中的關(guān)鍵統(tǒng)計(jì)指標(biāo)。

    使用AStalavista(軟件版本 4.0) (AStalavista,RRID:SCR 001815)在默認(rèn)參數(shù)下鑒別可變剪接基因和可變剪接事件。隨后,分析比較了主要的4類可變剪接IR (intron retention)、A3SS (alternative 3′splice site donor)、 A5SS (alternative 5′splice site acceptor) 和 ES (exon skipping)。 使用 Ontologizer(軟件版本 2.0) 和從Gene Ontology數(shù)據(jù)庫(kù) (GO,RRID:SCR_002811)得到的注釋信息對(duì)基因進(jìn)行富集分析。最后,基于表達(dá)水平 (每百萬個(gè)比對(duì)上的序列里面每千個(gè)堿基中的片段總數(shù))計(jì)算了基因在每個(gè)樣本中的組織特異性值 (Ts)。

    1.4 Iso-Seq文庫(kù)的構(gòu)建、測(cè)序與數(shù)據(jù)分析

    按照測(cè)序儀制造商 (PacBio)提供的方法,對(duì)Iso-Seq文庫(kù)進(jìn)行構(gòu)建和測(cè)序。首先根據(jù)生物信息學(xué)預(yù)測(cè)得到的轉(zhuǎn)錄本長(zhǎng)度分布 (原文附表S22)構(gòu)建了3類SMRTBell文庫(kù) (轉(zhuǎn)錄本長(zhǎng)度在1 Kb到2 Kb之間的3個(gè)文庫(kù),轉(zhuǎn)錄本長(zhǎng)度為2 Kb到3 Kb之間的2個(gè)文庫(kù),轉(zhuǎn)錄本長(zhǎng)度大于3 kb的4個(gè)文庫(kù)),然后使用PacBio測(cè)序平臺(tái)對(duì)這9個(gè)SMRT文庫(kù)進(jìn)行測(cè)序。

    通過分析PacBio RSII的測(cè)序數(shù)據(jù),得到了合理的全長(zhǎng)轉(zhuǎn)錄本異構(gòu)體,隨后合并不同文庫(kù)得到的轉(zhuǎn)錄本異構(gòu)體,并在對(duì)序列進(jìn)行插入、分類、聚類等處理后,去除冗余的數(shù)據(jù),最終得到了合理的轉(zhuǎn)錄本。組裝得到的轉(zhuǎn)錄本使用PASA(軟件版本2.0.2) (PASA,RRID:SCR_014656) 在默認(rèn)參數(shù)下比對(duì)到參考基因組上。最后,與短片段序列一樣,使用AStalavista在默認(rèn)參數(shù)下對(duì)得到的gtf文件進(jìn)行分析,從而鑒別可變剪接。

    1.5 進(jìn)化分析

    依據(jù)前人的研究方法,進(jìn)行基因家族的鑒定、系統(tǒng)發(fā)育樹的構(gòu)建和分化時(shí)間的預(yù)測(cè)。在原文附件和protocol.io提供了詳細(xì)的信息。

    1.6 木質(zhì)素生物合成相關(guān)基因的全基因組鑒定與分析

    從ENSEMBL數(shù)據(jù)庫(kù) (Ensembl,RRID:SCR_002344)下載了擬南芥 (基因組版本 10)、小麥(基因組版本3.1)、水稻 (基因組版本 7.0)、毛果楊 (基因組版本JGI 2.0.31)和高粱 (基因組版本3.1)5個(gè)物種基因組的序列,收集得到了經(jīng)過實(shí)驗(yàn)驗(yàn)證的140個(gè)木質(zhì)素生物合成相關(guān)的基因 (原文附表S28)。通過局部比對(duì)搜索 (BLAST) (美國(guó)國(guó)家生物技術(shù)信息中心 [NCBI]BLAST,RRID:SCR_004870)和結(jié)構(gòu)分析2種方法鑒定毛竹中木質(zhì)素生物合成相關(guān)的基因。BLAST搜索閾值設(shè)置為:E-value<1e-10,identity>40%,coverage>95%。然后使用 hmm-search(軟件版本 3.1b2)和 Pfam-A.hmm數(shù)據(jù)庫(kù) (2017.5.31發(fā)布的版本)對(duì)所有由BLAST搜索得到的序列進(jìn)行分析和人工矯正后,進(jìn)行系統(tǒng)發(fā)育分析。最后,使用yn00(PAML軟件包中預(yù)測(cè)同義替換率和非同義替換率的子軟件包)對(duì)木質(zhì)素生物合成相關(guān)的13個(gè)基因家族進(jìn)行同義替換率的計(jì)算。Ks率通過公式 T=Ks/2r(r=6.5×10-9)計(jì)算分化時(shí)間。

    對(duì)木質(zhì)素生物合成相關(guān)基因的編碼序列正選擇作用分析:首先使用PROBCONS(軟件版本1.12)(ProbCons,RRID:SCR_011813)對(duì)每個(gè)家族的蛋白質(zhì)序列進(jìn)行比對(duì) (迭代次數(shù)設(shè)為1 000,其他參數(shù)默認(rèn))。隨后將比對(duì)后的蛋白質(zhì)翻譯為相應(yīng)的編碼序列。接著使用Gblocks(軟件版本0.91b)軟件從比對(duì)的結(jié)果提取保守序列塊,并使用jModelTest(軟件版本2.1.6)依據(jù)貝葉斯信息準(zhǔn)則分析保守域得到最佳的核酸替換矩陣。然后,使用PhyML(軟件版本 3.0) (PhyML, RRID: SCR_014629) (1 000次檢驗(yàn))選擇最佳的核酸替換矩陣重新構(gòu)建了系統(tǒng)發(fā)育樹。最后,使用 PAML(軟件版本 4.8) 的branch-site模型對(duì)系統(tǒng)發(fā)育樹中的一些分支進(jìn)行正選擇作用分析 (詳細(xì)的信息見protocols.io)。

    2 結(jié)果與分析

    2.1 數(shù)據(jù)描述

    經(jīng)不同測(cè)序方法產(chǎn)生用于毛竹基因組組裝的基因組數(shù)據(jù)總共約有603.3 Gb。首先,全基因組序列(Whole-genome sequence,WGS)由新得到的約154 Gb數(shù)據(jù)和已有的約220 Gb數(shù)據(jù)進(jìn)行組裝。之后使用由157 Gb原始數(shù)據(jù)經(jīng)質(zhì)量控制得到的約17.58 Gb Hi-C有效序列進(jìn)行Hi-C輔助組裝 (原文附表S1)。此外,由Illumina和PacBio平臺(tái)分別產(chǎn)生的約379 Gb和5 Gb數(shù)據(jù)用于轉(zhuǎn)錄組分析 (原文附表S2-S7)。本研究通過使用染色體水平的參考基因組和大量轉(zhuǎn)錄組數(shù)據(jù),在毛竹中鑒定出了25 225個(gè)可變剪接基因和266 711個(gè)可變剪接事件。

    2.2 數(shù)據(jù)分析

    2.2.1 毛竹染色體水平基因組的組裝和基因注釋

    為提高毛竹基因組的質(zhì)量,本研究對(duì)61個(gè)文庫(kù)進(jìn)行測(cè)序 (原文附表S1),總共得到讀長(zhǎng)在76 bp到250 bp之間的約603.3 Gb的基因組數(shù)據(jù)。隨后為獲得高質(zhì)量的基因組,使用了不同的組裝策略。首先,WGS組裝得到的基因組大小達(dá)到1.91 Gb,其Contig N50和Scaffold N50的長(zhǎng)度分別達(dá)到55 Kb和894 Kb(原文附表S8)。由WGS組裝得到的新基因組的質(zhì)量和各統(tǒng)計(jì)量與之前的版本相比都有明顯的提高(原文附表 S9,S10),比如 Contig N50和 Scaffold N50分別增加172%和358%,且未知序列減少43%。其次,使用Hi-C數(shù)據(jù)進(jìn)行輔助組裝,得到了總長(zhǎng)度為1.91 Gb的基因組,其Contig N50和Scaffold N50長(zhǎng)度分別達(dá)到53.29 Kb和79.90 Mb(原文圖1A,1B)。由WGS組裝得到的約93.17%Scaffold可以錨定到24條染色體上 (原文附表S10),而且Scaffold N50的長(zhǎng)度增加了約89倍 (原文表1)。根據(jù)關(guān)聯(lián)圖和組裝結(jié)果,可以明顯區(qū)分開24條染色體。將毛竹的染色體比對(duì)到水稻基因組上面,發(fā)現(xiàn)平均覆蓋率約為59.77% (原文附圖S2,附表S11)。此外,使用人工細(xì)菌染色體 (BAC)序列、全長(zhǎng)cDNA和一些已知的毛竹基因序列對(duì)毛竹染色體水平的基因組進(jìn)行評(píng)估 (原文附圖S3,附表S12-14),發(fā)現(xiàn)染色體水平基因組較第一版基因組的覆蓋范圍更廣,且準(zhǔn)確性更高。

    在注釋重復(fù)序列后,染色體水平基因組將更有利于后續(xù)的基因注釋 (原文附表S15)?;诖罅康霓D(zhuǎn)錄組 (原文附表S16)、全長(zhǎng)cDNA和同源蛋白數(shù)據(jù),預(yù)測(cè)出了51 074個(gè)具有完整結(jié)構(gòu)的蛋白編碼基因 (原文附表S17),這些基因的內(nèi)含子和外顯子平均長(zhǎng)度分別為668 bp和284 bp(原文附表S18)。通過單分子實(shí)時(shí)測(cè)序數(shù)據(jù)和人工校驗(yàn)對(duì)不合理的注釋進(jìn)行校正,通過添加非編碼區(qū)域 (UTR)注釋和內(nèi)部結(jié)構(gòu)調(diào)整,改進(jìn)了約17%的基因模型 (原文附表S19)。基因組注釋完整性評(píng)估 (單拷貝同源物分析)顯示:毛竹 (95.2%)的注釋比玉米(92.2%)更完整,與水稻 (95.6%)接近 (原文圖1D,附表S20)。與前一版本的注釋相比,本版本中97.23%的基因模型在公共數(shù)據(jù)庫(kù)被鑒定,這有助于精確檢測(cè)可變剪接事件 (原文附表S21)。原文附表S22-S24和原文附圖S3-S9提供了基因模型預(yù)測(cè)和基因組進(jìn)化研究相關(guān)的詳細(xì)信息。此外,最新版本的基因組和基因注釋已經(jīng)通過GigaScience的GiGaDB資源庫(kù)發(fā)布。包括新發(fā)布的基因組序列、基因集、重復(fù)元件、tRNA、miRNA和基因簇等最新毛竹基因組的數(shù)據(jù)為基因組學(xué)研究、遺傳學(xué)研究、分子生物學(xué)實(shí)驗(yàn)等提供了可靠的數(shù)據(jù)資源。

    2.2.2 轉(zhuǎn)錄組數(shù)據(jù)分析

    為促進(jìn)毛竹全基因組水平的可變剪接圖譜研究,分析轉(zhuǎn)錄后水平影響可變剪接的因素,本研究利用Illumina HiSeq-4000平臺(tái)進(jìn)行高通量轉(zhuǎn)錄組測(cè)序(RNA sequencing,RNA-seq)??偣矊?duì) 26個(gè)獨(dú)立RNA樣本進(jìn)行雙端測(cè)序 (讀長(zhǎng)為150 bp,原文附表S2,附圖S10-11),測(cè)得每個(gè)樣本高質(zhì)序列平均約為9 000萬條 (約為 13.6 Gb),占原始序列的92.78%。高質(zhì)量序列中約80.57%可以比對(duì)到參考基因組唯一的位置上,并將其標(biāo)記為唯一序列(Unique reads,原文附表S3-4)。大部分序列比對(duì)到外顯子區(qū)域,外顯子的平均比對(duì)率為81.94%。余下的序列中,有8.46%比對(duì)到內(nèi)含子區(qū)域,有9.6%比對(duì)到基因間區(qū)域 (原文附表S5,附圖S12-13)。每個(gè)樣本的平均外顯子覆蓋率約為2 521倍 (原文附圖S14)。因此,大規(guī)模、高深度、高質(zhì)量的轉(zhuǎn)錄組測(cè)序和染色體水平的參考基因組,將有助于精準(zhǔn)識(shí)別基因組中的可變剪接。

    為準(zhǔn)確鑒定全長(zhǎng)的剪接體,使用PacBio平臺(tái)對(duì)由毛竹26個(gè)樣本混合得到的RNA進(jìn)行全長(zhǎng)可變剪接體 (FL-cDNA sequencing of alternatively spliced isoforms,Iso-Seq)的測(cè)序。根據(jù)全部樣本的轉(zhuǎn)錄本長(zhǎng)度分布,由混合樣本構(gòu)建了三類單分子實(shí)時(shí)Bell文庫(kù) (3個(gè)轉(zhuǎn)錄本長(zhǎng)度在1 Kb到2 Kb之間的文庫(kù),2個(gè)轉(zhuǎn)錄本長(zhǎng)度在2 Kb到3 Kb之間的文庫(kù),4個(gè)轉(zhuǎn)錄本長(zhǎng)度大于3 kb的文庫(kù)),并對(duì)9個(gè)文庫(kù)進(jìn)行測(cè)序,共產(chǎn)生約5 Gb的原始數(shù)據(jù)和214 372條序列(read-of-insert, ROI)。 在所有的ROI中,有133 599個(gè)是全長(zhǎng)ROI(包含了5′、3′和 poly(A) ),剩下的ROI則是非全長(zhǎng)序列 (原文附表S7,附圖S15)。本研究通過將ROI比對(duì)到新基因組進(jìn)行精確度評(píng)估,評(píng)估顯示每一個(gè)核苷酸的錯(cuò)誤率大約是2.05%,其中0.32%的錯(cuò)配,0.98%的插入,0.75%的刪除。

    2.2.3 毛竹中有大量的基因存在可變剪接

    借助高質(zhì)量基因組和高通量轉(zhuǎn)錄組數(shù)據(jù),使用前人的分析流程在全基因組范圍內(nèi)鑒定毛竹的可變剪接。結(jié)果顯示,總共在25 225個(gè)可變剪接基因中鑒定了266 711個(gè)特異的可變剪接事件。在所有鑒定的可變剪接基因中,有12 653個(gè)基因在基因注釋環(huán)節(jié)中被注釋為可變剪接基因,剩下的12 572個(gè)基因被認(rèn)定為新鑒定的可變剪接基因 (原文附圖S16)。

    作為Illumina RNA-Seq分析的平行實(shí)驗(yàn),本研究還用Iso-Seq數(shù)據(jù)在相同的流程下檢測(cè)可變剪接。結(jié)果顯示,共鑒定出2 218個(gè)可變剪接基因和4 246個(gè)可變剪接事件。對(duì)兩種數(shù)據(jù)得到的可變剪接結(jié)果進(jìn)行的重疊分析 (評(píng)估可變剪接基因預(yù)測(cè)結(jié)果的可靠程度)顯示,由Iso-Seq數(shù)據(jù)鑒定得到的81.21%可變剪接事件和97.34%可變剪接基因分別與由RNA-seq數(shù)據(jù)鑒定得到的完全重疊。隨后,對(duì)可變剪接進(jìn)行了分類,其中主要的可變剪接類型分別是IR、A3SS、A5SS和ES。平均80.37%的可變剪接事件和95.59%的可變剪接基因?qū)儆谶@主要的4類可變剪接類型 (原文附圖S17)。由PacBio與Illumina平臺(tái)數(shù)據(jù)分別鑒定得到的兩組可變剪接基因之間的覆蓋率很高,證明由計(jì)算機(jī)預(yù)測(cè)得到的可變剪接具有較高的可靠性。

    可變剪接事件數(shù)量與可變剪接基因數(shù)量、主要的4類可變剪接的基因數(shù)量之間存在強(qiáng)烈的正相關(guān)(Mann-Whitney U 檢驗(yàn), R2>0.91, P<0.05) (原文圖2C)。毛竹中主要的4類可變剪接由經(jīng)典的剪接模型 (GT-AG、GC-AC和AT-AC剪接位點(diǎn))從可變剪接事件中鑒定得到。如原文圖2B所示,在主要的4類可變剪接事件中,IR(38.22%)數(shù)量最多,接下來分別是 A3SS(20.20%)和 A5SS(10.48%),最少的是ES(2.92%)。

    基因功能富集分析顯示,所有樣本共有的885個(gè)可變剪接基因顯著富集于RNA代謝、mRNA加工、RNA修飾和RNA剪接等功能 (原文附表S25)??勺兗艚佑兄@著的組織和發(fā)育特異性,鑒定出了181 105個(gè)組織特異的可變剪接事件,是所有可變剪接事件的2/3(67.57%,標(biāo)識(shí)為 “組織間可變剪接”)。剩下的1/3的可變剪接事件通過比較單個(gè)組織內(nèi)的異構(gòu)轉(zhuǎn)錄本得到 (標(biāo)識(shí)為 “組織內(nèi)可變剪接” )(原文附圖S18)。

    轉(zhuǎn)座子 (TE)分析顯示,在26 366個(gè)基因 (占所有基因的51.62%)中存在著轉(zhuǎn)座子的插入,這些插入轉(zhuǎn)座子的總長(zhǎng)度約為46 Mb。通過分析轉(zhuǎn)座子插入內(nèi)含子的位置,發(fā)現(xiàn)插入轉(zhuǎn)座子的內(nèi)含子主要集中于基因的頭尾部分 (原文附圖S19)。此外,使用率最高和分布最多的剪接位點(diǎn)是GT-AG(占所有可變剪接事件的 97.31%),接下來的是 GC-AG(2.33%)和 GT-AT(0.32%)位點(diǎn) (原文附圖S20)。除經(jīng)典的剪接位點(diǎn) (GT-AG、GC-AG和ATAC)之外,將剩下的2 406個(gè)剪接位點(diǎn)標(biāo)記為非經(jīng)典剪接位點(diǎn),這些非經(jīng)典剪接位點(diǎn)包括2 373個(gè)GTAT類型的剪接位點(diǎn)和33個(gè)其他類型的剪接位點(diǎn)。

    2.2.4 毛竹中可變剪接的進(jìn)化分析

    通過在8個(gè)物種 (無油樟、擬南芥、油棕、二穗短柄草、水稻、浮萍、高粱和毛竹)中鑒定全基因組同源基因和系統(tǒng)發(fā)育樹的構(gòu)建 (原文圖3A,3B),根據(jù)同源基因起源的時(shí)間不同,定義了8個(gè)特異的同源基因數(shù)據(jù)集 (D8-D1)。比如,同源基因數(shù)據(jù)集7(D7)只包含在1.649億年 (Mya)到2.136億年前之間起源的同源基因 (原文圖3A)。此外,分別從上述的8個(gè)數(shù)據(jù)集中提取單拷貝基因集,標(biāo)記為D8s-D1s。毛竹特有的基因 (D1,包含4 023個(gè)同源基因)較不保守,而8個(gè)物種中都存在的基因 (D8,包含18 997個(gè)同源基因)則高度保守,基因保守程度由D8到D1單調(diào)下降。在所有數(shù)據(jù)集中都檢測(cè)到了可變剪接,但可變剪接基因的比例由D8到D1逐漸降低 (Mann-Whitney U檢驗(yàn),P<0.05)。同樣的趨勢(shì)也出現(xiàn)在單拷貝基因集 (D8s-D1s)中。因此,毛竹保守基因集中包含更多的可變剪接基因。

    通過對(duì)每個(gè)數(shù)據(jù)集中主要的4類可變剪接分布的研究,發(fā)現(xiàn)不同類型可變剪接的分布趨勢(shì)相同(原文圖3C),但不同類型可變剪接在不同數(shù)據(jù)集中所占的比例不同 (IR >A3SS>A5SS>ES,Chisquare檢驗(yàn),P>0.86)。D8中IR所占的比例為57.76%,約是D1中IR (16.95%)的3.4倍。其他類型可變剪接的比例隨著基因保守度的降低而增高。在兩類數(shù)據(jù)集中,可變剪接事件的數(shù)量從D8到D1、D8s到D1s逐漸降低 (原文圖3C)。此外,對(duì)組織表達(dá)特異性不同的基因 (maxTs,maxTs=1和maxTs=0分別表示基本表達(dá)和組織特異性表達(dá),詳細(xì)見方法部分)的可變剪接事件比較發(fā)現(xiàn),maxTs與D8-D1這8個(gè)數(shù)據(jù)集中基因的起源時(shí)間存在著負(fù)相關(guān)關(guān)系 (R2>0.86,P<0.01),即組織特異性隨著基因保守程度的降低而增大 (原文圖3D)。綜上所述,保守基因集傾向于包含更多的可變剪接基因和事件,且具備更低的組織特異性。

    為從整體了解可變剪接,認(rèn)識(shí)可變剪接與基因特征之間的關(guān)系和分析影響可變剪接的因素,研究了不同數(shù)據(jù)集中可變剪接的分布與基因特征之間的關(guān)系 (原文附圖S21)。分析表明不同數(shù)據(jù)集 (D8-D1)的基因與基因長(zhǎng)度、對(duì)應(yīng)編碼序列大小、內(nèi)含子大小、外顯子數(shù)量呈正相關(guān) (R2>0.9,P<0.05),與外顯子長(zhǎng)度、內(nèi)含子長(zhǎng)度呈負(fù)相關(guān) (R2>0.81,P<0.05)。除此之外,同樣分析了轉(zhuǎn)座子基因在8個(gè)數(shù)據(jù)集 (D8-D1)中的分布,大體上呈負(fù)相關(guān) (R2>0.77,P<0.05),表明越保守的基因存在越多的轉(zhuǎn)座子插入。

    2.2.5 木質(zhì)素生物合成相關(guān)基因家族的擴(kuò)張和其對(duì)功能多樣性的影響

    在擬南芥、二穗短柄草、水稻、楊樹、毛竹和高粱的基因組序列中系統(tǒng)地鑒定出了13個(gè)與木質(zhì)素生物合成相關(guān)的基因家族,這13個(gè)基因家族中的大多數(shù)都存在擴(kuò)張 (原文附表S26)。在毛竹基因組中每個(gè)基因都有著多個(gè)拷貝,其中木質(zhì)素生物合成相關(guān)基因家族的規(guī)模是最大的 (每個(gè)家族平均有19個(gè)成員)。過氧化物酶 (包含77個(gè)成員)和香豆酸-3-羥化酶 (包含3個(gè)成員)基因家族分別是毛竹中成員數(shù)量最多和最少的基因家族。此外,參與木質(zhì)素生物合成的基因在500萬年到1 600百萬年前之間分化,這與毛竹在700萬年到1 200萬年前之間發(fā)生的全基因組復(fù)制 (WGD)事件相對(duì)應(yīng)。

    之后,對(duì)木質(zhì)素生物合成相關(guān)基因的可變剪接進(jìn)行了分析 (原文圖4)。除了阿魏酸-5-羥化酶(F5H)基因家族中的可變剪接比例過低,查耳酮合成 (CHS))基因家族和咖啡酸鄰O-甲基轉(zhuǎn)移酶(COMT)基因家族沒有檢測(cè)到可變剪接基因之外,總計(jì)13個(gè)基因家族中的10個(gè) (超過總數(shù)的一半)存在可變剪接基因。在肉桂酸4-羥基肉桂酰輔酶A連接酶 (4CL)、羥基肉桂酰轉(zhuǎn)移酶 (HCT)、肉桂醇脫氫酶 (CAD)基因家族中觀察到高比例 (>75%)的可變剪接事件。最后,使用branch-site模型對(duì)參與木質(zhì)素生物合成的相關(guān)基因家族進(jìn)行正選擇作用分析表明,HCT和CAD這兩個(gè)基因家族中的一些基因檢測(cè)到正選擇作用。使用的最適模型得到的系統(tǒng)發(fā)育關(guān)系和對(duì)數(shù)似然比 (log likelihood ratio)等信息詳見原文附表S27。

    3 討論

    3.1 毛竹基因組的完善與質(zhì)量提升

    目前,隨著新技術(shù)的發(fā)展和更多高通量數(shù)據(jù)的出現(xiàn),高通量基因組測(cè)序和改進(jìn)的組裝技術(shù)普遍應(yīng)用于植物基因組研究中。在2013年,通過對(duì)毛竹基因組的初步分析,得到了毛竹的基因組、基因結(jié)構(gòu)、關(guān)鍵功能基因等數(shù)據(jù),還了解了毛竹中發(fā)生的WGD事件。通過本次研究,提高了毛竹基因組的精確度和完整性,更新了基因組注釋,且通過不同物種之間的比較研究精準(zhǔn)定位了毛竹在進(jìn)化上的位置。此外,利用從最新版本毛竹基因組得到的信息對(duì)毛竹不同的生物特性進(jìn)行了詳細(xì)的研究。染色體水平的參考基因組和準(zhǔn)確的注釋將有利于今后毛竹和其他近緣物種的基因組學(xué)研究。

    3.2 毛竹中可變剪接具有普遍性,且在不同組織中表現(xiàn)出多樣性

    RNA-seq和Iso-Seq高通量數(shù)據(jù)有助于精準(zhǔn)檢測(cè)低表達(dá)水平的轉(zhuǎn)錄本、鑒定完整的基因結(jié)構(gòu) (特別是在可變剪接的研究中),通過分析獲得了毛竹全基因組范圍的可變剪接圖譜。通過對(duì)毛竹中可變剪接的研究,加深了對(duì)轉(zhuǎn)錄后調(diào)控層面上的可變剪接的認(rèn)知,包括可變剪接基因與事件的鑒別、不同種類可變剪接的分布、不同種類剪接位點(diǎn)的使用率、可變剪接體外顯子的長(zhǎng)度分布等??勺兗艚颖徽J(rèn)為是基因數(shù)量有限的生物產(chǎn)生多樣性的主要機(jī)制。例如,果蠅中細(xì)胞黏附分子基因通過組合分別含有12、48、33和2個(gè)外顯子的4個(gè)可變剪接基因,最多可產(chǎn)生38 016(12×48×33×2) 種蛋白質(zhì)異構(gòu)體。在毛竹中,共鑒定出266 711個(gè)特異的可變剪接事件和25 225個(gè)可變剪接基因,平均每個(gè)樣本檢測(cè)到15 071個(gè)可變剪接事件和9 080個(gè)可變剪接基因。因此,可變剪接可能存在組織特異性,而且毛竹中實(shí)際的可變剪接比例可能被低估了。隨著測(cè)序深度的增加,可以從表達(dá)水平低的轉(zhuǎn)錄本中檢測(cè)出更多的可變剪接事件。

    根據(jù)觀察,毛竹中地下莖組織比根組織包含更多的可變剪接事件,這可能是因?yàn)檫@兩個(gè)組織在毛竹的生長(zhǎng)過程中扮演的角色不同。處于高速生長(zhǎng)期的竹筍沒有葉子不能進(jìn)行光合作用,所以竹筍生長(zhǎng)所需的大量營(yíng)養(yǎng)物質(zhì)和能量需要通過地下莖從成熟竹子運(yùn)輸?shù)缴L(zhǎng)部位。因此,作為地下莖植物,毛竹的地下莖在營(yíng)養(yǎng)物質(zhì)和能量運(yùn)輸中扮演著核心角色,這可能解釋了為什么會(huì)在地下莖組織樣本中檢測(cè)到更多的可變剪接事件。

    3.3 毛竹一些生物學(xué)特性可能與不同類型可變剪接的比例有關(guān)

    在全基因組范圍上鑒定可變剪接時(shí),經(jīng)常需要分析不同種類可變剪接之間的差異,因?yàn)椴煌N類可變剪接的頻率或比例的差異可能反映了pre-mRNA剪接的差異。不同種類可變剪接的分布結(jié)果表明:IR類型的可變剪接在毛竹中數(shù)量最多,且IR普遍存在于其他處在不同進(jìn)化位置上的植物中,這可能與IR的重要性有關(guān)。然而,與其他植物相比,毛竹中IR類型 (38.22%)和其他類型 (28.18%)的可變剪接比例更高。這可能與毛竹的特性或轉(zhuǎn)錄組測(cè)序深度有關(guān),也可能是兩者共同造成的,未來可以通過比較分析進(jìn)行辨別。此外,毛竹中主要的4類可變剪接的分布與擬南芥、大豆、玉米中的一致。然而,在動(dòng)物和酵母中不同種類可變剪接的分布與植物中的不同。在動(dòng)物和酵母中,可變剪接事件數(shù)量最多的一類是ES,接下來分別是A3SS和A5SS,最少的是IR??梢娭参锖蛣?dòng)物在不同種類可變剪接的分布上存在差異,這暗示了植物和動(dòng)物的基因組結(jié)構(gòu)和剪接位點(diǎn)識(shí)別機(jī)制存在差異。此外,對(duì)單一基因剪接位點(diǎn)的鑒定提供了理解可變剪接和剪接體結(jié)構(gòu)的重要數(shù)據(jù)。與各類可變剪接的分布情況不同的是,毛竹中主要的3種可變剪接位點(diǎn) (GT-AG,GC-AG和AT-AC)分布和之前在動(dòng)物和其他植物中觀察到的規(guī)律是一致的。

    3.4 包含更多可變剪接事件的保守基因可能在進(jìn)化和功能中扮演著重要的作用

    目前,基因保守性和可變剪接之間的關(guān)系仍不清楚。為了探究這個(gè)問題,在基因組層面分析兩類同源基因集中的可變剪接,這2類基因集按保守程度的不同又各自細(xì)分出8個(gè)子集,分別是D8-D1和D8s-D1s。分析表明可變剪接基因傾向于集中在高保守基因集中,且高保守基因集的可變剪接基因包含著更多的可變剪接事件。本研究在2類同源基因集中發(fā)現(xiàn)了相同的趨勢(shì),說明這個(gè)發(fā)現(xiàn)是可靠的。有研究表明復(fù)制是功能分化和新基因產(chǎn)生的主要原因之一,且保守基因在基因互作網(wǎng)絡(luò)中有著更高連接度的傾向。新基因因?yàn)樾录尤氲交蚧プ骶W(wǎng)絡(luò),所以其基因連接度較低,但隨著時(shí)間推移,新基因的連接度和重要程度逐漸增加。本研究表明,保守度高的基因比保守度低的基因傾向于擁有更多的可變剪接事件,這與保守度高的基因在基因互作網(wǎng)絡(luò)中連接度更高的趨勢(shì)是一致的。因此,推測(cè)在進(jìn)化過程中基因連接度的增加可能與可變剪接有關(guān)。此外,與保守度低的基因集相比,保守度高的可變剪接基因集的組織特異性低,表明在基因互作網(wǎng)絡(luò)中有著更高的連接度的基因可能在基礎(chǔ)功能中發(fā)揮著重要作用。因此,推測(cè)一些基因之所以重要可能是因?yàn)樗鼈儼嗟目勺兗艚邮录???勺兗艚幼鳛橐粋€(gè)重要的生物學(xué)過程,在生物功能多樣化的過程中扮演著重要的角色,這可能解釋了保守度高的可變剪接基因?yàn)槭裁窗嗟目勺兗艚邮录?。推測(cè)這個(gè)現(xiàn)象可能不僅存在于毛竹中,而且還存在于其他植物甚至是動(dòng)物中。

    3.5 毛竹中木質(zhì)素生物合成相關(guān)基因家族擴(kuò)增和可變剪接可能與WGD事件有關(guān)

    木質(zhì)素是一種由木質(zhì)素單體組成的復(fù)雜芳香族聚合物,它與纖維素和半纖維素相互作用共同構(gòu)成次生細(xì)胞壁。木質(zhì)素約占毛竹干重的25%。通過可變剪接分析和進(jìn)化分析,檢測(cè)到了木質(zhì)素生物合成相關(guān)基因家族的擴(kuò)增。結(jié)合木質(zhì)素生物合成相關(guān)基因的分化時(shí)間研究結(jié)果和我們前期的研究,估測(cè)毛竹基因組在前700萬年到1 200萬年前之間發(fā)生了一次WGD,這暗示在其進(jìn)化過程中可能存在一次四倍體事件。之后,四倍體祖先進(jìn)化成為現(xiàn)在的二倍體毛竹。WGD可以提供更多的基因拷貝,并通過產(chǎn)生新功能而加速基因的進(jìn)化。因此,毛竹中發(fā)生的WGD事件可能導(dǎo)致了木質(zhì)素生物合成相關(guān)基因的擴(kuò)增。此外,HCT和CAD這兩個(gè)基因家族被檢測(cè)到正選擇作用,且包含著更多的可變剪接事件。HCT介導(dǎo)對(duì)-香豆酰輔酶A(也被CHS介導(dǎo)產(chǎn)生類黃酮)產(chǎn)生木質(zhì)素。HCT與CHS都可以結(jié)合對(duì)-香豆酰輔酶A,它們之間存在相互競(jìng)爭(zhēng)作用。在毛竹中,HCT基因家族與CHS基因家族相比,有著更多的基因成員和可變剪接事件,說明HCT基因家族可能在與CHS基因家族對(duì)-香豆酰輔酶A競(jìng)爭(zhēng)結(jié)合中占據(jù)著主導(dǎo)地位。CAD可以催化不同的底物合成不同的木質(zhì)素。在毛竹中木質(zhì)素芳香聚合物由以下三個(gè)單體組成:對(duì)羥基苯基丙烷 (H)、香草醛 (G)、丁香醛(S)。前人的研究表明,在毛竹中G、S木質(zhì)素的含量高,而H木質(zhì)素的含量低。毛竹中CAD基因家族的擴(kuò)增和正選擇作用可能解釋了毛竹中不同木質(zhì)素單體比例差異帶來的基質(zhì)偏好性。毛竹木質(zhì)素生物合成相關(guān)基因中存在的大量可變剪接事件、基因擴(kuò)增和正選擇作用,這與毛竹強(qiáng)大的木質(zhì)素合成能力一致。

    4 結(jié)論

    為從進(jìn)化角度深入探究毛竹可變剪接情況,組裝出了染色體水平的基因組,并重新進(jìn)行了基因組注釋。通過分析染色體水平的基因組和大量轉(zhuǎn)錄組數(shù)據(jù),得到了毛竹中完整的可變剪接圖譜,鑒定出了25 225個(gè)可變剪接基因中的266 711個(gè)特異的可變剪接事件。此外,毛竹中可變剪接結(jié)果的綜合分析和8個(gè)代表性植物物種的比較分析表明,保守基因有著高表達(dá)和低組織特異性的趨勢(shì)。對(duì)毛竹中木質(zhì)素生物合成相關(guān)的基因進(jìn)行了可變剪接分析和進(jìn)化分析,觀察到木質(zhì)素生物合成相關(guān)的基因家族發(fā)生了擴(kuò)增,這些基因家族包含了大量的可變剪接和正選擇作用。綜上所述,本研究為毛竹獨(dú)特的材性研究和從進(jìn)化角度探索可變剪接提供了重要的資源。

    猜你喜歡
    毛竹附表木質(zhì)素
    附表
    附表
    節(jié)骨草和毛竹苗
    木質(zhì)素增強(qiáng)生物塑料的研究進(jìn)展
    上海包裝(2019年8期)2019-11-11 12:16:14
    附表3 湖南省2018年飼料添加劑、添加劑預(yù)混合飼料產(chǎn)品批準(zhǔn)文號(hào)名單(二)
    湖南飼料(2018年3期)2018-07-03 06:17:58
    湖南省2017年飼料添加劑產(chǎn)品批準(zhǔn)文號(hào)名單
    湖南飼料(2017年3期)2017-07-18 11:15:37
    一種改性木質(zhì)素基分散劑及其制備工藝
    天津造紙(2016年1期)2017-01-15 14:03:29
    一種新型酚化木質(zhì)素胺乳化劑的合成及其性能
    毛竹護(hù)筍養(yǎng)竹技術(shù)探討
    毛竹幼林四季管理技術(shù)
    麻豆一二三区av精品| 日韩av在线大香蕉| 他把我摸到了高潮在线观看| ponron亚洲| 三级男女做爰猛烈吃奶摸视频| 成年版毛片免费区| 亚洲国产日韩欧美精品在线观看| 欧美高清成人免费视频www| 少妇的逼水好多| 国产aⅴ精品一区二区三区波| 99久久精品热视频| 色哟哟哟哟哟哟| 少妇高潮的动态图| 内射极品少妇av片p| 搡老熟女国产l中国老女人| 校园春色视频在线观看| a级毛片a级免费在线| 亚洲av成人av| 黄色视频,在线免费观看| 日韩强制内射视频| 亚洲无线在线观看| 成人午夜高清在线视频| 免费搜索国产男女视频| 97热精品久久久久久| 欧美色视频一区免费| 欧美另类亚洲清纯唯美| 成人精品一区二区免费| 波野结衣二区三区在线| 中文亚洲av片在线观看爽| 欧美+日韩+精品| 国产女主播在线喷水免费视频网站 | 欧美最黄视频在线播放免费| 欧美+日韩+精品| 国产精品国产高清国产av| av.在线天堂| 99久久久亚洲精品蜜臀av| 国产伦精品一区二区三区视频9| avwww免费| 欧美一区二区国产精品久久精品| а√天堂www在线а√下载| 欧美激情在线99| 亚洲av免费高清在线观看| 无人区码免费观看不卡| av黄色大香蕉| 天堂网av新在线| 特级一级黄色大片| 小说图片视频综合网站| 美女大奶头视频| 我的女老师完整版在线观看| 亚洲熟妇中文字幕五十中出| 欧美性猛交╳xxx乱大交人| 国产高清有码在线观看视频| 国产精品美女特级片免费视频播放器| 国产精华一区二区三区| 97超级碰碰碰精品色视频在线观看| 真人一进一出gif抽搐免费| 蜜桃亚洲精品一区二区三区| 内地一区二区视频在线| 又爽又黄无遮挡网站| 亚洲avbb在线观看| 22中文网久久字幕| 久久久精品大字幕| 干丝袜人妻中文字幕| 麻豆av噜噜一区二区三区| 欧美xxxx黑人xx丫x性爽| 日本欧美国产在线视频| 精品久久久久久成人av| 老司机午夜福利在线观看视频| 色哟哟·www| 俄罗斯特黄特色一大片| 一边摸一边抽搐一进一小说| 久久午夜亚洲精品久久| 99久国产av精品| 成人一区二区视频在线观看| 国语自产精品视频在线第100页| 少妇人妻精品综合一区二区 | 不卡视频在线观看欧美| 色综合色国产| 国产爱豆传媒在线观看| 色综合亚洲欧美另类图片| 一个人免费在线观看电影| 国产淫片久久久久久久久| 最新在线观看一区二区三区| 精品国内亚洲2022精品成人| 亚洲在线观看片| 麻豆成人av在线观看| 欧美日本亚洲视频在线播放| 免费不卡的大黄色大毛片视频在线观看 | 国内久久婷婷六月综合欲色啪| 亚洲不卡免费看| 国产精品野战在线观看| 在线国产一区二区在线| 亚洲人成网站在线播放欧美日韩| 最好的美女福利视频网| 亚洲av美国av| а√天堂www在线а√下载| 久久九九热精品免费| 给我免费播放毛片高清在线观看| 97超级碰碰碰精品色视频在线观看| 特大巨黑吊av在线直播| 99精品在免费线老司机午夜| 国产探花极品一区二区| 亚洲中文字幕一区二区三区有码在线看| 97人妻精品一区二区三区麻豆| 久久久成人免费电影| 欧美性猛交╳xxx乱大交人| 99久久精品国产国产毛片| bbb黄色大片| 亚洲精品乱码久久久v下载方式| 真实男女啪啪啪动态图| 人妻制服诱惑在线中文字幕| 成人毛片a级毛片在线播放| 国产伦精品一区二区三区视频9| 日韩强制内射视频| 男人的好看免费观看在线视频| 国产日本99.免费观看| 99久国产av精品| 国产高清三级在线| 婷婷亚洲欧美| 一级a爱片免费观看的视频| 国产精品,欧美在线| 国产一区二区三区在线臀色熟女| 简卡轻食公司| 亚洲人成网站在线播| 真人做人爱边吃奶动态| 国产高潮美女av| 深夜精品福利| 天天躁日日操中文字幕| 国产人妻一区二区三区在| 亚洲第一电影网av| netflix在线观看网站| 欧美又色又爽又黄视频| 国产亚洲91精品色在线| 欧美xxxx黑人xx丫x性爽| 日本撒尿小便嘘嘘汇集6| 精品一区二区三区av网在线观看| 国产精品国产高清国产av| 欧美又色又爽又黄视频| 国产黄a三级三级三级人| 久久6这里有精品| 久久99热6这里只有精品| 国产aⅴ精品一区二区三区波| 亚洲va日本ⅴa欧美va伊人久久| 免费观看在线日韩| 男人舔女人下体高潮全视频| 国产黄色小视频在线观看| 国模一区二区三区四区视频| 韩国av在线不卡| 99久久成人亚洲精品观看| 美女被艹到高潮喷水动态| 91狼人影院| 国产女主播在线喷水免费视频网站 | 国产日本99.免费观看| 日本三级黄在线观看| 国产精品不卡视频一区二区| 色哟哟哟哟哟哟| 在线免费观看的www视频| 免费人成视频x8x8入口观看| 日本黄色片子视频| 中文字幕高清在线视频| 人妻夜夜爽99麻豆av| 97超视频在线观看视频| 欧美黑人巨大hd| 成年女人毛片免费观看观看9| 小蜜桃在线观看免费完整版高清| 身体一侧抽搐| 美女高潮喷水抽搐中文字幕| 俺也久久电影网| 国产激情偷乱视频一区二区| 成人国产麻豆网| 欧美色视频一区免费| 国产精品亚洲美女久久久| 久久99热这里只有精品18| 在现免费观看毛片| 亚洲国产精品合色在线| 综合色av麻豆| 高清在线国产一区| 无遮挡黄片免费观看| 一区二区三区免费毛片| 亚洲,欧美,日韩| 久久精品国产清高在天天线| 亚洲人与动物交配视频| 熟女电影av网| 狂野欧美激情性xxxx在线观看| 18禁黄网站禁片午夜丰满| 97人妻精品一区二区三区麻豆| 亚洲国产精品sss在线观看| 亚洲18禁久久av| 99在线人妻在线中文字幕| 欧美性猛交╳xxx乱大交人| 毛片女人毛片| 亚洲图色成人| 看免费成人av毛片| 美女免费视频网站| 一个人免费在线观看电影| 麻豆av噜噜一区二区三区| 久久久久免费精品人妻一区二区| 成人国产一区最新在线观看| 欧美精品国产亚洲| 午夜福利视频1000在线观看| 国产高清三级在线| 亚洲第一区二区三区不卡| 高清日韩中文字幕在线| 一本久久中文字幕| 欧美国产日韩亚洲一区| 成年版毛片免费区| 中文字幕av成人在线电影| 亚洲avbb在线观看| 国产av麻豆久久久久久久| 国产老妇女一区| 男女下面进入的视频免费午夜| 欧美zozozo另类| 欧美日韩亚洲国产一区二区在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲乱码一区二区免费版| 我要看日韩黄色一级片| 久久久久性生活片| 亚洲无线观看免费| 亚洲乱码一区二区免费版| 免费看美女性在线毛片视频| 免费看日本二区| 国产精品一区二区三区四区免费观看 | 国产高清三级在线| 身体一侧抽搐| 国产亚洲91精品色在线| 简卡轻食公司| 中文字幕av在线有码专区| av在线观看视频网站免费| 免费观看人在逋| 两人在一起打扑克的视频| 精品不卡国产一区二区三区| 亚洲av一区综合| 国内精品久久久久精免费| 国产一区二区在线观看日韩| 在线免费观看不下载黄p国产 | 乱系列少妇在线播放| 国产毛片a区久久久久| 国产精品日韩av在线免费观看| 久久精品国产亚洲av天美| 成人国产一区最新在线观看| 日本一二三区视频观看| 国产精品98久久久久久宅男小说| 亚洲无线在线观看| 啦啦啦韩国在线观看视频| 日韩精品有码人妻一区| 久久精品国产99精品国产亚洲性色| 国内精品美女久久久久久| 丝袜美腿在线中文| 不卡视频在线观看欧美| 极品教师在线免费播放| 成人欧美大片| 国产亚洲91精品色在线| av在线老鸭窝| 久久精品国产亚洲av香蕉五月| 欧美色欧美亚洲另类二区| 中文亚洲av片在线观看爽| 搡老妇女老女人老熟妇| 国产亚洲精品av在线| 五月玫瑰六月丁香| 黄色丝袜av网址大全| 亚洲乱码一区二区免费版| 久久久国产成人免费| 一个人看的www免费观看视频| 嫩草影院精品99| 久久久久久久亚洲中文字幕| 国产蜜桃级精品一区二区三区| 免费看av在线观看网站| av天堂在线播放| 亚洲最大成人手机在线| 国产精品久久久久久久电影| 亚洲精华国产精华精| 国产免费男女视频| 欧美性猛交╳xxx乱大交人| 亚洲三级黄色毛片| 日韩欧美三级三区| 日韩欧美国产在线观看| 久久久成人免费电影| 久久草成人影院| 成人国产一区最新在线观看| 可以在线观看的亚洲视频| 熟女电影av网| 日韩高清综合在线| av天堂中文字幕网| 毛片一级片免费看久久久久 | 色综合站精品国产| 国产探花在线观看一区二区| 亚洲自偷自拍三级| 久久久久久久亚洲中文字幕| 免费不卡的大黄色大毛片视频在线观看 | 男女那种视频在线观看| 日韩欧美一区二区三区在线观看| 欧洲精品卡2卡3卡4卡5卡区| 国产一区二区在线观看日韩| 欧美色欧美亚洲另类二区| 成人毛片a级毛片在线播放| 一本久久中文字幕| 国产精品98久久久久久宅男小说| 男女视频在线观看网站免费| 婷婷亚洲欧美| 亚洲av第一区精品v没综合| 中文字幕高清在线视频| 亚洲精品粉嫩美女一区| 又黄又爽又刺激的免费视频.| 热99在线观看视频| 极品教师在线视频| 97人妻精品一区二区三区麻豆| 一级黄片播放器| 国产高清视频在线播放一区| 成人性生交大片免费视频hd| 午夜福利高清视频| 别揉我奶头~嗯~啊~动态视频| 日本成人三级电影网站| 亚洲精华国产精华精| 看免费成人av毛片| 亚洲aⅴ乱码一区二区在线播放| 韩国av一区二区三区四区| 国产一区二区在线av高清观看| 男女那种视频在线观看| 久久久久久久久大av| 亚洲人成网站在线播| 成人三级黄色视频| 亚洲欧美精品综合久久99| 久久人人精品亚洲av| 18禁黄网站禁片午夜丰满| 麻豆久久精品国产亚洲av| 亚洲va在线va天堂va国产| 国产乱人伦免费视频| 日日干狠狠操夜夜爽| 国产精品99久久久久久久久| 男女之事视频高清在线观看| 不卡视频在线观看欧美| 成人午夜高清在线视频| 伦理电影大哥的女人| 日本免费a在线| 欧美高清成人免费视频www| 久久久久久国产a免费观看| 免费观看精品视频网站| 亚洲内射少妇av| 亚洲欧美日韩高清在线视频| 波多野结衣高清作品| 最近最新中文字幕大全电影3| 久久婷婷人人爽人人干人人爱| 中文字幕免费在线视频6| 国产 一区精品| 成人鲁丝片一二三区免费| 成年免费大片在线观看| 国产高清视频在线播放一区| 欧美xxxx性猛交bbbb| 国产高清视频在线播放一区| bbb黄色大片| 在线观看av片永久免费下载| 99热6这里只有精品| 欧美+亚洲+日韩+国产| 18禁黄网站禁片午夜丰满| 一进一出抽搐gif免费好疼| 久久久国产成人免费| 三级国产精品欧美在线观看| 别揉我奶头~嗯~啊~动态视频| 亚洲成人久久性| 精品欧美国产一区二区三| 国产极品精品免费视频能看的| 高清毛片免费观看视频网站| 永久网站在线| 可以在线观看的亚洲视频| 久久久久久九九精品二区国产| 国产精品嫩草影院av在线观看 | 亚洲黑人精品在线| 在线a可以看的网站| 色精品久久人妻99蜜桃| 全区人妻精品视频| 日本欧美国产在线视频| 国内少妇人妻偷人精品xxx网站| 国产人妻一区二区三区在| 22中文网久久字幕| 97超级碰碰碰精品色视频在线观看| av天堂在线播放| 97超级碰碰碰精品色视频在线观看| 国产av在哪里看| 麻豆av噜噜一区二区三区| 无遮挡黄片免费观看| 成年女人毛片免费观看观看9| 美女大奶头视频| 精品午夜福利视频在线观看一区| 亚洲美女视频黄频| 十八禁国产超污无遮挡网站| 欧美精品国产亚洲| 欧美潮喷喷水| 极品教师在线免费播放| 国内精品一区二区在线观看| 欧美xxxx黑人xx丫x性爽| 九九在线视频观看精品| a在线观看视频网站| 亚洲电影在线观看av| 在现免费观看毛片| 日韩精品青青久久久久久| xxxwww97欧美| 日韩中字成人| 高清在线国产一区| 国产亚洲91精品色在线| 国产精品人妻久久久久久| 欧美性猛交╳xxx乱大交人| 特级一级黄色大片| 亚洲五月天丁香| 亚洲成人免费电影在线观看| 1024手机看黄色片| 午夜福利在线在线| 人妻久久中文字幕网| 久久这里只有精品中国| 精品久久久久久久久久免费视频| 变态另类丝袜制服| 亚洲无线观看免费| 99久久无色码亚洲精品果冻| av.在线天堂| 国产真实乱freesex| 国产精品久久久久久亚洲av鲁大| 蜜桃亚洲精品一区二区三区| 色综合亚洲欧美另类图片| 男女那种视频在线观看| 亚洲av日韩精品久久久久久密| 韩国av一区二区三区四区| 一夜夜www| 99久久九九国产精品国产免费| 亚洲精华国产精华精| 日本黄色视频三级网站网址| 精品人妻熟女av久视频| 国产午夜精品久久久久久一区二区三区 | 一进一出抽搐动态| 国产午夜精品久久久久久一区二区三区 | 无人区码免费观看不卡| 亚洲精品日韩av片在线观看| 黄色丝袜av网址大全| 欧美黑人欧美精品刺激| 欧美高清成人免费视频www| 观看美女的网站| 免费无遮挡裸体视频| 老熟妇仑乱视频hdxx| 熟女人妻精品中文字幕| 久久人人精品亚洲av| 一卡2卡三卡四卡精品乱码亚洲| www.色视频.com| 欧美一级a爱片免费观看看| 精品久久久久久久久av| 午夜视频国产福利| 搡老妇女老女人老熟妇| 亚洲专区国产一区二区| 日韩亚洲欧美综合| 婷婷精品国产亚洲av在线| 国产v大片淫在线免费观看| 黄片wwwwww| 亚洲精品456在线播放app | 99国产极品粉嫩在线观看| 黄色视频,在线免费观看| 淫妇啪啪啪对白视频| 在线观看一区二区三区| 在线免费观看不下载黄p国产 | 欧美日韩乱码在线| 亚洲中文日韩欧美视频| 国产精品一区二区三区四区久久| 亚洲人成网站高清观看| 亚洲精品国产成人久久av| 亚洲va日本ⅴa欧美va伊人久久| 色精品久久人妻99蜜桃| 12—13女人毛片做爰片一| 久久99热6这里只有精品| 成人国产麻豆网| 久久久精品大字幕| 欧美性猛交╳xxx乱大交人| 老熟妇乱子伦视频在线观看| 午夜精品在线福利| 日韩精品青青久久久久久| 91麻豆av在线| 99热只有精品国产| 国产精品人妻久久久久久| 女人十人毛片免费观看3o分钟| 欧美区成人在线视频| 简卡轻食公司| 一级a爱片免费观看的视频| 精品人妻一区二区三区麻豆 | 岛国在线免费视频观看| 国产精品久久久久久久电影| 亚洲人成伊人成综合网2020| 日韩精品有码人妻一区| 国产探花在线观看一区二区| 制服丝袜大香蕉在线| 亚洲专区国产一区二区| 能在线免费观看的黄片| www日本黄色视频网| 九色国产91popny在线| 久久久久久大精品| 国产白丝娇喘喷水9色精品| 在现免费观看毛片| 亚洲自拍偷在线| 国产在视频线在精品| 色综合站精品国产| 99久久九九国产精品国产免费| 亚洲av免费高清在线观看| 国产三级在线视频| 尤物成人国产欧美一区二区三区| 丰满乱子伦码专区| 国产乱人伦免费视频| 久久久国产成人免费| 国产 一区精品| 久久国产精品人妻蜜桃| 欧美性猛交╳xxx乱大交人| 在线观看66精品国产| 乱系列少妇在线播放| 久久精品国产自在天天线| 精品久久久久久久久久免费视频| 免费大片18禁| 日本黄色视频三级网站网址| 日本黄色片子视频| 丰满人妻一区二区三区视频av| 亚洲精华国产精华液的使用体验 | 亚洲av一区综合| 久久香蕉精品热| 亚洲欧美激情综合另类| 免费av毛片视频| 国产精品人妻久久久久久| 最近最新免费中文字幕在线| 日本免费一区二区三区高清不卡| 午夜福利成人在线免费观看| 亚洲内射少妇av| 欧美激情久久久久久爽电影| 99精品久久久久人妻精品| 日韩中字成人| 窝窝影院91人妻| 99国产精品一区二区蜜桃av| 日韩欧美 国产精品| 老司机午夜福利在线观看视频| 88av欧美| 日本熟妇午夜| av女优亚洲男人天堂| 国产一区二区在线观看日韩| 成人一区二区视频在线观看| 亚洲国产精品成人综合色| 国产免费一级a男人的天堂| 免费看日本二区| 天堂动漫精品| 国内精品一区二区在线观看| 乱系列少妇在线播放| 精品国产三级普通话版| 麻豆一二三区av精品| 最近最新免费中文字幕在线| 99国产极品粉嫩在线观看| 亚洲天堂国产精品一区在线| 免费看美女性在线毛片视频| 精品一区二区三区视频在线观看免费| 桃色一区二区三区在线观看| 亚洲精品色激情综合| 91麻豆av在线| АⅤ资源中文在线天堂| 人妻制服诱惑在线中文字幕| 乱系列少妇在线播放| 国产精品av视频在线免费观看| 人人妻人人看人人澡| 亚洲精品久久国产高清桃花| 亚洲乱码一区二区免费版| 亚洲在线自拍视频| 欧美zozozo另类| 日韩一本色道免费dvd| 国产av一区在线观看免费| 俄罗斯特黄特色一大片| 国产一区二区三区视频了| 非洲黑人性xxxx精品又粗又长| 亚洲欧美日韩无卡精品| 亚洲男人的天堂狠狠| 他把我摸到了高潮在线观看| 亚洲欧美日韩高清在线视频| 免费高清视频大片| 欧美高清成人免费视频www| 偷拍熟女少妇极品色| 精品福利观看| 国产欧美日韩一区二区精品| 人人妻,人人澡人人爽秒播| 成人亚洲精品av一区二区| 亚洲精品色激情综合| 国产精品人妻久久久久久| 欧美激情在线99| 999久久久精品免费观看国产| 国产日本99.免费观看| 男女下面进入的视频免费午夜| 亚洲国产精品久久男人天堂| 他把我摸到了高潮在线观看| 久久精品国产鲁丝片午夜精品 | 中文字幕免费在线视频6| 国产伦精品一区二区三区视频9| 一卡2卡三卡四卡精品乱码亚洲| 亚洲美女视频黄频| 日本一本二区三区精品| 国产一区二区亚洲精品在线观看| 国产伦一二天堂av在线观看| 中亚洲国语对白在线视频| 两个人的视频大全免费| 免费看av在线观看网站| 男女做爰动态图高潮gif福利片| 成人综合一区亚洲| 欧美性感艳星| av天堂中文字幕网| 熟女人妻精品中文字幕| 18禁黄网站禁片免费观看直播| 老师上课跳d突然被开到最大视频| 麻豆国产97在线/欧美| 伊人久久精品亚洲午夜| videossex国产| 国产高清视频在线播放一区| 国产单亲对白刺激| 国产私拍福利视频在线观看| 18+在线观看网站| 天天躁日日操中文字幕| 亚洲成人久久性| 蜜桃亚洲精品一区二区三区| 免费不卡的大黄色大毛片视频在线观看 | 两个人视频免费观看高清| 又黄又爽又刺激的免费视频.| 国产成人福利小说|